A Flux-Coordinate Independent (FCI) approach for anisotropic systems, not based on magnetic flux coordinates, has been introduced in Hariri and 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 can efficiently deal with 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 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 Refs. 1–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 Ref. 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 developed in Ref. 8 and elaborated in Ref. 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 Refs. 8 and 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 and allows for the simulation of complex magnetic geometries. 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 Flux 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 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. First, 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 investigating turbulence in X-point geometries is crucial for predicting the L-H power threshold. Second, 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 Refs. 17–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. In fact, the expression of the parallel gradient operator, which the approach developed in Ref. 8 relies on, is valid everywhere including the X- and O-points. For this purpose, we use the code FENICIA presented in Refs. 8 and 9, where the Flux-Coordinate Independent field-aligned 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 Ref. 8. 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. A brief review of the FCI approach is given in Sec. II. Analytic solutions of a sound-wave propagation model for an island magnetic equilibrium are presented in Sec. III. Tests are performed at the exterior and interior of the island to show the convergence of the numerical solution to the analytic one in Sec. IV. Across the separatrix, where there are not yet known analytic solutions, two conclusive numerical convergence tests are presented in Sec. V, allowing us to show the robustness of the FCI approach in X-point magnetic configurations.

## II. REVIEW: THE FCI APPROACH

In order to fully describe a given field in a tokamak plasma, one needs high resolution in any given poloidal plane, but only a small number of these planes. Field-aligned coordinates allow for a sparse mesh along the field lines. They can be, however, constructed in a way that avoids the use of flux coordinates to discretize the fields in the poloidal plane. A new FCI approach, not based on flux variables to discretize the fields in the poloidal plane, has been developed in Ref. 8. The method exploits the elongated nature of microinstability driven turbulence whose scales perpendicular to the guiding magnetic field are typically of the order of the ion gyroradius, while parallel scales are of the order of the machine size. This paper extends the results of Ref. 8 to validate the FCI approach in an X-point configuration. The approach can be easily generalized to non-axisymmetric and stochastic magnetic field configurations. Details are beyond the scope of this paper and shall be discussed in an upcoming work. In the following, we consider a class of static low-*β* equilibria, such that the suitably normalized axisymmetric magnetic field can be written in the form

where we employ a three dimensional Cartesian reference system (*x*, *y*, *z*), such that $z\u0302$ 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 $\psi (x)$ such that

It follows that $bx=\u2202y\psi $ and $by=\u2212\u2202x\psi $. Magnetic surfaces can be labeled by the value of *ψ*. A point where *ψ* has a saddle point is an *X-point* and a point where *ψ* has a local extremum is an *O-point*. Both X- and O-points are such that $\u2207\psi =0$. Passing through the saddle are two separatrix legs, on which *ψ* is constant. Both closed and open field lines can be treated. The parallel derivative operator is given by

where the Poisson bracket operation is defined by $[\psi ,\xb7\u2009]=\u2202x\psi \u2202y\u2212\u2202y\psi \u2202x.$

The central idea is to work with a grid independent of magnetic flux surfaces, not using the magnetic coordinates, but computing the parallel derivatives directly along the magnetic field lines. One then looks for a change of coordinates from the original (*x*, *y*, *z*) to a new set $(\xi \alpha ,s)$, such that *s* can be treated as a slowly varying coordinate and only the two $\xi \alpha (\alpha =1,2)$ carry the information on the small scales. Taking advantage of what was learned in Ref. 6, we divide the domain in a certain number of sectors centered around *z _{k}*, and extending to the boundary in the (

*x*,

*y*) directions, with

*k*labeling a given sector. This leads to a set of transformations of the form

where *V ^{α}*(

**x**) and

*C*(

^{α}**x**) are yet unknown functions. In terms of the new variables the parallel derivative is given by

The FCI approach relies on expressing the parallel gradient operator as

thus eliminating the fast-varying derivatives. Here, *A* is a local function. One can show that it is always possible to construct such a coordinate system $(\xi \alpha ,s).$ The sufficient condition for the existence of the above transformation breaks at the extrema and saddle points of *ψ*, that is, at X- and O-points, respectively, where relation (6) is not valid. However, the key point for the FCI approach is that, precisely at the X- and O- points, the parallel gradient can still be expressed as a partial derivative along a single coordinate, namely *z*: $\u2207\u2225|X\u2212point\u223c\u2202z$. This, in turn, allows to construct the FCI approach everywhere, including the X and O-points, based on a field line map in the original Cartesian coordinates. In practice, considering a second-order finite difference representation, the expression of the discrete parallel gradient operator in (*x*, *y*, *z*) can be written as

So, a same point in space can be labeled by $(\xi 1,\xi 2,s+\Delta s)\u21d4(x+\Delta x,y+\Delta y,z+\Delta z)$, where (Δ*x*,Δ*y*) are obtained by integrating along the field lines

from (*x, y, z*) = (*x*_{0}, *y*_{0}, *z _{k}*) to $(x,y,z)=(x0+\Delta x,y0+\Delta y,zk+\Delta z)$ for a given displacement Δ

*z*= Δ

*s*.

More generally, the direct calculation of $\u2207\u2225$ at a given point is done by considering an arc of field-line passing through a point (*x*, *y*, *z*). One can show that there always exists a parametrization *x*(*τ*), *y*(*τ*), *z*(*τ*) to get the following system:

For any differentiable function *f*(*x, y, z*), the derivative along the field-line with respect to *τ* is computed as follows:

Right at the X-point, since *b ^{x}*

^{ }= 0 and

*b*

^{y}^{ }= 0, it trivially follows from Eq. (10) that $\u2207\u2225f$ reduces to $\u2202zf$. As a result, Eq. (10) is valid everywhere (provided that z can be used to parametrize unambiguously the field line, which is the case if

*B*does not change sign). Its finite difference form writes

_{z}where Δ*τ* = Δ*z*. Similarly, if toroidal coordinates are employed in the poloidal plane, $(R,\phi ,Z)$, the infinitesimal displacements become: $(dR,Rd\phi ,dZ)$. These must be proportional to the components of **B**. Thus,

where *τ* is a parameter for the position along the field line. It is convenient to fix *τ* such that $d\phi /d\tau =1$ (case for toroidal sectors). The following field line equations result

Then, for a function $f(R,\phi ,Z)$,

Again, right at the X-point, $B\xb7\u2207f=B\u2009\u2207\u2225f$ trivially reduces to $B\phi \u2202\phi f$. At a given point $\phi $, this leads to the following expression in finite difference form:

where $f(\phi \xb1\Delta \phi )$ corresponds to the value of $f(R,Z,\phi )$ at points $(R\xb1,Z\xb1,\phi \xb1\Delta \phi )$. Here, $R\xb1$ and $Z\xb1$ are obtained by integrating the field-line Eq. (13).

We conclude that the construction of the above scheme based on (6) is well-behaved everywhere, including the X-point. The construction of such a system allows decoupling of the geometric description in the poloidal plane from the treatment of the parallel derivatives. Any coordinate system can, therefore, be used in the poloidal plane. The parallel derivative is then computed in one go, by using values at the end points of arcs of magnetic field lines. This generally requires interpolating in the poloidal plane since, usually, end points are not grid nodes. Accuracy of the interpolation operation would constrain the resolution in the poloidal planes to be adequately high, but high resolution is needed anyway, in order to keep the necessary information on the fine structure and to carry out the operations in the poloidal plane to a satisfactory accuracy. Thus, the need to use interpolation for the parallel operations does not introduce substantial constraints as shown by using the Hermite cubic interpolation in Ref. 9. In light of these findings, a new modular code named FENICIA was developed.^{8,9} It uses the FCI approach and solves the general class of plasma models

where *S* is a structure of vectors representing the state of the system. *E*(*S*) is a nonlinear operator that can be treated explicitly without much penalty. *L* and *I* represent linear operators such that the reduced problem obtained by setting *E*(*S*) = 0 could be treated implicitly. The splitting of the right hand side (r.h.s.) between *E* and *I* is by no means unique and depends on the physics to be studied. As a general rule, only the physics occurring at the timescale of interest for the specific problem are treated explicitly. The main constraint with respect to a generic r.h.s. is that the implicit problem be linear. We further assume that *L* and *I* are time-independent. And boundary conditions are built into *L* beforehand.

## III. SOUND WAVE MODEL IN A MAGNETIC ISLAND GEOMETRY

As extensively discussed in Ref. 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–14} The nonlinear growth of these modes was also investigated in Ref. 15 and extended to cases where current nonlinearities are important as shown in Ref. 16.

The results presented in this paper are obtained using the code FENICIA. Consider a magnetic equilibrium characterized by a magnetic island whose half radial width is equal to $\delta x=2A$, with *A* being a parameter. Such an equilibrium corresponds to a magnetic field given by Eqs. (1) and (2), with

The considered domain, for an island centered at *x* = 1, is $xmin\u2264x\u2264xmax$ (in practice, *x _{max}* = 1 + Δ and

*x*= 1 − Δ, with 0 < Δ < 1) and $\u2212\pi \u2264y\u2264\pi $, with $0\u2264z\u22642\pi $.

_{min}Let us consider the following sound-wave model, belonging to Ref. 16, that propagates parallel to the magnetic field

$\varphi $ 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

*ρ*being the ion sound Larmor radius; and

_{s}*τ*= 1.

The aim of the present section is to construct analytic solutions of model (18) 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 $\u2207\u2225\psi =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 Appendix, 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. We 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 the right side (+) or the left side (−) of the exterior region is considered (cf. the Appendix). *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 Figs. 1(a) and 1(b), 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:

We consider wave-like solutions of model (18) 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. (18) leads to the following system:

The eigenfrequency solution of the dispersion relation depends on *ρ*

The eigenvectors are then

For a given expression of $\varphi 0(\rho )$, *u*_{0} can be calculated by using Eq. (28), with $\omega mn\xb1(\rho )$ given by Eq. (27). Analytic solutions of the model (18), provided by Eq. (25), will be compared to their numerical counterparts in Sec. IV. While the above solutions can be derived for *ρ* < 1 and *ρ* > 1, there is no obvious analytic solution to this problem at $\rho =1$. This issue is left for future investigation.

## IV. 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 to solve the system of Eq. (18) with a version embedding an island magnetic equilibrium. With respect to the model presented in Refs. 8, 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 $4A$ with *A* = 10^{−3} and a mode (*m*, *n*) resonant at *ρ* = *ρ _{mn}*, i.e., $m/n=1/(g(\rho mn)\u20092\u2009A)$. This means that $k\u2225=0$ at

*ρ*. For this value of

_{mn}*A*, we show $1/(g(\rho mn)\u20092\u2009A)$ as a function of

*ρ*for both the exterior and the interior of the island as plotted in Fig. 2. For rational values of $1/(g(\rho )\u20092\u2009A)$, there exists a resonant (

*m*,

*n*). Accordingly, the position around

*ρ*= 1 is directly constrained by high

*m*or high

*m*/

*n*values, especially when approaching the separatrix starting from the interior of the island. We initialize a perturbation of the form

where $\varphi 0(\rho )$ is a Gaussian structure centered around *ρ _{mn}* and having the following form:

with $\Delta \rho =\rho mn/m$ and *ρ _{bd}* is the value of

*ρ*at the boundaries $x={xmin,xmax}$ for

*y*= 0. The small radial width of the envelope $\varphi 0$, around the resonant surface

*ρ*, ensures that $k\u2225\u226a1$. Two cases are to be examined: the exterior of the island (

_{mn}*ρ*> 1), and the interior $(\rho <1)$. The only difference between the two cases comes from the two different expressions of $g(\rho )$ given by Eqs. (21) and (23). For the exterior region of the island, we consider a perturbation $\varphi 0(\rho )$ centered around

*ρ*= 1.25 for (

_{mn}*m,n*) = (24,1). The initial condition for a simulation of size $(nx,ny,nz)=(800\xd7800\xd720)$ is shown in Fig. 3(a). Note that time

*t*is normalized to the Bohm timescale (which is long by a factor $1/\rho *$ compared to the significant turbulence time) and the time step considered for the following simulations is $\Delta t=10\u22123$. In Fig. 3(b), 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. 3(c) and 3(d). The emphasis is now on showing that the exact slab solution given by (25) is recovered at the exterior of the island, where $\rho >1$. In Fig. 4, we plot the relative error

between the exact and the numerical solution (ϕ_{exact} and ϕ_{num}, respectively) as a function of time for increasing spatial resolutions $(nx,ny,nz)={(400,400,20);(600,600,20);(800,800,20)}$, where $\u27e8\xb7\u27e9=\u222b02\pi dy\u2009\u222b02\pi dz\u2009\u222b1xmaxdx$. The volume integral is defined over half of the domain (right side of the island). The calculations give similar results if considering the left side of the island. The numerical results show faster than linear convergence to the analytic solution, as seen in Fig. 4. The same tests are performed at the interior of the island where we consider a perturbation $\varphi 0(\rho )$ 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. 5(a), and the final solution at

*t*= 1 is given in Fig. 5(b). A zoomed view of the initial condition and the final solution is also shown in Figs. 5(c) and 5(d). The numerical solution is again compared to the analytic solution given by (25) for

*ρ*< 1. This is illustrated in Fig. 6, 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. The relative error for the resolutions considered here is yet 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.

## V. TESTS ACROSS THE SEPARATRIX

At the separatrix (*ρ* = 1), the analytic solutions (25) are no longer valid. As stated in Sec. III, 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 *rho* = 1, we instead consider a perturbation crossing the separatrix of the form

where $N0(\rho )$ 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 run with Δ = 0.1, *A* = 10^{−3} and (*m,n*) = (5,1). Figures 7 and 8 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). As it is apparent from Fig. 7, modes across the separatrix gradually evolve and shear. Eddies leaving one side of the box reenter on the other side due to the periodicity of the simulation domain. The elongated nature of structures along the magnetic field lines is identified through 3D illustrations of the same simulation, shown in Fig. 8. It becomes evident that the direction parallel to the magnetic field does not coincide with the direction of symmetry (*z*). Because our ultimate 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 (18). For this purpose, we perform a scan in *n _{x}* and in

*n*, the number of points in the

_{z}*x*and

*z*directions, respectively (in all simulations we considered

*n*=

_{x}*n*). For a fixed number of points in the

_{y}*z*direction,

*n*= 20, a scan over the following pairs is done in the

_{z}*x*and

*y*directions: $(nx,ny)={(100,100);(200,200);(400,400);(600,600)}$ until

*t*= 1. The energy-like quantity $E\u2261\u222b(\varphi 2+(1+1/\tau )\u22121\u2009u2)/2\u2009dxdydz$ is conserved by Eq. (18). It is plotted in Fig. 9(a) as a function of time. Notice that

*t*= 1 is a Bohm time that is long compared to turbulence time scales, which typically range between $1/\rho \u22c6$ and $\u223c10/\rho \u22c6$ with $\rho \u22c6$ being the ratio between the ion larmor radius and the system size. Figure 9(b) shows that the conservation of

*E*improves as the grid spacing Δ

*x*decreases. For a resolution (

*n*,

_{x}*n*) = (100,100), the relative change in Energy is equal to 4.5 × 10

_{y}^{−3}and in the well-resolved case, where (

*n*,

_{x}*n*)= (600,600), the relative change decreases to 0.5 × 10

_{y}^{−3}. When simulating turbulence, errors of the same order of magnitude are expected. For the special case where the box size is (

*n*,

_{x}*n*,

_{y}*n*) = (200,200,20) shown in the 2D and 3D snapshots above, we performed the run up to

_{z}*t*= 12. Though there are fine structures of the grid size scale appearing after

*t*= 6, the plot in Fig. 10 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 improves with increasing resolution in the

*z*direction. A box of size (

*n*,

_{x}*n*,

_{y}*n*) = (400, 400,

_{z}*n*) is considered where $nz={20,40,60,80,100}$. In Fig. 11(a), the energy variation converges towards zero as Δ

_{z}*z*decreases. Furthermore, its relative change is in the worst considered case equal to 5 × 10

^{−4}as shown in Fig. 11(b).

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*) = (400, 400) and a scan over

_{y}*n*= {5, 10, 20, 40, 60, 80, 100}. We start by considering the two solutions given by the simulations with

_{z}*n*= 5 and

_{z}*n*= 10. The idea is to calculate the difference between them and repeat the process over each pair of the entire set of solutions $\varphi $ given by all the simulations. This is called the

_{z}*moving difference*. In Fig. 12, the moving difference of the solution $\varphi $ for

*n*= {5, 10, 20, 40} at the first poloidal plane

_{z}*iz*= 1 is plotted in 2D at time

*t*= 1. At this stage, for the three pairs given by ${\varphi (nz=5),\varphi (nz=10)},{\varphi (nz=10),\varphi (nz=20)},{\varphi (nz=20),\varphi (nz=40)}$, we qualitatively obtain an error that is of the order of the solution. This means that convergence is not yet reached. In Fig. 13, the moving difference is calculated for pairs of solutions given by ${\varphi (nz=40),\varphi (nz=60)},{\varphi (nz=60),\varphi (nz=80)},{\varphi (nz=80),\varphi (nz=100)}.$ We observe that beyond

*n*= 40, the difference between each pair of solutions is nearly zero. This is an indication that convergence is reached, and a resolution of

_{z}*n*= 40 is sufficient to resolve the physics problem considered. This constitutes the main strength of the field-aligned FCI system of coordinates and validates its application to X-point configurations with an exponential convergence rate shown in Fig. 14. 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}*z*direction decreases. In fact, at

*n*= 40, the norm is approximately 10

_{z}^{−2}which is consistent with the fact that we used second order centered finite differences to compute the parallel gradient operator $\u2207\u2225$. 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*) = (400, 400) and

_{y}*n*= {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

_{z}*n*= 100 (Root Mean Squared of the difference). The resulting Fig. 15(b) shows that as

_{z}*n*tends to 100, the error tends to 0 with an estimation of the order of convergence given by the loglog plot in Fig. 15(b). The convergence is indeed fast, of the order of

_{z}*a*= 2.6, the corresponding slope of the loglog plot. This is once again in agreement with the use of second order finite differences to compute the parallel gradient operator. In summary, even with a highly elongated island in the perpendicular plane, which make the requirements for resolution particularly difficult, 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.

## VI. CONCLUSIONS

Employing coordinate systems which take advantage of the physical characteristics of magnetized plasma, namely the strong anisotropy of spatial scales between the parallel and the transverse directions (with respect to the equilibrium guiding magnetic field **B**) allow the simulation of more complex systems with today's computational resources. A new three-dimensional Flux-Coordinate Independent system, referred to as FCI in Refs. 8 and 9, is shown to have a number of advantages over earlier approaches. It first permits more flexible coding by decoupling the magnetic geometry from the computational grid. Second, it allows for coarse grids since it is field-aligned (the idea is to follow the field lines to calculate the parallel gradient operator). Third, 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, we show that the FCI approach can be extended to deal with X-point configurations such as magnetic islands. For verification of the approach, a study of sound-wave propagation across the X-point of a magnetic island was carried out using the code FENICIA (cf. Refs. 8 and 9). Analytic solutions to the sound-wave problem were constructed both inside and outside the island. The numerical results show that the FCI approach converges to the derived analytic solution. A test case is considered with an initial perturbation that straddles the magnetic separatrix including the X-point. In the absence of an analytic solution, conservation properties were verified and conclusive tests performed to show the convergence of the numerical solution with respect to the number of points in the *z* direction (direction of symmetry). This demonstrates the capability of the approach to study X-point physics with reduced numerical resolution in the parallel direction. If the magnetic island is highly elongated in the perpendicular plane, the requirements for the numerical resolution increase. Results show the robustness of the FCI approach also in this stressed test. The FCI approach allows for a coarser grid by exploiting the anisotropy of the system and for flux-coordinate independent operations in the poloidal plane. This opens the way to the implementation of complex magnetic geometries, without difficulties encountered with X- and O-points configurations. The flexible nature of FENICIA presented in Refs. 8 and 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.

## ACKNOWLEDGMENTS

F.H. acknowledges discussions with L. Villard, Y. Peysson, V. Naulin, and S. Brunner. F.H. would, in particular, like to express her appreciation to S. Cowley for his encouragement and valuable feedback on the FCI approach. Further thanks to X. Litaudon, E. Joffrin, S.-I. Itoh, and K. Itoh for advice and support benefiting this work.

This work, supported by the European Communities under the contract of Association between EURATOM and CEA, was carried out within the framework of the European Fusion Development Agreement. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Part of this work was also supported by the ANR Contract No. E2T2.

### APPENDIX: DERIVATION OF COORDINATES FOR AN ISLAND GEOMETRY

We seek to define proper coordinates associated to the field lines in the magnetic island equilibrium. For that purpose, let us consider a magnetic field given by Eqs. (1) and (2), with

The parallel gradient operators is expressed 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)\u2009\u2192\u2009(\rho \u2032,y\u2032)$ such that $\rho \u2032$ is a magnetic surface label, i.e., $\u2207\u2225\rho \u2032=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\u2009y\u223c1\u2212y2/2\u2009\u2192\u2009\rho \u2032\u2243(x\u22121)22A\u22121+y2/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\u2032(\rho \u2032,y\u2032)=\xb1(2A)\u2009(\rho \u2032+cos\u2009y\u2032)1/2$. At this point, to ensure the separation of variables, we define a new system: $(\rho \u2032,y\u2032)\u2009\u2192\u2009(\rho ,\eta )$ such that *ρ* = ρ′ and

We conclude from (A9) that

So, the new variable *η* can be expressed as

with *g*(*ρ*) chosen such that for all *x′*:

This gives

So

The function *η* is regular everywhere, but delicate near the X-point where *ρ* = 1 and $\u222b0ydy\u2032/(\rho \u2032+cos\u2009y\u2032)1/2$ diverges. When *ρ* < 1, the argument of the integral involves the term $(cos\u2009y\u2032+\rho )1/2$ which becomes negative at $cos\u2009y\u2032=\u2212\rho $. We choose to limit the integration to the upper bound *y _{M}* such that $cos\u2009yM=\u2212\rho $. The situation is described through a schematic of the island in Fig. 16. The explicit expression of

*g*(

*ρ*) at the exterior of the island (i.e.,

*ρ*> 1) can be derived as follows:

Use definitions of elliptic integrals of first kind

to get

Note that $K\u2009(2/(1+\rho ))$ diverges logarithmically as $\rho \u21921+$. Finally, for $\rho >1,\u2009g(\rho )$ can be written as

At the interior of the island (i.e., *ρ* < 1), however, the limit of integration is *y _{M}* such that $cos\u2009yM=\u2212\rho $. Thus

but now

In order to express *g*(*ρ*) as a function of *K*, the elliptic integral, we proceed by first rewriting the expression $I=\u222b0yM\u2009dy\u2032\u2009(cos\u2009y\u2032+\rho )\u22121/2$ by change of variables. We indeed have

where $0\u2264\alpha \u2264\pi /2$. Thus

Moreover, by differentiation, we get

Also,

yielding

and

Hence,

Finally, for *ρ* < 1, we have

For *x′*^{ }< 0, as shown in Fig. 16, *η′* needs to be completed. We define

so that $\pi /2\u2264\eta \u2032\u22643\pi /2$ when *x′* < 0. And

for *x′*^{ }> 0, but *y* < 0 if one wants $3\pi /2\u2264\eta \u2032\u22642\pi $ in the fourth quadrant.