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 λqXGC 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 λqXGC 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.

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 Bpol,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 Goldston3 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 IP = 15 MA with q95 = 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

h(RmpRmp,sep)=0.5h0exp[(0.5S/λq)2(RmpRmp,sep)/λq]Erfc[0.5S/λq(RmpRmp,sep)/S]+hBG,

where Rmp is the major radius along the outboard midplane, Rmp,sep is Rmp on the outboard separatrix surface, h(Rmp-Rmp,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), h0 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 hBG 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 nsep/nGW > 0.6 for λq ≈ 1 mm,4 where nGW 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 (λqXGC) 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 JINTRAC8 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 λqXGC 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 λqExp 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 λqXGC from XGC1 also follows the Eich scaling. This yields double-valued solutions for λqXGC between the high-current C-Mod plasma and the full-current ITER plasma if Bpol,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 λqXGC 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 λqXGC 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 effect12 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 result15 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 m2/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 λqXGC. 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.

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 Bpol,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 experiments11 with Bpol,MP reaching the ITER full-current case did not exist.

TABLE I.

Experimental discharges from three U.S. tokamaks that were used by XGC in Ref. 6 and that were part of the original “Eich-regression” exercise.1,2 The last row represents a full-current, full-power ITER scenario plasma. BT is the toroidal magnetic field strength at machine axis, IP is the plasma current, and Bpol,MP is the poloidal magnetic field strength at outboard midplane on the separatrix surface.

ShotTime (ms)BT (T)IP (MA)Bpol,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 
ShotTime (ms)BT (T)IP (MA)Bpol,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 λqExp 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 λqXGC = 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 λqEich(14)], represented by the solid line, together with the regression error represented by the two dashed lines. Here, we use the Eich formula #14 [λqEich(14) 0.63Bpol−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 λqXGC in the full-current Q = 10 ITER scenario plasma (15 MA, Bpol,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 λqEich(14).

FIG. 1.

XGC that has predicted the λq values in agreement with the Eich scaling formula in all three major U.S. tokamaks predicts λq = 5.9 mm in a model ITER plasma edge at IP = 15 MA, as shown in Ref. 5. After the publication of Ref. 5, two new simulation points are added: JET at 4.5 MA (Bpol,MP = 0.89 T, open red circle) and C-Mod at 1.4 MA (Bpol,MP = 1.11 T, open black star) at far-right bottom. They both follow the Eich scaling formula, approximately. A couple of new high field C-Mod experimental points are also shown together (at Bpol,MP = 1.1 T).

FIG. 1.

XGC that has predicted the λq values in agreement with the Eich scaling formula in all three major U.S. tokamaks predicts λq = 5.9 mm in a model ITER plasma edge at IP = 15 MA, as shown in Ref. 5. After the publication of Ref. 5, two new simulation points are added: JET at 4.5 MA (Bpol,MP = 0.89 T, open red circle) and C-Mod at 1.4 MA (Bpol,MP = 1.11 T, open black star) at far-right bottom. They both follow the Eich scaling formula, approximately. A couple of new high field C-Mod experimental points are also shown together (at Bpol,MP = 1.1 T).

Close modal

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 “ρi,pol is the poloidal ion Larmor radius at the outboard midplane separatrix point and “a” is the plasma minor radius. In the full-current ITER, ρi,pol/a 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,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 region17 

u||+(Ti/neBp)dn/dr=(1/eBp)[(k1)dTi/dredΦ/dr],
(1)

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 regime17 (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

uE/vi,polu||/viρi,pol/αa,
(2)

where uE = Er/B is the E × B flow speed, vi is the ion thermal speed, vi,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, uE 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 γΕ/ω* = duE/(ω* dr) which scales as, using Eq. (2) and the relations ω*∼kviρi/αa, duE/dr ∼ uE/LE, and assuming kρi,pol ∼ 1

γΕ/ω*(duE/vidr)(ρi,pol/ρi)αavi,pol(ρi,pol/αaLE)(ρi,pol/ρi)αaρi,pol/LE.
(3)

This relationship shows that γΕ/ω* decreases with ρi,pol/LE (ρip/a). Value of the XGC-found γΕ/ω* for representative device cases will be shown in Sec. VI, where detailed physics is discussed.

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 Bpol,MP value and the physical size in deuterium plasma. For this purpose, a JET 4.5 MA discharge19 is chosen that has the highest Bpol,MP (0.89 T) at the time of simulation (unfortunately, an experimental λqExp measurement does not exist on JET at this high value Bpol,MP). To be more specific, Bpol,MP for this JET plasma is only 36% lower than the full-current ITER plasma, and its linear size is a factor of 2 smaller than ITER. In this plasma, XGC1 finds λqXGC of about 0.64 mm, which is within the regression error bar from the Eich(14) value λqEich 0.72 mm (open red circle in Fig. 1). Thus, XGC1 indicates that there may be either a bifurcation of λqXGC between Bpol,MP = 0.89 T of JET and 1.21 T of ITER, or there is something other than the value of Bpol,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 Bpol,MP values up to 1.3 T (Ref. 11) which somewhat exceeds the full-current ITER value and found that the experimental λqExp still follows λqEich 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 Bpol,MP = 1.11 T. At this high value of Bpol,MP, though, we find λqXGC 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 λqEich(14) = 0.56 mm. As a result, XGC1's solution becomes double valued around the maximal C-Mod Bpol,MP values if Bpol,MP is used as the sole parameter and suggests the existence of hidden parameter(s) that was missed in Eich's regression process.

FIG. 2.

Footprint of the parallel heat-flux on the outer divertor plates, mapped to the outboard midplane, in one of the highest current C-Mod discharges (discharge number #1160930033).

FIG. 2.

Footprint of the parallel heat-flux on the outer divertor plates, mapped to the outboard midplane, in one of the highest current C-Mod discharges (discharge number #1160930033).

Close modal

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 DE to represent a set of λqExp data found from the present laboratory experimental measurements, DSE for a set of λq data found through high-fidelity simulation of the existing experiments, and DSF for a set of λq data found through high-fidelity simulation of future experiments. We use M to denote the machine-learning operation, FE for the modeling formula found by the operation M on the present experimental dataset DE, FSE for the modeling formula found by the operation M on DSE, and FP for the predictive modeling formula found by the operation M on all the datasets including DE, DSE, and DSF. DE and DSE do not need to have one-to-one correspondents.

For the validated high-fidelity simulations, we assume FEFSE as a pre-requisite condition, which is satisfied by XGC1 as shown in Secs. II and III. Thus, we have M(DE) FE and M(DSE) FE, with some allowance for error. We can then write down the following relations:

M(DEDSE)FEand
(4)
M(DEDSEDSF)FP[FE].
(5)

Here, FP[FE] means that the machine-learned formula FP reduces to FE 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 FP, within the simulation accuracy, by combining DE and DSE with DSF.

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 FP of Eq. (4). We present the result first: Fig. 3 depicts the simplest FP result from Eureka, as to be elaborated soon later in this section. Figure 3 contains the selected experimental dataset DE from NSTX, DIII-D, and C-Mod (marked with + symbols) as presented in Sec. II, and the corresponding simulation dataset DSE. 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, λqEich(14) = 0.63 Bpol,MP−1.19. The simple extrapolation to the future experiments from the present-day experimental dataset is represented by the solid black horizontal line.

FIG. 3.

A new λq formula (purple equation and dashed line) from simulation-anchored machine learning: Among the color-filled circles, NSTX (dark blue), DIII-D (orange) and C-Mod (gray) represent the XGC1 predicted dataset DSE for the experiments that has experimental dataset DE, marked with + sign with the same coloring scheme. The JET and the 15 MA ITER data points are the simulation dataset DSF that do not have experimental measurements. The solid black arrow shows a simple extrapolation from the present experiments to the unknown physics regime.

FIG. 3.

A new λq formula (purple equation and dashed line) from simulation-anchored machine learning: Among the color-filled circles, NSTX (dark blue), DIII-D (orange) and C-Mod (gray) represent the XGC1 predicted dataset DSE for the experiments that has experimental dataset DE, marked with + sign with the same coloring scheme. The JET and the 15 MA ITER data points are the simulation dataset DSF that do not have experimental measurements. The solid black arrow shows a simple extrapolation from the present experiments to the unknown physics regime.

Close modal

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: Btor (the toroidal magnetic field strength), Bpol,MP, q95 (the safety factor at the 95% poloidal-flux surface), PSOL (the power flow from core into the SOL), Rgeo (the geometric major radius), a (the plasma minor radius), Ip (the plasma current), and n/nGW (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 Bpol,MP by targeting all the present tokamaks, denoted here as the Eich regression number #14, with the squared correlation coefficient being R2 = 0.86. Our machine learning operation utilizes λqEich(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 size6,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 Er × 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: Bpol,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 λqEich formula for the present attached divertor experiments, i.e., FP[FE], 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 DEDSEDSF 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 Bpol,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

λq/λqEich(14)=f(Bpol,MP,a/ρi,pol)=f(Bpol,MP,a/ρi,pol,Bpol,MPa/ρi,pol),and=f(Bpol,MPa/ρi,pol).

Among the simulation-anchored formulas found by Eureqa

λqML=0.63Bpol,MP1.19[1.0+1.08×1010(Bpol,MPa/ρi,pol)4]withRMSerror=18.7%
(6)

is the simplest and lowest order expression for the heat flux width λqML 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 λqXGC for the full-current ITER plasma and reproduces λqEich(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 Bpol,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 Bpol,MP and a/ρi,pol, and that yield low RMSE include

λqML=0.63Bpol,MP1.19[1.0+1.961×1016(a/ρi,pol)6],RMSE=17.9%,
(7)
=0.63Bpol,MP1.19[1.0+1.68×1018Bpol,MP(a/ρi,pol)7],RMSE=17.0%
(8)
=0.63Bpol,MP1.19[1.0+1.91×1013(Bpol,MPa/ρi,pol)5],RMSE=17.0%
(9)
=0.63Bpol,MP1.19[1.0+3.46×104[4.04×105Bpol,MP(a/ρi,pol)3]]Bpol,MP,RMSE=16.1%.
(10)

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, Bpol,MP, a/ρi; mathematical operations; and variables), the evolutionary model search process in Eureqa, and the resulting λqML formulas (only one of them is shown).

FIG. 4.

Workflow used in the present supervised machine-learning study: The evolutionary model search is performed by the Eureqa program.20,21

FIG. 4.

Workflow used in the present supervised machine-learning study: The evolutionary model search is performed by the Eureqa program.20,21

Close modal

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.

FIG. 5.

Addition of three new ITER test cases to Fig. 3. The errorbar is relatively large at large Bpol,MP a/ρi,pol (20%, at the ITER 12.5 MA point) compared to the small Bpol,MP a/ρi,pol cases (6.3%, at the ITER 10 MA point).

FIG. 5.

Addition of three new ITER test cases to Fig. 3. The errorbar is relatively large at large Bpol,MP a/ρi,pol (20%, at the ITER 12.5 MA point) compared to the small Bpol,MP a/ρi,pol cases (6.3%, at the ITER 10 MA point).

Close modal

We note here that an extension of the high current (4.5 MA) JET plasma that is modeled toward the Bpol,MPa/ρ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 33- and 93-times greater. The first H-mode plasma scenario that will be explored in the initial ITER experimental phases with Ip = 5 MA25 is an excellent case to answer this question: It has Bpol,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 Bpol,MP a/ρi,pol for 5 MA ITER falling well within the present device range (see Fig. 5). For a quantitative comparison, the Bpol,MP a/ρi,pol value for the 5 MA ITER case is as small as 87, with Bpol,MP a/ρi,pol for all the present tokamak experiments falling between about 10 and 200. The test XGC1 simulation finds λqXGC = 2.2 mm, which satisfies the Eich formula value λqEich(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 λqXGC 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 Btor = 5.3 T and fusion gain of Q = 5 (Ref. 26) is an interesting case. Its toroidal magnetic field strength at the machine axis Btor = 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 Bpol,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 Bpol,MPa/ρ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 λqXGC by the large Bpol,MPa/ρ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 λqXGC 6.9 mm for this 12.5 MA ITER model plasma, as depicted in Fig. 5. This value is indeed somewhat greater than λqXGC 5.9 mm found on the full-current ITER model plasma, consistently with a slightly greater Bpol,MPa/ρ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 λqXGC 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 λqXGC 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 λqXGC (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 λqXGC 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 λqXGC value measurement.

FIG. 6.

Eich formula fitting over the raw XGC1 data, mapped to the outer midplane, for the 12.5 MA ITER case to obtain λqXGC = 5.5 mm. To emphasize the heat-load at the peak, the fit is performed using data between −5 and +6 mm. The blue fitting curve is cropped to show the fitting data range. The long tail into far SOL is ignored.

FIG. 6.

Eich formula fitting over the raw XGC1 data, mapped to the outer midplane, for the 12.5 MA ITER case to obtain λqXGC = 5.5 mm. To emphasize the heat-load at the peak, the fit is performed using data between −5 and +6 mm. The blue fitting curve is cropped to show the fitting data range. The long tail into far SOL is ignored.

Close modal

In the 12.5 MA ITER case, the raw data give the narrowest λqXGC(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 λqXGC, we try to emphasize the peak heat-load density around the separatrix leg. We find λqXGC(min) = 5.5 mm. We then smooth the footprint data until all the noisy fluctuations disappear before estimating the widest possible λqXGC(max). Here, we apply a nine-point (Δr ∼ 0.8 mm) moving-averaging in the radial direction and obtain λqXGC(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 ± 20% calculated from the maximal and minimal λqXGC 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

FIG. 7.

Eich formula fitting over the nine-point (Δr ∼ 0.8 mm) moving averaged XGC1 data along the divertor plates for the 12.5 MA ITER case. The same fitting range is used as in Fig. 6. λqXGC = 8.2 mm is obtained from this fitting. Again, the blue curve is cropped to show range of the data used for fitting.

FIG. 7.

Eich formula fitting over the nine-point (Δr ∼ 0.8 mm) moving averaged XGC1 data along the divertor plates for the 12.5 MA ITER case. The same fitting range is used as in Fig. 6. λqXGC = 8.2 mm is obtained from this fitting. Again, the blue curve is cropped to show range of the data used for fitting.

Close modal

There is a wide gap in the new parameter (Bpol,MPa/ρ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 λqML 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 λqXGC 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 (± 6%) are marked in Fig. 5. The difference in the λqXGC 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 λqXGC = 2.65 mm value found from Eq. (6). Table II summarizes the comparison. For reference, λqEich(14) = 0.53 mm. It can be seen that the simplest formula, Eq. (6), is the most consistent one with the XGC-found λqXGC value for the 10 MA ITER case.

TABLE II.

Ratio of the λqML values from various ML-found formulas to the XGC-found λqXGC(= 2.65 mm) for the 10 MA ITER scenario plasma where a/ρi,pol = 384 and Bpol,MP a/ρi,pol = 446 (T).

Formula no.λqML from various formulasRatio to λqXGC = 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 
Formula no.λqML from various formulasRatio to λqXGC = 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 

As explained in Sec. II, the new parameter a/ρi,pol, representing the ratio between the device size and the ion poloidal gyroradius ( 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 

FIG. 8.

Comparison of the edge turbulence pattern in δn/n obtained from XGC1 between 5 MA ITER and 15 MA ITER across the separatrix surface. Isolated blob structure can be seen in the 5 MA case, with the relative density fluctuation amplitude becoming large across the separatrix. On the other hand, turbulence becomes a connected streamer type in the 15 MA case.

FIG. 8.

Comparison of the edge turbulence pattern in δn/n obtained from XGC1 between 5 MA ITER and 15 MA ITER across the separatrix surface. Isolated blob structure can be seen in the 5 MA case, with the relative density fluctuation amplitude becoming large across the separatrix. On the other hand, turbulence becomes a connected streamer type in the 15 MA case.

Close modal

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 λqXGC 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.

FIG. 9.

Comparison of the phase correlation between the electron density fluctuation δn and the electrostatic potential fluctuation δΦ around the outboard midplane edge. It can be noticed that the electrons are nearly adiabatic in the near-SOL of the 5 MA ITER edge, while they are strongly non-adiabatic in the 15 MA ITER edge.

FIG. 9.

Comparison of the phase correlation between the electron density fluctuation δn and the electrostatic potential fluctuation δΦ around the outboard midplane edge. It can be noticed that the electrons are nearly adiabatic in the near-SOL of the 5 MA ITER edge, while they are strongly non-adiabatic in the 15 MA ITER edge.

Close modal

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 Spark29 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 (v2+v2)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 

FIG. 10.

Electron response correlation in velocity space to the edge turbulence just outside of the separatrix surface, found from unsupervised machine-learning. K-means clustering method, specifying only six groups, is used using APACHE Spark. It can be seen that the trapped electrons around the energy band [(v2+v2)/ve2]1/2 ∼ 2 respond to the edge turbulence in a correlated way. Reproduced with permission from Churchill et al., IEEE Trans. Plasma Sci. 48, 2661 (2020). Copyright 2020 IEEE.

FIG. 10.

Electron response correlation in velocity space to the edge turbulence just outside of the separatrix surface, found from unsupervised machine-learning. K-means clustering method, specifying only six groups, is used using APACHE Spark. It can be seen that the trapped electrons around the energy band [(v2+v2)/ve2]1/2 ∼ 2 respond to the edge turbulence in a correlated way. Reproduced with permission from Churchill et al., IEEE Trans. Plasma Sci. 48, 2661 (2020). Copyright 2020 IEEE.

Close modal

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, ω* ∼ vth(ρ/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 ω* ∼ Uprecess ∼ v(ρ/R)(B0/BP)—around the magnetic separatrix if the effective electron collision frequency is low νe* 1 and the local E × B-flow shearing rate is low. Using the XGC1 simulation parameters, we find νe*99, q95) 0.9 for the ITER 12.5 MA edge and νe*99, q95) 0.95 for the ITER 15 MA edge, where νe*99, q95) 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, q95) 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 λqXGC, 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=1/ρ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.

FIG. 11.

Structure of the mean, flux-surface-average electrostatic potential e ⟨Φ⟩ across the pedestal and magnetic separatrix from XGC1 as a function of the normalized radial coordinate ψΝ (poloidal magnetic flux) for (a) the 15 MA ITER edge where a relatively small variation across the separatrix surface can be seen that leads to a weak E × B flow and hence its shearing rate and (b) the 4.5 MA JET edge where the relatively strong e ⟨Φ⟩ variation across the separatrix is a typical phenomenon seen in present-day tokamaks.

FIG. 11.

Structure of the mean, flux-surface-average electrostatic potential e ⟨Φ⟩ across the pedestal and magnetic separatrix from XGC1 as a function of the normalized radial coordinate ψΝ (poloidal magnetic flux) for (a) the 15 MA ITER edge where a relatively small variation across the separatrix surface can be seen that leads to a weak E × B flow and hence its shearing rate and (b) the 4.5 MA JET edge where the relatively strong e ⟨Φ⟩ variation across the separatrix is a typical phenomenon seen in present-day tokamaks.

Close modal
FIG. 12.

Comparison of E × B flow shearing rate γE normalized to the electron diamagnetic frequency at ΨΝ = 1 among 1.5 MA DIII-D, 4.5 MA JET and 15 MA ITER, obtained from XGC1 (see Table I). The diamagnetic frequency is measured at the wavelength using Bpol,MP. a/ρi,pol values are 81, 231, and 473, respectively, for 1.5 MA DIII-D, 4.5 MA JET and 15 MA ITER.

FIG. 12.

Comparison of E × B flow shearing rate γE normalized to the electron diamagnetic frequency at ΨΝ = 1 among 1.5 MA DIII-D, 4.5 MA JET and 15 MA ITER, obtained from XGC1 (see Table I). The diamagnetic frequency is measured at the wavelength using Bpol,MP. a/ρi,pol values are 81, 231, and 473, respectively, for 1.5 MA DIII-D, 4.5 MA JET and 15 MA ITER.

Close modal

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 (ne) and temperature (Te) input profiles initially tried in XGC1, supplied from JINTRAC integrated modeling of a 15 MA ITER deuterium-plasma. The modeled ion temperature (Ti) profile is not shown, but is similar to Te, with its value somewhat higher (lower) than Ti 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 ne and Te,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 ne pedestal shape input (red line) as shown in Fig. 13(a), and the Te and Ti pedestal shapes plotted in Fig. 13(b) in red and yellow lines with an approximate power balance between the power crossing the separatrix ( 100 MW) and the total power deposited onto the divertor plates ( 90 MW).

FIG. 13.

XGC1-adjusted input pedestal-SOL density in (a) and temperatures in (b), which do not evolve noticeably further at the end of the simulation. The core profiles are from artificial models and not to be trusted.

FIG. 13.

XGC1-adjusted input pedestal-SOL density in (a) and temperatures in (b), which do not evolve noticeably further at the end of the simulation. The core profiles are from artificial models and not to be trusted.

Close modal

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 code16) 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 Er and toroidal torque.31 The neoclassical dielectric/polarization effect31,32 and the collisional damping of poloidal plasma rotation in a tokamak plasma17,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.

FIG. 14.

Two representative terms across the separatrix surface in the fluid radial force balance equation obtained from XGC in the 15 MA ITER simulation.

FIG. 14.

Two representative terms across the separatrix surface in the fluid radial force balance equation obtained from XGC in the 15 MA ITER simulation.

Close modal

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 λqXGC 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 λqXGC 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 λqXGC 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.

FIG. 15.

Turbulence structure across the outboard midplane separatrix (dashed line) of the 10 MA ITER plasma. (a) A mixture of blobs (dark blue and red at high amplitude) and streamers (most noticeable in the pale greenish yellow color) can be seen. (b) For an enhanced coloring of the mixture between isolated blobs and connected streamers, blob tops at |eδΦ/Te| 0.1 are colored white and the negligibly small fluctuation amplitude at |eδΦ/Te| 0.003 is colored black. Pale green and blue colors highlight the streamer range amplitude, showing connected structures across the separatrix surface.

FIG. 15.

Turbulence structure across the outboard midplane separatrix (dashed line) of the 10 MA ITER plasma. (a) A mixture of blobs (dark blue and red at high amplitude) and streamers (most noticeable in the pale greenish yellow color) can be seen. (b) For an enhanced coloring of the mixture between isolated blobs and connected streamers, blob tops at |eδΦ/Te| 0.1 are colored white and the negligibly small fluctuation amplitude at |eδΦ/Te| 0.003 is colored black. Pale green and blue colors highlight the streamer range amplitude, showing connected structures across the separatrix surface.

Close modal

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 Bpol,MP a/ρi,p used in the machine learning approach. As shown in Fig. 16(a), the out/in power ratio decreases as Bpol,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.

FIG. 16.

(a) Out/in divertor power deposition asymmetry factor Pasym decreases linearly with the new parameter Bpol,MP(T) a/ρi,p. The star mark represents the peculiar 12.5 MA case. (b) Pasym exhibits a similar decrease trend with a/ρi,p.

FIG. 16.

(a) Out/in divertor power deposition asymmetry factor Pasym decreases linearly with the new parameter Bpol,MP(T) a/ρi,p. The star mark represents the peculiar 12.5 MA case. (b) Pasym exhibits a similar decrease trend with a/ρi,p.

Close modal

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 Bpol,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.

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, λqXGC 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 λqXGC on the highest current JET plasma, which has Bpol,MP(=0.89 T) only 26% lower than the Bpol,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 Bpo,MP values, as can be seen in Fig. 1 or in the Eich-scaling reports1,2 (red circular dots). Were these old JET data from inaccurate experimental measurement or from real physics? Our simulation predicts that the highest-Bpol,MP JET discharge has λqXGC, that is consistent with the Eich scaling (red open circle in Fig. 1). This result could suggest a possibility for a λqXGC bifurcation between Bpol,MP = 0.89 T of JET and the Bpol,MP = 1.21 T of the 15 MA ITER discharge.

A more significant question then arises. In a C-Mod experiment, Bpol,MP was raised to the level of the full-current ITER and it was found that λqExp 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 Bpol,MP is the sole parameter in λqEich(14). This questioned the existence of a bifurcation of λqXGC with Bpol,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 λqXGC data points (together with the corresponding experimental data points λqExp) 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 Bpol,MP. The result, shown in Fig. 3, is a new simple formula for λqXGC that reduces to λqEich(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 Bpo,MP (a/ρi,pol), combination of the neoclassical E × B-flow shearing rate parameter ρi,pol/a and the ion orbit width parameter 1/Bpo,MP.

Tests of the new formula are performed using a 5 MA H-mode ITER plasma which has a Bpo,MP (a/ρi,pol) value similar to that in existing tokamaks, a 12.5 MA Q = 5 long pulse ITER plasma with Bpo,MP (a/ρi,pol) slightly greater than the full-current ITER plasma, and a 10 MA Q = 5 steady-state ITER plasma which has Bpo,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 λqXGC 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 Bpo,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 m2/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  0.2 m2/s (an averaged value across the separatrix surface 0.98  ΨN1.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 λqXGC 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 λqXGC 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* 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 metric36 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.

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 paper6 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.

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

1.
T.
Eich
,
A. W.
Leonard
,
R. A.
Pitts
,
W.
Fundamenski
,
R. J.
Goldston
,
T. K.
Gray
,
A.
Herrmann
,
A.
Kirk
,
A.
Kallenbach
,
O.
Kardaun
 et al,
Nucl. Fusion
53
,
093031
(
2013
).
2.
T.
Eich
,
B.
Sieglin
,
A.
Scarabosio
,
W.
Fundamenski
,
R. J.
Goldston
,
A.
Herrmann
, and
ASDEX Upgrade Team
,
Phys. Rev. Lett.
107
,
215001
(
2011
).
3.
R. J.
Goldston
,
Nucl. Fusion
52
,
013009
(
2012
).
4.
A.
Kukushkin
,
H. D.
Pacher
,
G. W.
Pacher
,
V.
Kotov
,
R. A.
Pitts
, and
D.
Reiter
,
J. Nucl. Mater.
438
,
S203
(
2013
);
I.
Veselova
,
E.
Kaveeva
,
V.
Rozhansky
,
I.
Senichenkov
,
A.
Poletaeva
,
R. A.
Pitts
, and
X.
Bonnin
,
Nucl. Mater. Energy
26
,
100870
(
2021
).
5.
M.
Greenwald
,
Plasma Phys. Controlled Fusion
44
,
R27
(
2002
).
6.
C. S.
Chang
,
S.
Ku
,
A.
Loarte
,
V.
Parail
,
F.
Köchl
,
M.
Romanelli
,
R.
Maingi
,
J.-W.
Ahn
,
T.
Gray
,
J.
Hughes
,
B.
LaBombard
,
T.
Leonard
,
M.
Makowski
, and
J.
Terry
,
Nucl. Fusion
57
,
116023
(
2017
).
7.
S.
Ku
,
C. S.
Chang
,
R.
Hager
,
R. M.
Churchill
,
G. R.
Tynan
,
I.
Cziegler
,
M.
Greenwald
,
J.
Hughes
,
S. E.
Parker
,
M. F.
Adams
,
E.
D'Azevedo
, and
P.
Worley
,
Phys. Plasmas
25
,
056107
(
2018
).
8.
M.
Romanelli
,
G.
Corrigan
,
V.
Parail
,
S.
Wiesen
,
R.
Ambrosino
,
P. D.
Silva Aresta Belo
,
L.
Garzotti
,
D.
Harting
,
F.
Koechl
,
T.
Koskela
,
L.
Lauro-Taroni
,
C.
Marchetto
,
M.
Mattei
 et al,
Plasma Fusion Res.
9
,
3403023
(
2014
).
9.
D. A.
D'Ippolito
,
J. R.
Myra
, and
S. J.
Zweben
,
Phys. Plasmas
18
,
060501
(
2011
).
10.
Y.
Xiao
and
Z.
Lin
,
Phys. Rev. Lett.
103
,
085004
(
2009
).
11.
D.
Brunner
,
B.
LaBombard
,
A. Q.
Kuang
, and
J. L.
Terry
,
Nucl. Fusion
58
,
094002
(
2018
).
12.
T.
Eich
,
P.
Manz
,
R. J.
Goldston
,
P.
Hennequin
,
P.
David
,
M.
Faitsch
,
B.
Kurzan
,
B.
Sieglin
,
E.
Wolfrum
,
and ASDEX Upgrade Team and EUROfusion MST Team,
Nucl. Fusion
60
,
056016
(
2020
).
13.
B.
Scott
,
New J. Phys.
4
,
52
(
2002
).
14.
B.
Rogers
,
J.
Drake
, and
A.
Zeiler
,
Phys. Rev. Lett.
81
,
4396
(
1998
).
15.
X. Q.
Xu
,
N. M.
Li
,
Z. Y.
Li
,
B.
Chen
,
T. Y.
Xia
,
T. F.
Tang
,
B.
Zhu
, and
V. S.
Chan
,
Nucl. Fusion
59
,
126039
(
2019
).
16.
E.
Kaveeva
,
V.
Rozhansky
,
I.
Senichenkov
,
E.
Sytova
,
I.
Veselova
,
S.
Voskoboynikov
,
X.
Bonnin
,
R. A.
Pitts
,
A. S.
Kukushkin
,
S.
Wiesen
, and
D.
Coster
,
Nucl. Fusion
60
,
046019
(
2020
).
17.
F. L.
Hinton
and
R. D.
Hazeltine
,
Rev. Mod. Phys.
48
,
239
(
1976
).
18.
D.
Stotler
,
C. S.
Chang
,
S. H.
Ku
,
J.
Lang
, and
G.
Park
,
J. Nucl. Mater.
438
,
1275
(
2013
).
19.
I.
Nunes
and
on behalf of JET Contributors,
Nucl. Fusion
53
,
073020
(
2013
).
20.
See https://www.nutonian.com/products/eureqa/ for an AI-powered modeling engine by Nutonian; now acquired by DataRobot.20 
22.
M.
Schmidt
and
H.
Lipson
,
Science
324
,
81
(
2009
).
23.
D. R.
Hatch
,
M.
Kotschenreuther
,
S.
Mahajan
,
P.
Valanju
, and
X.
Liu
,
Nucl. Fusion
57
,
036020
(
2017
).
24.
P. H.
Diamond
,
S.-I.
Itoh
,
K.
Itoh
, and
T. S.
Hahm
,
Plasma Phys. Controlled Fusion
47
,
R35
(
2005
).
25.
M.
Schneider
,
A. R.
Polevoi
,
S. H.
Kim
,
A.
Loarte
,
S. D.
Pinches
,
J.-F.
Artaud
,
E.
Militello-ASP
,
B.
Beaumont
,
R.
Bilato
,
D.
Boilson
,
D. J.
Campbell
,
P.
Dumortier
,
D.
Farina
,
L.
Figini
,
Y.
Gribov
,
M.
Henderson
,
R. R.
Khayrutdinov
,
A. A.
Kavin
,
F.
Köchl
,
T.
Kurki-Suonio
,
A.
Kuyanov
,
P.
Lamalle
,
E.
Lerche
,
V. E.
Lukash
,
A.
Messiaen
,
V.
Parail
,
K.
Särkimäki
,
A.
Snicker
,
D.
Van Eester
, and
ITPA Topical Groups on Energetic Particle Physics and Integrated Operation and Scenarios,
Nucl. Fusion
59
,
126014
(
2019
).
26.
S. H.
Kim
,
R. H.
Bulmer
,
D. J.
Campbell
,
T. A.
Casper
,
L. L.
LoDestro
,
W. H.
Meyer
,
L. D.
Pearlstein
, and
J. A.
Snipes
,
Nucl. Fusion
56
,
126002
(
2016
).
27.
A. R.
Polevoi
,
A. A.
Ivanov
,
S. Yu.
Medvedev
,
G. T. A.
Huijsmans
,
S. H.
Kim
,
A.
Loarte
,
E.
Fable
, and
A. Y.
Kuyanov
,
Nucl. Fusion
60
,
096024
(
2020
).
28.
R. M.
Churchill
,
C. S.
Chang
, and
S.
Ku
,
IEEE Trans. Plasma Sci.
48
,
2661
2664
(
2020
).
29.
See https://spark.apache.org for description of Apache Spark, a unified analytics engine for large-scale data processing.
30.
C. S.
Chang
,
S.
Ku
,
P. H.
Diamond
,
Z.
Lin
,
S.
Parker
,
T. S.
Hahm
, and
N.
Samatova
,
Phys. Plasmas
16
,
056108
(
2009
).
31.
C. S.
Chang
,
S.
Ku
, and
H.
Weitzner
,
Phys. Plasmas
11
,
2649
(
2004
).
32.
L.
Wang
and
T. S.
Hahm
,
Phys. Plasmas
16
,
062309
(
2009
).
33.
S. V.
Novakovskii
,
C. S.
Liu
,
R. Z.
Sagdeev
, and
M. N.
Rosenbluth
,
Phys. Plasmas
4
,
4272
(
1997
).
34.
J.
Seo
,
C. S.
Chang
,
S.
Ku
,
J. M.
Kwon
, and
W.
Choe
,
Phys. Plasmas
21
,
092501
(
2014
).
35.
G.
Park
and
C. S.
Chang
, “
Numerical study of collisional ripple diffusion in a tokamak plasma
,”
Phys. Plasmas
10
,
4004
(
2003
).
36.
M.
Greenwald
,
Phys. Plasmas
17
,
058101
(
2010
).