Plasma turbulence simulations in a diverted tokamak with applied resonant magnetic perturbations

The first results of three-dimensional, flux-driven, electrostatic, global, two-fluid turbulence simulations of a diverted tokamak configuration with applied resonant magnetic perturbations generated by a set of saddle coils are presented. The simulations of an L-mode plasma show that the heat flux pattern on the divertor targets is affected by the resonant magnetic perturbations, as a result of the interplay between turbulent cross field transport and parallel flows. The simulation results reveal the potential of resonant magnetic perturbations to reduce the heat flux to the wall. In fact, the peak of the toroidally- and time-averaged heat flux as well as its value integrated over the divertor decrease as the amplitude of the magnetic perturbation increases, while the plasma sources are held constant.


Introduction
Resonant magnetic perturbations (RMPs) are threedimensional, low-amplitude perturbations applied to a tokamak equilibrium magnetic field, usually generated by a set of saddle coils placed around the vessel. By creating a region of chaotic field lines at the edge of the confined region, RMPs are frequently used in tokamaks to control or even suppress edge localized modes (ELMs) [1][2][3], while keeping the core relatively unperturbed. At the same time, RMPs can enhance the radial transport of particles and heat in the scrape-off layer * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Conclusive experimental observations of the impact of RMPs on the divertor walls are still missing, calling for a dedicated simulation effort. On the one hand, an increase of the stationary heat flux to the divertor during the ELM suppression phase was observed [8], as well as total absence of differences in the peak heat flux to the divertor and power SOL falloff length with respect to the unperturbed situation [9]. On the other hand, experiments in different machines have shown the benefits of RMPs for handling the heat exhaust: when RMPs are applied, the broadening of the heat flux deposition profile at the targets has been reported [10], as well as the increase of the toroidally-averaged SOL decay length [11], thus providing hints for potential reduction of the heat flux on the divertor targets. In fact, while periodic coil configurations in the toroidal direction might create local maxima and minima of the heat flux deposition on the divertor targets, the effect can be toroidally averaged by applying alternating current in the coils [12], with a period several orders of magnitude larger than the turbulent time scale.
The impact of RMPs on the heat exhaust has been also investigated by using simulations. Turbulence flux-driven simulations in the presence of magnetic perturbations were performed considering a circular tokamak geometry. By using a cold-ion, isothermal, dissipative drift-Alfvén model, [13] reports on a strong suppression of density fluctuations in the presence of RMPs, attributed to a change of the equilibrium profiles of the density and electric field. By means of a twofluid, isothermal, drift-reduced Braginskii model, [14] reports on the decrease of the fluctuation levels, suggesting a stabilizing effect of the RMPs on turbulent transport. In [15], a nonisothermal gyrofluid model that takes into account the plasma screening of the applied RMPs is used. The simulations show that the radial transport induced by parallel motion along radially perturbed magnetic field lines is not increased when RMPs are applied, while RMPs lead to a decrease of the turbulent fluctuations.
In the present paper, we study the influence of RMPs on L-mode plasma turbulence by presenting the first results of global, three-dimensional, two-fluid, flux-driven turbulence simulations of a tokamak with applied RMPs in a diverted configuration using the GBS code [16][17][18]. For the first time, a diverted configuration is considered, and our simulations use a coil configuration around the tokamak vessel to emulate the RMPs. The size and parameters of the simulations presented here reflect those of a tokamak of, approximately, one-third the size of the TCV tokamak [19].
Our simulations show that RMPs have the potential to reduce the toroidally averaged heat flux to the divertor targets. While the reduction of the heat flux we observe is relatively modest, our results call for further experimental and theoretical investigations to optimise the RMPs in order to significantly reduce the heat flux.
The paper is organized as follows. Section 2 describes the physical model implemented by GBS. In section 3, we detail the magnetic field considered for our simulations. Section 4 presents the simulation results. Finally, we draw our conclusions in section 5.

GBS code
The simulations are performed with the GBS code, developed in the past decade to simulate turbulence in the tokamak boundary [16][17][18], and more recently extended to allow for the simulation of three-dimensional equilibrium magnetic fields, such as in stellarators [20]. GBS solves the drift-reduced Braginskii equations [21], valid in the high collisionality regime, an assumption often justified in the plasma boundary of magnetic fusion devices, as well as in low-temperature devices (e.g. TORPEX [22]). In GBS, all quantities are evolved in time without separation between equilibrium and fluctuating components. The plasma dynamics results from the interplay between plasma sources, turbulent perpendicular transport, and parallel flows that eventually lead to losses at the vessel walls.
For our simulations, we consider an electrostatic model and apply the Boussinesq approximation [16], both shown not to have an impact on the simulation of L-mode plasmas [23]. We also neglect the coupling to the neutral dynamics, although these effects are taken into account in the most recent version of the GBS code for tokamak axisymmetric simulations [17,18]. Within these approximations, the drift-reduced fluid model evolved by GBS is: which are closed by In equations (1)- (7), and in the rest of the paper, all quantities are normalized to typical SOL reference values. Density n, electron temperature T e and ion temperature T i are normalized to n 0 , T e0 and T i0 ; 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 constant in the present work, 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, where λ denotes the Coulomb logarithm [24]. Small numerical diffusion terms such as D n ∇ 2 ⊥ n and D ∥ n ∇ 2 ∥ n (and similar for the other fields) are introduced to improve the numerical stability of the simulations (the simulation results show that they have a negligible effect on turbulence since they lead to significantly lower perpendicular transport than turbulence). The terms S Te , S Ti and S n denote the electron temperature, ion temperature and density sources, respectively.
The geometrical operators appearing in equations (1)- is the unit vector of the magnetic field. They are expanded as described in [20].
The model in equations (1)-(7) is solved in a cylindrical coordinates system, (R, ϕ, Z), with R the radial coordinate, ϕ the toroidal angle and Z the vertical coordinate. The simulation domain is a torus of radius R 0 with a rectangular crosssection (grey box in figure 1). Equations (1)-(6) are advanced in time with a standard explicit Runge-Kutta fourth-order scheme, while spatial derivatives are computed with a fourthorder finite difference scheme [24,25]. We apply magnetic pre-sheath boundary conditions at the bottom (i.e. the divertor targets) and top boundaries of the simulation domain, as described in [26,27], except for Φ on the divertor plates where we use an insulating boundary condition. A summary of the boundary conditions is given in table 1.
The following parameters are used in our simulations. The extension of the rectangular cross section is 200 < R/ρ s0 < 740 and 0 < Z/ρ s0 < 800, The grid resolution is n R = 216, n Z = 320 and n ϕ = 80. The time step is of the order of 10 −5 R 0 /c s0 . We note that the values of χ ∥e,i are underestimated by approximately one order of magnitude to make the computational cost of the simulations manageable. However, since we assume that we are in the sheath-limited regime, where convection dominates over conduction, i.e. nT e v ∥e ≫ χ ∥e ∇ ∥ T e and temperature gradients along the parallel direction are small, underestimating the parallel heat diffusion coefficients does not have a relevant impact on the simulation results. Moreover, we note that the ratio m i /m e is set to a small value also for computational reasons. The m i /m e ratio determines if turbulence is driven by the inertial or the resistive branch of the ballooning instability. Indeed, turbulence in the L-mode regime is typically driven by the resistive branch of the ballooning instability, as shown in [23]. The resistivity used in the present simulations is sufficiently large that it compensates the small m i /m e used (i.e. (m e /m i )γ < ν with γ the growth rate of the resistive ballooning mode), placing the present simulations in the correct parameter regime.

Equilibrium and RMP magnetic field
The magnetic field used in the present simulations, B = B 0 + δB, is the superposition of the equilibrium field of a diverted tokamak configuration, B 0 , and a RMP generated by a set of saddle coils around the vessel, δB. The poloidal component of the diverted magnetic equilibrium, considered also in previous work [24], is generated by a current with a Gaussian profile centered at the location (R 0 , Z 0 ) inside the vessel, mimicking the effect of the plasma current, and a current filament located at (R 0 , Z 1 ) outside the vessel, mimicking the effect of a poloidal field coil that diverts the magnetic field lines to the bottom wall. The associated flux function is where R 0 = 500ρ s0 , Z 0 = 500ρ s0 , Z 1 = −140ρ s0 , σ = 100ρ s0 and α = 20ρ s0 . In equation (9), E 1 is the exponential integral defined as E 1 (x) =´∞ x e −t /t dt. The toroidal field is given by The perturbed magnetic field δB is computed numerically using the MAKEGRID code from the STELLOPT package [28], which uses the Biot-Savart law to determine the magnetic field at a specified location. The RMP coil configuration, displayed in figure 1, consists of six pairs of coils located around the vessel, with current alternating in the toroidal direction, thus creating a perturbation with a toroidal mode number n = 3. The coils are placed 5ρ s0 away from the outer wall of the device. They are 639ρ s0 long in the toroidal direction, covering an angle of approximately 0.858 rad, and are evenly spaced. This results in a gap between them of 141ρ s0 , corresponding to an angle of approximately 0.189 rad. The lower coils are placed at 100 ⩽ Z/ρ s0 ⩽ 300 and the upper ones at 500 ⩽ Z/ρ s0 ⩽ 700. This specific choice of coils is the result of an exhaustive investigation to find a configuration that provides a large resonant perturbation in the edge and a smaller resonant perturbation in the core, with many field lines in the SOL performing more than one poloidal turn when going from one divertor plate to the other (the high-field (HFS) and low-field side (LFS) divertor plates in figure 3). The amplitude of the perturbation is controlled by the current inside the coils, and is expressed by using the parameter ε = max (δB LCFS )/B 0 , where max (δB LCFS ) is the maximal value of |δB| at the last closed flux surface (LCFS). In order to identify the resonating flux surfaces for the imposed n = 3 perturbation, we compute the perturbation spectrum (figure 2), which reveal the flux surfaces that are expected to be more affected by the applied perturbations. The spectrum is given by the Fourier coefficients b mn (ψ) of the normalized amplitude of the perpendicular perturbation on a flux surface ψ, where θ * is a straight-field-line angle for the equilibrium magnetic field, q(ψ) the flux surface averaged safety factor of the equilibrium magnetic field and b mn = √ (b c mn ) 2 + (b s mn ) 2 . In figure 2, the safety factor is also displayed, highlighting the q = m/n values for n = 3. On a rational flux surface ψ r , where q(ψ r ) = m/n, the width of the islands generated by the perturbation is proportional to √ |b mn (ψ r )| [29]. When the width of the islands is larger than the distance between two rational flux surfaces, a region of chaotic field lines is created [30]. This is the case in our simulations. The perturbation spectrum for the considered set of coils shows that the islands are wider and more densely packed close to the edge than inside the core. Indeed, the Poincaré sections displayed in figure 3 show a chaotic field line region in the edge, which penetrates deeper into the core as ε increases. For simplicity, the response of the plasma to the externally applied RMPs is not considered in this work. Taking this response into account is a complex and a long term goal, where simplified models, for example based on introducing screening currents on rational surfaces [31,32], could be considered.  In figure 4, the footprints of the field lines on the divertor target are displayed for different perturbation amplitudes. These are computed by considering a set of field lines regularly spaced along the R and ϕ directions that start on the HFS divertor plate (Z = 0 and R < R 0 ). These field lines are followed and their landing point on the LFS divertor plate (Z = 0 and R > R 0 ), which we consider as the field line footprint, is evaluated. The number of poloidal turns, N p , described by a field line before landing on the LFS divertor is monitored and color-coded in figure 4. High values of N p are expected to yield large heat flux, if one assumes that convection in the parallel direction dominates over perpendicular transport. In fact, a field line with N p > 1 enters the hot plasma region, not accessible to any field line ending on the divertor plates in the unperturbed case, and allows particles and heat, otherwise confined, to flow along the field lines to the divertor plates. The regions of N p > 1 are lobe structures with a n = 3 periodicity.
The radial extension of the lobes on the LFS divertor increases linearly with the amplitude of the perturbation.

Simulation results
A simulation with ε = 0 is started from a noisy initial state and, after a transient, it reaches a quasi-steady state where sources, parallel and perpendicular transport, and losses at vessel balance each other. Starting from this quasi steady-state, we consider five simulations with ε = 2.12%, 3.71%, 4.24%, 4.91% and 5.83%, consecutively, and evolve them until a quasi-steady state is reached. The time-step used for these simulations is, respectively, 2 × 10 −5 , 2 × 10 −5 , 2 × 10 −5 , 1.7 × 10 −5 , 1.3 × 10 −5 and 1 × 10 −5 R 0 /c s0 . A convergence test was performed with a grid 1.5 times finer in all directions for the ε = 4.91% simulation, obtaining similar results to the simulation presented here. The temperature sources S Te = S Ti have a constant value inside the LCFS and vanish outside it, while the density source S n is a radially localized Gaussian around a closed flux surface inside the core that is close to the LCFS. All sources are independent of the toroidal angle ϕ.
The RMP degrades the confinement in the core, with the pressure dropping by 35% in the case of highest ε, with respect to the unperturbed case. We note, however, that our simulations overestimate the effect of the magnetic perturbations on core transport. In fact, the MHD response of the plasma, which is not included in the present simulations, tends to shield the equilibrium flux surfaces from the magnetic perturbations in the low resistivity region of the core (ν ∼ T −3/2 e ), hence potentially maintaining its good confinement. On the other hand, this shielding is less efficient in the edge and SOL regions because of the high plasma resistivity, thus the effect of the RMPs is expected to persist there. Figure 5 displays typical snapshots of the relative density fluctuations for different values of ε, showing fluctuations reaching order unity, as expected in the boundary of tokamaks, and the presence of coherent structures, denoted as blobs [33], particularly in the far SOL. In order to investigate the nature of the turbulent fluctuations, we carry tests where we toroidally average the curvature term in the vorticity equation, equation (6). This yields a suppression of the turbulent fluctuations and steepening of the gradients for all the values of ε considered, showing that turbulent transport is mostly driven by a ballooning instability [24,34]. An analysis of the number of blobs and their size is carried out using the analysis technique described in [35]. It consists of tracking a structure of enhanced density (at least 2.5 times above the standard deviation normalized to the time averaged values) that moves coherently (i.e. it exists for at least ∆t > 0.2). The blob velocity is then determined as the velocity of the center of mass of the identified structure. The blob analysis show very similar results for all simulations. In particular, the blob contribution to radial transport is found to be around 30% ± 10% in the far SOL region for all simulations, similarly to [36]. On the other hand, the contribution of the parallel flux to radial transport increases with ε, reaching 20% of its total value in the case of the highest perturbation.
Toroidally and time-averaged pressure profiles in the SOL are displayed in figure 6 for R > R LCFS at the midplane (Z = 484ρ s0 ) and at the divertor target (Z = 0). First, we note that the pressure decreases both at the midplane and at the divertor as the amplitude of the perturbation increases. Second, two decay lengths (the near and the far SOL) can be distinguished in the ε = 0 case at the midplane (dashed lines in the plot), as observed in other tokamak simulations [36,37]. These decay lengths are remapped at the divertor into a single decay length with a value that is in between the two decay lengths at the midplane. In the presence of the RMP, because of the chaotic nature of the field lines in the edge and SOL regions, the midplane and divertor pressure profiles are decoupled, in the sense that they cannot be directly related to each other by a local remapping of the pressure along each flux tube. At the midplane, the two decay lengths present in the ε = 0 merge into a single one and a region of constant pressure appears in the very far SOL, close to the wall. The origin of this flat density region might be related to the boundary condition applied to the outer wall and will be the subject of future investigations. The near SOL decay length of the ε = 0 simulation and the decay length of the other simulations agree well with the edge equilibrium pressure gradient length for ballooning-driven turbulence, derived in [24], within an error of 30%, hence confirming the fact that the nature of turbulence is not affected by the RMPs.
The equilibrium radial electric field, ⟨E r ⟩ t,ϕ , at the outboard midplane is shown in figure 7 (left). The radial coordinate r is defined with respect to the unperturbed equilibrium, i.e. E r = −∇Φ · ∇ψ/||∇ψ||. For ε = 0, we observe the formation of a well in the electric field, typical of the tokamak plasma edge, with E r < 0 in the core and E r > 0 in the far SOL. As the perturbation increases, the E r profile flattens, an  effect observed in previous simulations [14], but also experimentally in different tokamaks [38,39], where in some situations the sign of E r even reverses [40,41]. The peak value of the E × B shear in the highest ε simulation is reduced by 50% with respect to the ε = 0 case. However, the effect of the E × B shear on turbulence is expected to be small since the linear growth rate of the driving instability is larger than the shearing frequency in the present simulations. As a matter of fact, the turbulence fluctuation level, defined here as σ n /⟨n⟩ t , where σ n is the density standard deviation, decreases  with increasing perturbation, as seen in figure 7 (right). Indeed, it should be noted that the fluctuation level of the highest perturbation simulation is decreased by a factor larger than 2 with respect to the ε = 0 simulation. Since the input sources are the same for all values of ε and the parallel contribution to the radial transport increases with increasing perturbation, the radial turbulent E × B flux must decrease, and that is accomplished by the observed reduction in the fluctuation levels.
We now analyze the effect of the perturbation on the heat flux deposition on the divertor targets (corresponding to the Z = 0 plane). Since the fluxes to the divertor induced by the E × B and diamagnetic drifts are small when compared to the parallel flow (even when projected in the direction normal to the wall), we construct a measure of the heat flux to the wall as q = nT e v ∥e b, also neglecting the ion contribution since this is smaller than the electron one. We remark that the heat flux is evaluated at the boundary of the simulation domain, which corresponds to the magnetic pre-sheath entrance. The flux to the wall is then affected by the processes taking place inside the sheath, where a transfer of energy from the electrons to the ions takes place, and the total power being conserved. We project the parallel flux along the normal direction pointing outwards from the simulation domain, and therefore we compute the divertor heat flux as q = −nT e v ∥e B Z /B. The time averaged value of the divertor heat flux, ⟨q⟩ t , on the LFS divertor is displayed as a function of R and ϕ in figure 8.
In the unperturbed case, the heat flux is axisymmetric and localized around the leg of the separatrix displayed with a dashed black line (small deviations from axisymmetry are ascribed to the averaging over a finite time window). As the amplitude of the perturbation increases, three hot spots appear around the line where the separatrix is located in the ε = 0 case. Between these hot spots, three regions with low heat flux are observed. This number corresponds to the toroidal periodicity of the perturbation. It can also be observed that the heat flux increases with respect to the ε = 0 case in regions radially far from the separatrix. Figure 8 shows that the heat flux pattern does not follow exactly the pattern set by the field line footprints because of turbulence. In fact, when only the parallel transport is considered and turbulence is neglected, as in previous work [42], the heat flux pattern on the divertor matches that of the lobes generated by the footprints of the field lines. However, when turbulence is considered, the pattern is smoothed out.
To allow a comparison between heat flux deposition profiles for different RMP amplitudes, the heat flux is averaged over the toroidal direction. The result, ⟨q⟩ t,ϕ , is displayed as a function of R, for different perturbation amplitudes, in figure 9 (left). A remarkable result is the decrease of the peak heat flux as ε increases, which can be observed more clearly in figure 9 (right). The full analysis of the integrated heat flux over the walls shows that, as ε increases, the integrated heat flux on the bottom wall (where the divertor plates are) decreases by 29% for the highest value of ε (with respect to ε = 0 case), while it increases by 62% on the outer vertical wall, where the magnetic field is stronger due to the RMP coils and a parallel heat flux from the core to the vertical wall is thus allowed. These two effects compensate each other within ∼15%, ensuring the global balance in the simulations within the limitation of the physical and numerical scheme adopted in the present simulations [43].

Conclusions
In this paper, GBS simulations reveal the mechanisms behind plasma transport in the boundary of a diverted tokamak with applied RMPs. The perturbation is chosen to create chaotic regions in the plasma edge. Results show that the perturbation impacts mostly the time-averaged quantities, since the nature of turbulence and the transport due to blobs remain the same as the RMP amplitude is increased. It is observed that the near and far SOL decay lengths merge into a single one in the presence of the perturbation, and that the heat flux pattern on the divertor does not follow exactly the footprints created by the field lines, but instead is spread around them. As future improvements of the simulation model presented here, we mention the use of a set of more realistic boundary conditions on the outer wall, which might affect the plasma profiles in its proximity, as well as of a model to take into account the self-consistent plasma response to the externally applied magnetic perturbations. The peak heat flux on the divertor plates is effectively reduced as the perturbation is increased. This is achieved through the reduction of the total parallel flux towards the diverted targets (bottom wall of the simulation), compensated by an increase of the flux on the outer vertical wall. This confirms the potential of RMPs to help mitigating the power exhaust issue in future fusion reactors. Finally, it is seen that the well of the radial electric field on the plasma edge is flattened, effectively reducing the shearing frequency. Although this effect does not have an impact on the nature of the turbulent fluctuations, we note that it is a phenomenon observed in previous experiments and simulations.