Field reversed configurations (FRC) are characterized by a magnetic field topology, which exhibits the inversion of the external magnetic field through plasma sustained current and the subsequent presence of a null field surface. A monotonical radial decrease in the longitudinal magnetic field leads to the potential presence of harmonics of the ion cyclotron frequency of any order in the region included between the outer wall and the null field surface. What is the effective hot-plasma dispersion relation obtained through the convolution of a large ensemble of high harmonics fast waves (HHFW) confined in a finite radial region represents an open question that we attempt to address here. In particular, we discuss a combination of analytical modeling and numerical treatment, which allows us to retrieve the resulting high harmonic fast wave complex wavevector for any radial location of any FRC radial profile. Moreover, we show how the obtained hot-plasma HHFW wavevector relates to the cold-plasma solution, and how it depends on the plasma parameters.

## I. INTRODUCTION

In the energy balance of any nuclear fusion machine, plasma heating plays a fundamental role as it represents a process required for reaching fusion temperatures as well as for compensating energy losses. Historically, plasma heating has been carried out mainly through the injection of high energy neutral particles that are rapidly ionized in the plasma and transfer their energy to thermal ions first and then to electrons through Coulomb collisions. In addition to the injection of energetic particles, another mechanism that has been widely investigated for heating is the launch into the plasma of electromagnetic waves, which through different channels deposit energy in specific plasma regions. Plasma heating through radio frequency (RF) represents for present nuclear fusion experiments a rich field that includes several physical regimes with source frequency values that vary from below ion cyclotron to well beyond electron cyclotron. The main motivations that drive significant interest into the different RF schemes available are the flexibility of the mechanism in terms of operation points and the supposed efficiency in playing an auxiliary role with respect to neutral beam injection (NBI) to sustain both plasma temperature and non-inductive current drive.^{1}

TAE Technologies (TAE, hereinafter) is pursuing its scientific and technology program toward a nuclear fusion prototype through the design and construction of a series of field reversed configurations (FRCs) devices that so far have relied on NBI only to sustain the experimental plasma conditions.^{2–6} RF heating has not been applied in TAE's machines, yet it is under active study for plasma heating in future machines.^{7} For this reason, TAE has supported the development of a set of numerical tools, dubbed RFPisa, that are specifically designed for the FRC geometry and field topology to investigate possible RF scenarios. In fact, FRCs are so-called high-*β* devices, where *β* is the ratio of the plasma pressure to the magnetic pressure, and this implies that the wave accessibility constraints may differ significantly from what is typical in more standard low-*β* devices, leading to the impossibility to port directly to FRCs RF techniques and regimes that are well proved and established for other geometries. The RFPisa package includes a 1D full wave finite Larmor radius (FLR) solver, a 1D root finder for dispersion relation equations and an antenna Fourier decomposer, which allows to estimate the 3D field profiles generated by a 3D antenna model.

Among all possible RF scenarios that have been investigated in the past decades, a frequency regime that has attracted attention for applications in high-*β* devices is the intermediate frequency range between ion cyclotron and lower hybrid (LH), i.e., $ \Omega i \u226a \omega \u226a \omega LH$, which is usually referred to as a high-harmonic fast wave (HHFW) regime. In fact, in high-*β* devices, HHFW is expected to exhibit strong single pass absorption and, therefore, provide local power deposition.^{8} Modeling the aforementioned frequency regime requires a full wave approach, and while an extended literature on HHFW studies in tokamaks of different aspect ratios is available,^{9–14} there is not to the best of our knowledge previous full wave work that addresses HHFW in FRCs. Comparisons with other high-*β* devices such as spherical tokamaks are of marginal help because of the very different topology of the FRC magnetic field, and comparisons with low-*β* fusion plasmas are also infeasible. Nevertheless, there is some qualitative similarity with the reversed field pinch^{15} (RFP), and the methods that we have developed may be of use in that application also.

As a first step toward applying HHFW to FRCs, an antenna and power supply have been installed on the UCLA LAPD linear device, and an efficient RF-plasma coupling was found.^{7,16}

FRC's magnetic field radial profile is characterized by an inversion, and the subsequent presence of a null field point leads to a local cyclotron frequency that radially decreases first to zero and then increases again in magnitude although with negative sign. What is the resulting radial profile of the wavevector of a fast wave launched from the FRC high field outer boundary represents an intuitively puzzling question. Obtaining an analytical or numerical description of the HHFW wavevector is a necessary condition for any additional investigation related to wave-plasma coupling and energy deposition as it is a required quantity to compute the hot-plasma dielectric tensor and from there the wave fields. What can be easily envisioned is that given a source frequency value, a first harmonic whose order is mainly dictated by the external magnetic field (plasma current is negligible in the FRC outer part) is encountered at certain radial distance. The region between such a first harmonic and the null field point hosts all the other subsequent higher order resonances. The physical separation between the different resonances is given by the specific properties of the plasma equilibrium but because of the presence of a (theoretically) unlimited number of resonances confined in a finite region, and the distance between neighboring resonances reduces inevitably to zero with increasing order leading to a convolution of resonances. The novelty of the topic generates several open questions that include what are the effective dispersion function and propagation mechanisms in a field topology that is characterized by such a peculiar sequence of resonances. Answering these questions represents the first step required for investigating the efficiency and practicability of a possible HHFW system applied to FRC devices.

The outline of this paper is as follows: in Sec. II, we summarize the modeling that we adopt and that is based on the “quasi-local approximation” (QLA) wave equations derived by Marco Brambilla for the toroidal geometry;^{17} in Sec. III, we discuss the numerical technique that we apply to solve the high-order implicit equation obtained through the modeling; in Sec. IV, we show the FRC prototype equilibrium that has been selected to carry out our investigation, and in Sec. V, we present a series of results which are finalized to show how the obtained HHFW wavevectors scale vs both plasma and wave parameters. Finally, in Sec. VI, we draw the conclusions and provide an outlook for the next steps.

## II. MODEL EQUATIONS

In the following, we summarize the essential key steps of the quasi-local approximation model.^{17} The equations shown here include all the simplifications that apply and are written for the 1D cylindrical geometry taking into account periodic boundary conditions in the toroidal and longitudinal directions. The summary presented is consciously largely incomplete, but in our opinion, it is, nonetheless, sufficient to grasp the fundamental physics that is included in the effective equations that we solve and to develop a correct understanding of the results we present. We invite readers interested in a comprehensive description of the QLA model derivation to refer directly to Ref. 17.

*m*and

*n*represent the toroidal mode and longitudinal wavenumbers with

*n*as the multiplier of the minimum parallel wavevector $ k \u2225 0$. The current term $ J m , n ( r )$ can be written in terms of the kernel $ \Sigma \xaf ( m , n , r , r \u2032 )$ and the electric field $ E m , n ( r )$ as

*r*and

_{w}*r*is an integral function over the whole plasma domain. Implementing directly Eqs. (4) and (5) into Eq. (3) would lead to a tremendous set of integro-differential equations whose solution would be hardly feasible for any realistic configuration. A significant simplification can be achieved introducing the so-called eikonal form

^{18}for the radial dependence of the electric field, i.e.,

*ε*is the ratio of the wavelength to the equilibrium gradient lengths and $ \kappa m , n$ is the radial derivative of the eikonal function $ S m , n ( \epsilon r )$, i.e., $ \kappa m , n = \u2202 S m , n ( \epsilon r ) / \u2202 r$.

^{17,18}lead to

^{18}Expressing the curl with the aforementioned approximations, it is possible to recast the wave equation of Eq. (3) to the lowest order in

*ε*as

Once an explicit expression for the conductivity tensor is provided, the roots of the aforementioned Eq. (10) provide the desired dispersion relations.

^{19}as

*δ*is the angle between $ k \u22a5$ given by $ k \u22a5 = \kappa r \u0302 + k \theta \theta \u0302$ and the direction perpendicular to the magnetic field surface (which, in our 1D configuration, corresponds to the radial direction), and $ P \xaf p , \alpha $ is the

_{k}*p*th order tensor in local orthogonal coordinates whose elements are given by

*Z*is the plasma dispersion function, $ Z \u2032$ is its derivative,

*I*is the modified Bessel function of order

_{p}*p*, $ I p \u2032$ is its first derivative, and, finally, $ \nu \alpha , \u2009 x p \alpha $, and $ x p \alpha D$ are

*v*is the drift velocity and $ \Omega c \alpha $ and $ v th \alpha $ are the cyclotron frequency and thermal velocity of species

_{D}*α*, respectively. Each element of the dielectric tensor is given by a double summation over all plasma species and over all Bessel functions taken into account. The exact number of Bessel functions required for reaching convergence and a correct description of the region near the null field point depend on the plasma profile and the wave properties, but for the standard cases discussed here, we can consider that about 100 Bessel functions can be sufficient. Therefore, the final equation obtained through the determinant of the wave equations set represents a formidable algebraic expression. In fact, it is given by a summation of six triple products where each term in each triple product is, in turn, given by a summation over about one hundreds complex Bessel functions for each plasma species. It is worth noting that the unknown that we want to solve for is the radial complex component of the wavevector, i.e.,

*k*, and that in the above set of equations is “hidden” first in $ k \u22a5$ and then in the arguments of all the Bessel functions. The other two components, $ k \theta $ and $ k \u2225$ (i.e.,

_{r}*k*), are automatically given through the selection of the single (

_{z}*m*,

*n*) mode that we want to compute the wave equation for, i.e., $ k \theta = m / r$ and $ k \u2225 = n k \u2225 0$.

The validity of the solutions obtained from the afore described derivation requires the satisfaction of two main conditions,^{17} the radial wavelength is shorter than the equilibrium gradient scale length, and the integration domain should exclude mode-conversion layers. How these two conditions apply to our specific cases will be discussed in Sec. V together with the presentation of the results obtained.

## III. NUMERICAL APPROACH

Plasma equilibrium and source wave properties define uniquely all the coefficients and functions needed to define Eq. (10) but for *k _{r}*, which remains as the root that needs to be computed. Once the extended but straigthforward algebra is carried out, it is possible to proceed with a root finding technique to retrieve the admitted values of

*k*for any radial location. Although different numerical techniques and routines can be utilized, we found particularly convenient and efficient the so-called “downhill method,”

_{r}^{20}which has been rewritten from the original source and included in our numerical tools.

In addition to roots representing the HHFW wavevector, Eq. (10) admits several other solutions, and it is not, in general, possible to initialize the root finding process with seeds that will certainly converge to an HHFW value. In fact, in the hot-plasma regions, Bernstein waves are also present, and they have in the $ ( r , k r )$ space narrow intersections with HHFW waves. In the outer region of the FRC where temperature is significantly reduced and plasma current is negligible, one can expect that the HHFW solution does not differ significantly from the cold-plasma solution, and the latter can be utilized to benchmark and verify the root finding process.

We implement a two-step approach, first we retrieve all roots within a certain volume of the $ [ \mathbb{R} ( k r ) , I ( k r ) ]$ space, and then we proceed to separate the HHFW component through a tracer. The computation of the roots is based on the definition of a 3D lattice in the space $ [ r , \mathbb{R} ( k r ) , I ( k r ) ]$ with the initial seeds passed to the root finder distributed on the nodes of the lattice. Because of the reduced spatial distance between different resonances and the presence of several Bernstein waves distributed over different wavevector values, the size of the 3D lattice amounts up to thousands of nodes in the *r* direction and up to few tens in each *k _{r}* component direction. This makes the overall process for a standard case quantitatively cpu intensive as the same procedure must be repeated for each node.

The question “What is the highest order of the Bessel functions included in the expansion of the dieletric tensor” determines the extension of the region that can be treated, and, in particular, it defines the closest location to the null field point that can be properly solved. Here, we will show results obtained in the positive field region, but the negative field region can be treated with the same numerical approach.

## IV. GEOMETRY AND EQUILIBRIUM

Our study is one-dimensional in the radial direction, i.e., we assume a plasma column that is infinite in the longitudinal direction and uniform in the toroidal direction. These conditions are suitable to mimic the FRC midplane field topology and study RF scenarios that involve injection in the FRC central region. In Fig. 1, we show a typical two-dimensional FRC equilibrium, which exhibits a “quasi-one-dimensional” structure at the midplane, and the field inversion clearly visible through the direction of the field lines. Different injection schemes that include the proximity of the X points, and the mirror regions require an extension of the theoretical and numerical apparatus to the two-dimensional (*r*, *z*) space.

In order to carry out our numerical investigation, a 1D prototype FRC equilibrium, shown in Figs. 2 and 3, has been computed through TAE's code LR_eqMI.^{21} Although this particular equilibrium is not intended to reproduce a particular experimental case, it includes all relevant FRC midplane properties that are significant for the purpose of the present study, and its bulk parameters do not differ significantly from actual experimental values. In particular, we have background magnetic field $ B z \u223c 1000 G$, and the ion temperature and density at the null field point are $ T i \u223c 350 \u2009 \u2009 eV$ and $ n i \u223c 5.0 \xd7 10 13 \u2009 cm \u2212 3$, respectively. Therefore, the highest cyclotron frequency, located in the FRC outer region, is about 1.52 MHz, and in the following, we will utilize a source frequency value equal to 7.0 MHz or higher.

A key feature that affects the wave propagation and the HHFW wavevector limit value at the outer boundary is the presence of a low-density and low-temperature outer region. Plasma parameters of our prototype equilibrium satisfy in this outer region the requirements needed to apply the cold-plasma model. Moving inward from the outer boundary, we expect that our results match the analytical cold-plasma solution and then depart from that once the plasma temperature is high enough to break the cold-plasma assumptions. Our prototype equilibrium considers specific plasma parameter values, but moderate variations with respect to the equilibrium features considered are not expected to lead to significant changes in the obtained qualitative results.

## V. NUMERICAL RESULTS

Plasma data from the numerical equilibrium shown in Figs. 2 and 3 have been utilized to define all parameters required by the hot-plasma dielectric tensor as given by Eqs. (11) and (13). Although a lower number could still be sufficient, we have taken into account 150 Bessel functions and defined through Eq. (10) the algebraic equation that we need to solve. In Fig. 4, we show the collection of roots obtained through the deployment of TAE's numerical tool RFPisa for a given mode and frequency. For each radial location, we find multiple roots where each root corresponds to a particular wave type as Bernstein, fast, or slow. At a qualitative level, the fast wave can be easily distinguished from the other wave types as it exhibits all the expected resonances, and in low-temperature regions, the HHFW wavevector value is expected to replicate or be close to the value predicted by the cold-plasma model.

Through a specifically developed tracer, we can filter out and select the HHFW wavevector data from all others. In order to be effective, such a tracing process requires the relative incremental changes in the wavevector to be small, and therefore, a fine radial mesh is required to retrieve a 1D continuous function. The results shown here have been obtained utilizing between several tens and a few hundreds radial points per cm. Obviously, mesh requirements are, in principle, not the same over the whole radial range. Very high harmonics are closer and require higher resolution to be resolved.

In Fig. 5, we show the real and imaginary components $ \mathbb{R} ( k r )$ and $ I ( k r )$ of the HHFW vector as they have been selected by the tracer from Fig. 4, and we compare them with the corresponding cold-plasma analytical solution. Some observations are straightforwardly possible. In the outer cold-plasma region, there is a perfect agreement between the real component of the obtained numerical solution and the cold-plasma analytical solution, but we observe a dramatic departure between the two curves once temperature increases, and the high resonances are encountered. The exact match between our numerical result and the analytical cold-plasma solution plays also as a check of the validity and correctness of both the analytical equation we solve and the numerical technique we apply. In addition to the real component of *k _{r}*, the numerical solution provides also the imaginary component that is responsible for energy absorption and the possible evanescence of the wave. It is worth noting that for the first harmonics, an imaginary component different from zero is present in proximity of the resonances only, but moving inward, i.e., toward harmonics of higher order, $ I ( k r )$ shows oscillations with respect to non-zero values. This means that the wave is purely propagative in regions between the first harmonics only, and then, it is damped continuously. After a first increase in the resonance amplitude, both the width and the amplitude decrease with increasing order, and quite interestingly, the average value of $ \mathbb{R} ( k r )$ does not change significantly over the whole sequence of resonances. Moreover, it is possible to identify an asymptotic value that the HHFW wavevector converges to. This behavior represents a big difference with respect to the cold-plasma solution. The latter, in fact, diverges very rapidly, and the comparison between the two curves shows that the cold-plasma model cannot by any means be considered a valid approximation for this regime.

The specific behavior of the HHFW vector and its asymptotic behavior vs the different parameters cannot be predicted *a priori* from a simple analysis of Eq. (10). Therefore, we have carried out a series of simulations to verify how the single-mode result shown in Fig. 5 varies for other modes and frequencies. We start changing the value of the parallel component $ k \u2225$ and keeping constant the toroidal mode and the frequency. Results obtained for three different $ k \u2225$ values are shown in Fig. 6. We find that the behavior observed earlier at the higher harmonics (reduced resonance amplitude and asymptotic *k _{r}* value or trend) is anticipated at lower harmonics if the value of $ k \u2225$ increases. The differences in the average value of

*k*for different $ k \u2225$ values are marginal. Moreover, as a side note, it is worth noting that the radial distance between the expected resonances is well below the local wavelength. Therefore, the average value of

_{r}*k*is the main player in terms of expected physical effects.

_{r}In addition to verifying the dependence of the HHFW wavevector on $ k \u2225$, we have checked its dependence on the toroidal mode *m*. This is particularly relevant for experimental applications because the limited length of the antenna straps leads to a broad spectrum in the toroidal modes with higher toroidal modes still playing a significant role. We find that as long as the wave is propagative in the outer region, the differences in the results obtained for different *m* modes are marginal and do not play a significant role. At higher *m* values, an evanescent outer region could be present, see Fig. 7, and that could largely reduce the wave-plasma coupling efficacy as the propagative part of the wave that starts at a much inner location (at $ r \u223c 60 \u2009 cm$ for the case shown in Fig. 7). This is clearly a key feature to take into account for the design and location of a possible antenna system.

The radial profiles of Figs. 2 and 3 and the results shown in Figs. 5 and 6 allow us to verify the consistency of the obtained solutions with respect to the conditions required by the QLA approach. The average value of the magnetic field gradient scale length and the average value of the radial wavevector show that the requirement for the radial wavelength to be shorter than the equilibrium typical gradient scale length is satisfied for the large part of the radial range. In fact, such a condition is not met in proximity of the null-field point only where, however, the asymptotic behavior of the wavelength average value is already well established, and it appears not to be affected by variations of the accuracy level applied to the treatment of the narrow region around the null-field point. Similarly, expected mode-conversion layers are relatively sparse in the lower-order side of the sequence of high harmonics, and they become very dense for higher values of the harmonic order. In the latter case, although the density of conversion layers is significant, there are two effects that must be taken into account for an appropriate evaluation of the physical consequences of the integration over this region. The relative variation in $ k \u22a5$ becomes smaller, and the width of the layers becomes much narrower than the average wavelength. Such conditions lead to an integral over a region characterized by rapid oscillations that cannot be singularly resolved and contribute through their average value only. A more detailed discussion of a specific treatment for the proximity of the null-field point will be discussed elsewhere.^{22}

Finally, we have explored what is the dependence of the HHFW wavevector on the source frequency. Results are summarized in Fig. 8, where we show the same mode computed for three different frequency values. While the qualitative behavior is similar for the three cases, we observe, as expected, a different distribution of the resonances, but, contrary to what has been found for the dependences on *m* and $ k \u2225$, we note that the average and asymptotic values of *k _{r}* increase with the increase in the source frequency. Comparing the changes in the

*k*asymptotic value to the source frequency changes, we see that the two quantities appear to be strictly proportional, suggesting a constant limit phase velocity of the HHFW wave. Such a feature will be further investigated in future studies in combination with the transition from HHFW to helicon obtained through the increase in the source frequency.

_{r}## VI. CONCLUSIONS AND OUTLOOK

FRC devices are qualitatively distinct from low-*β* fusion plasma configurations in that the plasma current creates a locally reversed magnetic field, which is stronger than the background field provided by the coils. As a result, the magnetic field shows a null-field surface and an inversion in the inner part of the FRC. This makes direct comparisons with low-*β* fusion plasmas infeasible, and therefore, we have developed a methodology to describe the effects of reversed magnetic field profiles on the HHFW regime. Investigating such a regime appears as a challenging theoretical task because of the expected unlimited sequences of harmonics confined in the two regions on each side of the null-field surface. In order to study this regime, the quasi-local approximation^{17} has been reformulated to adapt it to the 1D cylindrical geometry, and the newly obtained QLA formulation has been utilized by TAE's numerical package to compute the effective dispersion relation. Such a combination between the QLA model and the RFPisa code has generated the capability to compute the HHFW wavevector straightforwardly for any mode. Knowing what are the profile and the average value of the HHFW wavevector is a strict requirement for any future investigation related to energy deposition and current drive. The results presented here have shown the inadequacy of the cold-plasma approximation to treat this regime. In particular, our results show the presence of an asymptotic value for *k _{r}* that depends mostly on the frequency and marginally on the toroidal and longitudinal mode numbers. The role played by the individual resonances is found to decrease with increasing values of $ k \u2225$ with the asymptotic wavevector value playing the effective role in the final physical results. Our next step is given by the utilization of the capability developed and presented here to compute the high-order terms of the hot-plasma dielectric tensor and fill the wave equation (details will be presented elsewhere

^{22}) that can include the response of both the thermal plasma and the fast ions. The exact modeling of the latter is still to be determined, but as a first approximation, fast ions can be included taking the equivalent Maxwellian of the neutral beam slowing-down distribution. The solution of the wave equation will provide all the quantities that are necessary to compute what is the energy deposited into each plasma and its radial distribution paving the way to the design of a possible RF system to be deployed in an FRC machine.

## ACKNOWLEDGMENTS

F.C. and L.G. are very grateful to Marco Brambilla for the enlightening discussions they had in the past and for having provided valuable scientific materials. Initial work on the HHFW tracer by James W. S. Cook is gratefully acknowledged.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Francesco Ceccherini:** Conceptualization (equal); Formal analysis (supporting); Investigation (equal); Methodology (equal); Software (lead); Validation (equal); Writing – original draft (lead). **Laura Galeotti:** Conceptualization (equal); Formal analysis (lead); Investigation (equal); Methodology (equal); Software (supporting); Validation (equal); Writing – review & editing (equal). **Daniel C. Barnes:** Writing – review & editing (equal). **Sean Dettrick:** Project administration (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.