A Singularity-Free Coordinate System for X-point Geometries

A Flux Coordinate Independent (FCI) approach for anisotropic systems, not based on magnetic flux coordinates has been introduced in [F. Hariri and M. Ottaviani, Comput. Phys. Commun., 184,2419 (2013)]. In this paper, we show that the approach can tackle magnetic configurations including X-points. Using the code FENICIA, an equilibrium with a magnetic island has been used to show the robustness of the FCI approach to cases in which a magnetic separatrix is present in the system, either by design or as a consequence of instabilities. Numerical results are in good agreement with the analytic solutions of the sound-wave propagation problem. Conservation properties are verified. Finally, the critical gain of the FCI approach in situations including the magnetic separatrix with an X-point is demonstrated by a fast convergence of the code with the numerical resolution in the direction of symmetry. The results highlighted in this paper show that the FCI approach should be able to address turbulent transport problems in X-point geometries.


I. INTRODUCTION
The presence of a magnetic field in plasmas is known to introduce a strong anisotropy in the system. This is is commonly met in fusion and most astrophysical plasmas. Developing numerical codes that take into account the strong plasma anisotropy and make efficient use of computational resources is crucial, for instance, to simulate ITER-like tokamak plasmas. Efforts have been made in this direction that lead to introducing optimized coordinate systems, so-called field-aligned coordinates, that take advantage of this anisotropy. Two-dimensional field-aligned coordinates, based on magnetic flux variables, were first presented in [1][2][3][4][5][6]. They have a fundamental drawback of not handling situations including the magnetic separatrix, due to the singularity of the field-aligned metric. The approach proposed later in [7] avoids this shortcoming, but still relies on flux coordinates, for instance (r, θ). A generalized 3-dimensional coordinate system, referred to as the Flux Coordinate Independent (FCI) system, not based on magnetic flux coordinates, is introduced in [8] and elaborated in [9]. The FCI coordinate system is constructed in a way that avoids the use of flux variables to discretize the fields in the poloidal plane. As discussed in [8,9], it allows the parallel derivative to be constructed by tracing the magnetic field lines from one perpendicular plane to the next, and interpolating to find the desired quantity. This frees us from using flux coordinates in the perpendicular plane, thus allowing for complex magnetic geometries free of singularities. The aim of the present work is to highlight the advantage underlying the use of the FCI field-aligned system for X-point magnetic configurations. Examples of such configurations are magnetic islands (driven by magnetohydrodynamics modes or embedded in the magnetic equilibrium such as in stellerators) and axi-symmetric X-points on the last closed field surface (LCFS) of tokamak plasmas.
These configurations, found in fusion as well as astrophysical plasmas, have a singularity in the metric of a coordinate system that uses the magnetic flux as a coordinate. This difficulty originates in the nature of the flux function, which has a saddle point (known as X-point) in one or more locations. Using the flux coordinate independent (FCI) approach avoids this problem and opens the way for the numerical simulations of plasma turbulence in a tokamak geometry encompassing both closed and open magnetic field line regions. Firstly, this is a necessary feature of codes addressing the question of the transition from the low (L) to the high (H) confinement regime in X-point geometry. It is worth reminding that ITER would operate in the H confinement mode to achieve the expected thermonuclear fusion power. Understandably, research in this direction is a hot topic and understanding turbulence in X-point geometries is crucial for predicting the L-H power threshold. Secondly, studies of the interaction between micro-turbulence and magnetic islands exhibit a significant interplay affecting their dynamics. The most recent investigations on this topic are reported in [17][18][19].
The ultimate result of this paper is showing that the FCI approach allows, in particular, not only a more natural treatment of the operations in the poloidal plane, but it also easily deals with X-point configurations and with O-points such as the magnetic axis, since it is constructed on coordinate systems with non-singular metric. For this purpose, we use the FENICIA code presented in [8,9], where the Flux indepENdent fIeld-aligned CoordInate Approach has been implemented. FENICIA is based on a modular numerical scheme specially designed to address the anisotropic transport of any set of equations belonging to the general class of models (1) in [8,9]. Since this code is easily extensible to different meshes and coordinate systems, we use it here to demonstrate the application of the FCI coordinate system to a magnetic island in a slab version of the code. The outline of this paper is as follows. First of all, analytic solutions of a sound-wave propagation model for an island magnetic equilibrium are presented in Section II. Then numerical tests are performed at the exterior and at the interior of the island to show the convergence of the numerical solution to the analytic one in Section III. Across the separatrix, where there are not yet known analytic solutions, two conclusive numerical convergence tests are presented in Section IV allowing us to show the robustness of the FCI approach in X-point and O-point magnetic configurations.

II. SOUND WAVE MODEL IN A MAGNETIC ISLAND GEOMETRY
Consider a class of static low-β equilibria, such that the suitably normalized axisymmetric magnetic field is given by where one employs a three dimensional Cartesian reference system (x, y, z) such thatẑ is the direction of the magnetic axis, the main magnetic field along z is constant and normalized to unity, and b(x) is the poloidal magnetic field in the poloidal plane (x, y). The two-dimensional vector x indicates the position in this plane. The poloidal field can be written in terms of a flux function ψ(x) such that Magnetic surfaces can be labeled by the value of ψ. Both closed and open field lines can be treated. As extensively discussed in [10], slowly growing tearing instabilities reconnect magnetic flux surfaces to form magnetic islands. Analytical methods were derived to examine the linear stability and radial distribution of tearing modes [11][12][13][14].
The nonlinear growth of these modes was also investigated in [15] and extended to cases where current nonlinearities are important as shown in [16].
Consider a magnetic equilibrium characterized by a magnetic island whose half radial width is equal to δ x = 2 √ A, with A being a parameter. Such an equilibrium corresponds to a magnetic field given by Eqs. (1)-(2), with: The considered domain, for an island centered at x = 1, is x min ≤ x ≤ x max (in practice, x max = 1 + ∆ and x min = 1 − ∆, with 0 < ∆ < 1) and −π ≤ y ≤ π, with 0 ≤ z ≤ 2π.
Let us consider the following sound-wave model that propagate parallel to the magnetic field φ is the electrostatic potential and u is the wave's parallel velocity. We define two dimensionless parameters: C = a/R × 1/ρ * where a is the tokamak minor radius, R is the tokamak major radius, ρ * = ρ s /a is the reduced gyro-radius with ρ s being the ion sound Larmor radius; and τ = 1. The parallel gradient operator is defined as follows: and the Poisson bracket is defined by [ψ, ·] = ∂ x ψ∂ y − ∂ y ψ∂x. The aim of the present section is to construct analytic solutions of model (4) in the magnetic equilibrium defined above, i.e. in the presence of a magnetic island. To do so, it is convenient to employ flux coordinates. It readily appears that ψ is a label of magnetic surfaces, since ∇ ψ = 0.
Let us introduce the following new set of coordinates: The coordinate η is a straight-field-line angle-like variable and has still to be calculated. As detailed in A, for η to be an angle ranging from −π to π, g(ρ) has to be given different expressions depending on whether ρ is bigger or smaller than 1, ρ = 1 being the radius at the separatrix of the island. There are indeed three regions, each with its own specific expression of g and η. The interior region of the island corresponds to values of ρ < 1, while ρ > 1 characterizes the exterior region of the island. One can actually distinguish two exterior regions: the left side and the right side of the island. g(ρ) in both of the exterior regions has the following expression: where the sign ± depends on whether one considers the right side (+) or the left side (-) of the exterior region (cf. A). K is the elliptic integral of the first kind: Conversely, g(ρ) can be shown to be as follows for the interior region ρ < 1: To illustrate this, contour plots of ρ and η are shown in figures ??-?? where ρ defines the island geometry and η defines the inner and outer regions of the island. It can be shown that the parallel gradient takes the following expression in terms of this new set of coordinates: Let us then look for wave-like solutions of model (4) of the form: with (m, n) standing for the wave numbers in η and z respectively, and ω(ρ) being the mode frequency. Injecting these expressions in Eq. (4) leads to the following system: The eigenfrequency solution of the dispersion relation depends on ρ: The eigenvectors are then For a given expression of φ 0 (ρ), u 0 can be calculated by using Eq. (15), with ω ± mn (ρ) given by Eq. (14). Then, Eq. (12) provides analytic solutions of the model (4), which will be compared to their numerical counterparts in the next section. While the above solutions can be derived for ρ < 1 and ρ > 1, there is no obvious analytic solution to this problem at ρ = 1. This issue is left for future investigation.

III. NUMERICAL TESTS AT THE EXTERIOR AND INTERIOR OF THE ISLAND
The first test aims at verifying that the numerical simulation of sound-wave propagation outside the X-point gives results in good agreement with the analytic solutions of this problem up to the order of accuracy of the chosen numerical scheme. FENICIA is used with a version embedding an island magnetic equilibrium. With respect to the model presented in [8,9], Dirichlet boundary conditions are still used in the x direction, but periodic boundary conditions are now used in the y and z directions. The results of our investigation are presented into the influence of a static axi-symmetric magnetic island, including the X-and O-points, on sound-waves. We consider a box of size 2∆ = 0.2 (normalized to the radial position of the separatrix), an island of size 4 √ A with A = 10 −3 and a mode (m, n) resonant at ρ = ρ mn , i.e: m/n = 1/(g(ρ mn ) √ 2 A). This means that k = 0 at ρ mn . For this value of A, we show 1/(g(ρ mn ) √ 2 A) as a function of ρ for both the exterior and the interior of the island as shown in Fig. ??. For rational values of 1/(g(ρ) √ 2 A), there exists a resonant (m, n). One sees that the position around ρ = 1 is directly constrained by high m or high m/n values, especially if one wishes to get closer to the separatrix starting from the interior of the island. We initialize a perturbation with the form where φ 0 (ρ) is a Gaussian structure centered around ρ mn and having the following form with ∆ρ = ρ mn /m and ρ bd is the value of ρ at the boundaries x = {x min , x max } for y = 0. The small radial width of the envelope φ 0 , around the resonant surface ρ mn , ensures that k ≪ 1. Two cases are to be examined: the exterior of the island (ρ > 1), and the interior (ρ < 1). The only difference between the two cases comes from the two different expressions of g(ρ) given by Eq. (10) in the previous section. For the former, we consider a perturbation φ 0 (ρ) centered around ρ mn = 1.25 for (m, n) = (24, 1). The initial condition for a simulation of size (n x , n y , n z ) = (800 × 800 × 20) is shown in Fig. ??. Note that time t is normalized to the Bohm timescale (which is long by a factor 1/ρ * compared to the significant turbulence time) and the time step considered for the following simulations is ∆t = 10 −3 . In Fig ??, we show the solution at time t = 1 . For clarity, a zoomed view of the initial condition and the final solution are also shown in Figs. ??-??. The emphasis is now on showing that the exact slab solution given by (12) is recovered at the exterior of the island where ρ > 1. In Fig. ??, we plot the relative error dx. The volume integral is defined over half of the domain (right side of the island). The calculations give similar results for the left side of the island. From Fig. ??, we see that the numerical results converge quickly to the analytic solution. The same tests are performed at the interior of the island where we consider a perturbation φ 0 (ρ) centered around ρ mn = 0.58 for (m, n) = (45, 1). The initial condition for a simulation of size (800 × 800 × 20) is shown in Fig. ?? and the final solution at t = 1 is given in Fig ??. A zoomed view of the initial condition and the final solution are also shown in Figs. ??-??. The numerical solution is again compared to the analytic solution given by (12) for ρ < 1 . This is illustrated in Fig. ?? where we plot the relative error as a function of time. We observe that the numerical solution to the sound-wave problem converges to the one derived analytically as the grid size ∆x decreases. However, the relative error for the resolutions considered here is still high (∼ 30%). We interpret this to indicate that the poloidal resolution is not sufficient to be able to resolve the m = 45 mode at the interior of the island.

IV. TESTS ACROSS THE SEPARATRIX
At the separatrix (ρ = 1), the analytic solutions (12) are no longer valid. As stated in section II, we do not have an obvious analytic solution at the separatrix. Given this limitation, rather than considering an initial condition as a function of the variable η, which diverges at ρ = 1, we instead consider a perturbation crossing the separatrix of the form where N 0 (ρ) is a Gaussian structure given by where ∆ρ is a parameter set to ∆ρ = 0.5 here for the perturbation to cover both the left and right-hand sides of the separatrix. A series of simulations were then run with ∆ = 0.1, A = 10 −3 and (m, n) = (5, 1). Figures (??-1) show the evolution of the electrostatic potential as a function of time (normalized to Bohm time). The time step for this simulation is ∆t = 10 −3 and the box size is (200 × 200 × 20). From Fig. ?? one sees that modes across the separatrix gradually evolve and shear. Due to the periodicity of the simulation box, we also observe eddies leaving one side of the box and reentering on the other side of the box. 3D illustrations of the same simulation are shown in Fig. 1 allowing us to see the elongated nature of structures along the magnetic field lines. It is also evident from that figure that the direction parallel to the magnetic field does not coincide with the direction of symmetry (z ). Because our ultimate (e) t=4 (f) t=5 FIG. 1: 3D snapshots of potential fluctuations at different simulation times goal was to prove the validity of the FCI approach for X-point geometries, three main tests are to be considered hereafter. The first test consists in showing the conservation of an energy-like quadratic quantity of the sound-wave model (4). For this purpose, we perform a scan in n x and in n z , the number of points in the x and z directions respectively (in all simulations we considered n x = n y ). For a fixed number of points in the z direction, n z = 20, a scan over the following pairs is done in the x and y directions: (n x , n y ) = {(100, 100); (200, 200); (400, 400); (600, 600)} until t = 1. The energy-like quantity E ≡ (φ 2 + (1 + 1/τ ) −1 u 2 )/2 dxdydz is conserved by Eqs. (4). It is plotted in Fig. ?? as a function of time. Notice that t = 1 is a Bohm time which is long compared to turbulence time scales which typically range between 1/ρ ⋆ and ∼ 10/ρ ⋆ with ρ ⋆ being the ratio between the ion larmor radius and the system size. The conservation of E is guaranteed as the grid spacing ∆x decreases. The line plot on graph ?? shows that it is effectively the case since in the worst case where (n x , n y ) = (100, 100) the relative change is equal to 4.5 × 10 −3 and in the well-resolved case where (n x , n y ) = (600, 600) the relative change is almost zero. Accumulating an error of a few percent in a few Bohm times would thus be tolerable when simulating turbulence. For the special case where the box size is (n x , n y , n z ) = (200,200,20) shown in the 2D and 3D snapshots above, we performed the run up to t = 12. Though there are fine structures of the grid size scale appearing after t = 6, the plot in graph ?? shows a relatively good conservation of the energy. The relative variation of energy as a function of time reaches at most 5%. Similarly, the conservation of the energy is well verified when scanning the number of points in the z direction. A box of size (n x , n y , n z ) = (400, 400, n z ) is considered where n z = {20, 40, 60, 80, 100}. In graph ??, the energy variation converges towards zero as ∆z decreases. Furthermore, its relative change is in the worst case equal to 5 × 10 −4 as shown in Fig. ??.
The second target test intends for the demonstration of the convergence of the numerical solution when crossing the separatrix region. For this purpose, a series of simulations were performed with a box of size (n x , n y ) = (400, 400) and a scan over n z = {5, 10, 20, 40, 60, 80, 100}. We start by considering the first two solutions given by the simulations having n z = 5 and n z = 10. The idea is to calculate the difference between the them and repeat the process over each pair of the entire set of solutions φ given by all the simulations. This is called the moving difference. In Fig. ??, the moving difference of the solution φ for n z = {5, 10, 20, 40} at the first poloidal plane iz = 1 is plotted in 2D at time t = 1. At this stage, for the three pairs given by {φ(n z = 5), φ(n z = 10)}, {φ(n z = 10), φ(n z = 20)}, {φ(n z = 20), φ(n z = 40)}, we qualitatively obtain an error that is of the order of the solution. This means that convergence is not yet reached. In Fig. ??, the moving difference is calculated for pairs of solutions given by {φ(n z = 40), φ(n z = 60)}, {φ(n z = 60), φ(n z = 80)}, {φ(n z = 80), φ(n z = 100)}. We observe that beyond n z = 40 the difference between each pair of solutions is nearly zero which is a good indication that convergence is reached. In other words, for this simulation, one does not need to go beyond n z = 40 to study the same physics. This constitutes the main strength of the field-aligned FCI system of coordinates and validates its application to X-point configurations with an exponential decay rate shown in graph ??. The graph of Fig. ?? indeed shows a plot line connecting all the fixed averages. It is called the moving average. More specifically, the average is calculated here by dividing the sum of the difference between each pair of solutions by the total number of solutions. It is thus a quantitative mean that allows us to conclude that the numerical solution converges as the resolution in the z direction decreases. In fact, at n z = 40, the norm is approximately 10 −2 which is consistent with the fact that we used second order centered finite differences to compute the parallel gradient operator ∇ . With this result, we conclude that the numerical solution of the sound-wave propagation problem given an equilibrium with a magnetic island has converged with the expected convergence order. To finish, we do a last test showing the order of convergence of the numerical solution. In fact, for the same set of simulations, i.e: (n x , n y ) = (400, 400) and n z = {5, 10, 20, 40, 60, 80, 100}, the test consists of choosing a reference case supposed to be the closest possible case to the analytic solution of the sound-wave problem. With this hypothesis, it is legitimate to calculate the average of the difference between all the simulations and the reference simulation chosen to be that with n z = 100 (Root Mean Squared of the difference). The resulting graph ?? shows that as n z tends to 100, the error tends to 0 with an estimation of the order of convergence given by the loglog plot in Fig. ??. The convergence is indeed fast, of the order of a = 2.6, the corresponding slope of the loglog plot. This is once again in good agreement with the use of second order finite differences to compute the parallel gradient operator. To sum up, even with a very elongated island in the perpendicular plane, and a high safety factor, which make the requirements for resolution particularly tough, the latter results show the robustness of the FCI system. At this stage, we can finally conclude that the FCI system of coordinates is quantitatively validated for X-point geometries.

V. CONCLUSIONS
Many computer resources can be spared by employing coordinate systems which take advantage of the physical characteristics of magnetized plasmas, namely the strong anisotropy of spatial scales between the parallel and the transverse directions (with respect to the equilibrium guiding magnetic field B). A new three-dimensional Flux-Coordinate Independent system, referred to as FCI in [8,9], is shown to have a number of advantages over earlier approaches. It firstly permits more flexible coding by decoupling the magnetic geometry from the meshgrid. Secondly, it allows for coarse grids since it is field-aligned (one follows the field lines to calculate the parallel gradient operator). Thirdly, the FCI system copes with X-point magnetic geometries, a situation where previous approaches based on magnetic flux coordinates have a fundamental problem due to the singularity of the field-aligned metric. In the present work, the focus is to show that the FCI approach can be extended to deal with X-point configurations such as magnetic islands. For this purpose, a study of sound-wave propagation across the X-point, in the presence of a magnetic island, was carried out using the code FENICIA (cf. [8,9]). Analytic solutions to the sound-wave problem are constructed both inside and outside the island. The numerical results show the ability of the FCI approach to converge to the derived analytic solution . Then, a test case is considered with an initial perturbation that straddles the magnetic separatrix including the X-point. Since there is no analytic solution in this case, conservation properties were first verified, then three conclusive tests are performed to show the convergence of the numerical solution with respect to the number of points in the z direction (direction of symmetry) allowing one to consider a much thinner grid to study the same Physics. In the latter test, the magnetic island is taken to be very elongated in the perpendicular plane, with a very high safety factor, making the requirements for resolution particularly tough. Even in this very demanding situation, results show the robustness of the FCI approach. We finally conclude this paper by the following: the FCI approach allows not only for a coarser grid by exploiting the anisotropy of the system, but it also allows for flux-coordinate independent operations in the poloidal plane which opens the way to the implementation of complex magnetic geometries free of singularities. It deals without difficulty with X-point configurations and with O-points such as the magnetic axis. The flexible nature of FENICIA presented in [8,9] allowed us to demonstrate the application of this coordinate system to a magnetic island in a slab, thus the validity of its application to X-point configurations. The FCI approach can thus be implemented in existing modular codes having pitfalls in dealing with the X-point geometry.

Acknowledgments
To be added later Appendix A: We seek to define proper coordinates associated to the field lines in the magnetic island equilibrium. For that purpose, let's consider a magnetic field given by Eqs. (1) and (2), with The parallel gradient operators can then be written as: In order to construct analytic solutions to the above model in the presence of a magnetic island, the idea is to define a coordinate system: (x, y) −→ (ρ ′ , y ′ ) such thatρ ′ is a magnetic surface label, i.e: ∇ ρ ′ = 0. One would write where ρ ′ is a normalization of ψ so that ρ ′ = 1 refers to the separatrix. Note that in the neighborhood of y = 0 (O-point), cos y ∼ 1 − y 2 /2 −→ ρ ′ ≃ (x−1) 2

2A
− 1 + y 2 /2 which defines the equation of an ellipse. With this set of coordinates, the spatial derivatives in the x and y directions write: where x ′ = x − 1. This leads to a parallel operator of the form The parallel gradient operator can then be written as where x ′ (ρ ′ , y ′ ) = ± (2A) (ρ ′ + cos y ′ ) 1/2 . At this point, to ensure the separation of variables, one defines a new system: (ρ ′ , y ′ ) −→ (ρ, η) such that ρ = ρ ′ and We conclude from (A9) that so, the new variable η can be expressed as with g(ρ) chosen such that, for all x ′ : The function η is regular everywhere, but delicate near the X-point where ρ = 1 and y 0 dy ′ /(ρ ′ + cos y ′ ) 1/2 diverges. When ρ < 1, the argument of the integral involves the term (cos y ′ + ρ) 1/2 which becomes negative at cos y ′ = −ρ. We choose to limit the integration to the upper bound y M such that cos y M = −ρ. The situation is described through a schematic of the island in Fig. ??. The explicit expression of g(ρ) at the exterior of the island (i.e: ρ > 1) can be derived as follows: to get π 0 dy ′ (ρ + cos y ′ ) −1/2 = 2 (1 + ρ) 1/2 K 2 1 + ρ (A17) Note that K (2/(1 + ρ)) diverges logarithmically as ρ −→ 1 + . Finally, for ρ > 1, g(ρ) can be written as At the interior of the island (i.e: ρ < 1), however, the limit of integration is y M such that cos y M = −ρ. Thus Finally for ρ < 1 we have, For x ′ < 0, as shown in Fig. ??, η ′ needs to be completed. We define for x ′ > 0, but y < 0 if one wants 3π/2 ≤ η ′ ≤ 2π in the fourth quadrant.