Coupled KIPP-EDGE2D modeling of parallel transport in the SOL and divertor plasma for the ITER baseline scenario

The KInetic code for Plasma Periphery (KIPP) models the parallel (along magnetic field lines) propagation of charged particles in the scrape-off layer (SOL) and divertor of tokamaks. An iterative coupling between KIPP and a 2D edge fluid code, EDGE2D, which in turn is coupled to the Monte Carlo solver for neutrals, EIRENE, was used to achieve a converged KIPP-EDGE2D-EIRENE solution. In the iterative coupling algorithm, KIPP transfers kinetic parallel electron and ion heat conductivities to EDGE2D, whereas EDGE2D returns to KIPP 2D distributions of macroscopic plasma parameters across the computations grid. The initial EDGE2D-EIRENE solution simulated the SOL and divertor plasma of the ITER inter-edge localized mode (ELM) baseline scenario. This work employs the same methodology as an earlier study based on the analysis of JET high radiative H-mode conditions, with strong nitrogen injection leading to partial detachment at divertor targets (Chankin et al 2022 Plasma Phys. Control. Fusion 64 095007). The results are qualitatively similar to those of the JET study in demonstrating the strong heat flux limiting effect, together with the rise in electron and ion temperatures in the main SOL. At the same time, the kinetic effects of the parallel propagation of charged particles are not expected to drastically change the target power deposition at divertor targets calculated by EDGE2D-EIRENE alone.


JET Contributors
* Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
The KInetic code for Plasma Periphery (KIPP) is a continuum Vlasov-Fokker-Planck code for parallel plasma transport (along magnetic field lines) in the scrape-off layer (SOL) and divertor. The main equations, including the normalization of parameters, are described in [1]. The code combines an implicit second order scheme for a full nonlinear Coulomb collision operator with an explicit second order scheme for the parallel free-streaming. Results of the code benchmarking can be found in [2] and references therein. In the present paper, the code was run in the 1D2V mode, with one spatial coordinate (along several magnetic field lines) and two velocity variables: parallel and gyro-averaged perpendicular velocities.
The present paper is a continuation of an earlier study on the iterative coupling of KIPP with EDGE2D-EIRENE based on the JET high radiative H-mode inter-edge localized mode (ELM) plasma conditions with 8 MW of input power into the computational grid and strong nitrogen (N 2 ) injection leading to the partial detachment at divertor targets [3]. EDGE2D-EIRENE is a code package combining coupled plasma 2D edge fluid code EDGE2D and neutral Monte Carlo solver EIRENE [4][5][6].
Owing to high electron and ion temperatures (T e and T i ) upstream (along field lines in the direction away from the targets) generating high energy super-thermal charged particles with low collisionality, especially in the near SOL close to the separatrix, an impact of kinetics on plasma conditions at divertor targets could be expected. This effect was confirmed in the KIPP-EDGE2D modeling, but it was not strong enough to drastically change target plasma conditions. The dominant contribution to the outer target (OT) power flux (OT receives higher power flux than the inner one) was found to come from the total, electron plus ion, convection, which is adequately described by the (fluid) EDGE2D code, typically coupled with EIRENE (below, reference to EIRENE will be dropped when describing data exchange between KIPP and EDGE2D). Electron conductive power flux, which contains super-thermal electrons' power flux attributed to nonlocal effects of high energy electron propagation along field lines, was found to be the second largest. At the same time, the ion conductive power flux was found to be much lower, as is typically the case in EDGE2D modeling. The total target power flux density when including KIPP kinetic results increased only by a factor of ∼1.5 compared to that calculated by EDGE2D-EIRENE runs, across the whole target. It is noteworthy that kinetic effects of parallel charged particle propagation in the near SOL, with the highest power flux into the divertor, were found not to trigger a 'burnthrough' leading to a substantial increase in the target power load at locations where EDGE2D-EIRENE predicts partial divertor detachment. This was attributed to the presence of high density (and low temperature) plasma in the divertor at these locations, causing significant Maxwellization/thermalization of super-thermal charged particles by Coulomb collisions (see [3] and references therein).
In this paper, edge plasma parameters expected in the ITER baseline scenario are described in section 2. The setup of the EDGE2D-EIRENE case relevant to the ITER baseline scenario is described in section 3. Limitations of running EDGE2D-EIRENE cases and of the iterative KIPP-EDGE2D coupling algorithm are highlighted in section 4. Iterative KIPP-EDGE2D coupling results, including parallel KIPP and EDGE2D electron and ion power fluxes, are presented in section 5. Global power balance in iterative KIPP-EDGE2D cases is analyzed in section 6. The results of this work are compared with the previous study of the JET plasma in section 7. A summary and future plans are presented in section 8.

ITER parameters at the plasma edge used in the iterative KIPP-EDGE2D coupling
As a starting point of this work, the EDGE2D-EIRENE solution for the edge plasma covering the outer core, SOL and divertor regions was obtained. The 'standard' ITER parameters were adopted, also referred to as the 'ITER baseline scenario' [7], with D-T reactions generating ∼500 MW of fusion power (Q = 10) in 15 MA/5.3 T (q 95 ≈ 3) discharges, with power flux into the SOL P SOL ∼ 100 MW [8]. It was assumed that in order to avoid the destruction of divertor target plates by excessive heat loads, a large fraction of P SOL had to be radiated by injecting light impurities into the plasma. The physics basis for edge ITER parameters, focusing on divertor physics, was presented in [8]. A conservative approach with a safety margin on the maximum power load to divertor targets by limiting the power to ≈5 MW m −2 was adopted. This was to be achieved by injecting/puffing Ne, preferably, into the plasma. To predict divertor parameters in ITER, an output from the 2D fluid code SOLPS-4.3 coupled to EIRENE was widely used in [8].
SOLPS-4.3 simulations assumed, for simplicity, beryllium (Be) wall and divertor targets, instead of the more realistic ITER-like wall (ILW) in the present study, where the main chamber wall is covered with Be, while divertor targetswith tungsten (W). Ion species included in the plasma mix were deuterium, representing both deuterium (D) and tritium (T) isotopes, and helium (He). P SOL was fixed at 100 MW, assumed to be equally split between ion and electron channels. Fixed ion fluxes into the computational grid from its inner core boundary were imposed: 9.1 × 10 21 s −1 for D + ions and 2 × 10 20 s −1 for He 2+ ions, consistent with core fuelling and fusion power production at Q DT = 10. Spatially constant transport coefficients for plasma diffusion, D ⊥ = 0.3 m 2 s −1 , and ion and electron heat conductivities, χ i = χ e = 1 m 2 s −1 , were used. At the outer edge of the plasma grid, radial decay lengths of 3 cm for densities and temperatures of all charged species were set as boundary conditions at all poloidal locations.
The output from the SOLPS-4.3 code and similar 2D fluid codes, including key output parameters such as neutral particle recycling in the divertor and power load to divertor targets, is known to depend strongly on the neutral density, or pressure, in the divertor. This can be controlled by varying the D puff introduced as D 2 molecules at half the prescribed rate counted in the number of D atoms per second. In this study, we were aiming to approximately replicate SOLPS-4.3 plasma conditions obtained by the D puff = 1.85 ×10 23 s −1 and Ne puff = 2.1 ×10 20 s −1 . Both puffs were from the top of the machine. At such puff levels, the maximum power flux density onto divertor plates was slightly below 5 MW m −2 .
With the parameters specified above, SOLPS-4.3 simulations yielded a total radiated power from the computations grid equal to 56.6 MW, consisting of 41.3 MW Ne radiation and 15.3 MW D radiation. Electron density n e at the poloidal ring of the grid just outside of the separatrix was n e,sep = 4.5 ×10 19 m −3 . The ratio of the Ne impurity density n Ne , counting all charge states and averaged poloidally along the same poloidal ring from the X-point to X-point, by the separatrix electron density, was C Ne = n Ne /n e = 0.008, or 0.8%. We will compare the output from EDGE2D-EIRENE cases with that from SOLPS-4.3 in the following sections.

Setup of the initial EDGE2D-EIRENE case and its output
Before KIPP-EDGE2D coupling iterations begin, a well converged reference EDGE2D-EIRENE case, also referred to as the iteration 0 case, must be established. The computational grid for running EDGE2D cases adopted here, together with pumping surfaces in the divertor, was taken from earlier EDGE2D cases for modeling ITER plasmas in pure hydrogen, by Harting [9]. The grid, combining the EDGE2D grid for the plasma and the triangular EIRENE grid for neutrals, is shown in figure 1. The position of the pumping surface is indicated in red. This grid is similar to the one used in figure 3 of [8]. The number of cells was, however, different. In the grid used for EDGE2D cases there were ten poloidal rings in the core and in the private flux region (PFR), and 32 poloidal rings in the SOL. The number of cells/rows in the poloidal direction in the SOL, from one divertor plate (inner) to the other (outer), is 104. Anomalous transport coefficients for perpendicular/radial transport in the present EDGE2D modeling are as in the SOLPS-4.3 cases, as stated in the previous section: 0.3 m 2 s −1 for particle diffusion and 1 m 2 s −1 for ion and electron heat conductivities, constant across the computational domain.
The EIRENE model processes used by EDGE2D were a subset of the Kotov 2008 model used by SOLPS-4.3 (socalled NIMBUS-like neutral model). The main difference is that the Kotov model includes additional molecular processes which become important in detached conditions. The detailed description of the model, including its benefits in speeding up the runs, can be found in section B.4 of [10].
Similarly to the SOLS-4.3 modeling, drifts were not included in the present EDGE2D runs. The inclusion of drifts can change target asymmetries. However, in the limit of a highly radiating divertor with target power loads below 5 MW m −2 these changes are not large (see [11] with results of the ITER baseline scenario modeling with drifts using SOLPS-ITER).
No heat flux limiters were specified in the initial EDGE2D-EIRENE case described in this section. Thereafter, only heat flux corrections following from KIPP runs and imported into EDGE2D as multipliers for transport coefficients were used. Unlike in the SOLPS-4.3 cases described in [8], where the tungsten (W) wall was assumed, the true ILW was specified in the EDGE2D-EIRENE modeling, with Be covering the main wall and W the target plates. The plasma ion mix, at the same time, was different from that in the SOLPS-4.3 cases in that it does not include He, leaving only D out of the main ion species, but containing Be. The impact of Be on all plasma parameters, however, was found to be insignificant.
The absence of He from the plasma mix in the EDGE2D-EIRENE cases was partly motivated by the desire to reduce the central processing unit (CPU) time consumption. KIPP calculations used 6048 cores (processors) of the Marconi high performance computer (HPC). Since the most CPU consuming part of KIPP runs are Coulomb collision calculations, each spatial cell of the EDGE2D computational grid in the SOL and divertor regions was allocated two cores, one dedicated to Coulomb collisions for electrons and one dedicated to D + ions. Since Ne ions either have too low concentrations (for very low and very high charge states), or are strongly collisional (for ions with charges 4-6), Ne could be treated according to the fluid approach. Its densities were imported from EDGE2D-EIRENE solutions, while their impact on electronion (e-i) and ion-ion (i-i) collisions (for D + ions) was accounted for via the dependence of Coulomb collisions on plasma Z eff (for e-i collisions) or Ne (as well as Be) ion contributions to Rosenbluth potentials used to calculate transport coefficients in the velocity space for D + ions. The distribution functions of the impurity ions species were assumed to be drifting Maxwellians. The reduction of the number of species by a factor of three (from using D + , He + , H +2 species) to using only D + reduced the total CPU consumption by a factor of 2.5, when electrons are taken into account. This helped to satisfy limitations imposed on the CPU consumption, in node-hours, for using Marconi in KIPP calculations.
The input power into the computational grid from the inner core boundary was 100 MW, equally split between ion and electron channels, as in the SOLPS-4.3 cases. In the initial EDGE2D-EIRENE case discussed here ('iteration 0') high recycling conditions in the divertor and acceptable maximum power loads on divertor targets, below 5 MW m −2 , were achieved by the 1.2 ×10 23 s −1 D puff from the top of the machine and 5.5 ×10 20 s −1 Ne puff from the position below the outer midplane (OMP) (in the SOLPS-4.3 case referred to in the previous section, D puff = 1.85 ×10 23 s −1 is higher, while Ne puff = 2.1 ×10 20 s −1 -lower than in the EDGE2D case). In steady state conditions of the EDGE2D-EIRENE run both D input into the grid, which also included 1 ×10 22 s −1 D flux through the inner core boundary responsible for core fuelling, and the Ne puff, were matched by fluxes to the pump. Electron and D + densities at OMP were 3.98 ×10 19 m −3 and 3.44 ×10 19 m −3 , respectively. Radiated power on Ne was 49.5 MW, and on D it was 8.5 MW, giving the total radiated power of 58 MW in the whole computational grid. Finally, the impurity (mostly Ne) relative concentration C Ne = n Ne /n e , an important parameter for coupling SOL and divertor solutions to core codes and defined in the previous section, was 1.14%, higher than in [8]. Figure 2 shows radial profiles of T e , T i and electron density n e across cells at the OMP position for the 0th (initial) and second iteration EDGE2D-EIRENE cases. See section 5 for numbering of the iterations. Higher T i in the core region is attributed to lower ion density owing to the presence of ion Ne impurity species, since for the same input power into the grid through the ion channel as through the electron channel the ion power flux is carried by fewer ions than electrons. In the SOL, a lower T e is caused mostly by a much higher parallel electron heat conduction (towards the target plate) than ion.
Target plate profiles of the total power flux density, plasma temperatures T e and T i and electron density n e , covering the  For poloidal rings outside of these regions parallel power flux corrections attributed to kinetic effects were not introduced into EDGE2D, hence original EDGE2D coefficients were used in EDGE2D-EIRENE runs. The reason for such a restriction will be explained in the next section.
At the OT, which receives higher power flux, the maximum target power flux density did not exceed 4.5 MW m −2 , being substantially below this level at the strike point, indicating partial detachment at the target. T e and T i at both inner and outer strike points did not exceed 1 eV.

Limitations of running EDGE2D-EIRENE cases and of the iterative KIPP-EDGE2D coupling algorithm
Owing to numerical instabilities occurring in detached plasmas with a large number of radiating impurities in EDGE2D-EIRENE runs, in all cases a small time step of 5 × 10 −7 s had to be used, making the total duration of each EDGE2D-EIRENE run longer than one month, until the case reaches steady state conditions. The duration of the whole modeling work was mainly limited by the high real time consumption of EDGE2D-EIRENE cases, whereas KIPP calculations took much less time.
Regarding KIPP runs, the main limitation was on the allocated number of node hours for running the code on the Marconi HPC. As was pointed out in section 3, this in particular was the motivation for not using He in the ion mix in both EDGE2D-EIRENE and KIPP runs. Dimensionless collisionality of He + and He 2+ is not much different from that of D + ions, hence these ions would need to be treated kinetically by KIPP, doubling the number of cores (processors) requested for Coulomb collision calculations. Parallel free-streaming of charged particles in KIPP takes much less CPU time than Coulomb collisions. In order to further reduce the CPU time consumption, KIPP calculations were not performed to model kinetic effects on outer poloidal rings of the EDGE2D grid, where power flux density to divertor targets was negligible. Out of 32 available poloidal rings of the EDGE2D grid, KIPP calculations were performed only on the first 21 rings counting from the separatrix position.
Unlike in earlier JET cases [3], where the number of cells along the parallel/poloidal direction in KIPP calculations had to be doubled to obtained smooth, less oscillatory solutions, in this study KIPP used the same grid as EDGE2D for the placement of cell centers and cell faces along field lines. This was due to smoother parallel macroscopic plasma parameter profiles obtained in EDGE2D-EIRENE solutions for ITER compared to in [3]. And this prompted renewed attempts to extend the iterative coupling onto divertor regions, unlike in [3], where such a coupling was applied only to the main SOL plasma, outside of divertors. Such attempts were, however, unsuccessful due to large jumps in derivatives of plasma parameters. The problem may be attributed to an inconsistency between KIPP and EDGE2D grids, leading to an inconsistency between calculated fluxes. In KIPP, calculated parallel fluxes are assumed to be along flux tubes with their axes passing through EDGE2D cell faces. At the same time, the EDGE2D grid is 2D in the poloidal cross section, consisting of quadrilaterals which may, especially near the X-point position, strongly deviate in shape from rectangles. Since in EDGE2D cases cell face fluxes are integrated along cell face boundaries, different algorithms are used to calculateKIPP and EDGE2D fluxes through cell faces. The more nonrectangular shapes of EDGE2D cells are, the larger the difference between KIPP and EDGE2D calculated fluxes may become. The difference is particularly large not only at cells close to the X-point, as pointed out above, but also at some cells in the divertor, especially near targets where EDGE2D cells may be strongly nonrectangular. It has to be pointed out that the iterative coupling algorithm with KIPP and SOLPS codes was successfully tested using a 1D grid in [12].
Apart from the possible causes of inconsistencies between KIPP and EDGE2D fluxes mentioned above, there may also be other inconsistencies related to EDGE2D fluxes themselves. In figure 6 of [3] and the surrounding text, an example of unphysical features of the total parallel power flux from EDGE2D is presented, with the nonmonotonic power flux near the entrance to the inner divertor and an unexpected reduction of the derivative of the power flux near the entrance to the outer divertor.
As a result of all the above problems, the KIPP-EDGE2D iteration algorithm in this study was only applied to the main SOL region, outside of divertors, as in the earlier JET study [3]. For the divertor plasma, original EDGE2D transport coefficients for parallel heat fluxes, without corrections for kinetic effects, were used.

Iterative KIPP-EDGE2D coupling results
In this study, KIPP uses the spatial grid generated for EDGE2D with the same number of cells. As in the earlier study [3], in the parallel/poloidal direction KIPP cells are counted from the inner target (IT) to the outer target (OT). The same cell numbering is used when presenting EDGE2D-EIRENE output, despite the fact that, internally, EDGE2D uses an opposite cell numbering, from OT to IT. Poloidal rings of the EDGE2D grid will be referred to according to their names in EDGE2D, beginning with 's01 ′ for the ring just outside of the separatrix, up to the ring 's21' deep in the SOL (as was pointed out earlier, KIPP calculations do not cover the whole EDGE2D grid in the SOL, which consists of 32 rings).
As pointed out above, the initial EDGE2D-EIRENE case here is referred to as the 0th iteration EDGE2D (or EDGE2D-EIRENE) case. The following KIPP case uses macroscopic plasma parameters from this EDGE2D case maintained by the thermal source and after reaching steady state its output is referred to as the 0th iteration KIPP case. After that, effective electron and ion parallel heat conductivities are extracted from the KIPP solution in the main SOL, divided by those from EDGE2D, and two 2D matrices of the ratios, for χ e and χ i , are transferred to EDGE2D before the start of the run to serve as multipliers for EDGE2D heat transport coefficients. The next EDGE2D-EIRENE case, based on the KIPP-derived multipliers for heat conductivities, after reaching the steady state, is referred to as the first EDGE2D case. It is followed by the next, first iteration KIPP case, and so on.
Zeroth iteration EDGE2D-EIRENE and KIPP cases differ in parallel conductive power flux profiles for electrons and ions in the main SOL owing to the kinetic effect of 'heat flux limiting', with kinetic conductive power fluxes being lower than the ones following from EDGE2D, which are close to Braginskii's heat fluxes analytically derived for strongly collisional plasma conditions. This is a well known effect encountered when kinetic code modelers compare their results with those of 'fluid' (using transport coefficients derived for strongly collisional plasmas) codes, see e.g. review paper [13]. This effect is caused by insufficient collisionality of SOL plasmas where super-thermal charged particles responsible for the bulk of the heat flux escape 'hot' plasma regions by moving along field lines without being sufficiently thermalized, leaving a 'hole' in the high energy part of the distribution function resulting in lower heat fluxes. For the same input power into the computational grid lower parallel power fluxes result in higher peak T e and T i along field lines.
It should be noted that, since only D + ions were modeled kinetically in KIPP, whenever parallel or perpendicular (onto the target) 'ion' power fluxes from KIPP are quoted or presented in figures, by the 'ion' fluxes one should understand D + fluxes only. Consequently, in order to compare similar quantities between EDGE2D-EIRENE and KIPP solutions, D + only power fluxes will always be implied as 'ion' fluxes extracted from EDGE2D, unless otherwise stated.
Both heat flux limiting and rising peak T e,i values were obtained in the earlier study [3]. Zeroth iteration parallel EDGE2D-EIRENE heat fluxes were substantially above those from 0th iteration KIPP cases. The fluxes became closer to each other in subsequent iterations. Only two KIPP-EDGE2D iterations were found to be sufficient to achieve a good match between KIPP and EDGE2D heat fluxes in [3] (see e.g. figures 3 and 4 in this reference). The same is found in the present, ITER related study. However, in this study the third iteration EDGE2D-EIRENE case underwent a 'thermal collapse' in the outer divertor, apparently caused by insufficient power flux to the divertor, where lower divertor T e leads to higher n e (as plasma pressure tends to reach equilibrium along field lines), increasing radiative power losses and resulting in turn in even lower T e , and so on. Such transitions were not mentioned in the SOLPS4.2 modeling of ITER plasmas in [8] despite a large number of cases with large variations of input parameters were run leading to a wide range of n e,sep . Plasmas with an enhanced radiation from divertors may, however, also be reactor relevant. Highly radiative SOLPS-ITER runs were recently carried out in the study [14] aiming at exploring ITER operation at a higher degree of detachment. In this study, by increasing D and Ne puff rates proportionally, the peak divertor energy flux decreased from 5 to 3 MW m −2 , with the position of the impurity radiation peak shifting towards the X-point. At the highest neon puff levels with steady state solutions, T e is reduced below 1 eV across 50 cm of each divertor target.
Parallel profiles of T e , n e and electron conductive power fluxes q cond e for iterations 0-3 of EDGE2D-EIRENE cases for the poloidal ring just outside of the separatrix, ring s01, are plotted in figure 4. The 'thermal collapse' affected global plasma parameters such as total radiated power (increased) and total power to divertor targets (decreased). The fourth iteration of the KIPP-EDGE2D coupling was also performed, but it produced almost the same output as the third iteration, indicating that the new plasma state is stable. Its results are not plotted in figure 4. Changes in the profiles near the entrance to the outer divertor and in the divertor itself between the second and third iteration cases shown in figure 4 make these profiles qualitatively similar to those seen at corresponding positions at the inner divertor side. The plasma in the inner divertor is usually 'colder' (lower T e,i ) and denser (higher n e ) than in the outer divertor, with a higher degree of detachment. Analysis of the profiles at other poloidal rings shows that the transition to the new plasma state in the outer divertor is localized to near SOL rings, closer to the separatrix position, and the profiles' changes seen in figure 4 almost completely disappear at ring s10. Maximum T e (eV) values along poloidal ring s01 are 154.1, 166.9, 171.1, 170.9 and 172.4 for iteration numbers 0-4, confirming that the transition in the outer divertor state is a localized effect that does not result in significant changes of upstream parameters.
T e and q cond e profiles for the third iteration case in figure 4 deep in the main SOL are almost superimposed onto the second iteration profiles. Since plasma conditions strongly changed between the second and third iteration cases, this alone cannot serve as a proof that the quasi steady state plasma conditions were reached already at the second iteration, so that it can be used as a case characteristic of the converged KIPP-EDGE2D iterative solution. There is, however, strong evidence that such conditions were indeed achieved already after only two iterations, as was the case in the early study of JET plasmas [3]. This evidence is considered next.
As was pointed out earlier in this section, there are two closely related features of parallel plasma parameter profiles caused by the kinetic effect of the heat flux limiting. One is the reduction of parallel heat fluxes of electrons, q cond e , and ions, q cond i , and the other is the rise in peak T e,i values along field lines, with the temperatures' rise being the consequence of the reduction in heat conductivities. Parallel conductive ion power fluxes, together with ion temperature, from 0th and second iteration EDGE2D-EIRENE and KIPP cases are shown in figures 5(a) and (b), respectively, along field lines belonging to the poloidal ring s01. The choice of ion, rather than electron power fluxes, is motivated by much lower dimensionless ion collisionality in the SOL compared to electron, due to higher upstream ion temperature, see figure 2, resulting in a substantially stronger heat flux limiting effect for ions. As can be seen in figure 5(a), in regions of large T i gradients, ion conductive power flux calculated by KIPP, q cond,KIPP i , is substantially lower than that calculated by EDGE2D, q cond,EDGE2D i . The difference in the shapes of the two power profiles is attributed to nonlocal kinetic effects of the ion heat flux transport in KIPP which is less dependent on local ion temperature gradients ∇ || T i . After two KIPP-EDGE2D iterations, the difference between q cond,KIPP i and q cond,EDGE2D i in the main SOL region, outside of divertors, is almost eliminated, as can be seen in figure 5(b). This good match is obtained by using KIPP/EDGE2D multipliers for χ i down to ∼0.25 in the high ∇ || T i zone close to the inner divertor, and ∼0.2-close to the outer divertor. Hence, two iterations are sufficient to achieve a good match between q cond,KIPP i and q cond,EDGE2D i , indicating that a quasi steady state plasma condition was reached. Note that the good match between the two power fluxes was achieved mainly by a significant drop in q cond,EDGE2D i , while q cond,KIPP i has not changed much. The much lower q cond,EDGE2D i in figure 5(b) than that shown in figure 5(a) is partly compensated by the higher convective power flux. Note also that q cond,EDGE2D i is much less than q cond,EDGE2D e , see below. A good match between conductive power flux density profiles in KIPP and EDGE2D was achieved at all poloidal rings for both electrons and ions.
A figure similar to figure 4 in this paper, but for ion, rather than electron, temperature and conductive power flux, can be found in [3] (figure 4), where gradual changes in the parameters occurred in iterations 0-5. From this figure, it follows that the largest changes occur between iterations 0 and 1. The changes can be best seen and resolved when one analyses parallel T i profiles. The T i rise between iterations 1 and 2 is by a factor of 3-4 smaller than between iterations 0 and 1. The last three iterations, from the second to the fifth, cause T i changes which are comparable to the changes between the first and second iterations, indicating a sharp slowdown in T i changes as the iteration number is increased. The earlier, JET results presented in [3] are relevant to the discussion about whether enough iterations were made and whether the coupled quasi steady state conditions were reached since, as will be discussed below in section 7, ion and electron dimensionless collisionalities at the separatrix position are very close to each other in JET and ITER studies.
It must be pointed out that the improvement in the understanding of the physics of the kinetic effects of super-thermal charged particle propagation and Coulomb collisions that can be achieved in the earlier [3] and present studies is based on running EDGE2D-EIRENE cases with fixed perpendicular (radial) ad hoc anomalous transport coefficients. Therefore, some potentially important features related to the impact of super-thermal charged particles on perpendicular transport coefficients may be missed, since they can only be described by turbulence codes with the inclusion of kinetic effects.
Next, we will analyze T e , T i , n e and electron and ion parallel power fluxes in the 2nd iteration EDGE2D-EIRENE and KIPP cases for three positions in the SOL: at ring s01, with the highest peak T e , T i values in the main SOL, ring s07 with the highest power flux to the target, and ring s15 located deep in the SOL, with relatively low plasma temperature, density and target power fluxes. Results on target power fluxes will be presented only for the OT divertor, which receives higher power flux than the IT. Also, the plasma in the inner divertor has a higher degree of detachment, and parallel T e,i profiles from EDGE2D-EIRENE solutions exhibit large, order of magnitude cell-to-cell drops across the entrance to the inner divertor, which makes both EDGE2D-EIRENE and KIPP results less reliable.  i profiles in most of the main SOL region can be seen. Unlike for electrons, q conv i is large and has an opposite direction to conductive power fluxes in most of the main SOL, which is caused by ion flows driven by neutral ionization sources. A similar 'flow reversal' can be seen in figure 13 of the JET study [3], except that there this feature was seen deep in the SOL at ring s11. Ion power fluxes are much lower than electron, so the total power flux in the main SOL is in the expected direction-towards the closer target. Power fluxes at the OT cannot be well resolved in figures 6(a) and (b). They will be plotted separately in section 5.4 for both iterations 0 and 2 for all positions (poloidal rings) along OT.

Profiles along ring s07
Figure 7(a) shows the same quantities as in figure 6(a), but for ring s07. A very good match between q cond,KIPP e and q cond,EDGE2D e in the main SOL and parts of the divertor regions can be seen. The match between the power fluxes in divertors may, however, be incidental since KIPP/EDGE2D multipliers were not applied to EDGE2D cells in the divertor, as was pointed out in the previous section. Figure 7(b) shows the same quantities as in figure 7(a), but for ions. As at the ring s01 ( figure 6(b)), ion conductive power flux is directed towards the inner target almost across the whole flux surface, except for positions in the outer divertor close to the target.
In figure 7(b), as well as in other figures beginning with figure 6(a), a certain mismatch between conductive power fluxes from EDGE2D-EIRENE and KIPP in the middle of the main SOL can be seen. It occurs near the stagnation point with maximum T e,i . where EDGE2D conductive power fluxes change sign, since they scale approximately as -T 5/2 ∇ || T. As was explained in [3], the difference between the nature of KIPP (nonlocal) and EDGE2D (local) conductive power fluxes often leads to different signs of the two fluxes near the stagnation point for T e,i . In order to match the two power fluxes, negative KIPP/EDGE2D multipliers must be transferred to EDGE2D for the calculation of coefficients χ e and χ i . EDGE2D, however, does not allow positive coefficients before ∇ || T for conductive power fluxes. For this reason negative multipliers were not applied to such cells, that is, they were set to 1, impeding convergence between the two power fluxes in these cells.  One feature of T i profiles shown in figures 5-8(b) is their strong variation in the main SOL, which is normally not seen in EDGE2D-EIRENE cases. This is explained by the stronger ion heat flux limiting effect at the outer (low toroidal field) side of flux surfaces with lower multiplication factors for q cond,EDGE2D i compared to standard (without the use of KIPP/EDGE2D multipliers) χ i values calculated internally by EDGE2D. Since the ion parallel conductive power flux is relatively small compared to that for electrons, and of the same order of magnitude as its convective flux, spatial variation of heat flux limiting factors may cause significant variation of T i along field lines. Note that the degree of T e variations in the main SOL, as seen in figures 6-8(a), is much lower than that of T i .

Power fluxes at OT
Contributions of electron and main ion conductive and convective parallel power fluxes at the OT, normal to the target surface, with the exception of the ion conductive power flux from EDGE2D, q cond,EDGE2D i , which is zero, for the 0th iteration of both EDGE2D-EIRENE and KIPP runs, are plotted in figure 9(a) Ion conductive power flux from KIPP, q cond,KIPP i , is not much lower than q cond,KIPP e , except for positions near the strike point. At the strike point position, at ring s01, the two power fluxes are almost equal to each other. This can be attributed to much higher peak T i than T e upstream in the main SOL, with T i,max = 281 eV and T e,max = 154 eV, and a strong dependence of nonlocal transport of super-thermal charged particles on their energy.
Even though q cond,KIPP e seems relatively low at the strike point, it is a factor ≈70 higher than the corresponding EDGE2D-EIRENE power flux, q cond,EDGE2D e . Compared to the q cond,KIPP e maximum at ring s07, such a low power flux at the strike point should be attributed to the effect of strong thermalization (or Maxwellization) of super-thermal electrons by Coulomb collisions at ring s01 due to much higher n e and longer connection length (caused, in particular, by a high magnetic flux expansion near the X-point) along this ring. Calculations of sums of products of cells' n e values by cell sizes along field lines in the outer divertor show that this quantity is larger by a factor of six at ring s01 than at ring s07, making Coulomb collisions' thermalization of high energy electrons much more efficient along ring s01.
The same quantities as shown in figure 9(a) are plotted in figure 9(b) for the second iteration KIPP and EDGE2D-EIRENE cases. Almost all power fluxes, except for the ion conductive power flux from KIPP, q cond,KIPP i , are lower than in figure 9(a), as expected from the reduced parallel power flux to the outer divertor caused by the effect of higher flux limiting. The higher peak value of q cond,KIPP i compared to q cond,EDGE2D e in the second iteration case around ring s15 may be explained by a particularly large T i rise between these two iterations, with the peak T i value in the main SOL increasing from T i,max = 134.6-172.1 eV (28% increase), compared to a much more moderate, only by 7.5%, rise in peak T e values, from 61.7 to 66.3 eV.
At more outer SOL positions, beginning approximately with ring s18, the recycling level in the divertor drops causing a substantial target n e drop, see e.g. the n e profile in figures 3, 3 points outward from the s15 ring, causing in turn the rise in the target T e . As a result, EDGE2D electron conductive power flux, q cond,EDGE2D e , increases and becomes comparable to the KIPP flux q cond,KIPP e . Radially (in the direction towards the wall) at the OMP position T 2 e drops faster than n e , resulting in the increase in dimensionless electron collisionality which scales as n e T -2 e against electron density and temperature. The plasma becomes more collisional, and KIPP and EDGE2D target electron power fluxes become closer to each other.
One comment on the difference between KIPP and EDGE2D conductive power fluxes and how it may impact the power balance should be made. KIPP power fluxes in the divertor may be considered as results of only one iteration in assessing the effect of kinetics, applied to EDGE2D-EIRENE runs of different iterations. KIPP results in the divertor region are not transferred back to EDGE2D for each following EDGE2D-EIRENE iteration. Owing to the absence of this feedback effect, KIPP power fluxes do not have to comply with the global power balance calculations which are here based on EDGE2D-EIRENE results (see next section). Therefore, the KIPP results should rather be considered as an indication of the potential strength of kinetic effects in modifying target power fluxes. It is important in this respect that KIPP calculations predict a fairly small impact of kinetic effects on the total target power load for positions close to the strike point, where upstream T e,i in the main SOL are the highest. At the same time, further out in the SOL, at positions near the maximum target power load, a possible increase of the total power load by factor of ∼1.5 attributed to electron kinetic effects follows from figure 9(b).
It must be pointed out that, similarly to the earlier study related to the JET plasma [3], in the iterative KIPP-EDGE2D-EIRENE modeling described here, non-Maxwellian features of the electron distribution function in the calculation of rate constants for atomic physics processes contributing to particle and power sources were not taken into account. As described in section 5 of [3], the presence of high energy super-thermal electrons in the divertor causes enhanced excitation and ionization rates due to electron collisions with both neutrals and ions, leading to enhanced volumetric power losses and a higher degree of detachment [15].

Global power balance in iterative EDGE2D-EIRENE cases
Sums of electron and ion (D + ), conductive and convective, power fluxes at the OT for all EDGDE2D-EIRENE iteration cases are plotted in figure 10. As one can see, the maximum power flux in the 0th iteration EDGE2D-EIRENE case is substantially below the maximum of the total power flux at the OT shown in figure 3. The difference is mostly attributed to the inclusion of the recombination power flux in figure 3. The absence of power fluxes by impurities is less important, since near the target impurity densities in EDGE2D-EIRENE cases are low due to the effect of the main ion-impurity thermoforce driving impurity ions away from the target.
The impact of inclusion of kinetic effects on the reduction of the target power load due to the iterative KIPP-EDGE2D coupling in the main SOL is well seen in figure 10. The third iteration EDGE2D-EIRENE solution is quite stable, as the fourth iteration does not result in any appreciable change in the power flux profiles. For the same input power into the numerical grid, the reduction of the target power flux is balanced by higher volume radiation in the divertor. The total number of electrons in the outer divertor, calculated by summing up multiples of n e by cell sizes (see section 5.4), is increased by 16% from the 0th to second iteration, and further increases by 34% from the second to the third iteration. The total number of electrons in other regions: inner divertor and main SOL, does not show such a steady rise as in the outer divertor: from the 0th to first iteration, these quantities increase by 12% and 3.1%, respectively, showing small drops from the first to second iteration, and more significant drops, by 5.9% and 3.5%, for the inner divertor and main SOL regions, respectively. The total electron content in the whole SOL region, including the divertor, is increased by 5.7% from the 0th to the second iteration, and then by 1.8% from the second to the third iteration.
Volumetric power losses, as well as target power loads, can be extracted from 'print' files formed after the completion of each EDGE2D-EIRENE run. The results are presented for 'marco-zones': core, main SOL, outer and inner SOLdivertor (divertor regions excluding PFR), and outer and inner PFR regions, for volumetric power losses, and power fluxes at boundaries between them, including the separatrix, wall, target plates and boundaries between the main SOL and inner and outer divertors. The reduction of total power fluxes to OT with an increase in the iteration number, seen in figure 10, is matched by a steady increase in the total volumetric radiation in the outer SOL-divertor region from 27.3 to 32.6 and 40.3 MW for iterations 0, 2 and 3, respectively. Volumetric power losses include Ne and 'hydrogen' (D) radiation, atomic and molecular dissociation, and charge-exchange power losses. Similarly to the situation with the number of electrons in the main SOL and SOL-divertor regions discussed above, volumetric power losses in the main SOL and inner SOL-divertor regions do not exhibit a consistent trend with iteration numbers, except for a ≈13% increase between the second and third iterations.
Finally, the figure of merit of the relative impurity density at the separatrix, parameter C Ne discussed in section 2, is 1.14%, 1.07%, 1.03%, 1.28% and 1.28% in iterations 0, 1, 2, 3 and 4, respectively. These values are higher than in the earlier SOLPS4.3 modeling of the ITER baseline scenario, ∼0.8%, for comparable plasma conditions in the SOL used in this study.

Comparison with earlier results on the JET plasma
The EDGE2D-EIRENE cases in the earlier JET study [3] and the present ITER study have different radiating impurities (N in JET and Ne in ITER) and very different input powers into the EDGE2D grid. In [3], P input was 8 MW, while in this study it is 100 MW. Such a large difference in P input makes JET cases underpowered compared to ITER cases, regardless of possible scalings (P input /R or P input /R 2 , where R is the major radius; in ITER the major radius is a factor of two larger than in JET) that may be used to compare expected dimensionless plasma parameters. It may therefore seem surprising that, in many respects, the results obtained in [3] and in this study are qualitatively similar.
The two studies, however, have very similar dimensionless collisionalities upstream in the main SOL plasma. In the JET study [3], separatrix T e and n e at the OMP position along poloidal ring s01 for the 0th iteration EDGE2D-EIRENE cases were 88.8 eV and 2.44 × 10 19 m −3 , while in ITERthey were 153 eV and 3.98 × 10 19 m −3 . The field line lengths from inner to OT are L || = 106.4 m and 185.4 m for JET and ITER, respectively. Using the approximate scaling for dimensionless electron collisionalities ν * e ∝ n e L || /T 2 e gives for ITER and JET ν * e almost the same values, with the ν * e,ITER /ν * e,JET ratio being ≈0.96. For ions, using T i instead of T e (160 eV for JET and 277 eV for ITER) gives ν * i,ITER /ν * i,JET = 0.95. Electron collisionality also depends on electron-impurity collisions via Z eff (linear dependence). At OMP positions, Z eff in JET and ITER cases are 1.93 and 2.48 for JET and ITER EDGE2D-EIRENE cases at rings s01, respectively. If these numbers were taken into account (more precise numbers would need to take into account Z eff distribution along field lines) the ν * i,ITER /ν * i,JET would rise to 1.22, making the ITER case even more collisional than the JET one. Since for plasma kinetic effects involving parallel propagation and Coulomb collisions of charged particles ν * e values are of prime importance, given the above estimates for dimensionless collisionalities, similarities between the results of [3] and the present study can indeed be expected.
Owing to the same technique being applied to the iterative KIPP-EDGE2D coupling algorithm, in both JET and ITER two iterations were sufficient to reach a quasi steady state solution with a good match between parallel conductive power fluxes q cond,KIPP e,i and q cond,EDGE2D e,i . In the case of JET, further iterations, from third to fifth, did not significantly change the plasma state. In contrast to JET cases, after the second iteration the ITER EDGE2D-EIRENE case has undergone a transition to a new plasma state with higher density, lower temperature and higher volume radiation in the outer divertor covering poloidal rings closer to the separatrix (up to ring s10, approximately), as was pointed out in section 5. The impact of this transition on the OT power flux profile can be seen in figure 10. The formation of a higher n e , lower T e,i and higher radiation zone in the outer divertor has also reduced T e in cells adjacent to the entrance to the outer divertor, see the third iteration T e profile for the ring s01 in figure 4 as an example. The transition made plasma parameter profiles at entrances to divertors and inside divertors more in-out symmetric. The new plasma state is characterized by increased impurity content in the plasma, which is evidenced by higher C Ne at the separatrix, as pointed out in section 6.
In both JET and ITER studies, super-thermal electrons and ions originating from regions with highest T e,i upstream in the main SOL were largerly thermalized/Maxwellized by the time they reached OT on poloidal rings closer to the separatrix. This can be seen in figures 9(a) and (b) where the target q cond,KIPP e is much lower than its maximum value some distance away from the separatrix. Also, in both JET and ITER studies further out in the SOL, where a maximum in target q cond,EDGE2D e is reached, q cond,KIPP e significantly exceeds q cond,EDGE2D e . The discrepancy between these two power fluxes is a result of the absence of a feedback loop between KIPP and EDGE2D-EIRENE iterations in the divertor, leading to violations in conservation laws. Similarly to large discrepancies between parallel fluxes q cond,KIPP e,i and q cond,EDGE2D e,i in the main SOL between 0th iteration KIPP and EDGE2D-EIRENE cases, the target power flux discrepancies are indicative of a possible potential contribution of kinetic nonlocal transport of high energy charged particles (mostly electrons) to the increase in the total target power load. In both JET and ITER studies, however, the potential contribution of the non-local kinetic transport to the target power load calculated by EDGE2D-EIRENE can increase the power load by a factor of less than two, once large contributions of convective electron and ion power fluxes are taken into account. And this is probably a high-end estimate, since in both JET and ITER studies power losses by superthermal energetic electrons and ions on inelastic collisions with neutrals and impurities were not taken into account.
Further out in the SOL, in the region of low OT power load, both JET and ITER studies predict small differences between q cond,KIPP e and q cond,EDGE2D e . This can be explained by smaller T e excursions along field lines together with the impact of the heat flux limiting effect. As for the conductive power flux q cond,KIPP i , it cannot make any significant contribution to the OT power load, even though it is by orders of magnitude larger than q cond,EDGE2D i and may become comparable to q cond,EDGE2D e at some locations along the target.

Summary and future plans
An iterative coupling between KIPP and EDGE2D-EIRENE codes, earlier applied to JET high radiative H-mode conditions with strong nitrogen injection [3], was extended to the ITER baseline scenario, with Ne injection reducing the target power load. These studies provide scientifically-based kinetic downward corrections to parallel heat transport coefficients (transferred from KIPP to EDGE2D), replacing rather arbitrarily defined 'heat flux limiters' commonly used in 2D edge fluid codes in SOL regions with high T e,i where fluid heat transport coefficients, derived in the approximation of a strongly collisional plasma, may lead to unrealistically high power fluxes. Despite being successfully tested in 1D geometry [12], due to the presence of strongly nonrectangular cells in 2D EDGE2D grids based on realistic magnetic equilibrium, and large cellto-cell plasma parameter variations close to the X-point, the iterative coupling algorithm could not be applied to positions at entrances to and inside divertors, neither in JET nor in ITER cases. As a result, kinetic corrections were not applied in these regions. One step in improving the situation in this respect will be the use of a much finer spatial grid in divertor regions, eliminating at least the problem of large cell-to-cell variations of plasma parameters in the poloidal direction. The maximum allowable number of cells in the poloidal direction in the SOL, including divertors, in EDGE2D is limited to 154, while the number of cells in the SOL and divertor in the poloidal direction in the grid used in the present study is 105, including 27 × 2 = 54 cells in the two divertors. There is therefore some room for improvement of the spatial resolution for coupled KIPP-EDGE2D studies, especially if the generated grid will have more cells around the X-point and e.g. only in the outer divertor where the plasma is typically less collisional than in the inner divertor.
One other planned extension of the present work is to focus on the heat flux limiting effects in the main SOL region, analyzing results of the iterative KIPP-EDGE2D coupling of a number of cases with different upstream dimensionless collisionalities and parameterizing the results in terms of the dependence of flux limiters on e.g. upstream dimensionless collisionalities. This could provide scientifically-based scalings for ion and electron heat flux limiters to be incorporated into today's 2D edge fluid codes. The focus on the 'hotter' main SOL regions, while using generic EDGE2D parallel heat flux conductivities in the divertor, as was done in the present study, will allow the use of substantially larger time steps in KIPP calculations owing to much lower Coulomb collision frequencies in the main SOL compared to the divertor plasma, releasing node-hours of the allocated CPU time to cover a larger number of coupled KIPP-EDGE2D solutions.
Despite the input power into the EDGE2D grid in the JET being a factor of 12.5 lower compared to the ITER case, upstream collisionalities just outside the separatrix were extremely close to each other. Similarities between the results of the earlier JET and the present ITER studies confirm that the upstream dimensionless collisionality is, as expected, a figure of merit to ascertain the role of kinetic effects in the SOL. In both studies, only two iterations (two transfers of KIPP corrections to EDGE2D) were sufficient to reach EDGE2D-EIRENE solutions, which were close to steady state conditions. In these plasma conditions, a good match between KIPP and EDGE2D conductive parallel power flux profiles in the main SOL, q cond,KIPP e,i and q cond,EDGE2D e,i , were achieved. In the present ITER study, progressive iterations led to the reduction of q cond,EDGE2D e,i in the main SOL and in the divertor regions. In the case of the outer divertor, which receives a higher power flux than the inner divertor, q cond,EDGE2D e,i reduction is accompanied by an increase in n e , a reduction in T e,i and an increase in volumetric power losses, resulting in a lower target power load.
In the present study, unlike in the earlier JET study, during the third iteration the outer divertor plasma has undergone a transition into a higher density, lower temperature and higher volumetric power loss state with substantially lower target power load. The new plasma state was stable and unaffected by the fourth iteration. The transition was confined to the part of the divertor within ∼1 cm (when mapped to the OMP) distance outside the separatrix. In the new plasma state, the total Ne impurity ion concentration (the sum of all Ne +Z ionization states) at the separatrix is higher.
KIPP cases were run after each EDGE2D-EIRENE case across the SOL and divertor regions. But, since kinetic corrections were not passed onto EDGE2D in the divertor, the iterative coupling did not take place there. KIPP runs were performed using plasma parameter profiles from EDGE2D-EIRNE cases maintained by thermal sources. The absence of a feedback from KIPP to EDGE2D meant that parallel power fluxes extracted from KIPP could violate the global power balance. The fact that the plasma in the divertor is almost always in the high collisionality regime does not justify the noninclusion of KIPP results into the iterative feedback loop because of strongly nonlocal transport of high energy super-thermal charged particles originating in the main SOL and reaching divertor targets. KIPP results in the divertor should therefore be considered only as indicative of the potential strength of kinetic effects obtained in stand-alone cases which trace the behavior of super-thermal charged particles against the background of the thermal plasma.
KIPP results for the target power deposition in the outer divertor indicate that, due to high density in the divertor at field lines near the separatrix, thermalization/Maxwellization of high energy super-thermal electrons reduces nonlocal kinetic power fluxes to levels below convective power fluxes. The present results confirm the conclusion of the JET study that in the near SOL kinetic effects are unlikely to cause a power flux 'burnthrough' leading to a substantial increase in the target power load at locations where EDGE2D-EIRENE predicts partial divertor detachment. On the contrary, the impact of a minority of super-thermal electrons on atomic rate coefficients, which were not taken into account, can increase the volumetic power loss in the divertor, leading to stronger detachment.
While the KIPP conductive ion power flux was found not to significantly influence OT power deposition profiles, the electron conductive power flux made a large contribution to the total target power flux at locations where EDGDE2D-EIRENE predicts the highest target power load. But even at these locations the total target power flux, with the replacement of q cond,EDGE2D e by q cond,KIPP e , could only increase the power load by a factor of less than two due to the presence of relatively large electron and ion convective power flux contributions. Further out in the SOL, where target power fluxes are strongly reduced, OT kinetic convective power flux q cond,KIPP e is found to be not much different from q cond,EDGE2D e . Overall, the conclusion of this study confirms that of the earlier JET study; that under ITER-relevant conditions with high radiation from the divertor reducing the target power load down to levels acceptable for the operation of a future tokamak reactor from the engineering point of view, kinetic effects of the parallel propagation and Coulomb collisions of charged particles do not seem to be critically important in calculating the target power load.

Data availability statement
The data that support the findings of this study are available upon request from the authors.
All data that support the findings of this study are included within the article (and any supplementary files).