D full-wave simulation of waves in space and tokamak plasmas

Simulation results using a 2D full-wave code (FW2D) for space and NSTX fusion plasmas are presented. The FW2D code solves the cold plasma wave equations using the finite element method. The wave code has been successfully applied to describe low frequency waves in planetary magnetospheres (i.e., dipole geometry) and the results include generation and propagation of externally driven ultra-low frequency waves via mode conversion at Mercury and mode coupling, refraction and reflection of internally driven field-aligned propagating left-handed electromagnetic ion cyclotron (EMIC) waves at Earth. In this paper, global structure of linearly polarized EMIC waves is examined and the result shows such resonant wave modes can be localized near the equatorial plane. We also adopt the FW2D code to tokamak geometry and examine radio frequency (RF) waves in the scape-off layer (SOL) of tokamaks. By adopting the rectangular and limiter boundary, we compare the results with existing AORSA simulations. The FW2D code results for the high harmonic fast wave heating case on NSTX with a rectangular vessel boundary shows excellent agreement with the AORSA code.


Introduction
In the solar system, many planets including Mercury, Earth, Jupiter, Saturn, Uranus, and Neptune, and their moons are magnetized.The solar wind acts to confine the planetary magnetic field in a region known as the magnetosphere and interaction of the solar wind with these magnetospheres leads to the formation of large scale field aligned currents along with various electrostatic and electromagnetic waves.
Ultra-low frequency (ULF) waves are one example of such magnetospheric waves that exist in the lowest observed frequency range (mHz-Hz).One very common ULF wave occurs near the ion cyclotron frequency (Ω ion ) and is called the ElectroMagnetic Ion Cyclotron (EMIC) wave.
ULF waves can lead to significant acceleration, transport, and loss of electrons during magnetic storms [1].They can be a source for plasma transport at the magnetopause and have been shown to be responsible for ion outflows [2][3].Such outflows may dominate the energy of the ring current during magnetic storms and lead to significant modifications of magnetospheric dynamics.EMIC waves are also considered to be a dominant radiation belt loss mechanism leading to rapid pitch angle scattering radiation belt electrons of almost all energies on the time-scale of seconds with subsequent loss to the ionosphere [4].In addition, these waves have also been used as a diagnostic tool to measure plasma densities [5][6].
There are many numerical efforts to examine ULF waves in planetary magnetospheres.To date, many models for waves in the magnetosphere are based primarily on WKB-based methods, such as local dispersion relation analysis or ray tracing models [7][8].While these models can follow the flow of energy due to a particular excited wave, this approach will fail if the wave enters a region in space in which it can couple power to other propagating modes via tunneling and /or mode conversion.One the other hand, time-dependent fluid wave models have been used for MHD waves [9] and EMIC waves [10].Although they successfully demonstrate mode conversion process for MHD and EMIC waves, since MHD wave models are only valid for waves with ω << Ω ion , where ω is angular frequency, and multi-ion fluid wave code adopted a slab coordinate so that the magnetic field curvature effect on waves could not be captured.
A recently developed 2D full-wave (FW2D) code using finite element method [11] overcomes these shortcomings using an approach that describes wave propagation, mode effects for arbitrary plasma and magnetic field configurations.The FW2D wave code has been successfully used to describe low frequency waves in Earth and Mercury's multi-ion magnetospheres (i.e., dipole geometry).The results include generation and propagation of externally driven ULF waves via mode conversion at Mercury [11].Ref. [11] showed that efficient mode conversion from the fast compressional waves to the ion-ion hybrid (IIH) resonance occurs at Mercury and such mode-converted waves globally oscillate similar to field-line resonance at Earth.Mode coupling, refraction and reflection of internally driven field-aligned propagating left-handed EMIC waves at EPJ Web of Conferences 157, 02005 (2017) Earth [12] has also been demonstrated.Ref. [12] suggested that polarization reversal and mode conversion occurs at the crossover frequency, where the left-hand (LHP) and right-hand (RHP) polarization modes cross each other in the dispersion relation, from LHP to RHP polarization and such RHP EMIC waves can propagate to the inner magnetosphere, which is consistent with 1D full-wave calculations [13,14].
Furthermore, we have adapted the FW2D code to tokamak geometry and have used it to examine radio frequency (RF) waves in the scape-off layer (SOL) of tokamaks, the region of the plasma between the last closed flux surface (LCFS) and the tokamak vessel.Taking into account the SOL region is an important aspect because the interaction between fast waves and the SOL region can be significant.More details, specifically, for National Spherical Torus eXperiment (NSTX) can be found in Refs.[15][16][17][18][19][20][21].This code is ideal for waves in SOL plasma, because realistic boundary shapes and arbitrary density profiles can be easily adopted in the code and the SOL plasma can be approximated as cold plasma.
In this paper, we use the FW2D code to examine EMIC waves at Earth's magnetosphere and RF waves in the SOL region of a tokamak plasma.It is found that (a) linearly polarized EMIC waves generated via mode conversion are localized near the magnetic equator; (2) FW2D and AORSA results show an excellent agreement in the rectangular boundary; and (3) FW2D results with a realistic limiter boundary are significantly different to results with the rectangular vessel boundary.

FW2D Model Description
We have developed a 2D full wave code, so called the FW2D code, appropriate for general geometries [7], which solves the cold plasma wave equations in twodimensions.Assuming time dependence, exp(−iωt), the linear, cold plasma wave equation takes the form: where E is the perturbed electric field, c is a speed of light, ε is the dielectric tensor, and j e is the (localized) external current source, which generates the waves within our model.In space plasma, the source could be instability in the equatorial region excited, for example, by free energy in velocity space, and its structure could be specified from kinetic simulations in the equatorial magnetosphere.Alternatively, the source could be external compressions in the outer magnetosphere that propagate energy into the inner magnetosphere.
The dielectric tensor ε is naturally expressed in coordinates aligned along and across the local ambient magnetic field direction b=B/|B|.For an axisymmetric magnetospheric and tokamak model, Eq. ( 1) is most naturally expressed in cylindrical (r, z, φ) coordinates.The field is represented in terms of its projections along (b) and perpendicular (ξ, η) to B, thus we pick a local and orthogonal basis (b, ξ, η=b×ξ).For magnetospheric plasma, we define ξ = φ with the assumption that B ⋅ φ=0 while for tokamak geometry, we define ξ = r×b.
The wave equation in Equation ( 1) has been solved on an unstructured triangular mesh using the finite element method.One advantage of using the finite element method is that the local basis functions that are employed can be readily adapted to boundary shapes and can be packed in such a way as to provide higher resolution in regions where solutions may exhibit singular behavior.

Global structure of EMIC waves at Earth
We demonstrate the linearly polarized EMIC waves at Earth using the FW2D code.EMIC waves are in the Pc 1-2 (0.1-5Hz) frequency range, which are excited below the proton gyrofrequency and commonly observed in the magnetosphere.The polarization is one of the key characteristics that can be used to identify wave properties of EMIC waves.Their polarization is generally reported to be LHP, however, RHP and linear polarization have also been reported.Recent statistical studies showed that the occurrence of linearly polarized waves is comparable to LHP waves even in the source region [22].Furthermore, these waves predominantly occur in the inner magnetosphere independent of magnetic local time (MLT) [23].
The dispersion relation of heavy ions indicates that any branch of the wave mode should be predominantly LHP or RHP.Since linear polarization occurs only near the crossover frequency or at oblique propagation near the multi-ion hybrid resonances, the linearly polarized EMIC waves are suggested to occur at such critical frequencies [6,[24][25].
Ref. [6,24] calculated the efficiency of mode conversion of a radially propagating compressional wave at the IIH resonance using a simplified 1D slab cold plasma model that captures the essential features of the IIH resonance, however, global structure at Earth's magnetosphere has not been investigated.
By adopting realistic plasma density model, we examine global structure of such resonant waves using the FW2D code.In the simulation, we adopt a dipole magnetic field.We also use the global core plasma model (GCPM) [26], which provides an empirically derived plasma density as a function of geomagnetic and solar conditions throughout the inner magnetosphere.Figure 1 shows the plasma densities calculated by the GCPM at the MLT=12, where is on the magnetic meridian (magnetic line of longitude) facing the Sun, on 1/22/2000 for Kp = 1.Here, r and z are the radial and south-north directions normalized to Earth's radius (R E ).We also show contour of the Alfvén speed (V A ) in Figure 2. Because the electron density profile has an asymmetry structure between northern and southern hemispheres, the Alfvén speed profile in Figure 2 also shows asymmetry.In the inner magnetosphere the minimum value of Alfvén speed occurs at the magnetic equator near r ~ 3.2 R E (i.e., in the plasmasphere).
The waves are initially launched by specification of j e = (0, 0, j ϕ ) at r = 6.4 R E near the magnetic equator as a compressional fast waves (FW) and the frequency is f = 0.4 Hz, which is between local oxygen (Ω O ) and helium (Ω He ) cyclotron frequencies.Reflecting boundary conditions are prescribed at the Earth's surface, while absorbing boundary conditions are used at all other boundaries.
Figure 3  The FW mode with a long wavelength along the z direction appears in the E ϕ component in Figure 3(b).For ω<<Ω e , ω pe , where ω pe and Ω e are electron plasma and cyclotron frequencies, dispersion relation of the FW mode can be written as where n || (n⊥) is refractive index parallel (perpendicular) to the ambient magnetic field, and ε R , ε L , and ε S are the components of the cold dielectric tensor [27].Equation (2) exhibits two cutoffs at n || 2 ~ ε R(L) .FW mode launched near r = 6.4 R E propagates toward Earth quasiperpendicular to the ambient magnetic field.These waves reach the cutoff location near r ~ 4 R E where ε L = 0 is satisfied and are partially reflected.The rest of the FW energy tunnels through the stopgap between the cutoff and resonance (n || 2 ~ ε S ~ 0) and finally encounters the ionosphere.The FW amplitude in the inner magnetosphere is concentrated/trapped near the equatorial plane where Alfvén speed has minimum value as shown in Figure 3.
The IIH resonant waves appear only in E η showing field-aligned structure (See Figure 3a).Because the waves are initially launched as a pure FW mode, the waves exhibiting field-aligned structure in E η provide evidence of mode conversion [11].The IIH resonant waves reach higher magnetic latitude, but the waves reflect before reaching the altitude where the frequency matches the oxygen gyrofrequency, which is unlike to the IIH resonant waves at Mercury [11].Such localization of the EMIC waves, which is predicted by Ref. [28], is similar to the LHP EMIC waves [12].Ref. [12] also showed that the EMIC waves generated in the LHP propagate along the field line and reach the heavier ion cyclotron frequency.However, the waves cannot tunnel through the cut-off and heavier ion cyclotron frequency even though the heavy ion concentration ratio is 5% (because the stop gap is too large compared with the wavelength).Because EMIC waves are simultaneously detected at the ground as well as space on conjugate field lines [29][30], it is particularly confounding that the full wave calculations show that the waves cannot reach higher magnetic latitude, and the resolution of this descrepancy likely involves kinetic effects.

HHFW propagation in the SOL of tokamak
Recent experimental studies of high harmonic fast wave (HHFW) heating on the NSTX [31] have demonstrated that up to 60% of the coupled HHFW power can be lost in the SOL (i.e., the region of the plasma between LCFS and tokamak vessel, where the magnetic field lines are open and the plasma density has a sharp profile) [15][16].
Numerical efforts to examine waves in the SOL have been performed using the AORSA code [32] and they showed a large amplitude electric field occurs in the SOL when edge densities are high enough that the fast waves can propagate in the SOL [20][21].However, they adopted a rectangular boundary, which differs from the vacuum vessel structures.For comparison we adopt (1) a rectangular boundary to benchmark with the AORSA results and (2) a limiter boundary to examine boundary effects on HHFW propagation in the SOL of NSTX.

Comparisoin between FW2D and AORSA
In order to verify the FW2D, we have performed a benchmark between the AORSA and FW2D codes for a rectangular boundary case as shown in Ref. [20].
Figure 4 show the wave total amplitude of the wave electric field (normalized to its maximum value) calculated by AORSA (Figs.In this figure, we define that r and z are radial and vertical directions of the tokamak, respectively, while ϕ is a toroidal direction.For the FW2D simulation, we adopt radial direction between 0.8m < r < 1.7m to focus on the SOL and reduce computational time.We also assumed a strong collisions at the inner boundary near r = 0.8m to avoid the wave reflection at left boundary.The wave frequency is f = 30 MHz and the antenna is assumed to be located at r = 1.575m.An electron density is assumed to have exponential profile as used in AORSA simulations [20][21] and electron density in front of the antenna is N ant = 1.5 × 10 18 m -3 in electrondeuterium plasma.
In Figure 4, AORSA and FW2D simulations show excellent agreement.Even in the core plasmas, although the cold plasma model in FW2D, the FW2D solution shows a wave field behavior similar to the AORSA results.Therefore, these results confirm that the FW2D code is successfully verified.

Limiter boundary vs. rectangular boundary
Although previous numerical studies using AORSA clearly demonstrated that substantial HHFW power loss occurs in the SOL, the rectangular boundary differs from physical vacuum vessel structure and could affect the wave solutions.
In order to overcome this limitation, we adopt a limiter boundary to the FW2D code similar to Ref. [33] for various N ant as shown in Figure 5.In this case, we assume that the antenna is located along the limiter EPJ Web of Conferences 157, 02005 (2017) boundary for -0.5m < z < 0.5m.Similar to rectangular boundary cases as shown in Figure 4, we adopt an artificial inner boundary at r = 0.7m with strong collision.However we assume the limiter to be a perfect conductor, thus waves totally reflect at the boundary.For comparison, wave solutions from the AORSA code are also presented in Figure 5(d)-(f).
For N ant = 1 × 10 m -3 in Figure 5(a) and 5(d), because the FW cut-off is 'closed', the FW field is strongly localized in front of antenna and waves are evanescent.On the other hand, for N ant = 2.0 × 10 18 m -3 and 4.0 × 10 18 m -3 , the electron density is sufficiently large to 'open' the FW cut-off and thus standing waves clearly appear in the SOL.
The FW2D simulation shows that the maximum amplitude of the wave electric field for N ant = 4.0 × 10 18 m -3 seems to be stronger than for N ant = 2.0 × 10 18 m -3 .On the other hand the AORSA shows the maximum amplitude of standing wave modes occur for N ant = 2.0 × 10 18 m -3 and the peak of the wave electric field in the SOL for the AORSA simulation occurs outside region of the limiter boundary.In addition, for N ant = 4 × 10 18 m -3 , standing wave amplitude for z > 0 is stronger than for z < 0, which is different with respect to the AORSA showing stronger wave electric field appear at z < 0.
Therefore, the preliminary results in this paper clearly show that the boundary shape can significantly affect the FW propagation in the SOL and it is definitely necessary to adopt a realistic boundary shape of tokamak to investigate RF wave propagation.

Summary
In this paper, we present wave simulations using 2D fullwave FW2D code in Earth's magnetosphere and in the NSTX spherical tokamak.For the Earth case, we show strong mode conversion from incoming FW to linearly polarized IIH resonant waves occurs and such modeconverted waves are localized near the equatorial plane.Preliminary simulations for HHFW in NSTX show that the wave electic field in the SOL plasma is sensitive to the boundary shape (rectangular vs limiter boundaries) and suggest that the realistic boundary should be considered to determine wave power loss in the SOL.
The work at the Princeton University was supported by NASA grants (NNH13AV37I and NNH14AY20I) and DOE contract DE-AC02-09CH11466.
shows two-dimensional wave solutions of the electric field components (a) |E η | and (b) |E ϕ |, respectively.In this figure, white dotted lines represent where the wave frequency matches Ω O or Ω He , and the wave amplitude of |E η | and |E ϕ | is normalized to the maximum value of the |E η |.
a and c) and FW2D (Figs. b and d) for a single toroidal mode n ϕ = -21 (Figs. a and b) and n ϕ =-12 (Figs.c and d) for NSTX shot 130608.

Fig. 2 .
Fig. 2. Alfven speed profile calculated using density profiles in Figure 1 and dipole magnetic field configuration.Yellow bar near r/R E = 6.4 is where the initial FW waves are launched.

Fig. 3 .
Fig. 3. Wave solutions using the FW2D; (a) |E η |, electric field amplitude normal direction to local magnetic field line, (b) |E φ |, electric field amplitude azimuthal (toroidal) component, respectively.Here each power has been normalized to the maximum value of |E η |.Bold white dotted lines indicate the oxygen (Ω O ) and helium (Ω He ) cyclotron freuqencies, and light dotted lines are magnetic field line at r equator /R E = 4.8, 5.3, and 5.8, respectively.

Fig. 4 .
Fig. 4. Two-dimensional simulation results of the total electric field amplitude using the AORSA and FW2D codes for (a-b) n ϕ = -21 and (c-d) n ϕ = -12 for N ant = 1.5 × 10 18 m -3 .Each electric field amplitude has been normalized to its maximum value.White and red solid lines are LCFS and FW cut-off, respectively.

Fig. 5 .
Fig. 5. Two-dimensional electric field amplitude from (a-c) FW2D simulation with a limiter boundary and (d-f) AORSA simulation assuming a rectangular boundary for N ant = 1, 2, and 4 × 10 18 m -3 , respectively.In this figure, white, black and light red solid lines are the LCFS, the limiter boundary, and the FW cut-off, respectively.