Global fluid turbulence simulations in the scrape-off layer of a stellarator island divertor

Isothermal fluid turbulence simulations have been performed in the edge and scrape-off layer (SOL) of an analytic stellarator configuration with an island divertor, thereby providing numerical insight into edge turbulence in regions around islands in a stellarator. The steady-state transport follows the a curvature drive that is inverse to the major radius ($1/R$) toward the outboard side, but large fluctuations are present throughout the island divertor region, with the average wavelength of similar size to the island width. The system exhibits a prominent $m=2$, $n=5$ mode, where m is the poloidal mode number and n is the toroidal mode number, although other modes are present. The amplitude and radial extent of the density fluctuations are similar throughout the edge and SOL, but can decrease near island O-points. The fluctuations exhibit a predominantly positive skewness on the outboard midplane, indicating blob-like perturbations for the transport into the outer SOL. It is determined that a point on the separatrix is generally more correlated with regions outside of the SOL than a nearby reference point which does not lie on the separatrix.


Introduction
Advances in stellarator optimization have led to unprecedented improvement in neoclassical transport in Wendelstein 7-X (W7-X) (Beidler et al. 2021) such that anomalous transport is now contributing a significant portion of the transport, especially in the scrape-off layer (SOL) (Pablant et al. 2018).The island divertor concept, originally proposed nearly five decades ago (Karger & Lackner 1977), utilizes low-order magnetic islands which intersect the wall to manage heat and particle exhaust (Feng et al. 2006;Feng & W7-X-team 2022).The island divertor is advantageous in that it potentially provides a higher whetted area, and has been able to enter very stable detachment regimes (Pedersen et al. 2019).Due to this island divertor, the SOL of W7-X provides a unique environment for SOL physics, as the toroidally discontinuous intersection of magnetic islands creates complicated topologies which can inhibit numerical investigation.As such, turbulence simulations in stellarators have focused on the core transport using either a gyrokinetic (Xanthopoulos et al. 2007;Bañón Navarro et al. 2020;Maurer et al. 2020;Singh et al. 2022) or, less commonly, a fluid approach (Kleiber & Scott 2005).Recent work using local fluid simulations has informed interpretation of experimental measurements (Shanahan, Dudson & Hill 2018;Killer et al. 2020;Shanahan et al. 2021;Huslage et al. 2023), without requiring modelling of the complex and numerically challenging SOL topology.Developments of numerical methods have provided the opportunity for global fluid turbulence simulations in stellarator geometries (Shanahan, Dudson & Hill 2019;Coelho et al. 2022), but a general understanding of global phenomena at the edge and SOL of stellarators is only in its infancy.In this work, we utilize an isothermal model in the BOUT++ (Dudson et al. 2009) framework to explore the nature of turbulence in the island divertor region at the edge of an analytic stellarator configuration (Coelho et al. 2022).
Section 2 details the methods used in this work, including the plasma model, geometry and computational grid.Section 3 discusses the simulation results, and is divided into two subsections; the steady-state transport properties of the system and the dynamics of the fluctuations in the island SOL.The implications of this work and its context within previous work are discussed in § 4.

Methods
The strength of BOUT++ (Dudson et al. 2009) is its flexibility: models and numerical methods can be easily changed to suit the problem to be addressed.Here, we use an isothermal plasma model in a complex numerical scheme, as the turbulent dynamics is available despite numerous assumptions, but the geometry is necessarily complicated.

Isothermal plasma model
In this work we exploit a reduced magnetohydrodynamic model from the Hermes family of models (Dudson & Leddy 2017;Leddy et al. 2017;Huslage et al. 2023), which do not separate fluctuations from the background.The simulations here are isothermal, and evolve number density n, vorticity ω, parallel ion momentum m i nv and Ohm's law where p e is the electron pressure, e is the electron charge, m i is the ion mass and φ is the plasma potential, for the parallel electron velocity v e .Quasineutrality is assumed, so that the electron and ion densities are equal: n e = n i = n.Here, ν = 1.96τ ei (m i /m e ) is the resistivity.The thin layer (Oberbeck-Boussinesq (Oberbeck 1879)) approximation is made in calculating the potential φ from vorticity ω such that where Ω = eB/m i is the ion cyclotron frequency, where B is the magnetic field strength and n 0 is a constant typical number density, in this work set to be n 0 = 1 × 10 18 m −3 .The resulting equations for (electron) number density n, vorticity ω and parallel momentum m i nv i are where S n is the density source, p e/i is the pressure for electrons and ions, respectively, and the drift terms for E × B, v E×B and ion and electron magnetic drifts, v mag,e/i , are defined as where E is the electric field, T e,i is the electron and ion temperature, respectively, b is a unit vector in the direction of the magnetic field, and f is an arbitrary function.Here, the notation is such that While the model allows terms for anomalous diffusion, these terms are not necessary for numerical stability and, as the nature of their origin is unknown, are neglected in this work.An imposed parallel diffusion of 4.2 m 2 s −1 is used in Ohm's law to suppress numerical instability.

BSTING and BOUT++
The BSTING project (Shanahan et al. 2019) has recently allowed for simulations of stellarator geometries using the BOUT++ framework (Dudson et al. 2009), in which exploitation of the flux-coordinate-independent (FCI) method for parallel derivatives (Hariri & Ottaviani 2013;Hariri et al. 2014;Hill, Hariri & Ottaviani 2015;Shanahan, Hill & Dudson 2016;Stegmeir et al. 2016;Hill, Shanahan & Dudson 2017) allows for simulation of complex magnetic geometries including magnetic islands and chaotic field regions.The implementation of three-dimensional metric tensor components in BOUT++ has allowed for poloidally curvilinear FCI grids (Shanahan et al. 2019) that do not include grid points in the core of the plasma -where the fluid approximation due to high collisionality can break down -and allow for the poloidal grid alignment to flux surfaces or plasma-facing components.As such, FCI simulations can be performed with resolutions comparable to field-aligned tokamak simulations.

Geometry and initial conditions
The simulation geometry used here is similar to that in Coelho et al. (2022).Namely, a stellarator magnetic field with an outer m = 9 island chain was created using the same Dommaschk potential (Dommaschk 1986) as in Coelho et al. (2022), where the magnetic field varies to lowest order by 1/R, where R is the major radius -similar to (for instance) a tokamak.The BSTING framework utilizes curvilinear grids, however, meaning that, in contrast to the work presented in Coelho et al. (2022), here, the core is neglected, and the outer poloidal surface is an ellipse with the same maximum dimensions of the outer boundary as used in Coelho et al. (2022).The inner surface is aligned to a vacuum flux surface, and the elliptical outer surface provides a plasma-facing component which intersects islands at discrete toroidal locations (for instance, figure 1b), as in an island divertor.
The numerical grid utilized in this work was of the size (radial, toroidal, poloidal) (nx, ny, nz) = (68, 128, 256) with an average resolution of (dx, dy, dz) ≈ (1, 75, 1)ρ i , where ρ i is the ion Larmor radius.The simulation domain spans the entire toroidal extent (0 to 2π).Despite the difference in numerical implementation, the magnetic field is identical to that in Coelho et al. (2022), including the outer magnetic islands, as shown in figure 1.The radial extent of the SOL can be as small as 10 Larmor radii, and perturbations often exhibit a similar perpendicular size.As such, fluctuations cannot be categorized as edge or SOL localized.A more detailed discussion of the fluctuations follows in § 3.2.The simulations are flux driven, with a particle source of the form where S 0 is a prescribed constant used to balance the sinks, r 0 is the radial location the peak and w is the Gaussian root-mean-square width.Here, r 0 and w are chosen such that the source peaks at the same radial location in the closed-field-line region as the source in Coelho et al. (2022).Since the particle source introduces an external drive into the system, the dynamics which originates farther inward (toward the core) must be carefully treated -and is therefore not included in the computational domain.The equilibrium parameters -the background density (n 0 = 1 × 10 18 m −3 ), the magnetic field (B 0 ≈ 0.3T) and the background temperature (T 0 = 10 eV) have been chosen to most closely align with Coelho et al. (2022).These parameters are similar to those found in university-scale experiments, such as TORPEX (Furno et al. 2008;Shanahan & Dudson 2016).Parallel sheath boundary conditions (Walkden et al. 2016) are imposed such that the species are accelerated to the local sound speed and implemented using the leg-boundary-fill method at the ends of magnetic field lines (Hill et al. 2017).Perpendicular boundary conditions, with the exception of plasma potential, are zero gradient, allowing flows into the domain from the inner boundary and out of the domain at the outer boundary.The plasma potential is calculated from vorticity using PETSc (Balay et al. 1997) with Dirichlet boundary conditions such that φ boundary = 2.8T e .The density is set to an approximate initial profile, and an initial perturbation in vorticity is added to the system; after a transient phase, the simulation settles into a state where the particle source is balanced by the sinks in the sheath, characterized by fluctuations throughout the domain.The simulation is then run such that the time scale exceeds 1 ms, which is calculated in approximately 30 000 core hours.

Results
The dynamics simulated here can be separated into two scales; large-scale, steady-state flows of the plasma and the perturbations on top of these macroscopic aspects.The following two subsections will look into this dynamics in more detail.To further examine the flows in the island SOL presented here, one can plot the Pearson correlation coefficient (Freedman, Pisani & Purves 2007) of density signals relative to a given point in the SOL, as is done in Shanahan & Dudson (2014).In this work, we correlate the time trace of every point in the domain to a reference point, defining the Pearson correlation coefficient P fg between a reference point f and a sample point g as  For this reason, we will look at the relative correlation when choosing two reference points; one correlation with a reference point at the separatrix P sep , and one with a reference point at an island O-point, P O .
Figure 5 shows the difference of two correlations given by (3.1); one where the reference point is in the separatrix, f sep , and one where the reference point is at an O-point, f O .The reference points are shown as white markers in figure 5(a).By plotting their difference,

Perturbations near the edge and SOL
Having examined the characteristics of the large-scale dynamics, attention is now turned to the dynamics of the perturbations.The system exhibits a strong m = 2, n = 5 mode, although other modes are present.The fluctuations are largek ⊥ ρ s is generally less than 1, as indicated in figure 6.As the magnetic islands are, at widest, only a couple of tens of ρ i the fluctuations can at times exceed the width of the island SOL.To determine the nature of the fluctuations, we examine central moments of the time signal, where the nth central moment of a quantity x is given by μ n = (x − x ) n , where the angled brackets ( ) indicate a time average.The standard deviation of x, σ x , is the square root of the second moment (n = 2) and the skewness, which will be used later, is proportional to the third moment (n = 3) normalized to the standard deviation.The standard deviation of the density signal, σ n is shown in figure 7. Figure 7 indicates that fluctuations are present throughout the domain.The vertical cross-section, figure 7(b), indicates a decrease in σ n outside the outer midplane island O-point relative to the neighbouring X-point regions, suggesting that fluctuations are reduced here.Figure 7 indicates that the fluctuations, despite their large scale relative to the system size, fall off radially.Using the correlation coefficient, (3.1), it is determined that the fluctuations have radial correlation lengths which vary from approximately 8ρ i on the outboard midplane of the verticalcross-section (figure 5b) to approximately 50ρ i on the inboard side of the horizontal cross-section (figure 5a).Since the islands in the SOL are wider at the locations of longer correlation lengths, this poloidal variation in correlation length could be related to the island width.Furthermore, the standard deviation of the vorticity perturbations, σ ω , also shows fluctuations present throughout the domain, figure 8. Figure 8(b) indicates again that the outboard midplane island O-point in the vertical cross-section exhibits reduced fluctuation amplitudes, whereas the neighbouring X-points show stronger fluctuations.Vorticity localized to islands can in turn lead to transport around the island separatrices.Island-localized potentials leading to flows around the island separatrix have been seen in W7-X (Killer et al. 2019).It is also possible to plot the radial flux of the perturbations, Γ⊥ = ñ ṽr , illustrated in figure 9.The radial perturbation flux can increase near X-points, and is indeed seen in several studies in other geometries (Poli, Bottino & Peeters 2009;Hill et al. 2015;Choi 2021), but this phenomenon is not universal in the simulations presented here and a concrete conclusion cannot be drawn.The radial perturbation flux shown in figure 9(b) indicates stronger outboard activity, with a change in the sign of the radial flux near the outboard midplane.
To determine the transport nature of the fluctuations, the skewness of the profiles is plotted in figure 10, where a positive skewness indicates blob-like transport (Walkden et al. 2017), and negative skewness indicates areas dominated by the propagation of negative perturbations (holes).From figure 10, it is apparent that the transport near the outer boundary, at the inboard and outboard midplane locations mostly include positive  perturbations, whereas in the middle of the domain, particularly the top and bottom of the configuration exhibit a predominantly negative skewness, where negative perturbations are more prevalent.

Summary and implications
In this work, a detailed analysis of fluid turbulence in a stellarator island divertor is presented.An isothermal fluid turbulence model was used to simulate turbulence in an analytic stellarator geometry, finding that the fluctuations are present throughout the island divertor region.It was determined that the steady-state dynamics is mostly consistent with the 1/R curvature drive, but the sheath connection in the cross-sections which intersect the boundary introduces a discontinuity.The plasma outside the SOL is more highly correlated with the point on the separatrix than the island O-point, indicating the separatrix as a transport channel into the private flux region.The fluctuations are exhibited throughout the domain, with a smaller amplitude near the outboard midplane island in the vertical cross-section.The radial perturbation flux can, but must not necessarily manifest near island X-points.Future work will look to relax the isothermal approximation and remove the islands from the geometry in order to more fully assess the impact of islands on the transport and perturbation dynamics, although this is a laborious task due to the non-intuitive nature of Dommaschk potentials.While it is true that by adjusting higher-order coefficients one can almost remove islands, this process is described by the author in Dommaschk (1986) as only achieved 'by trial.'It could be more interesting, therefore, to explore more realistic geometries, such as the various W7-X configurations which have been more extensively investigated.

FIGURE 1 .
FIGURE 1. Poincaré plots (black) of the magnetic field using parameters from Coelho et al.(2022), including the inner (red) and outer (blue) surfaces of the grid.

FIGURE 2 .
FIGURE 2. Time-averaged plasma potential; the contour lines indicate flows due to E × B advection.Transport towards the outboard side due to 1/R curvature drive is the dominant flow within the system.In this figure, the magnetic field is predominantly into the page.

FIGURE 3 .
FIGURE 3. Time-averaged radial flux, indicating regions of radial transport.A negative radial flux indicates radially inward flux, whereas a positive flux is outward.

FIGURE 4 .
FIGURE 4. The Pearson correlation coefficients when considering a reference point (a) in the separatrix, P sep , and (b) when choosing a reference point at a magnetic O-point, P O .These reference points are marked with a white X and a white dot, respectively.A positive value indicates positive correlation with the reference point.

FIGURE 5 .
FIGURE 5.The relative Pearson correlation coefficients when considering a reference point in the separatrix, P sep , and when choosing a reference point at a magnetic O-point, P O .These reference points are marked with a white X and and circle in figure 5(a).A positive value indicates stronger correlations with fluctuations at the separatrix, whereas a negative value indicates correlation with fluctuations at an island O-point.
FIGURE 6.(a) Average radial fluctuation size and (b) average poloidal fluctuation size on a flux surface indicated by the dashed surface in figure 1.

FIGURE 7 .
FIGURE 7. Standard deviation of the density; perturbations are apparent at all poloidal locations, with a slight increase on the outboard side.

FIGURE 8 .
FIGURE 8. Standard deviation of vorticity.The horizontal cross-section (a) indicates better poloidal symmetry compared with the vertical cross-section (b).

FIGURE 9 .
FIGURE 9. Time average of the radial perturbation flux.