Global fluid simulation of plasma turbulence in stellarators with the GBS code

The implementation of three-dimensional magnetic fields, such as the ones of stellarators, in the GBS code (Ricci et al 2012 Plasma Phys. Control. Fusion 54 124047; Giacomin et al 2022 J. Comput. Phys. 464 111294) is presented, and simulation results are discussed. The geometrical operators appearing in the drif-reduced Braginskii equations evolved by GBS are expanded considering the typical parameter ordering of stellarator configurations. It turns out that most of the operators have a similar structure as the one implemented in the tokamak axisymmetric version of the code. In particular, the perpendicular laplacian only acts on the poloidal plane, which avoids the need of a three-dimensional solver for the electrostatic potential. The simulation of an island divertor stellarator is then presented, showing the derivation of the magnetic equilibrium in detail and extending the results in (Coelho et al 2022 Nucl. Fusion 62 074004). Although the island magnetic field-lines divert the plasma towards the strike points of the walls, the islands do not seem to have a significant impact on the turbulence properties. The dominant mode, identified as interchange-driven, is field-aligned and breaks the stellarator toroidal symmetry. The radial and poloidal extensions of the mode are of the same order, in contrast to typical tokamak simulations. This has consequences on the poloidal dependence of turbulent transport.


Introduction
As stellarators assume an increasing interest as a viable alternative to the tokamak [1,2], the importance of improving the predictive capabilities of plasma turbulence simulations in three-dimensional magnetic fields grows.Significant effort radial and poloidal directions are similar.Although such level of coherency and size of the structures is not usual in tokamaks, in stellarators such as TJ-K, a large coherent mode with low poloidal mode number is observed (usually m = 3 or m = 4 in hydrogen and helium discharges) [11], a feature that GBS simulations of TJ-K were able to retrieve [12].
In our initial work reporting a simulation of a stellarator with an island divertor [9], the drive of the dominant mode is identified as a curvature-driven mode with a ballooning nature.The investigation of the nature of the mode is achieved by means of a simplified linear theory.The surprising result that the non-linear turbulent flux is larger on the high-field side (HFS) than on the low-field side (LFS) is not explained.
In this paper we extend the results presented in [9].We present a detailed description of the model used to perform stellarator simulations and its implementation in GBS.We analyse in more detail the simulation results and, in particular, the nature of the mode dominating the transport, and develop a more complete linear theory that takes into account the stellarator geometry, which is able to explain the peaking of the turbulent transport at the HFS.
The present paper is organized as follows.Section 2 describes the extension of the GBS code to 3D magnetic fields.In section 3, we present the Dommaschk potentials which are used to generate stellarator vacuum fields.In section 4, results of an island divertor stellarator simulation is presented.Finally, our conclusions are drawn in section 5.

Extension of the GBS code to 3D configurations
GBS [13][14][15] is a three-dimensional, global, two-fluid, flux-driven code that solves the drift-reduced Braginskii equations [16], valid in the high-collisionality regime that often characterizes the plasma boundary of magnetic fusion devices, as well as the core of low-temperature devices such as the TORPEX basic plasma physics experiment [17] or the TJ-K stellarator [12].GBS evolves all quantities in time, without separation between equilibrium and fluctuating components.In the version of GBS with non-axisymmetric magnetic fields, the electrostatic limit is taken, the Boussinesq approximation is employed [13] and the gyroviscous terms, as well as the coupling to the neutral dynamics, are neglected, although these are implemented in the most recent version of the GBS code for tokamak simulations [14].Within these approximations, the drift-reduced model evolved by GBS is: (2) (3) which are closed by In equations ( 1)- (7), and in the rest of the paper, density n, electron temperature T e and ion temperature T i are normalized to the reference values n 0 , T e0 and T i0 , respectively; electron parallel velocity V ∥e and ion parallel velocity V ∥i are both normalized to the sound speed c s0 = T e0 /m i ; vorticity ω and the electrostatic potential Φ are normalized to T e0 /(eρ 2 s0 ) and T e0 /e; time is normalized to R 0 /c s0 , where R 0 is the machine major radius; perpendicular and parallel lengths are normalized to the ion sound Larmor radius, ρ s0 = √ T e0 m i /(eB 0 ), and R 0 , respectively.The normalized parallel current is j ∥ = n(V ∥i − V ∥e ) and the magnetic field B is normalized to the magnitude of the field on axis, B 0 .
The dimensionless parameters appearing in equations ( 1)-( 7) are the normalized ion sound Larmor radius ρ * = ρ s0 /R 0 , the normalized electron and ion parallel heat diffusivities, χ ∥e and χ ∥i , considered constants here, the ion to electron temperature ratio τ = T i0 /T e0 , the normalized electron and ion viscosities, η 0e and η 0i , which we assume to have constant values, and the normalized Spitzer resistivity, ν = ν 0 T −3/2 e , with where λ denotes the Coulomb logarithm [18].Small numerical diffusion terms such as ∥ n (and similar for the other fields) are introduced to improve the numerical stability of the simulations (the simulation results show that they Boundary conditions applied at the top, bottom, inner and outer domain boundaries.The derivative ∂s = s • ∇ is along the direction normal to the surface, s.The upper signs apply if the magnetic field is directed towards the wall, while the lower signs apply in the opposite case.It is defined F T = 1 + τ T i /Te and Λ = 3 [21].

Top and bottom walls
Inner and outer walls have a negligible effect on turbulence since they lead to significantly lower perpendicular transport than turbulence).The terms S n , S Te and S T i denote the density, electron temperature and ion temperature sources, respectively.Magnetic presheath boundary conditions, described in [19,20], are applied to all quantities at the top and bottom boundaries of the simulation domain.For the stellarator simulations presented here, the density and vorticity at the boundary satisfy, respectively, ∂ s n = 0 and ω = 0, where ∂ s = s • ∇ is the derivative along the direction normal to the wall, s.This choice is due to stability reasons.A summary of the applied boundary conditions is given in table 1.
The GBS simulation domain is a torus of radius R 0 with a rectangular cross-section of size L R × L Z .The physical model in equations ( 1)-( 7) is discretized using a regular cylindrical grid (R, ϕ, Z), with R the radial coordinate, ϕ the toroidal angle and Z the vertical coordinate.Equations ( 1)-( 6) are advanced in time using an explicit Runge-Kutta fourth-order scheme, while spatial derivatives are computed with a fourth-order finite difference scheme.
The normalized geometrical operators appearing in equations ( 1)-( 7) are the parallel gradient For their numerical implementation in three-dimensional magnetic field configurations, we expand these operators in the following small parameters: where B R and B Z are the radial and vertical components of the magnetic field; the normalized mirror ratio, ∆ = (B max − B min )/B, where B is the toroidally averaged value of B along the magnetic axis and (B max − B min ) is the toroidal ripple amplitude; and the ratio between perpendicular and parallel turbulence length scales, σ = l ⊥ /l ∥ .We then retain only the leading order terms in these expansion parameters as shown in appendix A. For the stellarator configuration considered in this work δ ∼ 0.1, ∆ ≲ 0.1 and a posteriori we verify that σ ∼ 0.01, confirming the validity of our expansion.In the particular case of a vacuum field, the operators are, after expansion, Since the geometrical operators ( 9)-( 13) have the same structure as the ones implemented in the axisymmetric version of GBS [14], their implementation does not require major changes in the structure of the code thanks to the use of a non field-aligned spatial discretization algorithm [22,23].In fact, as in the previous versions of GBS, an Arakawa scheme is applied to implement the Poisson brackets [23], and the same solver used to invert equation (7) can be considered, since the perpendicular Laplacian only acts on the RZ-plane (to lowest order in the expansion parameters δ, ∆ and σ).

Stellarator vacuum magnetic field based on the Dommaschk potentials
The stellarator vacuum magnetic field generated by any set of external coils can be described using appropriate basis functions that completely span the space of field solutions.A vacuum field obeys ∇ × B = 0, meaning that B = ∇V, where V is denoted as the magnetic potential.Since ∇ • B = 0, V satisfies a Laplace equation, ∇ 2 V = 0.There are several sets of analytical functions that solve Laplace's equation in a torus, namely, toroidal field harmonics, spherical harmonics or the Vuillemin-Gourdon set [24,25].A set of functions that has advantages over the other sets is the one commonly known as Dommaschk potentials [26][27][28], which are given in cylindrical coordinates by where a m,l , b m,l , c m,l and d m,l are free real constants and D m,l , N m,l are explicit functions given in [26].The first term on the right-hand side of equation ( 14) gives the 1/R dependence of the toroidal component of the magnetic field.The Dommaschk potentials are advantageous over others sets of solutions because the V m,l functions consist of simple algebraic expressions (powers, logarithms, sines and cosines) whose numerical evaluation is computationally convenient.Furthermore, the indices m and l correspond, respectively, to the toroidal and poloidal periodicities of V m,l , a property that for example the spherical harmonics lack [26].As a consequence, if the field is generated with potentials restricted to m = kN fp and k ∈ N, then the stellarator has field period N fp .In addition, stellarator symmetry [29] requires a m,l = d m,l = 0 for even l, and b m,l = c m,l = 0 for odd l.
We have constructed a five-field period (N fp = 5) stellarator-symmetric configuration with a 5/9 chain of islands surrounding a closed flux surface region by taking b 5,2 = c 5,2 = 1.5, b 5,4 = 10, a 5,9 = −d 5,9 = −7.5 × 10 9 , and imposing all other coefficients to vanish.Figure 1 shows the corresponding Poincaré plot of the magnetic field at different toroidal angles.The potential V 5,9 is added in order to control the width of the 5/9 islands.In fact, this potential resonates with the 5/9 rational surface and can be used to increase the size of the islands.If V 5,9 becomes too large, a chaotic region is generated due to the overlap of different island chains.This is shown in figure 2 where the Poincaré sections with a 5,9 = −d 5,9 = −15 × 10 9 reveal a chaotic region outside the last closed flux surface (LCFS) instead of an island chain.
The rotational transform ι is entirely provided by the rotation of the flux surfaces since no torsion is present in the configuration considered here (the magnetic axis lies on a plane) nor current.We can estimate the expected rotational transform on axis, ι 0 , from the elongation of the ellipses [30,31], that is ι 0 = (N fp /2)(r max − r min ) 2 /(r 2 min + r 2 max ) ≈ 0.5, where r max and r min are the ellipses major and minor radius, which is in very good agreement with the ι evaluated numerically by field-line tracing and depicted in figure 3. The ι profile flattens at ι = 5/9 ≈ 0.555, where the island chain is present.We note that the magnetic shear is very small.The amplitude of B varies approximately as 1/R, as in a tokamak.Therefore, a LFS and a HFS can be identified, similarly to tokamaks.This can be seen in figure 4 where we show |B| on the LCFS.In fact, |B| peaks at θ = π and is minimized around θ = 0, where θ is the polar angle defined with respect to the magnetic axis that rotates anti-clockwise, with the outboard midplane corresponding to θ = 0 and the inboard midplane to θ = π.
The height and width of the rectangular computational domain cross-section can be taylored such that the magnetic field of the islands intercept the top and bottom of the simulation box, as shown in figure 5, emulating an island divertor stellarator.In such a configuration, the LCFS is not in contact with the wall.Heat and particles outflowing from the core reach the island region and are transported along the magnetic field of the islands, eventually striking the top and bottom walls at specific toroidal locations.The strike points of the island magnetic field lines at the top of the simulation box are indicated in figure 6(a).A similar structure is found on the bottom part of the box.On the other hand, the Dommaschk configuration without islands, shown in figure 2, yields a limited stellarator.In such configuration, the LCFS is in contact with the wall.
It is noteworthy to mention the similarities between our island divertor stellarator model and W7-X.Our configuration leads to large connection lengths in the island region (L c ∼ 100R 0 ) as seen in figure 5(right), similarly to W7-X where connection lengths of the order of hundreds of meters are found in the edge [32].In addition, in W7-X the magnetic shear is also small.Finally, a chain of islands surrounds the closed flux surface region in W7-X and the plasma is evacuated through the magnetic islands that strike nonaxisymmetric divertor targets, although in our case the divertor targets are axisymmetric.

Simulation parameters and convergence
We use the following parameters for the simulations discussed in this paper: points and a time-step of 2.9 × 10 −6 R 0 /c s0 .The convergence of the simulation is assessed by performing a simulation with grid N R × N Z × N ϕ = 250 × 150 × 200 and a second one with grid N R × N Z × N ϕ = 250 × 150 × 250, which is at the limit of our computational capabilities.In both cases the dynamics of the simulation is similar to the one presented here.The sources for density and temperature, S n = S Te = S T i , are localized around a magnetic surface near the LCFS, as shown in figure 7, mimicking the ionization of recycled neutrals.We have performed simulations with forward and reversed magnetic field by reversing the total B, i.e. the sign of equation (14).We denote these two cases with B ϕ > 0 and B ϕ < 0, respectively.In what follows, if not mentioned, we are referring to the simulation with forward magnetic field (B ϕ > 0), which corresponds to an anticlockwise toroidal magnetic field seen from above.

Quasi-steady state and balance between particle sources and fluxes
The simulations are started from a noisy initial state and, after a transient, reach a quasi-steady state where sources, parallel and perpendicular transport and, ultimately, losses at the vessel balance each other.If we volume-integrate equation (1) we obtain   The second term in equation ( 16) corresponds to the E × B flux, the third term to the electron diamagnetic flux, and the fourth term to the electron parallel flux.These terms are integrated over the plasma volume inside the LCFS, and the time evolution of the integral is shown in figure 8.The plot shows that the simulation reaches a quasi-steady state at t ≈  The sum of all terms on the left-hand side of ( 16) balances the source term within an error of 14%, reflecting the numerical error of the volume integration and the approximation used in evaluating the time derivative ∂ t ( ´ndV).In fact, while time derivatives are discretized with a fourthorder Runge-Kutta scheme in GBS, we have used a secondorder finite-differences scheme when computing the term ∂ t ( ´ndV).This difference, along with the low sampling frequency of the data used to evaluate the quantities in figure 8, leads to errors of the order of 10%.

Equilibrium profiles
In figure 9 we show the equilibrium (time-averaged over the quasi-steady state) density and electron temperature profiles, which are very similar.Density and temperature peak in the center and decay radially across the flux surfaces, confirming that the prescribed configuration confines particles and heat within the LCFS.We note that the gradients are steeper at the HFS than at the LFS, which has consequences on the position of the peaks of the turbulent flux, as discussed in section 4.5.As expected, most of the particles and heat is deposited where the island magnetic field lines strike the top and bottom of the simulation domain.This is shown in figure 6(b), where the equilibrium electron pressure, ⟨p e ⟩ t = ⟨nT e ⟩ t , on the top of the simulation box is presented.The pressure peaks at the footprints of the island field-lines on the vessel.
In figure 10 we show the equilibrium profiles of the electrostatic potential and parallel current.j ∥ has a sinusoidal-like profile along the poloidal direction in the island region.In fact, ⟨∇ • j⟩ t = 0 implies ∇ ∥ j ∥ ∼ |∇ ⊥ j dia |, and since j dia = ∇p × B/B 2 , the Pfirsch-Schlütter current modulation j ∥ ∼ cos(θ) is expected from the 1/R nature of the magnetic field.We note that the potential profile in the core region is quite flat leading to an almost-zero radial electric field, E r .In fact, in the core, where fluctuations have small amplitude, ⟨E r ⟩ t ≡ −∂ r ⟨Φ⟩ t ∼ ∂ r ⟨p i ⟩ t / ⟨n⟩ t (see derivation in [33]), as confirmed by figure 11 where this relation along R, at ϕ = 0 and Z = 0, is shown to be satisfied in the core (950 ρ s0 ≲ R ≲ 1050 ρ s0 ).Although not clear in the plot due to saturation of the colorbar scale, we note that strong parallel currents form at the magnetic presheath entrance in the region where the islands strike the simulation boundary, at the top and bottom walls.This is a result of the potential deviating substantially from ΛT e at the top and bottom walls.We also note that a potentially-nonphysical negative potential is seen in the top boundary.

Nature of turbulent fluctuations
Snapshots of the density and potential fluctuations, evaluated as n = n − ⟨n⟩ t and Φ = Φ − ⟨Φ⟩ t , are shown in figure 12.A mode with toroidal mode number n = 2 and poloidal mode number m = 4 (corresponding to k y ρ s0 ≈ 0.04, where y is the binormal direction, while x and z are the radial and parallel coordinates) dominates the global dynamics of the system.The mode is coherent and rotates in the ion diamagnetic direction.We remark that the radial and poloidal wavenumbers of the mode are of the same order, k x ∼ k y .This is in contrast with ballooning-driven tokamak edge turbulence, where k x ∼ k y /L p with L p the equilibrium pressure length scale [34], meaning that k x < k y .
In figure 13 we display a snapshot of the potential fluctuations on the LCFS and its Fourier spectrum in straight-fieldline coordinates (θ * ,ϕ), with θ * defined such that B•∇θ * B•∇ϕ = ι is  constant on a flux surface, showing that the dominant mode has m = 4 and n = 2 (Coelho et al [9] wrongly reports that the toroidal mode number is 5).Alongside the dominant mode, two sub-dominant modes corresponding to the side-bands m = 3 and m = 5 are present.We remark that the dominant and subdominant modes are global, i.e. they do not satisfy the 5-field periodicity of the stellarator.Such symmetry-breaking modes are reported to exist in the TJ-K stellarator experiment [11], as confirmed by GBS simulations of this machine.This underlines the importance of simulating the whole torus and not just one field-period.We note that the dominant mode is field-aligned.Indeed, its poloidal and toroidal mode numbers are linked to the rotational transform through ι = n/m.In fact, in a straightfield-line coordinate system (s, θ * , ϕ), where s is a flux surface coordinate, the parallel gradient can be written as |B| (ιm − n) in Fourier space, and therefore k ∥ → 0 implies ι = n/m.In the current configuration, ι ≈ 0.5 throughout the whole volume, yielding k ∥ = 0 for n/m = 2/4.In our previous Letter [9], we show, by means of a non-local linear theory, that the mode 2/4 is one of the few linear modes that can transport the E × B flux observed in the non-linear simulation.On the other hand, a simulation we performed with only one field-period shows a dominant mode with n/m = 5/9, a consequence of limiting the toroidal mode number to multiples of 5. A simulation with a configuration where the islands are removed (see figure 2) shows that the same dominant mode, n/m = 2/4, appears (see figure 14), leading us to conclude that the islands are not an important element in determining the properties of the dominant mode.Furthermore, the pressure profile is not flat inside the islands.Indeed, such flattening could only occur if the time-scale for the convection along the islands flux tubes, τ ∥ , is shorter than the rotation timescale of the mode, τ mode .Since the mode rotates with approximately the ion diamagnetic frequency, s T e /(Ω i a) is the ion diamagnetic drift velocity, Ω i the ion cyclotron frequency and a the plasma minor radius, and denoting L c as the connection length of the island field lines, we have that since L c /a ≳ 100R 0 /a and k y ρ s = 0.04.This shows that convection along the island magnetic field lines occurs on a longer time-scale than the rotation of the mode, thus preventing the islands to all have the same flat pressure profile along them.This can be seen in the plots of the equilibrium density and temperature (figure 9).Based on a reduced non-local linear theory, Coelho et al [9] points out the ballooning nature of the dominant mode [9].In addition to this study, the nature of the mode can also be determined by comparing the amplitude of density fluctuations with those of the potential.Indeed, for resistive ballooning instabilities, ∇ ∥ Φ ∼ ν j∥ .From equation ( 4) this implies Φ ≫ ñ.In our simulation, the amplitude of the potential fluctuations is about twice as large as the amplitude of density fluctuations, as can be seen in figure 12.While this is not a conclusive observation that precludes the presence of other instabilities (drift-waves, for example), it is compatible with the previous statement on the ballooning nature of the dominant mode.

The turbulent E × B flux
The time-averaged E × B flux is characterized by the stellarator toroidal periodicity (n = 5).In fact, the instantaneous radial E × B flux (assuming B purely toroidal) is given by in the coordinate system (s, θ, ϕ), where s is a flux surface label, and ξ(s, θ, ϕ) is a geometrical factor such that in a stellarator has the same toroidal periodicity of the magnetic field and therefore can be written as ξ(s, θ, ϕ) = MN |ξ s MN | e i(Mθ−NN fp ϕ) (in the case of a tokamak with circular flux surfaces it reduces to ξ = 1/r).Since the density and potential fluctuations can be written as n = mn |n m n | e i(mθ−nϕ+ωt) and Φ = , with ∆ the phase difference between density and potential and ω the frequency of the mode, the time-averaged flux is The spectrum in figure 13 shows that all the dominant modes have n = 2 and thus equation ( 19) simplifies to (20) showing that the toroidal periodicity of the time-averaged flux is dictated by the toroidal periodicity of the geometrical factor (which is a multiple of the field period).On the other hand, the poloidal dependence is an intricate convolution between the poloidal dependencies of the geometrical factor and the density and potential fluctuations.As seen in figure 15, the flux peaks on the HFS rather than on the LFS, in contrast to typical observations in tokamaks.In the case B ϕ > 0 the peaks are localized around θ = 3π/2, while in the B ϕ < 0 case around θ = π/2.
We argue that the shift of the peak of the flux along the poloidal angle is a result of the poloidal dependence of the linear instabilities, as well as their convection due to an equilibrium E × B drift.We consider a reduced model of equations ( 1)-( 7) that includes the main elements of the ballooning instability [35]: We linearize the model by assuming that the density can be written as n = n 0 (x, z) + n 1 (x, z)e i kyy+γt , where x is a local coordinate normal to the flux surface, z a coordinate along the magnetic field line and y the local binormal coordinate; n 0 is the background and n 1 the perturbation (n 1 ≪ n 0 ).A similar expression is used for T e , while for the potential and the electron parallel velocity we assume Φ 0 = 0 and V ∥e0 = 0. We retain the z dependence of the equilibrium density and electron temperature because the gradients depend on the poloidal angle (they are typically steeper on the HFS).In addition, we retain the x dependence of n 1 (and others), since we observe in the non-linear simulation ∂ x n 1 ∼ k y n 1 .This contrasts with previous local linear theories of the drift-reduced Braginskii equations, where variations along the x direction are neglected, i.e. k x ≪ k y is assumed [35][36][37].The present linearization leads to the following set of equations whose derivation is detailed in appendix B: where , and κ n and κ g are the normalized normal and geodesic field-line curvatures, i.e. the projection of the curvature vector κ = −b × [∇ × b] along the normal and binormal directions (normalized to R 0 ).Note that, when the field is reversed, κ g changes sign while κ n is not affected.By Fourier transforming the parallel and radial directions, ∂/∂z → ik ∥ and ∂/∂x → i k x , a quadratic equation for the growth-rate is obtained: While the k ∥ term is stabilizing, the drive of the ballooning instability appears on the right-hand side of equation ( 23) and it results from the density and temperature equilibrium gradients and the normal and geodesic components of the curvature.Note that sign(B ϕ ) cancels out since κ g changes sign when the field is reversed.This means that the dispersion relation in equation (23) depends on the direction of the magnetic field only through the equilibrium profiles (which, in the non-linear simulation, are slightly different due to the different frequency with which data is saved).
In an infinite aspect-ratio tokamak where k x ≪ k y , one deduces from equation ( 23) that γ 2 ∼ cos(θ), consistent with  the common observation that turbulence is more developed at the LFS (θ ≈ 0) than at the HFS (θ ≈ π).On the other hand, the geodesic curvature is important when k x ∼ k y , as it is observed in our non-linear simulations, thereby influencing the position of the eigenmodes maximum amplitude.
The solution of the 2D eigenvalue problem in equation ( 22), namely the product of the eigenmodes |n 1 Φ 1 | at the radial position where the density gradient has its maximum, is shown in figure 16 for k y ρ s0 = 0.04, considering B ϕ > 0 and B ϕ < 0. The eigenmodes peak at θ = π/2 and θ = 3π/2, respectively, where the normal curvature is negative (but not where it has its minimum value) and the geodesic curvature is positive, as seen in figure 17 for the case B ϕ > 0 (the same applies for B ϕ < 0).
The peak position of the eigenmodes differ from the peak positions of the E × B flux observed in the non-linear simulation.In fact, the turbulent fluctuations are transported along the poloidal direction by the equilibrium E × B drift, ⟨V E×B ⟩ t , which is predominantly in the poloidal direction as shown in figure 18.The distance travelled by a turbulent eddy within the instability time-scale, V θ E×B /γ, is of the order of π l/2, with l the arc length of the LCFS in a single toroidal cross-section, in both the forward and reversed cases.
We draw the attention to the fact that the relevance of the geodesic curvature is pointed out to explain observations B.8.1)in the TJ-K stellarator.Experiments in this machine systematically observe the peak of the level of the density and potential fluctuations where κ n is negative and κ g is positive [38,39].In the experiment, however, this region also coincides with the peaking of the E × B flux, whereas in our simulation the effect of the poloidal ⟨V E×B ⟩ t convection is important.

Conclusions
The present paper describes the non-axisymmetric version of the GBS code, recently extended to handle three-dimensional magnetic configurations [9].Using the expansion parameters δ, σ and ∆, the geometrical operators appearing in the driftreduced Braginskii equations solved by GBS are expanded up to leading order, resulting in a perpendicular Laplacian that acts on the poloidal RZ plane and Poisson brackets that can still be numerically implemented using the Arakawa algorithm.This makes GBS implementation of the three-dimensional configurations similar to the axisymmetric version of the code.
We use the Dommaschk potentials to generate two vacuum magnetic fields, namely a configuration where a chain of islands surrounds the closed flux region mimicking an island divertor stellarator, and a configuration where the islands are replaced by a chaotic region, yielding a limited stellarator configuration.With the parameters of our simulations, we observe that the presence of the islands has a minor impact on the turbulent dynamics of the plasma.
The identification of the nature of the dominant mode as a ballooning-driven instability, based on a simplified nonlocal linear theory [9], is supported here by noticing that the amplitude of the potential fluctuations is larger than that of the density.The mode is field-aligned, k ∥ → 0, and its toroidal and poloidal mode numbers obey ι = n/m since ι = 0.5 and (n, m) = (2, 4).We note that the dominant mode does not have the toroidal periodicity of the magnetic field.Symmetrybreaking modes are reported in the TJ-K stellarator experiment [11].This observation in experiments, and now in simulations, underscores the importance of simulating the whole torus and not just one field-period.
When the magnetic field is reversed, B → −B, the E × B flux peaks around θ = π/2, in contrast to the B ϕ > 0 case where it peaks around θ = 3π/2.A two-dimensional linear theory suggests that the instabilities develop at θ = π/2 when B ϕ > 0, and at θ = 3π/2 when B ϕ < 0 because the equilibrium gradients are poloidally asymmetric and tend to peak at the HFS and the radial and poloidal elongations of the mode are similar, k x ∼ k y .In these conditions, not only the normal curvature plays an important role, but also the geodesic one, effectively shifting the most destabilizing region away from θ = 0.The mode is then transported in the poloidal direction by the equilibrium E × B drift, providing a possible explanation for the location of the peak of the turbulent flux.
One of the surprising outcomes of the simulations presented here is that the mode driving transport occurs on a large scale, it is coherent and has k x ∼ k y .This is not the case in typical tokamak simulations and it is a subject currently under investigation.Preliminary results indicate that this difference is due to the small value of the magnetic shear and the ellipticity of the flux surfaces.
A natural next step in the study of island divertor configurations is the simulation of W7-AS and W7-X configurations.Although W7-AS is no longer operating, there is a large body of literature reporting on experimental measurements in different configurations [40], which can help with further validation of the GBS code.Moreover, the simulation of W7-AS constitutes a step forward towards the simulation of W7-X, a machine that shares many features with its predecessor but whose size is still prohibitively large for simulation purposes at the moment, in particular because of the time necessary to solve the Poisson equation for the electrostatic potential, the bottleneck of GBS [14].Furthermore, we believe that the role of the neutrals cannot be neglected if one is interested in simulating meaningful scenarios in such large machines.This is the reason why a kinetic neutrals model, currently implemented in the tokamak version of the code [41], is being ported to the stellarator version of GBS. for x and y small (in the general case, e 2 and e 3 also depend on x and y).The covariant basis is formulated as: Note that g(x = 0, y = 0) = g −1 (x = 0, y = 0) are equal to the 3 × 3 identity matrix.
We can now express the parallel gradient and thus As usual, we normalize the parallel gradient by R 0 : Since we are interested in discretizing the toroidal direction rather than the z direction, we make use of the following property: allowing us to write the normalized operator as Regarding the Poisson Brackets,  and we renormalize x/ρ s0 → x and y/ρ s0 → y.
We now consider the curvature operator, The operator is normalized by R 0 ρ s0 : where x/ρ s0 → x, y/ρ s0 → y, R 0 κ n → κ n , and R 0 κ g → κ g .Note that κ n does not change sign when the field is reversed, but κ g does.Finally, the perpendicular Laplacian can be expressed as The perpendicular laplacian is evaluated according to We obtain:  The normalization is meaning that the last two terms are R 0 /ρ s0 times smaller than the first two and therefore we neglect them.To derive a reduced set of equations able to describe the ballooning instability, we start from the model in equation (21).We assume that quantities vary as f = f 0 (x, z) + f 1 (x, z)e ikyy+γt .For the potential and the electron parallel velocity we assume f 0 = 0.The z-dependence of the equilibrium density and electron temperature is retained because the gradients depend on the poloidal angle (they are typically steeper on the HFS).These assumptions lead to system in equation (22).
To numerically solve the system in equation ( 22), we recast it in matrix form, γ Âx = Bx.We define a two-dimensional regular grid, x i = x 0 + i∆x and ϕ = j∆ϕ with i = 0, . . ., N x − 1 and j = 0, . . ., N ϕ − 1.The operators ∂ z , ∂ x and ∂ 2 x are implemented with finite differences.We neglect the magnetic shear and assume infinite aspect-ratio, which allows us to set the grid spacing in the flux surface direction, ∆x, to be the same at all z positions.This spacing is defined as the distance between the flux surfaces at the initial position where the field lines start to be followed (in our case we choose the starting position as R = x i , Z = 0 and ϕ = 0).Note that the equilibrium gradients, ∂n 0 /∂x and ∂T e0 /∂x, depend on x and ϕ and are computed from the non-linear simulations.The curvatures κ n and κ g depend on ϕ and very weakly on x.They are shown in figure B1 on the LCFS.For comparison, the reference values κ n = − cos(ιϕ) and κ g = sin(ιϕ) of an infinite aspect-ratio tokamak with circular flux surfaces are also shown.

Figure 4 .
Figure 4. Amplitude of the magnetic field on the LCFS as a function of the polar angle θ defined with respect to the magnetic axis and the toroidal angle ϕ.Only one field period is displayed.

Figure 5 .
Figure 5. (Left) 1/10 of the simulation box.Note how the height and width of the cross-section of the simulation domain are adjusted to create an island divertor stellarator.(Right) Connection length of the field lines normalized to the major radius, Lc/R 0 .

Figure 6 .
Figure 6.(a) Top view of the simulation box.Black dots correspond to field-lines of the island chain leaving the box (B Z > 0) and red dots to field-lines entering (B Z < 0).The contours correspond to the time-averaged electron pressure obtained from the GBS simulations on the target.(b) Time-averaged electron plasma pressure on the top of the simulation box, as obtained from the GBS non-linear simulation in steady-state.

Figure 7 .
Figure 7. Density, electron temperature and ion temperature sources at two different poloidal planes.

Figure 8 .
Figure 8. Volume integral of each term in the density equation.The system reaches a quasi-steady state at t ≈ 70 R 0 /c s0 .The mode m = 4 appears at t ≈ 35 R 0 /c s0 .

Figure 9 .
Figure 9. Equilibrium profiles of density (left) and electron temperature (right), obtained by time-averaging the simulation results.Top and bottom correspond to the toroidal planes ϕ = 0 and ϕ = 0.6, respectively.

Figure 10 .
Figure 10.Equilibrium profiles of electrostatic potential (left) and parallel current (right), obtained by time-averaging the simulation results.Top and bottom correspond to the toroidal planes ϕ = 0 and ϕ = 0.6, respectively.

Figure 12 .
Figure 12.Snapshot of fluctuations of density (left) and electrostatic potential (right).Top and bottom correspond to the toroidal planes ϕ = 0 and ϕ = 0.6, respectively.

Figure 13 .
Figure 13.Snapshot of potential fluctuations on the LCFS (left) and its Fourier spectrum (right).The poloidal angle θ * is such that field lines (black arrows) are straight in the (θ * , ϕ) plane.

Figure 14 .
Figure 14.Snapshot of density fluctuations at two different planes in the configuration without islands.

Figure 15 .
Figure 15.Time-averaged radial E × B flux on the LCFS for positive (left) and negative (right) magnetic field.

Figure 16 .
Figure 16.Product of the density and potential eigenmodes, |n 1 Φ 1 |, for the forward and reversed magnetic field along the toroidal direction.The vertical dashed lines corresponds to the positions θ = π/2 (red) and θ = 3π/2 (grey).These are the eigenmodes that correspond to the eigenvalue with largest real part.In both cases, γ ≈ 3.

Figure 17 .
Figure 17.Equilibrium gradients, retrieved from the non-linear GBS simulation (top) and the normal and geodesic curvatures (bottom) used in the model in equation (22).

Figure 18 .
Figure 18.Equilibrium E × B drift, ⟨V ExB ⟩ t , in the forward and reversed field simulations.

Figure B1 .
Figure B1.Normal, κn, and geodesic, κg, components of curvature κ along the field line that spans the LCFS of the Dommaschk equilibrium.For comparison, the curvatures of an infinite aspect-ratio tokamak with circular flux surfaces are shown in red, where ι = 0.547.