Understanding and predicting divertor heat-load width λ_{q} is a critically important problem for an easier and more robust operation of ITER with high fusion gain. Previous predictive simulation data for λ_{q} using the extreme-scale edge gyrokinetic code XGC1 [S. Ku *et al.*, Phys. Plasmas **25**, 056107 (2018)] in the electrostatic limit under attached divertor plasma conditions in three major US tokamaks [C. S. Chang *et al.*, Nucl. Fusion **57**, 116023 (2017)] reproduced the Eich and Goldston attached-divertor formula results [formula #14 in T. Eich *et al.*, Nucl. Fusion **53**, 093031 (2013) and R. J. Goldston, Nucl. Fusion **52**, 013009 (2012)] and furthermore predicted over six times wider λ_{q} than the maximal Eich and Goldston formula predictions on a full-power (Q = 10) scenario ITER plasma. After adding data from further predictive simulations on a highest current JET and highest-current Alcator C-Mod, a machine learning program is used to identify a new scaling formula for λ_{q} as a simple modification to the Eich formula #14, which reproduces the Eich scaling formula for the present tokamaks and which embraces the wide λ_{q}^{XGC} for the full-current Q = 10 ITER plasma. The new formula is then successfully tested on three more ITER plasmas: two corresponding to long burning scenarios with Q = 5 and one at low plasma current to be explored in the initial phases of ITER operation. The new physics that gives rise to the wider λ_{q}^{XGC} is identified to be the weakly collisional, trapped-electron-mode turbulence across the magnetic separatrix, which is known to be an efficient transporter of the electron heat and mass. Electromagnetic turbulence and high-collisionality effects on the new formula are the next study topics for XGC1.

## I. INTRODUCTION

A challenge for ITER operation is the ability of the divertor plates to withstand the steady plasma exhaust heat that will be deposited on the surface along a narrow toroidal strip. A simple data-based regression using macroscopic parameters from attached-divertor experiments on all the present devices (formula #14 in Refs. 1 and 2) shows that the heat-flux width follows a scaling 1/B^{γ}_{pol,MP} where B_{pol,MP} is the magnitude of the poloidal magnetic field on the outboard midplane separatrix surface and γ = 1.19. References 1 and 2 also present other possible regression formulas that are valid for certain chosen device sets. There has also been a heuristic model by Goldston^{3} based on the neoclassical orbit-driven ion losses for weakly collisional edge plasma, which resulted in a similar result to that in Refs. 1 and 2. For ITER H-mode operation at I_{P} = 15 MA with q_{95} = 3, these regression and heuristic formulas yield at most λ_{q} ≈ 1 mm for the divertor heat-flux width measured at outboard midplane after being mapped from the divertor plates along the magnetic field lines. Here, λ_{q} is defined in the following fitting formula:^{1,2}

where R_{mp} is the major radius along the outboard midplane, R_{mp,sep} is R_{mp} on the outboard separatrix surface, *h*(R_{mp}-R_{mp,sep}) is the input function to the fitting formula (namely, the divertor heat-flux profile data at outboard midplane after being mapped from the divertor plates along the magnetic field lines), *h*_{0} is the peak value of *h*, S is a spreading parameter which makes the heat flux profile deviate from an exponential decay, *Erfc* is the complementary error function, and *h*_{BG} is the background heat-flux.

For this range of λ_{q} in ITER, the peak divertor power fluxes in attached divertor conditions are beyond the design limits of the stationary heat loads of the ITER divertor target, thus requiring the divertor operation in deeply semi-detached or detached conditions in which the plasma power is dissipated over a larger area by atomic radiation from hydrogenic-isotope atoms and impurities in the divertor chamber. The operational range for such a deeply semi-detached or detached divertor operation decreases with smaller λ_{q}, and is restricted to very high plasma separatrix densities and radiative fractions, requiring n_{sep}/n_{GW} > 0.6 for λ_{q} ≈ 1 mm,^{4} where n_{GW} is the critical plasma density inside the pedestal top above which the plasma tends to have a deteriorated confinement and even disrupt.^{5} This raises concerns regarding their compatibility with the good H-mode energy confinement required to achieve Q = 10 operation in ITER and the increased probability for plasma disruption. In addition, such a small λ_{q} poses additional challenges for the control and sustainment of the semi-detached or detached divertor conditions since the power fluxes during transient re-attachment may significantly exceed the stationary heat flux design limits of the ITER divertor.

However, it is questionable if such a simple extrapolation from present experiments is valid as there may be differences in the fundamental edge physics between ITER and the present devices. Any extrapolation from present experiments to ITER may need to be on a more fundamental, first-principles-based kinetic physics. This was the purpose of the gyrokinetic study in Ref. 6, utilizing the edge gyrokinetic particle-in-cell code XGC1.^{7}

First, the heat-flux width (λ_{q}^{XGC}) predictions from the XGC1 gyrokinetic model reproduced the carefully chosen representative experimental data from three U.S. tokamaks within the regression error bar of the Eich scaling study.^{1,2} Total-f gyrokinetic simulations were performed until an approximate gyrokinetic power balance was achieved in XGC1 between separatrix surface and divertor plates at the level of core heating power. A minor adjustment by the total-f XGC1 code of the experimentally measured or model profiles across the magnetic separatrix was made before approximate power balance was achieved. Second, the same XGC1 code was used to predict the heat-flux width on the full-current (15 MA) Q = 10 ITER plasma, with the caveat that the initial ITER plasma input to XGC1 from the reduced model code JINTRAC^{8} may not be in agreement with the total-f gyrokinetic code XGC1. As a matter of fact, a significant adjustment from the initial JINTRAC edge plasma happened before XGC1's achievement of an approximate gyrokinetic power-balance between the power-crossing at separatrix and the heat load at divertor plates at the level of heat-source at the burning core.

Actual experimental plasma profiles that satisfy the Grad–Shafranov equilibrium relation required only a minor adjustment before a gyrokinetic quasi-equilibrium is reached in the total-f XGC1. However, the reduced-model predicted plasma profiles (such as those for ITER) often require a significant adjustment, in accordance with the radial plasma transport fluxes, before a gyrokinetic quasi-equilibrium is reached consistently with the magnetic equilibrium, as shown in Ref. 6 and later in the present report. There is an underlying assumption here that a deterministic gyrokinetic plasma profile state exists in accordance with external constraints when starting from different but nearby reduced-model predicted plasma profiles as long as the external heat source profiles, the wall recycling coefficients, and the boundary conditions are identical. The most interesting finding from the study was that the same gyrokinetic code that reproduced experimental λ_{q} in the present tokamak plasmas predicted that λ_{q} in the full-current ITER model-plasma in attached divertor condition would be over six-times wider than what could be maximally extrapolated from the various Eich scaling formulas and the Goldston formula. More details can be found in Ref. 6.

Understanding the physics cause behind such a significant broadening of λ_{q}^{XGC} in the full-current ITER Q = 10 edge plasma has remained as a critical research issue for the XGC group. A subsequent data analysis showed that the edge turbulence pattern across the magnetic separatrix changes from the space-time isolated “blobs”^{9} in all the present tokamaks to radially extended and connected “streamers”^{10} in the full-current ITER Q = 10 scenario that are typically seen in the ion-scale microturbulence such as the ion-temperature-gradient (ITG) driven turbulence and the trapped-electron-mode (TEM) turbulence. This gives us a strong hint that there is a fundamental physics change between the present tokamak edge plasma and the full-current ITER edge plasma in the XGC1 electrostatic simulation.

Another strong clue arises from the recent high-current experiments on Alcator C-Mod tokamak.^{11} With the poloidal magnetic field strength as strong as that of the ITER full-current Q = 10 plasma, experimental λ_{q}^{Exp} values in the Alcator C-Mod experiments still follow the Eich scaling. An XGC1 simulation has been performed on one of these high-current C-Mod plasmas and confirmed that the gyrokinetic λ_{q}^{XGC} from XGC1 also follows the Eich scaling. This yields double-valued solutions for λ_{q}^{XGC} between the high-current C-Mod plasma and the full-current ITER plasma if B_{pol,MP} (or the macroscopic parameters used in Eich *et al.*) is the sole independent parameter, indicating the existence of other hidden parameter(s).

It is the purpose of the present paper to conduct a systematic search for the hidden parameter(s) and the corresponding new physics by utilizing deeper data analyses, high-fidelity physics knowledge, and a convenient machine-learning tool in search of an improved λ_{q}^{XGC} scaling formula that can encompass not only all the present experimental results, but also the gyrokinetic predictions for the full-current (15 MA) ITER result. Three more simulations are performed on different ITER model plasmas to successfully test the new scaling formula. The present study opens up doors to several deeper edge-physics research topics, as will be pointed out in Secs. V–VII. Study of the electromagnetic and high-collisionality effects on λ_{q}^{XGC} is left for future work.

We note that there is recent empirical modeling showing some widening of the near-scrape-off layer (near-SOL) upstream power-width due to a high collisionality effect^{12} in present tokamaks that could represent the relative importance of the interchange effect on drift-wave turbulence,^{13,14} aiming for semi-detached or detached divertor plasmas. In this work, we confine our study to the low recycling, attached divertor plasma conditions and do not attempt to study the high collisionality effect of Ref. 12. There is a BOUT++ fluid turbulence simulation result^{15} which shows broadening of λ_{q} in the 15 MA Q = 10 ITER plasma. Since fluid modeling does not contain the kinetic physics that are essential in the present work, such as the finite ion orbit width and trapped electron modes, we do not attempt to compare the present work with Ref. 15. There is also a SOLPS-ITER, which is a new version of the SOLPS plasma boundary code package described in Ref. 16, transport modeling of the 15 MA ITER discharge, with an arbitrarily chosen radial diffusion coefficient that shows an anomalous electron thermal diffusivity at 1 m^{2}/s in the SOL could broaden λ_{q} to 3–4 mm.^{16}

The paper is organized as follows: In Sec. II, for the sake of completeness, we briefly summarize the previous results from Ref. 6. In Sec. III, we present new simulation results that answer some questions left by Ref. 6. In Sec. IV, we utilize a machine learning program to find a new scaling formula for λ_{q}^{XGC}. In Sec. VI, we test the new predictive formula by performing simulations on different ITER model plasmas. In Sec. V, we describe the new physics understanding in relation to the new scaling formula. We present summary and discussion in Sec. VII.

## II. A BRIEF SUMMARY OF THE PREVIOUS XGC1 SIMULATION RESULTS

In this section, for the sake of completeness, we briefly summarize the previous XGC1 simulation results reported in Ref. 6 as the basis for the discussions presented in this paper. Table I shows the seven simulation cases studied in Ref. 6, chosen in collaboration with three major U.S. tokamaks and the ITER Organization. The discharges were selected to cover a wide range of the then experimentally available B_{pol,MP}, the poloidal magnetic field magnitude at the outboard midplane on the magnetic separatrix surface. Discharges from three U.S. tokamaks were part of the discharge set used in the regression analysis in Eich *et al.*^{1,2} In all the discharges, the ion magnetic drift direction is toward the single magnetic X-point, and the inter-ELM (edge localized modes) divertor plasma is in the attached regime. It should be noted here that at that time at which the work in Ref. 6 was being conducted, the highest-field C-Mod experiments^{11} with B_{pol,MP} reaching the ITER full-current case did not exist.

Shot . | Time (ms) . | B_{T} (T)
. | I_{P} (MA)
. | B_{pol,MP} (T)
. |
---|---|---|---|---|

NSTX 132368 | 360 | 0.4 | 0.7 | 0.20 |

DIII-D 144977 | 3103 | 2.1 | 1.0 | 0.30 |

DIII-D 144981 | 3175 | 2.1 | 1.5 | 0.42 |

C-Mod 1100223026 | 1091 | 5.4 | 0.5 | 0.50 |

C-Mod 1100223012 | 1149 | 5.4 | 0.8 | 0.67 |

C-Mod 1100223023 | 1236 | 5.4 | 0.9 | 0.81 |

ITER full-current scenario | Flattop | 5.3 | 15 | 1.25 |

Shot . | Time (ms) . | B_{T} (T)
. | I_{P} (MA)
. | B_{pol,MP} (T)
. |
---|---|---|---|---|

NSTX 132368 | 360 | 0.4 | 0.7 | 0.20 |

DIII-D 144977 | 3103 | 2.1 | 1.0 | 0.30 |

DIII-D 144981 | 3175 | 2.1 | 1.5 | 0.42 |

C-Mod 1100223026 | 1091 | 5.4 | 0.5 | 0.50 |

C-Mod 1100223012 | 1149 | 5.4 | 0.8 | 0.67 |

C-Mod 1100223023 | 1236 | 5.4 | 0.9 | 0.81 |

ITER full-current scenario | Flattop | 5.3 | 15 | 1.25 |

Figure 1, without counting the 4.5 MA JET and the 1.4 MA C-Mod points that will be used in Sec. III, shows the simulation results for λ_{q} from XGC1 in comparison with the experimental results λ_{q}^{Exp} of Refs. 1 and 2, with the symbols from XGC1 improved from Fig. 16 of Ref. 6 to resemble the corresponding experimental symbol shapes. The inaccuracy in the ITER λ_{q}^{XGC} = 5.9 mm point position in Fig. 16 of Ref. 6 is corrected in Fig. 1. As can be seen from all the open symbols, the XGC1 predictions for the present tokamaks agree well with the Eich scaling for λ_{q} from formula #14 in Ref. 1 [hereafter referred to as λ_{q}^{Eich(14)}], represented by the solid line, together with the regression error represented by the two dashed lines. Here, we use the Eich formula #14 [λ_{q}^{Eich(14)} $\u2248$ 0.63B_{pol}^{−1.19 }mm] because it contains data from all the tokamaks. Turbulence across the magnetic separatrix and in SOL was always of blob type in the present devices in the XGC1 simulations, as measured in some experiments. A blob is a magnetic-field-aligned intermittent plasma structure which is considerably denser than the surrounding background plasma and highly isolated in the two directions perpendicular to the equilibrium magnetic field.^{9} However, the XGC1-predicted λ_{q}^{XGC} in the full-current Q = 10 ITER scenario plasma (15 MA, B_{pol,mp} = 1.21 T) is about six-times greater than what could be maximally predicted from various Eich formulas/Goldston formula or about 12 times greater than λ_{q}^{Eich(14)}.

In Ref. 6, a possibility for this large deviation for the full-current ITER was hypothesized to be from a much longer radial correlation length of the edge turbulence across the separatrix surface caused by the low neoclassical E × B shearing rate in the ITER full-current Q = 10 plasma. In Sec. VI, it will be shown that the turbulence with much longer radial-correlation length has a streamer structure, which is usually observed in ITG and TEM driven turbulence. This hypothesis was drawn from the fact that the neoclassical physics strength, and thus the neoclassical E × B flow shearing rate, becomes weaker as *ρ _{i,pol}/a* becomes smaller, where “

*ρ*is the poloidal ion Larmor radius at the outboard midplane separatrix point and “

_{i,pol}”*a”*is the plasma minor radius. In the full-current ITER,

*ρ*is an order of magnitude smaller than that in the highest-current C-Mod plasma. In the present tokamak devices, XGC1 found that the divertor heat-flux width physics is dominated by the ion neoclassical drift motions,

_{i,pol}/a^{6}in agreement with Ref. 3, in spite of the existence of large-amplitude blobby turbulence across the separatrix and in the SOL.

A quick demonstration of the neoclassical E × B dependence on ion banana width can be given by using the standard neoclassical radial force balance equation in the closed field-line region^{17}

where ⟨u_{ǁ}⟩ is the flux-surface averaged parallel fluid-flow velocity and k is a collisionality-dependent parameter, that is, 1.17 when ions are in the banana regime^{17} (ions near the magnetic separatrix in the full-current ITER edge are in this regime, but the value k = 1.17 may not be accurate in the edge plasma). Neglecting, for the sake of a simpler argument, the temperature gradient term, whose gradient and coefficient are significantly smaller than the density gradient term for k = 1.17, we can simplify and rearrange Eq. (1) to

where u_{E} = E_{r}/B is the E × B flow speed, v_{i} is the ion thermal speed, v_{i,pol} is the poloidal component of the parallel thermal speed, ρ_{i,pol} is the ion gyroradius in the poloidal magnetic field, and αa is an expression for the density gradient scale length expressed in terms of a parameter α and the plasma minor radius. For H-mode pedestals in the conventional aspect-ratio tokamak edge, α does not vary widely but stays around ∼0.05. It can be easily noticed from Eq. (2) that the plasma gradient term ρ_{i,pol}/αa is the driver for the radial electric field, or equivalently for the E × B flow, that is, mostly in the poloidal direction. As the device size becomes greater relative to the ion poloidal gyroradius, u_{E} becomes smaller in proportion. For the full current 15 MA ITER, ρ_{i,pol}/a is about six times smaller than the 1.5 MA DIII-D case of Table I. The physically meaningful E × B shearing parameter is $\gamma \Epsilon $_{/}$\omega *$ = du_{E}/($\omega *$ dr) which scales as, using Eq. (2) and the relations $\omega *$∼kv_{i}ρ_{i}/αa, du_{E}/dr ∼ u_{E}/L_{E}, and assuming kρ_{i,pol} ∼ 1

This relationship shows that $\gamma \Epsilon /\omega *$ decreases with ρ_{i,pol}/L_{E} ($\u221d$ρ_{ip}/a). Value of the XGC-found $\gamma \Epsilon /\omega *$ for representative device cases will be shown in Sec. VI, where detailed physics is discussed.

## III. NEW XGC1 SIMULATIONS

The XGC family codes are equipped with a built-in Monte Carlo neutral particle transport capability using ionization and charge exchange cross sections for neutral–plasma interaction. A recycling coefficient R = 0.99 is used for the divertor heat-load width simulations presented here for generation of neutral marker-particles at Frank–Condon energy (3 eV) in front of the material wall wherever the ions are absorbed. For a more detailed introduction, we refer the reader to Ref. 18. In addition to the built-in Monte Carlo neutral particle transport routine, XGC family codes can utilize the DEGAS2 Monte Carlo neutral particle code as a subroutine, which can start the neutral particle recycling process from molecular neutral birth, with volumetric and surface recombination. The latter features are not utilized in the present simulations; hence, our study is limited to the attached, low-recycling divertor regime. We also use a simple cooling profile in the divertor chamber to keep the electron temperature on the outboard separatrix surface close to the input value.

The first new XGC1 simulation is to test an existing experimental plasma, that is, closest to the full-current ITER in both the B_{pol,MP} value and the physical size in deuterium plasma. For this purpose, a JET 4.5 MA discharge^{19} is chosen that has the highest B_{pol,MP} (0.89 T) at the time of simulation (unfortunately, an experimental λ_{q}^{Exp} measurement does not exist on JET at this high value B_{pol,MP}). To be more specific, B_{pol,MP} for this JET plasma is only 36% lower than the full-current ITER plasma, and its linear size is a factor of $\u2248$ 2 smaller than ITER. In this plasma, XGC1 finds λ_{q}^{XGC} of about 0.64 mm, which is within the regression error bar from the Eich(14) value λ_{q}^{Eich} $\u2245$ 0.72 mm (open red circle in Fig. 1). Thus, XGC1 indicates that there may be either a bifurcation of λ_{q}^{XGC} between B_{pol,MP} = 0.89 T of JET and 1.21 T of ITER, or there is something other than the value of B_{pol,MP} which sets the full-current ITER case apart from the present experimental scaling.

Shortly after the JET simulation described above was performed, experiments at C-Mod raised B_{pol,MP} values up to 1.3 T (Ref. 11) which somewhat exceeds the full-current ITER value and found that the experimental λ_{q}^{Exp} still follows λ_{q}^{Eich} approximately. This was an excellent comparison case to be studied by XGC1. Accordingly, we chose the C-Mod discharge #1160930033 with 1.4 MA of plasma current and B_{pol,MP} = 1.11 T. At this high value of B_{pol,MP}, though, we find λ_{q}^{XGC}$\u2009\u2245$ 0.38 mm (see Fig. 2 and the open black star symbol at the far-right bottom of Fig. 1), which is even somewhat smaller than λ_{q}^{Eich(14)} = 0.56 mm. As a result, XGC1's solution becomes double valued around the maximal C-Mod B_{pol,MP} values if B_{pol,MP} is used as the sole parameter and suggests the existence of hidden parameter(s) that was missed in Eich's regression process.

## IV. A SIMULATION-ANCHORED, PREDICTIVE MACHINE LEARNING STUDY

In this section, we use a supervised machine learning program in search of the possible hidden parameter(s). A machine learning program is basically a systematic interpolation and regression technique utilizing mathematical tools. A machine-learning program can yield answers much more rapidly and systematically than human interaction with ordinary spreadsheets can. Any presently available dataset forms an underdetermined system, which is only a subset of all the possible datasets and which may not be good for extrapolation into a new regime where the governing physics phenomena may be different. An extrapolation path from the present data knowledge alone could lead us to a completely wrong direction. However, if a first-principles model can be used to study the new regime and make predictions in accordance with the new governing physics, the simulation results can “anchor” the machine learning into the new physics direction, at least as far as the simulation correctness in the specific target regime is concerned. The “anchoring” high-fidelity simulation points do not have to be many to lead the machine learning prediction into the intended direction, but the more the better for accuracy. Of course, the accuracy of the simulation-anchored predictive machine learning will only be as good as the anchoring high-fidelity model accuracy, which will improve as the computational power increases (or a high-fidelity analytic model). We caution here that the simulation must be well-validated on the present experimental data before adding the anchoring data. The extrapolated predictions must also be validated continuously against new experiments when available.

In this section, we use this “anchored machine learning” concept to search for a predictive analytic scaling formula by combining the experimental and predictive-simulation datasets for the divertor heat-flux width λ_{q}. We use the symbol **D**^{E} to represent a set of λ_{q}^{Exp} data found from the present laboratory experimental measurements, **D**^{SE} for a set of λ_{q} data found through high-fidelity simulation of the existing experiments, and **D**^{SF} for a set of λ_{q} data found through high-fidelity simulation of future experiments. We use **M** to denote the machine-learning operation, **F**^{E} for the modeling formula found by the operation **M** on the present experimental dataset **D**^{E}, **F**^{SE} for the modeling formula found by the operation **M** on **D**^{SE}, and **F**^{P} for the predictive modeling formula found by the operation **M** on all the datasets including **D**^{E}, **D**^{SE}, and **D**^{SF}. **D**^{E} and **D**^{SE} do not need to have one-to-one correspondents.

For the validated high-fidelity simulations, we assume **F**^{E} ≈ **F**^{SE} as a pre-requisite condition, which is satisfied by XGC1 as shown in Secs. II and III. Thus, we have **M**(**D**^{E}) $\u2192$ **F**^{E} and **M**(**D**^{SE}) $\u2192$ **F**^{E}, with some allowance for error. We can then write down the following relations:

Here, **F**^{P}[$\u2283$**F**^{E}] means that the machine-learned formula **F**^{P} reduces to **F**^{E} in the present-day experimental space. In other words, using predictions from simulation on the unexplored future experiments, the simulation-anchored machine-learning operation can be made to possess the predictive capability **F**^{P}, within the simulation accuracy, by combining **D**^{E} and **D**^{SE} with **D**^{SF}.

To achieve this goal, we use an artificial intelligence (AI)-based modeling engine Eureqa.^{20,21} Eureqa uses supervised machine learning techniques to conduct an evolutionary model search to find the best combination of the user-specified mathematical building blocks that fit labeled training data, not only equation parameters, but also the form of the symbolic equation which best fits the data.^{22} Starting with a series of random expressions, the algorithm combines the best-fitting expressions with each other until it gradually finds formulas which fit the data. Eureqa also applies a penalty in proportion to the complexity of the formula so as to prevent overfitting. While trial-and-error single fits could be performed using different forms of equations on combinations of parameters, using symbolic regression frees us from specifying the form of equations to fit the data, resulting in more generic equations.

Our attempt is to find a new predictive scaling formula **F**^{P} of Eq. (4). We present the result first: Fig. 3 depicts the simplest **F**^{P} result from Eureka, as to be elaborated soon later in this section. Figure 3 contains the selected experimental dataset **D**^{E} from NSTX, DIII-D, and C-Mod (marked with + symbols) as presented in Sec. II, and the corresponding simulation dataset **D**^{SE}. The purely predictive 4.5 MA JET and 15 MA ITER simulations, for which experimental measurements do not exist, are also contained Fig. 3. We have normalized all the λ_{q} values in **D** of Eqs. (4) and (5) to the Eich scaling formula #14, λ_{q}^{Eich(14)} = 0.63 B_{pol,MP}^{−1.19}. The simple extrapolation to the future experiments from the present-day experimental dataset is represented by the solid black horizontal line.

Observables in tokamak plasmas are functions of many variables and the machine learning can be a many-variable operation. Eich *et al.* used the nine well-known macroscopic variables for a thorough data regression,^{1,2} which spans the macroscopic plasma-operation space rather completely: *B*_{tor} (the toroidal magnetic field strength), *B*_{pol,MP}, *q*_{95} (the safety factor at the 95% poloidal-flux surface), *P*_{SOL} (the power flow from core into the SOL), *R*_{geo} (the geometric major radius), *a* (the plasma minor radius), *I*_{p} (the plasma current), and *n*/*n*_{GW} (the density ratio to the Greenwald density). Multiple possible formulations are found from data regression in Refs. 1 and 2 depending on the combination of the target tokamaks, but the main dependence of the divertor heat-flux width is found to be on *B*_{pol,MP} by targeting all the present tokamaks, denoted here as the Eich regression number #14, with the squared correlation coefficient being R^{2} = 0.86. Our machine learning operation utilizes λ_{q}^{Eich(14)} as the normalization factor.

We note that Refs. 1 and 2 did not consider microscopic kinetic parameters. Among the microscopic kinetic parameters, there is a dimensionless quantity that could be as important as the macroscopic parameters: the ratio between the ion banana width to the device size^{6,23} as elaborated at the end of Sec. II. The ratio between the ion banana width and the machine size determines the strength of neoclassical physics [see Eq. (2)], including the important background E_{r} × B-flow shearing rate [see Eq. (3)] which controls plasma turbulence.^{24} Plasma turbulence could then affect the cross field spread of the divertor heat-load (characterized by λ_{q}). For this reason, we introduce a new parameter “*a*/*ρ*_{i,pol}” to be used for a physics-based feature-engineering in the supervised machine-learning in Eureqa. Comparison of the normalized E × B-flow shearing rate for example tokamaks that have different *a*/*ρ*_{i,pol} values will be presented in Sec. VI.

Our first try in the present work is to accept the regression result of Refs. 1 and 2, thus accept that there is little dependence of λ_{q} on all other macroscopic parameters, and utilize only two parameters in the machine learning program Eureqa: B_{pol,MP} inherited from Refs. 1 and 2 and the kinetic parameter a/ρ_{i,pol}. If this simplified approach does not work satisfying our three conditions—to resolve the double valued solution issue, to agree with the well-validated λ_{q}^{Eich} formula for the present attached divertor experiments, i.e., **F**^{P}[$\u2283$**F**^{E}], and to encompass the full-current ITER Q = 10 result—then we will have to ignore the work done in Refs. 1 and 2 and perform a many variable machine learning study from scratch.

Application of the data sets **D**^{E} $\u222a$ **D**^{SE} $\u222a$ **D**^{SF} to Eureqa then gave us numerous possible predictive modeling formulas, most of which turn out to be some complicated and meaningless functional combinations of the input parameters B_{pol,MP} and a/ρ_{i,pol}. Three physics-based search-formulas are given to Eureqa to shorten the search time to one hour on a MacBook Pro equipped with a 2.6 GHz Intel Core i7 four-core processor

Among the simulation-anchored formulas found by Eureqa

is the simplest and lowest order expression for the heat flux width λ_{q}^{ML} derived by this machine learning approach with a reasonably low mean square error (RMS error = RMSE = 18.7%). Equation (6) is depicted in Fig. 3 using the dashed purple curve. A lower-order formula could not be picked because the mean square error jumps to above 50%. The formula agrees with λ_{q}^{XGC} for the full-current ITER plasma and reproduces λ_{q}^{Eich(14)} for all the present-day tokamak data. The predictive simulation on the 4.5 MA JET plasma (for which the experimental data does not yet exist) contributes valuably to the fourth power law in the B_{pol,MP} a/ρ_{i,pol} dependence. Notice here that in Fig. 1, the right-most data point used for the XGC1 simulation is from the high field C-Mod. In Fig. 3, however, the right-most data point became the JET simulation point indicating that the highest-field JET case is the closest present tokamak device to the full-current 15 MA ITER as far as λ_{q} is concerned in this parameter space.

Other example formulas that utilize combinations of B_{pol,MP} and a/ρ_{i,pol}, and that yield low RMSE include

All these formulas yield fitting curves that have similar levels of RMSE to Eq. (6), matching the λ_{q} values for the existing tokamaks and the “anchored” full-current ITER as well as Eq. (6) does. However, they have higher order and/or more complicated parameter dependencies, which could make the fitting curve behave differently in the gap region between the present tokamaks and the full-current ITER. In Sec. V, we test Eqs. (6)–(10) by performing XGC simulations on three more ITER model plasmas. The results do not suggest that we should switch Eq. (6) to a more complicated formula. Besides Eqs. (6)–(10), there are other highly complicated and non-smooth formulas Eureqa has produced that try to fit details of the noisy data with much lower mean-squared error (as low as RMSE ∼ 4.5%). However, theses formulas do not reproduce the smooth Eich experimental formula and do not satisfy the requirement to reproduce the Eich regression #14 formula.

A schematic diagram for the workflow used to find the above machine-learned formulas is depicted in Fig. 4, showing the inputs (labeled experimental and simulation data for λ_{q}, B_{pol,MP}, a/ρ_{i}; mathematical operations; and variables), the evolutionary model search process in Eureqa, and the resulting λ_{q}^{ML} formulas (only one of them is shown).

## V. TEST OF THE NEW FORMULA

The new ML-found formula is tested on three different ITER model plasmas: (i) the first H-mode plasma to be explored in the initial phases of ITER operation at 5 MA,^{25} (ii) an H-mode hybrid plasma at 12.5 MA providing long pulse operation with fusion yield Q = 5,^{26} and (iii) an H-mode steady-state plasma at 10 MA providing steady-state operation with Q = 5.^{27} These three ITER model plasmas have distinctively different values of the kinetic parameter *a*/*ρ*_{i,pol} at the outboard midplane edge. The 5 MA plasma has *a*/*ρ*_{i,pol} that is well within the present tokamak range, but its physical size is the same as the full-current ITER plasma; the 12.5 MA hybrid plasma has *a*/*ρ*_{i,pol} slightly above the 15 MA plasma and is thus a good test problem to confirm/refute the large *a*/*ρ*_{i,pol} effect found on the 15 MA plasma; and the 10 MA steady-state plasma has *a*/*ρ*_{i,pol} deep in the gap region between the high-field JET plasma and the 15 MA ITER plasma. As for the original 15 MA Q = 10 ITER discharge, all the new ITER points assume deuteron plasma only and do not include impurity species, but with realistic electron mass. For a visual introduction, results from three new cases are depicted in Fig. 5, as additions to Fig. 3, before being described below in more detail.

We note here that an extension of the high current (4.5 MA) JET plasma that is modeled toward the *B*_{pol,MP} *a*/*ρ*_{i,pol} value of the 15 MA ITER discharge could have been an option instead of the 10 MA ITER case. The plasma equilibrium has to be made up in both cases, which would certainly not be in gyrokinetic equilibrium and must be evolved significantly by XGC1 before power balance between the separatrix and the divertor plates is reached. We choose the 10 MA ITER case here because of the relevance of the 10 MA ITER H-mode scenario for steady-state demonstration at Q = 5. Our simulation can be taken as a gyrokinetic base for predictions of a future real experiment that is planned to be executed and that can be compared with future SOLPS-ITER simulations for these plasmas. A JET experiment at much higher plasma current than 4.5 MA in the present divertor geometry is beyond the capabilities of the device and thus cannot be realized (nor will it be simulated by fluid codes).

5 MA ITER case

After the previous XGC publication of the significantly enhanced divertor heat-flux width in the ITER full-current scenario plasma,^{6} a question naturally arose if the enhancement could simply be from the pure size-effect: ITER is about three-times as large as DIII-D and nine-times as large as Alcator C-Mod in linear size, with its plasma volume approximately 3^{3}- and 9^{3}-times greater. The first H-mode plasma scenario that will be explored in the initial ITER experimental phases with I_{p} = 5 MA^{25} is an excellent case to answer this question: It has B_{pol,MP} = 0.43 T, similar to a high-field DIII-D plasma and a low-field C-Mod plasma (see Table I and Fig. 1), while the plasma size essentially the same as the full-current ITER. The a/ρ_{i,pol} value of 201 is also similar to a typical JET plasma value, with our new parameter B_{pol,MP} a/ρ_{i,pol} for 5 MA ITER falling well within the present device range (see Fig. 5). For a quantitative comparison, the B_{pol,MP} a/ρ_{i,pol} value for the 5 MA ITER case is as small as 87, with B_{pol,MP} a/ρ_{i,pol} for all the present tokamak experiments falling between about 10 and 200. The test XGC1 simulation finds λ_{q}^{XGC} = 2.2 mm, which satisfies the Eich formula value λ_{q}^{Eich(14)} = 1.7 mm approximately within the regression error bar. This result thus excludes the pure size effect from the possible cause for the large λ_{q}^{XGC} found for the full-current 15 MA Q = 10 ITER plasma.

- (ii)
12.5 MA Q = 5 long pulse ITER hybrid scenario case

The 12.5 MA ITER hybrid scenario plasma with *B*_{tor} = 5.3 T and fusion gain of Q = 5 (Ref. 26) is an interesting case. Its toroidal magnetic field strength at the machine axis B_{tor} = 5.3 T is the same as the full current 15 MA case. However, because of the stronger Shafranov shift due to the higher beta and a somewhat smaller major radius of outer-midplane separatrix, the value of *B*_{pol,MP} (=1.22 T) for the 12.5 MA case is about the same as that (1.21 T) in the 15 MA discharge. Due to the smaller ion temperature at the edge (we use plasma values at ψ_{N} = 0.99), the new parameter *B*_{pol,MP} *a*/*ρ*_{i,pol} for the 12.5 MA case is actually slightly greater than the 15 MA case (592 T vs 572 T). This is an interesting case that may be at odds with conventional ITER H-mode plasmas between 5 and 15 MA (with similar beta and H98 = 1) in the *a*/*ρ*_{i,pol} kinetic parameter space, but an excellent second case for testing the broadening of λ_{q}^{XGC} by the large *B*_{pol,MP} *a*/*ρ*_{i,pol} effect. A peculiarity of this plasma scenario will appear again in the discussion on the in–out asymmetry of the divertor power load in Sec. VI. Our simulation shows that λ_{q}^{XGC} $\u2248$ 6.9 mm for this 12.5 MA ITER model plasma, as depicted in Fig. 5. This value is indeed somewhat greater than λ_{q}^{XGC} $\u2248$ 5.9 mm found on the full-current ITER model plasma, consistently with a slightly greater *B*_{pol,MP} *a*/*ρ*_{i,pol} value. Thus, our new formula passes this test, too.

At this point, we mention the error/uncertainty range in the Eich-formula fitting of the XGC1 data for the 12.5 MA ITER case. The Eich fitting formula, as described in Refs. 1 and 2, itself is well defined. The uncertainty range of λ_{q}^{XGC} fitting for the present devices was smaller than the Eich regression error range and was not discussed in Ref. 6 (the ITER 10 MA case can be used as an example to be presented later in this section). However, at such a large λ_{q}^{XGC} as in the 12.5 MA ITER case, we find that the noisy fluctuations in the heat-flux footprint are surfacing in the raw simulation data due to the smallness of the radial resolution compared with λ_{q}^{XGC} (see Fig. 6). This type of fluctuation in the XGC footprint is most likely from numerical noise due to particle noise and may not represent what is seen in the experiment. Possible difference between the numerical heat-flux measurement and the experimental thermal sensor measurement is the reason why we call *λ*_{q}^{XGC} the “heat-flux” width instead of the “heat-load” width. A long tail into the far scrape-off layer (SOL) can be noticed, which is unimportant for the peak divertor heat-load density. We can smooth out the footprint until the noisy fluctuation disappears. This introduces arbitrariness and uncertainty in the λ_{q}^{XGC} value measurement.

In the 12.5 MA ITER case, the raw data give the narrowest λ_{q}^{XGC}(min) fitting due to the sharp peak near the separatrix leg (see Fig. 6), caused by the parallel electron heat flow. In our Eich-formula fitting of λ_{q}^{XGC}, we try to emphasize the peak heat-load density around the separatrix leg. We find λ_{q}^{XGC}(min) = 5.5 mm. We then smooth the footprint data until all the noisy fluctuations disappear before estimating the widest possible λ_{q}^{XGC}(max). Here, we apply a nine-point (Δr ∼ 0.8 mm) moving-averaging in the radial direction and obtain λ_{q}^{XGC}(max) = 8.2 mm (see Fig. 7). The point depicted in Fig. 5 is the midpoint between these two values, with the error bar of about $\xb1$ 20% calculated from the maximal and minimal λ_{q}^{XGC} values. This type of uncertainty analysis was not performed on the 15 MA case in Ref. 6, but it can be assumed that a similar level of uncertainty exists.

- (iii)
10 MA Q = 5 steady-state ITER scenario case

There is a wide gap in the new parameter (*B*_{pol,MP} *a*/*ρ*_{i,pol}) space between the high-current JET plasma and the 15 MA ITER plasma. To check the validity and accuracy of the new machine-learned λ_{q}^{ML} formula, it is necessary to have at least one predictive simulation deep in the gap region as explained earlier. For this purpose, we pick the 10 MA Q = 5 ITER steady-state model plasma (see Fig. 5). XGC1 finds that λ_{q}^{XGC} from the raw footprint data are 2.5 mm and from the smoothed data are 2.8 mm. If we take 2.5 mm as the theoretical minimum value and 2.8 mm as the theoretical maximum value, the median value 2.65 mm and the error bar ($\xb1$ 6%) are marked in Fig. 5. The difference in the λ_{q}^{XGC} fitting between the raw data and the smoothed data are not as great as the 12.5 MA case since the finite radial grid size has already provided some smoothing (given that the spreading is lower than at 12.5 MA). As can be seen from Fig. 5, the validity of the new simple formula is remarkably good.

Since the 10 MA ITER case is located deep in the gap between 4.5 MA JET and 15 MA ITER, this is a good case to check the consistency of Eqs. (6)–(10) with the λ_{q}^{XGC} = 2.65 mm value found from Eq. (6). Table II summarizes the comparison. For reference, λ_{q}^{Eich(14)} = 0.53 mm. It can be seen that the simplest formula, Eq. (6), is the most consistent one with the XGC-found λ_{q}^{XGC} value for the 10 MA ITER case.

Formula no. . | λ_{q}^{ML} from various formulas
. | Ratio to λ_{q}^{XGC} = 2.65 mm
. |
---|---|---|

Eq. (6) | 2.77 mm | 1.05 |

Eq. (7) | 0.86 mm | 0.32 |

Eq. (8) | 1.79 mm | 0.68 |

Eq. (9) | 2.30 mm | 0.87 |

Eq. (10) | 2.24 mm | 0.84 |

## VI. NEW PHYSICS UNDERSTANDING AND ITS RELEVANCE TO THE PREDICTIVE FORMULA

As explained in Sec. II, the new parameter *a*/*ρ*_{i,pol}, representing the ratio between the device size and the ion poloidal gyroradius ($\u2248$ ion banana width in the edge plasma) comes from the important kinetic micro-physics that was not part of the macro-parameter set used in Refs. 1 and 2. This ratio determines the strength and weakness of the neoclassical effects, which include the background E × B-flow shearing rate [see Eq. (3)] that can control plasma turbulence.^{24} As the *a/ρ _{i,pol}* ratio becomes higher, the neoclassical E × B-flow shearing effect gets weaker, turbulence modes that were otherwise suppressed by a strong shear-flow could surface and, at the same time, the E × B-shear-flow driven turbulence can recede.

To investigate if there is a physics difference between the full-current ITER edge and the tokamak edge that follows the Eich/Goldston-scaling, we compare the turbulence property between the full-current 15 MA ITER edge that has much greater *a*/*ρ*_{i,pol} than in today's tokamaks than that of the 5 MA ITER edge. This choice removes the pure, absolute size effect in the comparison. Figure 8 depicts a snapshot pattern of the normalized electron density fluctuation δn/n obtained from the XGC1 simulations around the outboard midplane across the magnetic separatrix surface (vertical dashed line). It can be seen that across the outboard separatrix surface of the 5 MA ITER H-mode plasma, plasma turbulence is of the isolated blob type as seen in both XGC1 simulations and laboratory experiments on today's tokamaks.^{9} However, in the zoomed-in figure for the 15 MA full-current ITER, the turbulence is of radially extended/connected streamer type as usually seen in ITG and TEM turbulence.^{10}

For a deeper understanding of the turbulence modes, we study the phase correlation between the electron density fluctuation δn and the electrostatic potential fluctuation δΦ, and plot them in Fig. 9. When the electrons behave adiabatically, which is a typical signature of TEM modes, the phase correlation vanishes and the radial transport vanishes. It can be easily noticed that the electrons in the near-SOL have small phase correlation coefficient between δn and δΦ; hence, they are more adiabatic in the 5 MA ITER edge, which is the region where the λ_{q}^{XGC} footprint is measured, while they are strongly non-adiabatic in the near-SOL of the 15 MA ITER edge—actually, the strongly non-adiabatic region starts just inside the separatrix into the near-SOL. This is indication that the streamer type fluctuations seen in the 15 MA ITER have a strong TEM component. An ITG dominant turbulence has a stronger adiabatic electron response.

The third data analysis we performed is a simple unsupervised machine-learning analysis of the electron-response correlation to the edge turbulence just outside of the separatrix surface.^{28} The K-means clustering method in APACHE Spark^{29} is used to divide the electron response into six groups with each group represented by different colors. The result is depicted in Fig. 10 as a contour plot in two-dimensional velocity space (reprint from Fig. 3 of Ref. 28). It can be seen that the electrons are grouped mostly in energy—a sign of kinetic-energy dependent oscillations—except around $(v\u22252+v\u22a52)$^{1/2} ∼ 2 where there is a distinctively different response between the trapped and passing electrons. In this energy band, dark navy blue and medium sapphire blue are separated at the trapped-passing boundary. This is a sign of trapped electron mode driven turbulence. Different behavior around v_{ǁ} ∼ 0 in the trapped electron response band is not a surprise since the deeply trapped electrons around the outboard midplane do not experience much toroidal precession drift (TEMs are driven by resonance between toroidal precession drift of the trapped electrons and drift waves). Higher number of clustering groups could show a more detailed and gradual change. The vertical Landau resonance pattern in accordance with *k*_{ǁ}v_{ǁ} ∼ *ω* is not seen, indicating that the turbulence may not be from ITG modes. Besides, there is evidence in the literature that ITG modes cannot survive in the SOL.^{30}

All three pieces evidence (streamer-like structures, non-adiabatic electrons, and different response of trapped electrons from passing electrons at a specific energy band) suggest that the turbulence modes are TEMs. It is well-known that the streamer-type TEM turbulence is highly effective in transporting plasma energy along the radial streamers for electrostatic potential perturbations on the order 10^{−2} relative to the electron thermal energy.^{10} At the same time, evidence exists that blobby turbulence may not be effective in the radial transport of plasma energy and that the heat-flux spreading seen in present devices is mostly from the ion neoclassical orbit effect.^{3,6} Details of the electron and ion transport in blobby turbulence are the subject of an on-going study. We note here that due to the high drift frequency in the H-mode edge, *ω** ∼ v_{th}(*ρ*/*L*) with a short gradient scale length *L*, the weakly collisional trapped electron modes can easily be triggered at higher electron kinetic energies—according to the resonance relation *ω** ∼ *U*_{precess} ∼ v(*ρ*/*R*)(*B _{0}*/

*B*

_{P})—around the magnetic separatrix if the effective electron collision frequency is low ν

_{e*}$\u2272$ 1 and the local E × B-flow shearing rate is low. Using the XGC1 simulation parameters, we find ν

_{e*}(ψ

_{99}, q

_{95}) $\u2243\u2009$ 0.9 for the ITER 12.5 MA edge and ν

_{e*}(ψ

_{99,}q

_{95}) $\u2243$0.95 for the ITER 15 MA edge, where ν

_{e*}(ψ

_{99}, q

_{95}) is defined using the plasma density and temperature at ψ

_{99}, but the safety factor q is measured at ψ

_{95}. We also find that ν

_{e*}(ψ

_{99}, q

_{95}) for the ITER 5 MA edge is similarly low, indicating that the low electron collisionality is not a sufficient condition for the occurrence of a wide λ

_{q}

^{XGC}, but only a necessary condition (requiring a weak E × B-flow shearing rate also).

In fact, together with the low electron collisionality, a weak E × B-flow shearing rate across the separatrix surface in the high current ITER edge is observed in XGC1, while a strong E × B shearing rate is always observed in XGC1—and in the laboratory experiments—in the edge of present tokamaks. Figure 11(a) depicts the mean electrostatic potential profile in the pedestal and across the separatrix of the 15 MA ITER plasma and, for comparison, the equivalent for the JET 4.5 MA plasma in Fig. 11(b). Vertical axes are approximately scaled to be proportional to the pedestal temperature for each plasma: 5 keV for the 15 MA ITER pedestal and 1.75 keV for the 4.5 MA JET pedestal. A large difference in the E × B-flow shearing rates across the magnetic separatrix can be easily implied from these figures. The actual E × B-flow shearing rate across Ψ_{N} = 1 (normalized to diamagnetic frequency at $k\u22a5=1/\rho i,pol)$ is in fact compared in Fig. 12 for the JET 4.5 and 15 MA ITER discharges, together with the 1.5 MA DIII-D case. We comment here in passing that the zonal flow oscillations are more noticeable in the 15 MA ITER edge, which will be further subject for future study.

For reference, we show in Fig. 13 the plasma density and temperature profile inputs used in the XGC1 simulation of the 15 MA ITER plasma which produced Fig. 11(a). The blue lines represent the electron density (n_{e}) and temperature (T_{e}) input profiles initially tried in XGC1, supplied from JINTRAC integrated modeling of a 15 MA ITER deuterium-plasma. The modeled ion temperature (T_{i}) profile is not shown, but is similar to T_{e}, with its value somewhat higher (lower) than T_{i} in the core (pedestal) region. As explained in Ref. 6 and earlier in this paper, XGC1 found that the ion-scale turbulence level was too high to maintain the JINTRAC-modeled n_{e} and T_{e,i} profiles and, as a result, the plasma power flow across the separatrix and to the divertor plates was an order of magnitude higher than the edge power flow of 100 MW expected in a Q = 10 ITER burning plasma (50 MW additional heating, 100 MW alpha heating and 50 MW of core radiation). Following the direction of XGC1's pedestal profile relaxation, we ended up with the n_{e} pedestal shape input (red line) as shown in Fig. 13(a), and the T_{e} and T_{i} pedestal shapes plotted in Fig. 13(b) in red and yellow lines with an approximate power balance between the power crossing the separatrix ($\u2248$ 100 MW) and the total power deposited onto the divertor plates ($\u2248$ 90 MW).

The plasma profiles deep in the core region, manufactured to have similar electron and ion pressure as in the original JINTRAC model, are not to be trusted since the core turbulence had not yet been established by the time the XGC1 simulation was stopped. This is done to save computational time and is based on the criterion that the turbulence at the separatrix/SOL and the divertor heat flux footprint are saturated. The central plasma profiles still stay at the manufactured input level without being given a chance to evolve to a power balance. It will be an important future work to perform a much longer simulation, especially with electromagnetic turbulence, to find the self-organized plasma density and temperature values in the pedestal and central core of 15 MA ITER that are consistent with the 150 MW additional + alpha heating and turbulent/neoclassical transport. We also note here the following: (i) the outer divertor power-load was only ∼25% higher than that at the inner divertor in the 15 MA ITER plasma, unlike in the present tokamaks (and in fluid modeling of attached ITER burning plasmas with the SOLPS-ITER code^{16}) where XGC1 finds that the outer divertor power-load is almost twice higher; (ii) the divertor heat-flux width on outer divertor target is not well correlated with the plasma decay length in the near-SOL along the outer midplane (the so-called density SOL width). The cause of observation (i) is an equilibrated ion power deposition between inboard and outboard divertor plates, while the inboard electron power load is only about half of the outboard power load as observed in the present tokamak simulations. Preliminary results on the parameter dependence of the out/in divertor power deposition asymmetry will be presented later in this section. The observation (ii) indicates that the plasma energy crossing below the outboard midplane may be more important than the flux-tube connection effect between the outer divertor and outboard midplane. These topics are not well studied and are as yet inconclusive yet. They require more careful study in the future.

We caution here that the flux-surface-averaged mean electrostatic potential ⟨Φ⟩ in the far SOL shown in Fig. 11 may not be physically meaningful. Only the shape of ⟨Φ⟩ in the near-SOL—and radially inward—needs to be considered physical, with an unknown additive constant. First of all, what is solved in the gyrokinetic Poisson equation is not the absolute electrostatic potential value itself, but the first and the second derivatives of the electrostatic potential under a given boundary condition. Second, we use an artificial Dirichlet boundary condition (⟨Φ⟩ = 0) at the flux surface where the field lines connect to a material surface. In the case of the 15 MA ITER plasma, the contact of the plasma with the first wall occurs at the low field side. In other words, our axisymmetric electrostatic potential in the SOL is non-zero only in the region where the field lines intersect the inner and outer divertor plates without being intercepted by the first wall. Since the first wall surface touches the edge plasma only in certain small areas, large areas of the flux surface are filled with plasma which continues into the first wall shadow. In the real tokamak plasmas, this flux surface may have a mean positive ⟨Φ⟩ value relative to the limiter surface on the order of electron thermal energy. The reason for using an artificial ⟨Φ⟩ = 0 Dirichlet boundary condition before reaching the real material wall in these simulations is that when the particle number density becomes too low in the limiter/first-wall shadow, our axisymmetric Poisson solver sometimes does not give a converged solution. As a consequence of these assumptions in the far-SOL, we can only discuss the mean radial electric field and its shearing rate in the near-SOL, across the magnetic separatrix, and inward into pedestal in Fig. 11.

There could also be a question of how the steep H-mode pedestal gradient can be supported in the radial force balance equation at Ψ_{N} > 0.98 of the full-current ITER edge plasma where the radial electric field is small, as shown in Fig. 11(a). For the sake of argument, we use the radial force balance equation (1) derived for the closed flux surface, even though it may not be highly accurate across the separatrix surface. XGC1 finds that the plasma gradient across the magnetic separatrix (0.98 < Ψ_{Ν} < 1.01) is maintained by the local co-current parallel/toroidal flow across the magnetic separatrix [see Eq. (2) for a simpler equation]. We demonstrate this phenomenon in Fig. 14 by showing two representative forces across the Ψ_{Ν} = 1 surface: the radial force term from co-current toroidal flow (≈⟨v_{ǁ}⟩, green line) which is of the same order of magnitude and opposite to the radial density gradient force (dashed line). Other terms are less significant and are not shown in the figure. The physics origin for this phenomenon is the X-point orbit-loss driven E_{r} and toroidal torque.^{31} The neoclassical dielectric/polarization effect^{31,32} and the collisional damping of poloidal plasma rotation in a tokamak plasma^{17,33} can easily suppress the weak radial electric field, but the weak toroidal viscosity cannot easily suppress the toroidal rotation. Without the radial electric field opposing the X-point orbit-loss driven toroidal flow, the toroidal flow can replace the role of the radial electric field. A discussion of the physics of kinetic co-current edge momentum generation across the magnetic separatrix by X-point orbit loss torque can be found in Refs. 31 and 34.

The spatial turbulence pattern of the 10 MA steady-state ITER edge plasma is of special interest since it shows only a partial enhancement of λ_{q}^{XGC} compared with the expected experimental scaling value. It can be seen from Fig. 15 that the temperature-normalized electrostatic potential fluctuation across the outboard-midplane magnetic separatrix is a mixture of blobs (isolated structures at high amplitude, red and blue) and streamers (connected structures at low amplitude). The streamer feature has not been seen in the XGC1 simulations of present tokamaks, where only the blob feature has been observed. The partial enhancement of λ_{q}^{XGC} in the 10 MA ITER edge appears to be from the low amplitude streamers, which are known to be highly effective carriers of heat from core-region turbulence studies,^{10} as explained earlier. This is valuable information. The large enhancement of λ_{q}^{XGC} in the 15 MA or 12.5 MA ITER plasma is not from a sudden physics bifurcation, but is a gradual effect occurring as a result of the transition from blobs to streamer transport. An explicit transport mechanism study of kinetic electron and ion particles as they pass through the blobs and streamers in the open-field line region under parallel streaming and perpendicular drift motions is presently underway using an *in situ* data management technology. It will be reported in the near future.

Another noteworthy observation we have made from the gyrokinetic ITER simulations is the dependence of the power deposition ratio between the outer and inner divertor plates on the new scaling parameter B_{pol,MP} a/ρ_{i,p} used in the machine learning approach. As shown in Fig. 16(a), the out/in power ratio decreases as B_{pol,MP} a/ρ_{i,p} increases from the 5 MA plasma to the 10 and 15 MA plasmas. At 5 MA, the out/in ratio of ∼1.7 is similar to the present tokamak values. At 15 MA, the ratio decreases to 1.25. The peculiar 12.5 MA plasma (star mark), though, shows an irregular behavior compared to the other cases. This could mean that the reduction amount of the outer/inner power deposition ratio in the ITER 15 MA could be subject to some unknown effects that need to be studied.

Figure 16(b) depicts the same graph as in Fig. 16(a), but now as a function of a/ρ_{i,p}. The same trend is found, meaning that the out/in power deposition ratio behavior cannot be definitely identified due to the enhanced B_{pol,MP}(T) a/ρ_{i,p value} or the enhanced a/ρ_{i,p} value. It appears that the reduction in the inner/outer divertor power deposition ratio from 5 MA, to 10 MA, and to 15 MA is related to the co-current directional parallel plasma flow, thus positive poloidal flow, across the separatrix surface (see green line in Fig. 13) which could bring more plasma power to the inner divertor plates. In common with several other detailed phenomena observed from the simulations, further work is required to provide a more definitive answer to this question.

## VII. SUMMARY AND DISCUSSION

The XGC1 gyrokinetic particle-in-cell code in the electrostatic mode, with which the predictive divertor heat-flux width simulations have reproduced the experimentally measured λ_{q} from the three U.S. major tokamaks in the attached (inter-ELM H-mode) low recycling divertor regime, has reported a much wider divertor heat-flux width, λ_{q}^{XGC} for the full-current (15 MA) ITER model plasma than expected on the basis of the experimental scaling.^{6} Several new simulations are performed to answer some essential questions following the previous report. How would XGC1 predict λ_{q}^{XGC} on the highest current JET plasma, which has *B*_{pol,MP}(=0.89 T) only 26% lower than the *B*_{pol,MP}(=1.21 T) of the full current ITER? This question was especially worthy to answer because the old JET experimental data showed some broadening-like data points even at lower *B*_{po,MP} values, as can be seen in Fig. 1 or in the Eich-scaling reports^{1,2} (red circular dots). Were these old JET data from inaccurate experimental measurement or from real physics? Our simulation predicts that the highest-*B*_{pol,MP} JET discharge has λ_{q}^{XGC,} that is consistent with the Eich scaling (red open circle in Fig. 1). This result could suggest a possibility for a λ_{q}^{XGC} bifurcation between *B*_{pol,MP} = 0.89 T of JET and the *B*_{pol,MP} = 1.21 T of the 15 MA ITER discharge.

A more significant question then arises. In a C-Mod experiment, *B*_{pol,MP} was raised to the level of the full-current ITER and it was found that λ_{q}^{Exp} still follows the Eich formulas. An XGC1 simulation performed and agreed with the experimental finding (see the black open star symbol at the far-right bottom of Fig. 1), giving rise to double valued solutions if *B*_{pol,MP} is the sole parameter in *λ*_{q}^{Eich(14)}. This questioned the existence of a bifurcation of *λ*_{q}^{XGC} with *B*_{pol,MP} and suggested a hidden parameter outside of the macroscopic parameter set used in Refs. 1–3.

A supervised machine-learning tool is applied to all the λ_{q}^{XGC} data points (together with the corresponding experimental data points λ_{q}^{Exp}) obtained for the existing tokamaks and the full-current 15 MA Q = 10 ITER plasma, with a feature engineering of adding the physics-based kinetic parameter a/ρ_{i,pol} to *B*_{pol,MP}. The result, shown in Fig. 3, is a new simple formula for *λ*_{q}^{XGC} that reduces to *λ*_{q}^{Eich(14)} in the present tokamak regime including the highest current C-Mod case, that reproduces the full-current ITER result, and that is physically meaningful. The new additional simplest dependence parameter is found to be *B*_{po,MP} (a/ρ_{i,pol}), combination of the neoclassical E × B-flow shearing rate parameter ρ_{i,pol}/a and the ion orbit width parameter 1/*B*_{po,MP}.

Tests of the new formula are performed using a 5 MA H-mode ITER plasma which has a *B*_{po,MP} (a/ρ_{i,pol}) value similar to that in existing tokamaks, a 12.5 MA Q = 5 long pulse ITER plasma with *B*_{po,MP} (a/ρ_{i,pol}) slightly greater than the full-current ITER plasma, and a 10 MA Q = 5 steady-state ITER plasma which has *B*_{po,MP} (a/ρ_{i,pol}) in the gap between the highest current JET and the full current 15 MA ITER. The new simplest formula well survives against these tests, as depicted in Fig. 5. Other, more complicated formulas suggested by the machine learning program did not do well against the 10 MA ITER test, which lies deep in the gap region between the JET and the 15 MA ITER points in the new parameter space.

In an effort to study the new physics that leads to the *λ*_{q}^{XGC} broadening in the full-current ITER and that is consistent with the new parameter, three independent data analyses are performed. The study identifies the new physics to be weakly collisional, trapped-electron driven turbulence, gradually dominating over the blobby turbulence as the new parameter *B*_{po,MP} (a/ρ_{i,pol}) increases.

We comment here that the main difference in the present gyrokinetic simulation results from the low divertor pressure case in the recent 15 MA ITER result of Kaveeva *et al.,*^{16} which used SOLPS-ITER code with an assumed anomalous electron thermal diffusivity of 1 m^{2}/s in SOL are as follows: (i) much smaller value of the E × B-flow shear across the separatrix, (ii) ∼2× wider outer divertor heat-load width, (iii) weaker outboard/inboard power load ratio, and (iv) smaller effective heat diffusivity at $\u2248$ 0.2 m^{2}/s (an averaged value across the separatrix surface 0.98 $\u2264$ Ψ_{N} $\u22641.02)$. The physics relationship between the ∼2× wider outer divertor heat-load width and the eventual relaxation to ∼2× wider edge pedestal width has not been established from the present gyrokinetic simulations. As stated in our previous report,^{6} λ_{q}^{XGC} saturates before the ∼2× relaxation of the pedestal width is reached. The above-quoted effective radial diffusion coefficient is only a ball-park number. Radial plasma fluxes fluctuate significantly along the field line depending upon the space-time varying turbulence structure and, thus, a “flux-surface-averaging” is employed to obtain a statistically accurate value in a core-region plasma. In the open field region and across the separatrix surface, the survival time of an individual particle motion is short due to divertor-plate intersection and atomic physics, hence the flux-surface-averaging is limited and yields a higher statistical error. An advanced data analysis technique is under development to resolve this issue, by accurately following the individual particle motions in the turbulent field while obtaining statistical transport information, in a similar way to the transport measurement used in stochastic systems [see Eqs. (10) and (11) in Ref. 35 and the quoted references therein].

We note here that the present simulation is electrostatic. Even though the electrostatic XGC has reproduced *λ*_{q} in the present tokamaks, the effect of the electromagnetic turbulence on *λ*_{q}^{XGC} in the high Q ITER edge is of interest. Present studies are conducted under the low-recycling attached divertor conditions, corresponding to the condition relevant to Refs. 1–3. ITER will have to operate in the semi-detached or detached divertor regimes for high Q plasmas. These subjects, and others, are left for future study. In addition, a way to test the new formula in the present experiments is of interest. This may require finding or creating a plasma with ν_{e*} $\u2272$ 1 and a low-sheared E × B flow near the magnetic separatrix surface.

A shortfall not mentioned in the main text is the lack of a systematic validation metric^{36} from the XGC1 simulation results due to the small number of extreme-scale simulations and highly limited availability of the experimental primacy hierarchy data in the edge plasma. Systematic validation of limited number, extreme-scale simulations is an active research topic in the uncertainty quantification community.

## ACKNOWLEDGMENTS

We acknowledge helpful discussion with M. Romanelli, T. Eich and R. Goldston in the early phase of the study. We thank R. Maingi, J.-W. Ahn, T. Gray, B. LaBombard, T. Leonard, M. Makowski, and J. Terry for their contribution to the original paper^{6} with experimental data from three U.S. tokamaks and valuable discussions. Without their help, this work could not have been launched.

The research was mostly supported by the U.S. DOE Office of Science SciDAC program, both from ASCR and FES, through the award the SciDAC-4 Partnership Center for High-fidelity Boundary Plasma Simulation (HBPS) under the Contract No. DEAC02–09CH11466 with Princeton University for Princeton Plasma Physics Laboratory. Computing resource is provided by OLCF and ALCF via INCITE program. The simulations of ITER plasmas have been carried out as part of the coordinated work programme of the Pedestal Confinement and Stability Modelling ITER Scientist Fellows group.

ITER is the Nuclear Facility INB No. 174. The views and opinions expressed herein do not necessarily reflect those of the ITER Organization.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## References

^{20}