Progress towards modeling tokamak boundary plasma turbulence and understanding its role in setting divertor heat flux widths

The high confinement (H-mode) regime^1,2 is currently considered the most promising scenario to achieve magnetic confinement fusion.³ However, the periodic burst of ELMs (edge localized modes)^4,5 in the H-mode plasmas is a serious problem, which may cause excessive erosion and damage to the plasma facing components (PFCs) in the future tokamaks. Extensive work has been done to understand ELM-free operating conditions. Experimentally, there is a continuous edge fluctuation present in the enhanced D_α (EDA) H-mode on Alcator C-Mod, which is named a quasi-coherent mode (QCM), and it is found that the EDA H-mode is naturally stable to ELMs.^6–14 The QCMs have characteristics between coherent and broadband fluctuations as they oscillate around a given frequency but have a wide spectrum. Typically, the measured QCM frequencies are centered around 100 kHz and 1–2 cm^–1 for k_θ on C-Mod. As an important part of the electromagnetic turbulence, the QCMs are thought to be responsible for the enhancement of the perpendicular particle transport. With the existence of QCM in the plasma edge, the periodic bursts of ELMs can be eliminated. Similar fluctuations have been observed on DIII-D¹⁵ and in ELM free H-mode for both EAST¹⁶ and HL-2A.¹⁷ As a combination of the QCM and the broadband turbulence, this kind of electromagnetic turbulence is not only well characterized on multi-machines but also by multiple diagnostics on C-Mod, for example, the two-coil magnetic probe, mirror Langmuir probe,⁹ phase contrast imaging (PCI),¹² gas puffing imaging (GPI),^7,11 polarimetry,¹⁰ and reflectometry.¹⁴ In order to get better understandings of the QCMs, a series of analyses have been compared with the measurements: Initial modeling with BOUT tentatively identified the QCM as a resistive-ballooning/x-point mode.¹² Subsequently, it was pointed out by Myra et al. that the drift-Alfven wave may also be a candidate for the QCM.¹³ The physics of this mode is explored in this paper, using a model in which both the resistive-ballooning mode and the drift-Alfven wave are used.

It has long been known that the competition between the cross field and parallel transport in the scrape-off layer (SOL) determines the heat flux width λ_q. This is a critical parameter needed to characterize the particle and power exhaust on divertor in tokamaks, due to its being inversely proportional to the peak power density. For the ITER divertor design,¹⁸ the heat flux onto the target under H-mode steady state operation should not exceed the limit: 10 MW m^–2. Therefore, the prediction of the heat flux width is important for future devices. Experimentally, the inverse scaling dependence of λ_q with plasma current I_p (equivalently related to the poloidal magnetic field at the outer mid-plane $Bpol,MP$⁠, so called Eich scaling) is observed on multiple machines,^19–23 such as C-Mod, NSTX, DIII-D, JET, and EAST, in H-mode plasmas.²⁴ In order to find the physics behind this scaling, a heuristic drift-based (HD) model with origins in the neoclassical ion orbit width is proposed²⁵ and agrees well with the scaling law. In addition, a theoretical analysis has been done for DIII-D: the two-point theory within the sheath-limited model has a good correlation with the scale length at all collisionalities, while this theory only agrees at high collisionalities within the Spitzer resistivity model.²⁶ Lately, an inverse dependence of the width on the poloidal magnetic field is also obtained with electrostatic gyrokinetic simulations by using the XGC1 code.²⁷ 

A three dimensional electromagnetic fluid turbulence code BOUT++ provides a flexible framework to study the edge physics in tokamaks.^28,29 With implementation of different modules, such as six-field and Gyro-Landau-Fluid (GLF) 3+ 1 models, BOUT++ has been used to study the ELM physics and turbulence transport for both circular and X-point divertor geometry in the last few years.^30–35 They have shed valuable insight on the nature of edge turbulence and the nonlinear saturation mechanisms of turbulence¹² and other edge physics phenomena such as transport. The theoretical and simulation results of the GLF extension of the BOUT++ code are summarized in Ref. 32, which contributes to increasing the physics understanding of ELMs. Recent studies on edge plasma turbulence and divertor heat flux width^33–35 with the six-field two-fluid electromagnetic module show the characteristic of QCM-like modes and achieve good agreements of divertor heat flux widths to the experimental measurements within a factor of 2 for DIII-D, EAST, and C-Mod. Based on the success of BOUT++ simulations mentioned above, the six-field module is exploited in order to improve the understanding of the role of turbulent modes in controlling edge transport and resulting scaling of the heat flux width.

In general, the electromagnetic turbulence on C-Mod has been well characterized in the previous experiments^6–14 and BOUT++ fluid simulations.^12,33 This study extends the previous work to systematically investigate the electromagnetic turbulence and the heat flux width with the improved radial electric field E_r model. Instead of calculating E_r from the ion pressure gradient term (diamagnetic E_r), the E_r is calculated from the plasma transport equations with the sheath potential in the SOL and the plasma density and temperature profiles inside the separatrix from the experiment.

In this paper, the underling physics to be further explored are as follows: (1) how does the radial electric field E_r with sheath potential change the simulation result; (2) the impact of the E_r profile, magnitude or shear, on the electromagnetic turbulence and the heat flux. The first question will be answered by including the sheath potential in the E_r model when conducting the plasma current I_p scan, while the second question will be discussed through an E_r scan. The six-field two-fluid physics model and the transport model for the self-consistent E_r calculation as well as the detailed simulation inputs are described in Sec. II. The possible correction to the simulation result by E_r with the sheath potential in the SOL and the sensitivity study of the electromagnetic turbulence and the heat flux to E_r are introduced in Sec. III. In addition, the turbulence heat flux width is compared to the Eich scaling law and the comparison with the Goldston heuristic drift model²⁵ is given in this section as well. Section IV is the summary.

A set of four EDA H-mode discharges with different plasma current I_p used for the C-Mod simulations. The detailed main plasma parameters are shown in Table I with shot numbers #1100223012, #1100212023, #1100303017, and #1160729008. The first three equilibriums of the C-Mod EDA H-Mode discharges with the I_p from 0.8 to 1.0 MA are the same as the equilibriums used in the previous study in Ref. 33. The discharge parameters are B_T = 5.4 T for attached lower single-null divertor plasma operation (with the B × ∇B drift directed to the active X-point), at conventional aspect ratio of 3.2. Additionally, a really high field discharge #1160729008 of C-Mod EDA H-mode with toroidal magnetic field B_t = 7.8 T and poloidal magnetic field $Bpol,OM=1.11$ T will be introduced in this study, which is very important because some of the operation parameters are more like the future fusion reactors.

The plasma profiles of the four shots, as shown in Fig. 1, used in these simulations are taken from fits of a modified tanh function to experimental data, mapped onto a radial coordinate of normalized poloidal flux ψ. For the purposes of modelling, the profiles for normalized poloidal magnetic flux ψ  >  1.0 are assumed to decrease linearly and extend smoothly into the SOL region. It should be noted that the uncertainty in the separatrix position calculated by EFIT code (Equilibrium FITting code) is of order dψ ∼ 0.01 for normalized poloidal magnetic flux at separatrix ψ_sep = 1, leading to the uncertainty in measured separatrix density and temperature. Clear trends are shown in Fig. 1 for the pedestal pressure profiles: when the plasma current I_p increases, the pedestal height increases and/or the scale length decreases, which is consistent with empirical results in Ref. 36. The major simulation model and inputs used in this study are the same as previous settings in Ref. 33 which are summarized below.

The electromagnetic six-field two-fluid model evolves vorticity ϖ, ion density n_i, ion temperature T_i, electron temperature T_e, perturbed magnetic vector potential $A∥$⁠, and ion parallel velocity $V∥$⁠. This model includes physical phenomena such as Peeling-Ballooning modes, drift waves, ion acoustic waves, thermal conductivities, energy exchange, Hall MHD effects, toroidal compressibility, and electron-ion friction. This model is based on the full Braginskii equations, which are described explicitly as in Refs. 28–31.

The electric field profiles are not measured in the edge region across the separatrix in these discharges. In this paper, the transport module under the BOUT++ framework is used to calculate the electric field E_r from a set of transport equations, including a vorticity equation from a quasi-neutrality constraint. The BOUT++ transport module is separated from the six-field module. The model, equations, and the other simulation settings are described in Ref. 37, where the discharge #1100223012_1149 of C-Mod is used which is the same as the low current I_p = 0.8 MA discharge in this paper. In the transport module, transport coefficients are inverted from experimental profiles inside separatrix. The role of transport module in this work is to set up initial plasma profiles and E_r profile across the separatrix for turbulence simulations.

Based on the experimentally measured plasma density and temperature profiles inside the separatrix, the effective particle and heat diffusivities can be inverted from the set of plasma transport equations [Eqs. (1)–(3) in Ref. 37] with sources only from the inner radial boundary (flux driven sources). With the effective particle and heat diffusivities extended into the SOL and inner flux boundary conditions, we can calculate the plasma density and temperature profiles across the separatrix into the SOL. Inside the separatrix, the plasma profiles are exactly matched to the experimental profiles by design. While in the SOL, the plasma profiles strongly depend on the assumptions how the effective particle and heat diffusivities are extended into the SOL. Here, we assume an initial linear decay of the profiles in the SOL which smoothly connect the experimental profiles inside the separatrix. Then the electric field can be calculated from the vorticity equation [Eq. (5) in Ref. 37] with sheath boundary conditions.³⁷ The SOL plasma density and temperature profiles are also calculated from transport equations. It is worth noting that the cross-field drifts have been ignored. It remains as a research topic for transport calculations how to simultaneously match both plasma profiles and electric field profile in the boundary region across the separatrix. The generated E_r profiles are compared in Fig. 2. It can be seen clearly that the sheath potential modifies E_r in the SOL by comparing the dashed (⁠$Er∇Pi$⁠) to the solid (⁠$Er∇Pi+sh$⁠) curves. Based on the solved E_r, the impact of E_r with sheath potential on the electromagnetic turbulence and the heat flux width will be presented. The E_r profiles used in the simulations are shown in Fig. 3; here we will consider all the four C-Mod discharges with different plasma current.

The simulations are radially global with both edge and SOL regions. The typical simulation domain is shown in Fig. 4. The red curve is the last closed flux surface (LCFS), and the black outline is the vacuum vessel of the C-Mod device. The short black curves represent the equilibrium mesh, typically with a resolution of 260 radial × 64 poloidal grid points, which is generated by using the kinetic EFIT for magnetic equilibria. The background impurity is taken into account in order to satisfy quasi-neutral constraints. The flux-limiting parallel thermal conductivities are used to make up the kinetic correction of parallel transport from collisionless to collisional regimes. The flux limited parallel thermal conduction^38,39 is used by setting the effective parallel thermal conduction coefficients as

The parallel heat conduction coefficients $χ∥e, χ∥i$ are set in the form of Braginskii expressions

The heat flux can be limited by setting the “flux limiting coefficient” $αj$ to a value greater than zero. This calculates the contribution of the fraction of the free streaming to the heat conduction for each species as

where the kinetic effects are introduced through α_j. According to the kinetic calculations,^31,38 the parallel heat flux increases as α_j is increased from the sheath-limited (α_j = 0.05) to free streaming (α_j = 0.8–1.0). Here, the flux-limiting coefficient α_j is setting as 0.3,³¹ if not mentioned.

Dedicated efforts are made to maintain the experimentally measured plasma profiles by adding special sources inside the separatrix, while the SOL plasma profiles are allowed to freely evolve. Figure 5 shows an example of the typical evolution of the ion density profiles at the outside midplane at different times, while the density at the pedestal top is fixed.

Interaction with plasma sheath is a complex problem. Here, the relatively simple sheath boundary conditions are included by assuming that the ion velocity matches the sound speed at the sheath boundary⁴⁰ 

where j represents different ion species. The parallel current through the sheath is described in terms of electron and ion current

the sheath heat fluxes are

where γ_e = 4.8 and γ_i = 4.8 are electron and ion sheath heat transmission factors, respectively. The temperature gradient is set by using the sheath heat flux divided by the parallel conductivity. Ion density and vorticity are set as Neumann boundary condition in front of the target. The sheath boundary condition (SBC) changes the parallel heat flux from being flux limited, as shown in Eq. (2), to be determined by the SBC in front of divertor target.

The boundary plasma turbulence determines the anomalous cross field transport which sets the divertor heat flux widths. In this section, the characteristics of the electromagnetic turbulence are analyzed and compared with the experiment. After that, the comparison of the turbulence heat flux width with the experimental Eich scaling law and the Goldston heuristic drift model is given. In addition, the sensitivity study of the electromagnetic turbulence and the heat flux to E_r is discussed.

Based on the simulation models and inputs described in the previous sections, a dozen nonlinear simulations are performed. The characteristics of the electromagnetic turbulence are analyzed.

A typical frequency spectrum of the ion density fluctuation is shown in Fig. 6 from BOUT++ nonlinear simulations for the C-Mod EDA H-mode discharges, indicating a co-existence of both QCM and broadband turbulence. The broadband turbulence is represented by the red dashed curve, which is a common phenomenon observed on various fusion devices while the peak of the black solid curve between the two black dashed vertical lines is the QCM with a frequency centered around 100 kHz with a spectrum spreading from 80 kHz to 130 kHz. As a part of the electromagnetic turbulence, the QCMs have been preliminarily studied in the previous BOUT++ simulations with a diamagnetic E_r model³³ which shows: (1) the frequency and wave number of QCMs are in agreement with experiment; (2) the QCMs are dominated by resistive-ballooning modes and drift-Alfven wave instabilities; (3) the large SOL turbulence is originated from peak gradients in the pedestal, and not from local instabilities in the SOL; and (4) Blobby type transport which is evidenced by the non-Gaussian probability density function and statistical (skewness-kurtosis) analysis. Consistently, we find a clear blobby structure in the plot of the toroidal slice of the $nĩ$ fluctuation in Fig. 7 in the present research. And the contour plot of fluctuations versus radius and time at the outside midplane in Fig. 8 agrees with the previous results as well, which shows that the blobby turbulence originates in the pedestal peak pressure gradient region inside the magnetic separatrix and nonlinearly spreads across the separatrix into the SOL. As a complementary to the previous study, in this section, the further studies regarding the fluctuation level, locations, and mode spectra of the electromagnetic turbulence, as well as the effect of magnetic flutter on the turbulence spreading are conducted.

The fluctuation level of the main plasma parameters like temperature, density, and magnetic potential is used for the estimation of the turbulence transport experimentally. Langmuir probes are used for the fluctuation measurement of QCMs on C-Mod.⁹ The measured fluctuation level is 45% and 30% for temperature and density, respectively. Figure 8 shows the simulated fluctuation of electron temperature, ion temperature, and ion density at the outside midplane. The relative fluctuation is calculated and normalized to the local mean value for each variable. The simulated values near separatrix are similar to the experimental measurements which are around 50% for both electron and ion temperature, and 45% for the ion density.

The simulated radial and poloidal locations of the QCMs are compared with experiment. It is generally consistent with measurement of probes,⁹ reflectometry,¹⁴ and the previous study in Ref. 33: (1) radially, the plots of fluctuations versus radius and time in Fig. 8 show that the radial location of the mode is spreading across the separatrix, which originates at the peak pressure gradient position in the pedestal, not in the SOL. As shown in Fig. 9, the radial mode structure versus E_r from simulations indicates that the QCMs are spreading around the E_r well minimum near separatrix, which is consistent with the latest research by multiple diagnostics in Ref. 7, although the probe data show that the mode is near the separatrix;⁹ (2) poloidally, the poloidal slice of the nonlinear mode structure in Fig. 7(a) shows that the mode locates at the low field side with the characteristic of resistive-ballooning mode, even though it spreads across the whole cross section, but still dominant outside.

The simulated frequency and wave number spectra are the most important features to compare with experiment. In this section, the spectra of QCM from BOUT++ simulation are compared with PCI data. The PCI measures line integrated density fluctuations along vertical chords ranging from R = 0.63 to 0.76 m. Figure 10 shows an example of frequency versus radial wave number spectrum for QCMs from the PCI diagnostic on C-Mod EDA H-mode discharges, which exhibits clear spectra of QCM with frequency between 120 kHz and 180 kHz and radial wave number k_R around 4–8 cm⁻¹. The plot of the spectra has symmetrical positive and negative values of wave number due to the PCI chords going across the main plasma vertically from the top of the torus to the bottom.

FFT analysis of the electromagnetic fluctuations from BOUT++ simulations⁴¹ has a similar characteristic of the mode spectra as shown in Fig. 11. (The spectra have a polar symmetry in the first and third quadrants with the same sign of frequency and wave number, and so does the mode velocity. The comparison with PCI can just focus on the first quadrant.) The poloidal wave number k_θ in Table II is inferred from k_R by mapping to the outside midplane as illustrated in Fig. 15 in Ref. 42. The formula used in this conversion is briefly introduced here:¹² $kθ(θ1)/kθ(θ2)≈[Bθ(θ2)R22]/[Bθ(θ1)R12]$⁠, where θ is the poloidal angle, $Bθ$ is the poloidal magnetic field, R is the major radius, and the subscripts 1 and 2 represent two different poloidal locations.

The detailed comparison of the simulated mode frequency and $kθ$ with the PCI data is summarized in Table II. Generally, the simulated mode spectra are in agreement with the PCI data with the mode frequencies around 100 kHz, $kθ$ around 1–2 cm⁻¹. Both of the experiment and the simulations show that the QCM frequencies increase with increasing plasma current (with an outlier, $Ip=0.9$ MA), while the wave number remains almost the same. One should note that the simulated mode spectra for the high field discharge (⁠$Ip=1.4$ MA, $Bpol,OM=1.11$ T) are also in agreement with experiment. The frequency is around 169 kHz with corresponding $kθ$ around 2.1 cm^–1, consistent with experimental measurements where the frequency is around 130 kHz and k_θ is around 2 cm^–1. Basically, the spectra almost have no change in comparison to cases with diamagnetic E_r in the previous study.³³ So the E_r in the SOL has a little effect on the QCM spectra. However, the E_r inside the separatrix has a strong effect on the QCM spectra. We find that QCM frequency is proportional to E_r in the sensitivity study in Sec. III D.

One important electromagnetic effect is that the turbulence spreading across the separatrix is enhanced by the additional magnetic flutter in the parallel thermal conduction. Based on the magnetic probe measurements in the C-Mod EDA H-mode discharges, the experimental data show that the electromagnetic components are important part of the QCMs.⁹ The diagnostics show a typical amplitude of radial magnetic field B_r fluctuation to be about 0.3 mT in the QCM mode layer. So it is interesting to know the role of the magnetic flutter on the turbulent transport.

The thermal conduction can be written as Eq. (13) in Ref. 31 

The first term is that the parallel derivatives are only taken along the equilibrium magnetic field lines, where $κ∥$ is the parallel thermal conduction coefficient, the subscript 0 means the equilibrium variables, and j represents different species. The second, third, and fourth terms are generated due to the perturbed magnetic field. This magnetic perturbation is considered in all the other terms in Eqs. (1)–(5) in Ref. 31 through the perturbed magnetic field in the form of $∂∥=(b0+b̃)·∇$ and $b̃=∇A∥×b0/B$⁠, where $b̃$ is the normalized vector of the perturbed magnetic field, B is the total magnetic field, and $A∥$ is the perturbed parallel vector potential. Again, the magnetic flutter here, in this paper, only means the terms in the parallel thermal conduction.

By comparing the fluctuations of electron temperature without magnetic flutter in Fig. 12(b) to that with magnetic flutter in Fig. 8(b), we find that the T_e fluctuation across separatrix is enhanced by the magnetic flutter. However, the magnetic flutter almost has no effect on the ion density [Figs. 8(a) and 12(a)] and ion temperature fluctuations.

Since the electromagnetic turbulence can be well described in the BOUT++ simulations, then here comes the question: How does the electromagnetic fluctuation impact on the divertor heat flux width. So our focus in this section is using the BOUT++ code to study the role of electromagnetic turbulence on scaling of divertor heat flux width. First of all, the big picture is presented according to the schematic diagram in Fig. 7 to help understanding the physics which builds up the divertor heat flux width.

First, based on the section above, we find on the upstream that the turbulence shows the nature of blobby transport [Fig. 7(b)]. It can be seen in Fig. 8 that the electromagnetic turbulence spreads across the separatrix, and the poloidal slice of the nonlinear mode structure and the flow pattern of the heat flux in the plasma edge region in Fig. 7(a) show that the electromagnetic fluctuation causes particle and heat to be turbulently transported radially down their gradients across the separatrix into the SOL. Then in the SOL, the parallel transport will take the heat flow towards the divertor target. Finally, on the divertor footprint, the radial profiles of parallel heat fluxes are mapped to the outside midplane according to the Eich-fitting formula.²⁴ One example of the Eich-fitting of the divertor heat flux footprint width is shown in Fig. 7(c), where the red curve from the fitting matches well with the BOUT++ produced footprint (black dotted curve). The Eich-fitting formula is described in Ref. 24 and is given as follows:

where $s¯=s−s0=(Rsep−R)·fx$⁠, s is the radial path length along the divertor, q_BG is the background heat flux, f_x is the effective flux expansion, S is the divertor power spreading parameter, q₀ is the peak heat flux, s₀ is the radial path length at the strike point, R_sep is the major radius at the separatrix at outer midplane, R is the major radius at outer midplane, and λ_q is the SOL width. The fitting function is a convolution of a Gaussian with an exponential function. Here the exponential captures transport processes in the SOL. In order to compare with the experimental results, the formula we described here is the same fitting method as used by the experimentalists.

The simulated divertor heat flux widths (black dotted curves) for the selected C-Mod discharges are shown in Fig. 13. The radial profile mapped to the outside midplane shows a good agreement with the Eich-fitting formula. Figure 14 shows the results of the Eich-fitting λ_q versus poloidal magnetic field $Bpol,MP$ from BOUT++ simulations in comparison with the experimental Eich scaling law. The red circles are the simulation results with the self-consistently solved E_r (diamagnetic E_r + sheath potential in the SOL) calculated from the BOUT++ Transport module without magnetic drifts.³⁷ The solid curve is the result of a regression of the width, using the Eich-fitting formula to the ITPA multi-machine database, and the dashed curves show the error bars from the experimental regression.²⁴ The detailed comparison of simulated data with experiment is listed in Table III. Here the simulated λ_q and the amplitude of the heat fluxes are in reasonable agreement with C-Mod experimental data within a factor of 2. In Fig. 14, the previous results using the diamagnetic E_r³³ are plotted (green stars). By comparing the red circles and green stars, we find that the simulation results with $Er∇Pi+sh$ (diamagnetic E_r in core, sheath physics in SOL) get closer to the experiment; typically, the divertor heat flux width is suppressed by E_r with sheath potential. In parallel to this paper, the heat flux width studies on DIII-D³⁴ and EAST³⁵ by BOUT++ six-field model have been conducted and the results show that the three devices follow the experimental heat flux width scaling of the inverse dependence on the poloidal magnetic field with an outlier for DIII-D, while for the C-Mod high field discharge the heat flux width scatters in a relatively wide range.

On the theory side: the heuristic drift-based (HD) model²⁵ which has its origins in the neoclassical ion orbit width provides a model of the Eich scaling. With the hypothesis that the magnetic (⁠$∇B$ and curvature B) drifts against the near-sonic parallel flow and the anomalous radial electron transport provides the heat flux across the separatrix, this model gives a width $λq∼qρs$⁠, where ρ_s is the ion gyro-radius, and q is the safety factor. This model provides both a $1/Ip$ scaling and the absolute widths in reasonable agreement with experimental observations. In comparison with BOUT++ simulations, the following Goldston formula²⁵ is used in this paper:

where a is the minor radius, R is the major radius, $Z¯$ is the effective number charge, e is the electron charge, B_p is the poloidal magnetic field, $A¯$ is the effective atomic mass number, m_p is the proton mass, and T_sep is the separatrix temperature.

The corresponding parameters for the C-Mod deuterium plasma cases are: a = 0.22 m, R = 0.68 m, $Z¯$=1, and $A¯$=2. The following parameters are used to estimate the heat flux width based on the Goldston model: The poloidal magnetic field at the outside midplane $Bpol,MP$ (B_p) which are listed in Table I, and the electron temperature at the separatrix $Te,sep$ for I_p = 0.8–1.4 MA. The plasma profiles are plotted as a function of distance from the last closed flux surface (LCFS). However, there are uncertainties in the actual location of the LCFS, as the position reported by EFIT often does not make physical sense from the power balance analysis. The uncertainty in the actual location of the LCFS leads to the uncertainty in the value of $Te,sep$ from the measurements. Here we evaluate the uncertainties in the Goldston model by assuming a range of the separatrix position $dψ∼0.01$ for normalized poloidal magnetic flux ψ with ψ = 1 at the separatrix. Then the uncertainty in the value of the $Te,sep$ is determined and calculated for each cases. The data points of the heat flux width calculated by the HD model are added to the data base in Fig. 15. The blue circles are the data from the Goldston model and the error bars as the short vertical blue lines.

In the HD model, there is no electric drift. So a model including this drift is needed to investigate the electric drift effect. Actually, it is the radial electric field E_r which leads to the uncertainty of the Pfirsch-Schluter flow due to the lack of $Er×BT$ drift. So the improved E_r model described in Sec. II is used in the BOUT++ simulations which can generate a more accurate $Er×BT$ drift in the SOL. On the other hand, the anomalous radial transport can also be changed directly by the sheath potential due to the $Eθ×BT$ drift. To find out the E × B drift effects, full numerical calculations to estimate the anomalous radial transport and the parallel flows are conducted by using the BOUT++ turbulence code. The simulation results are compared with the Goldston HD model in Fig. 15. We find that the turbulence simulations show the similar trend to the HD model and the difference is within a factor of 2. So not only the Goldston model can provide an I_p dependence based on the neoclassical theory, but the turbulence can do the same as well.

By comparing the E_r profiles used in these simulations in Fig. 3, we find that the higher I_p cases have larger E_r (amplitude and/or shear) in general. Figure 16 shows the electron heat transport coefficients in the SOL decrease when the E_r increases (the larger I_p case has larger E_r for these cases) with an outlier (the $Ip=1.4$ MA case). Just like the poloidal ion gyro-radius from the magnetic drifts determines the heat flux width of ions, these results show that the electron divertor heat flux width is determined by turbulent transport and its inverse current I_p scaling comes from the impact of radial electric field E_r on the turbulence. Due to the different equilibriums and experimental profiles and many other different parameters are used in the I_p scan, the fact that E_r suppresses the radial transport in the SOL is not obvious, so it is necessary to study the effect of E_r on the electromagnetic turbulence and the heat flux width in the section below based on an E_r scan. In the scan only the E_r changes and keeps all the other parameters the same.

As we know the magnitude and shear of the electric field E_r always have some strong effects on the turbulence behavior, which may eventually change the overall behavior of the heat flux.^43–47 In this section, the impact of E_r on the electromagnetic turbulence and the divertor heat flux width will be presented via an E_r scan on the high field discharge (⁠$Ip=1.4$ MA). The E_r profiles at the outside midplane for the scan are shown in Fig. 17. Here we will consider four different E_r profiles: (1) simple diamagnetic E_r from the pressure gradient term: $Er∇Pi=1/(niZe)dp/dr$⁠; (2) calculated electric field from transport equations with sheath boundary conditions as described in Sec. II, $Er∇Pi+sh$⁠; (3) half of E_r in second case $0.5×Er∇Pi+sh$⁠, where both the magnitude and shear of E_r are decreased; and (4) twice of E_r in second case $2×Er∇Pi+sh$⁠, where both the magnitude and shear of E_r are increased.

One important finding of the physics behind the electromagnetic turbulence is that the mode frequency is proportional to the radial electric field E_r. The details of the comparison of the mode frequency and poloidal wave number are shown in Fig. 11 and Table IV. It is shown that the mode spectra have changed a lot due to the change of E_r, where the frequency becomes larger when the absolute value of the E_r increases. However, the poloidal wave number k_θ still remains almost the same. The behavior of the frequency versus E_r is consistent with previous BOUT simulations.¹² 

Here we also want to find out the impact of E_r on the divertor heat flux width. Based on the Eich-fitting method, the heat flux widths are calculated. The data points of the heat flux width from the E_r scan are shown in Fig. 14 and Table IV. We find that the heat flux width is suppressed from 1.21 mm to 0.81 mm when E_r increases from half to 2 times. To understand this, the radial distributions of the turbulent transport coefficients are calculated. Figure 18 shows that the radial heat transport in the SOL is inversely proportional to E_r. Generally speaking, the reduced radial electric field E_r leads to the enhanced radial SOL transport and larger heat flux width. So from the perspective of turbulence transport, the E_r plays a critical role to determine the electron heat flux width, which in turn provides an inversely I_p dependence as discussed in Sec. III C.

In summary, the improved BOUT++ model with self-consistent calculation of the radial electric field E_r across the separatrix into the SOL is utilized to understand the underlying physics mechanisms of the electromagnetic turbulence in setting divertor heat flux widths. The BOUT++ simulations show the follow results:

On the upstream, the blobby turbulence originates in the pedestal peak pressure gradient region inside the magnetic separatrix and nonlinearly spreads across the separatrix. The electromagnetic fluctuation provides anomalous radial transport which causes particle and heat to be turbulently transported radially down their gradients across the separatrix into the SOL. The electromagnetic fluctuations show the characteristics of both QCMs and broadband turbulence

• The frequencies of the QCMs are proportional to plasma current I_p, while the poloidal wave numbers are almost not changing;

• Location of the mode is generally consistent with experiment measurements: the GPI experimental data suggest that the QCMs are centered near the E_r well minimum, and the probe data show that the mode is near the separatrix. BOUT++ simulation data show that QCMs are oscillating around the E_r minimum near the separatrix.

• Poloidal mode spectra are in agreement with the Phase Contrast Imaging (PCI) data.

In the SOL, the electromagnetic fluctuations yield the radial transport, and then the parallel transport takes particles and heat flow towards the divertor target. On the divertor footprint, the heat fluxes follow the Eich scaling law.

The impact of the radial electric field E_r on the electromagnetic turbulence and the divertor heat fluxes is demonstrated as follows:

• the E_r changes the turbulent characteristics of the electromagnetic turbulence; from the simulations we find that the frequencies of the modes are proportional to E_r, which is consistent with previous BOUT simulations;¹² 

• the divertor heat flux widths are inversely proportional to E_r: E_r strongly impacts on transport, and the radial heat transport in the SOL is inversely proportional to E_r. Which means that the reduced radial electric field E_r leads to the enhanced radial SOL transport and larger heat flux width.

Comparison with the Goldston HD model, the turbulence simulations show a similar trend to this model.

Not only I_p but also E_r plays an important role on the heat flux width: The self-consistent boundary E_r including sheath potential is needed for the turbulence simulations to get the consistent results of the heat flux width with the experiment. This study suggests a new way to mitigate the heat load on the divertor targets both in current existing tokamaks and in future reactor design. The heat flux width can be broadened by reducing the edge rotation to decrease the edge E_r and/or its shear, therefore increasing the radial turbulent transport. Since this would presumably have a detrimental impact on pedestal confinement and fusion performance, the technical methods and scientific feasibility need to be verified.

The authors would thank Dr. N. Yan and Y. M. Wang for providing the analysis tools, M. V. Umansky, A. Dimits, I. Joseph, T. D. Rognlien, and T. F. Tang for useful discussions. This work was performed under the auspices of the CSC (Grant No. 201506340019) and supported by the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract Nos. DE-AC52-7NA27344, LLNL-JRNL-742019. This material is based upon the work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences. Experimental work was supported by U.S. Department of Energy Agreements DE-FC02-99ER54512, DE-SC0014264, using Alcator C-Mod, a DoE Office of Science user facility. This work was also supported by the National Natural Science Foundation of China under Grant No. 11375191.