The prediction of power fluxes and plasma-wall interactions impacted by MHD processes during ITER operation [disruption, Edge Localized Modes (ELMs), 3D magnetic fields applied for ELM control, etc.] requires models that include an accurate description of the MHD processes themselves, as well as of the edge plasma and plasma-wall interaction processes. In this paper, we report progress on improving the edge plasma physics models in the nonlinear extended MHD code JOREK, which has capabilities to simulate the MHD response of the plasma to the applied external 3D fields, disruptions and ELMs. The extended MHD model includes E × B drifts, diamagnetic drifts, and neoclassical flows. These drifts can have large influences, on e.g., divertor asymmetries. Realistic divertor conditions are important for impurity sputtering, transport, and their effect on the plasma. In this work, we implemented kinetic and fluid neutral physics modules, investigated the influence of poloidal flows under divertor conditions in the ITER PFPO-1 (1.8T/5MA) H-mode plasma scenario, and compared the divertor plasma conditions and heat flux to the wall for both the fluid and kinetic neutral model (in JOREK) to the well-established 2D boundary plasma simulation code suite SOLPS-ITER. As an application of the newly developed model, we investigated time-dependent divertor solutions and the transition from attached to partially detached plasmas. We present the formation of a high-field-side high-density-region and how it is driven by poloidal E × B drifts.

## I. INTRODUCTION

Simulation of transient phenomena driven by MHD in tokamak plasmas and their control have so far mostly focused on the MHD processes while the associated plasma-edge and plasma–wall interaction aspects of these processes have been addressed with specialized edge plasma codes in open-loop.^{2,3} This consists in taking the magnetic geometry as given by the MHD codes as fixed in time and using the edge-plasma codes to model the plasma in these conditions. This has two drawbacks: one is that transient modeling of edge-plasma and plasma wall interactions driven by MHD processes cannot be accurately done and the second is that the varying behavior of the edge plasma and plasma–wall interactions does not feedback on the MHD. The work presented here aims to address these two drawbacks by improving the edge plasma physics description in the non-linear extended MHD code JOREK, which has demonstrated capabilities to model the MHD and core plasma behavior during Edge Localized Mode (ELM) suppression by external 3D fields,^{4,5} ELMs^{6,7} and disruptions^{8} to cite few important examples.

JOREK is a fully implicit 3D non-linear extended MHD code for realistic tokamak X-point plasmas.^{1,9,10} The JOREK code has typically been developed and used for studying MHD instabilities such as ELMs and disruptions.^{6,11,12} The present version of JOREK includes some divertor physics in terms of sheath boundary conditions (Mach-one velocity and sheath transmission factors) and a simple neutral fluid model (reflection of outgoing ions as neutrals, ionization, recombination, and radiation). This model has been extensively applied for simulations of ELMs^{12} and ELM control by RMPs^{5} and by pellet triggering.^{13}

The aim of this paper is to describe, benchmark, and apply an improved physics model for the neutrals in the JOREK code. Previously, a kinetic particle framework^{14} has been added to the JOREK fluid MHD models. This particle model has been used to study the transport of Tungsten in the 3D ergodic fields due to ELMs.^{15} Kinetic particles have also been used to simulate the fast particle destabilization of TAE modes.^{16} This particle framework will be used to replace the neutral fluid model with a fully kinetic neutral model (not including molecular physics at this time).

In this paper, we first summarize the essentials of JOREK in Sec. II. Then, in Sec. III, the atomic neutral-plasma interaction physics of interest are derived and presented in Euler equations form (and the form as used for JOREK in Sec. III A) in a three-fluid model; both forms are particle, momentum, and energy conserving. Then, it is reduced with JOREK's assumptions to a two-fluid model. The focus here lies on volumetric ionization and recombination. Charge-exchange (CX) will only be used in the kinetic model. After the derivation, Sec. IV describes the implementation of the kinetic neutral model in JOREK. For the fluid neutral model, we refer to Fil *et al.*^{17} In Sec. V, we further present a benchmark with SOLPS-ITER simulations of a ITER PFPO-1 5MA H-mode scenario^{18,19} and the difference between utilizing fluid or kinetic neutral models within JOREK. We report on the energy conservation properties in Sec. V F. As a first application, the effect of poloidal *E* × *B* drifts on the heat and particle fluxes in this ITER scenario is investigated. This is investigated in the conduction-limited regime for both steady-state solutions in Sec. VI and in dynamic simulations in Sec. VII. In Sec. VIII, we present the dynamic transition toward plasma detachment and the formation of a High-Field-Side High-Density (HFSHD)-region driven by poloidal *E* × *B* drifts. Section IX is the conclusion.

The importance of *E* × *B* drifts in the SOL and divertor has been identified in Refs. 20–22. It is typically very difficult to measure the electric fields and the *E* × *B* drifts velocity directly in experiments. Hence, there is a need of understanding the effect of electric fields by means of modeling and simulations. Since the implementation of *E* × *B* drifts in codes such as SOLPS-ITER^{23} and SOLEDGE2D,^{24,25} the *E* × *B* drifts have attracted increasingly more attention in recent years; also in the MHD plasma modeling community.^{25–32} These edge plasma studies highlight the impact of drifts on divertor asymmetries and SOL flows in multiple machines. Changes to heat fluxes and to detachment behavior by drifts might impact the effective plasma operation space with acceptable power loads on the wall in future fusion devices. The interplay between *E* × *B* drifts and radiating impurities has been identified as a necessary condition for the onset of the experimental “detachment cliff” in DIII-D with a nonlinear dependence.^{33}^{,} *E* × *B* drifts further have been found to play a fundamental role in the development of the high-field-side high-density (HFSHD)-region in ASDEX-Upgrade.^{34,35} The HFSHD-region with radiating impurities can be a precursor for a X-point radiator plasma, which if not controlled can lead to disruptions. However, if controlled well can lead to ELM-suppressed regimes.^{36}

## II. SHORT SUMMARY OF JOREK

This section shortly summarizes the base MHD model in JOREK and the reduced MHD form to provide a background for the derivations of the plasma-neutral interactions derived in Sec. III. We start with the base MHD model, then state the reduced MHD assumptions and how to obtain the weak form of the reduced MHD equations. A more complete explanation of the model can be found in the study by Hoelzl *et al.*^{1}

### A. Full MHD

*ρ*) Eq. (3), and the total pressure (

*p*) Eq. (4). The normalization of the equations is performed such that parameters as the vacuum permeability

*μ*

_{0}and the Boltzmann constant

*k*do not appear in the equations,

_{b}**E**is defined by Ohm's law including resistivity,

*η*, and drift-ordered (diamagnetic) terms,

### B. Reduced MHD

^{1}The ansatz considered here assumes that the time dependent part of the magnetic potential is dominated by the toroidal component. The ansatz for the magnetic field is deduced by approximating $ B \varphi $ as the vacuum $ F 0 / R$ toroidal field. After some calculations, we can obtain an expression of the poloidal component of the velocity ( $ v pol = ( e \varphi \xd7 v ) \xd7 e \varphi $), which captures

*E*×

*B*effects and the ion diamagnetic drift velocity; and the parallel component of the velocity field $ v \u2225 = v \xb7 B$,

The velocity stream function *u* is defined as $ \Phi / F 0$. In this paper, the diamagnetic drift velocity is not used.

*ψ*, velocity stream function

*u*, toroidal current density

*j*, toroidal vorticity

*w*, mass density

*ρ*, temperature $ T = T i + T e$, and parallel velocity $ v \u2225$. There is also the possibility to run with a two-temperature model in JOREK, but is not used is this work. Any other time-evolved parameter (e.g., the neutral density for the fluid neutral model

^{17}) requires another equation with another unknown to solve.

### C. Finite elements and time-integration

The physics variables are Discretized using finite elements based on bicubic Bézier surfaces^{10} in the poloidal plane and by a Fourier series in the toroidal direction. These higher order elements provide $ G 1 \u2212$ continuity, i.e., continuity of the variables and their gradients in the R,Z plane. Even higher order continuity is available.^{37} The iso-parametric finite element grid is flux-aligned after running a Grad–Shafranov solver on the initial equilibrium. Furthermore, the grid can be fine-tuned to add resolution in regions of interest (e.g., the pedestal region and close to the strike-points).

For the time-evolution of the MHD fluid, the fully implicit time-evolution schemes, Crank–Nicholson or Gears^{38} (i.e., a BDF2 scheme), are utilized. These second-order accuracy schemes are unconditionally stable; therefore, the limitations of the time steps are determined by timescales of physical processes present in the simulation.

## III. NEUTRAL MODEL OVERVIEW

We describe two different neutral models used based on the same physics. First the fluid neutral model, and a new approach for the kinetic neutrals that are two-way coupled to the plasma fluid. For now, we are only looking at atomic species as a first step. The plasma-neutral reactions we take into account are electron-impact ionization, effective radiative recombination reactions, and exclusively for the kinetic neutrals, charge exchange (CX). For the neutral fluid model, neutral transport is described as a diffusive process with a spatially constant diffusivity.^{17}

We first derive a more general three-fluid model and thereafter reduce it with approximations to a two-fluid model and write it in the fluid equations form that JOREK uses. This is then the form implemented for the fluid neutral model. For the kinetic neutral model, this gives us a basis for what coupling terms are needed.

*γ*is the ratio of specific heats, $ \varphi ion$ is the ionization potential, and

*X*is an arbitrary initial energy higher than the ionization potential. The reaction equation plus the energy terms are shown in the following equation:

With ionization, the ion gains $ m i / m n$ times the neutral energy, the electron gains the rest of the energy but loses the ionization potential. In an electron fluid description, it does not matter how the energy is distributed between the two electrons. If the recombination results in the neutral being in an excited electronic state $ e \u2212 + i + \u2192 rec n * + photon$, called dielectronic recombination, the excited state will cascade down to the ground (or metastable) state. This will result in most of the ionization potential energy being converted into photons. Generally, we may assume that all the radiatively recombined electrons cascade to the ground state.^{39} For implementation in JOREK, we should take the plasma fluid macroscopic velocity and kinetic energy into account as well. This results in an extra $ 1 2 m v 2$ in the energy terms.

### A. Continuity equation

^{40}The reactions rates are then

*ion*and

*rec*or

*ir*to indicate whether the partial derivative is associated with ionization, recombination or both,

### B. Momentum equation

*i*,

*e*,

*n*denoting the property belonging to species ions, electrons, and neutrals.

### C. Energy equation

#### 1. Ionization

#### 2. Recombination

### D. Assumptions and overview

To use the model in JOREK, we can reduce the equations and remove some terms by implementing assumptions made in JOREK. We reduce the ion and electron continuity equation by regarding only one plasma fluid. Let the center-of-mass density be $ \rho \u2261 m i n i + m e n e$, and *n _{e}* =

*n*in the absence of impurities. We then assume the limit of $ m e \u2192 0$. Then, $ \rho = m i n i$ and $ m i / m n = 1$. This further results in the electron momentum equations being zero and all terms with electron mass disappearing.

_{i}*T*, and nothing for the energy to go to. This is the energy lost due to radiation. In the ADAS database the Continuum and line power driven by recombination is combined with Bremsstrahlung of dominant ions (PRB). We can remove $ \Gamma rec T e$ from electron energy equation if and only if we explicitly take radiation into account. We have to be careful here, as the PRB recombination is the sum of radiation from free-electron recombination, the dielectronic-cascade + Bremsstrahlung power coefficients. This includes the dielectronic-cascade power loss from an excited energy state to the ground state. This energy is not lost from the electron fluid energy equation. When we explicitly write the electron energy loss in terms of radiation power density, we should write $ P rad = ( P PRB \u2212 \Gamma n rec \varphi ion )$ to correct for the ionization potential energy. In terms of radiated power coefficients,

_{e}## IV. KINETIC NEUTRAL MODEL

The kinetic neutral model is part of the kinetic particle framework developed in JOREK.^{14} The particle framework is very flexible for the implementation of many other physical models. In this section, we first present the general way of coupling particles to the fluid description of the plasma, then focus on fully kinetic neutrals.

### A. Pusher and projections

^{43}in toroidal geometry in the time varying electric and magnetic fields of the MHD fluid. Particles are followed in both ( $ R , Z , \varphi $) space and in the local ( $ s , t , \varphi $) space of the cubic Bezier finite elements. The ionization of the neutrals gives rise to a source of ions, momentum, and energy in the fluid. To add these sources to the fluid equations, the discrete particle distribution needs to be converted to the continuous, in space, description of the finite element basis. For example, the ion source $ S \rho ( x ) = \u222b m f ( x , v ) d v$, where

*f*, the distribution of ionized neutrals is represented in the finite element $ H i j ( s , t )$ and Fourier basis: $ S \u0303 \rho ( x ) = \u2211 ijk p \rho , ijk H i j ( s , t ) H \varphi , k ( \varphi )$, with

*p*being the expansion coefficients. In the weak form, this yields a set of equations,

_{ijk}*w*. At every ionization event this weight is reduced until it reaches a critical lower value at which it is removed. For the conservation properties of the combined fluid–particle system, it has been found essential to collect the contributions to the right-hand side particle sums at every particle time step. To reduce the noise in the finite element representation (or as regularization), the following equation is solved instead of Eq. (22),

_{n}*λ*and

*ζ*are the regularization coefficients and $ v *$ the test function of the finite element method (weak form). The regularization coefficients can be chosen with different values for the

*n*= 0 harmonic (as opposed to higher harmonics) and in the perpendicular and parallel directions. This method of projection allows us to project an arbitrary number of any moments of the particles while only calculating one matrix factorization.

### B. Fluid and kinetic boundary interaction

^{14}based on the Eckstein sputtering coefficient.

^{44}Both the fluid and kinetic particles can reflect and sputter on the wall. For recycling (i.e., plasma hitting the boundary and being re-introduced in the domain as a neutral), the TRIM-code database

^{45}is used to determine the probability of fast reflection or thermal desorption based on the electron temperature and wall material. For the plasma fluid, the boundary values are known and we can integrate the total reflected or sputtered flux on the boundary,

As the incident plasma fluence is only defined based on the MHD fluid time stepping scheme, sputtering and reflection events are performed only once per fluid time step. This means, after a fluid time step, superparticles are introduced in a sputtering event, then in the kinetic particle loop, they are evolved per kinetic time step just like the other superparticles that already existed. Although for the kinetic particles this a pulsed particle source, as it is coupled back to the plasma fluid only every fluid time step, it does not result in oscillating or pulsed solutions.

### C. Moments and coupling scheme

*S*. This results in the following source terms in the JOREK form,

_{E}Although the moments of kinetic particles are coupled as explicit terms, they are multiplied by implicit MHD terms [e.g., **v** in Eq. (15)] to maintain the highest order accuracy. In these coupling schemes, when kinetic particles are the dominant species, the accuracy of the time-evolution scheme might be of a lower order. Since we have knowledge on the dependency of rate coefficients on the plasma density and temperature (*n _{e}*,

*T*), it is possible to improve the order of accuracy by making the coupling terms (semi-) implicit in future work.

_{e}### D. Plasma–neutral interactions

As the kinetic particle framework is a Monte Carlo model, we need to sample from probabilities to determine which interaction(s) will happen. Ionization, CX, and line radiation are dependent on the kinetic neutrals population and it should be calculated based on the kinetic time step. Phenomena that create neutrals from the fluid plasma are only well-defined over an entire fluid time step ( $ \Delta t fluid$) and are not an interaction that can be based on the existing superparticles in the plasma domain. We, thus, use separate routines to initialize the superparticles once per fluid time step for plasma–boundary interaction, puffing, and recombination. However, in principle, it is possible to initialize these phenomena more than once per *t _{fluid}*. Section IV D 1 describes the implementation of the plasma–neutral interactions.

#### 1. Ionization

#### 2. Charge-exchange

*T*and macroscopic drift velocity $ v drift$ and exchange its momentum with the superparticle,

_{e}#### 3. Line radiation

#### 4. Recombination

The location of the recombined superparticle is sampled in the poloidal plane from the local element space (s,t) with two uniform distributed random numbers $ U rng [ 0 , 1 ]$, and in the toroidal direction, $ \varphi = \varphi plane + 1 2 U rng [ 0 , 1 ] \Delta \varphi $ to place it in the sector associated with the poloidal plane with coordinate $ \varphi plane$ and distance between planes $ \Delta \varphi $.

## V. BENCHMARK WITH SOLPS-ITER

To test the implementation of the neutrals code development, JOREK simulations are compared to SOLPS-ITER simulations of an ITER (pre-fusion power operation) PFPO-1 H-mode scenario without impurities (first published in Park *et al.*^{18}) In this section, we compare the plasma power and peak heat flux to the divertor walls as a function of the upstream density (i.e., the outer mid-plane separatrix density). Here, the simulations are performed for a range of different neutral fueling rates in the divertor. The details of the scenario are $ B T = 1.85 T , \u2009 I p = 5 M A$, *P _{heating}* = 20 MW, $ D \u22a5 = 0.3$ m

^{2}s

^{−1}, $ \chi i , \u22a5 = \chi e , \u22a5 = 1.0$ m

^{2}s

^{−1}. Parallel heat transport is based on the Braginskii closure.

^{46,47}The fluid neutral diffusivity was chosen to match the SOLPS-ITER neutral density decay from the strike point following the field lines at low puffing rates. This resulted in $ D n \u2248 2 \xd7 10 2$ m

^{2}s

^{−1}. For the kinetic neutrals, such an assumption is not needed. To reach a quasi-steady-state, JOREK simulation were run between $ 11 \u2009 and \u2009 55$ ms depending on the neutral gas puffing rate.

Although SOLPS-ITER solves plasma transport differently and has hydrogen molecules included, we still expect to see similar behaviors with JOREK. It is known that with temperatures under 1 eV, molecules are important for the correct detachment behavior and the exact behavior of detachment front.^{18} However, to benchmark as a first step let us look at the power incident on the inner and outer divertor targets and the peak heat fluxes.

Although more advanced options are available in JOREK, to implement and study neutral physics we do not use everything at once. Here, the visco-resistive reduced MHD is used as model where the time-evolved parameters are $ \psi , j , u , w , \rho , T = T i + T e , v \u2225$, and *ρ _{n}*, the neutral density for the addition of fluid neutrals.

### A. Simulation grids and model reduction

Figure 1 shows the JOREK FE (finite element) grid extending up to the first wall/divertor and the SOLPS-ITER plasma and kinetic grids. SOLPS-ITER does not simulate the full core, and the plasma does not reach the wall except for the divertor. The red lines are possible puffing and pumping surfaces. The JOREK grid does include the core. SOLPS-ITER uses two separate grids, one for the plasma fluid, one for the kinetic particles. The JOREK kinetic particles exist in real space (R, Z, $\varphi $) and local finite element coordinates.

As the fueling rate scan in the study by Park *et al.*^{18} was performed with SOLPS-ITER without drifts, we have reduced the JOREK models here to obtain driftless behavior where also the magnetic field is frozen. For this, the time-evolution of the parameters $ \psi , j , u$, and *w* is removed. This fixes the magnetic field and removes the *E* × *B* flows.

For this benchmark, the neutral gas puffing location was chosen to be in the divertor. In steady-state and without impurities, the puff location has no significant impact in the near-SOL and divertor region.^{18} However, this could have an effect on time-dependent behavior and the effectiveness of the fueling the plasma.

The pumping surface in the JOREK grid is in the private region boundary. Here, diffusively transported neutrals are pumped out.

Figure 2 shows the resolution of the FE grid on top of the plasma density. To resolve the sharp (density) gradients close to the strike points, a resolution on the order of 1 mm is need in the direction along the field. Insufficient resolution will alter the divertor solution of the simulation. The figure shows the gradients to be well resolved. Monte Carlo noise of the projections of the kinetic particles depends on the number of superparticles per finite element. A finer grid, thus, requires more superparticles for similar noise levels.

### B. Benchmark with fluid neutrals

In this subsection, we present the results of the fueling scan done in JOREK with fluid neutrals compared to the SOLPS-ITER results. Only the results without drifts are directly comparable between JOREK and SOLPS-ITER. Here, only a few metrics are compared as a global overview.

First, the heat load to the wall as function of the neutral gas puffing rate was investigated in Ref. 48. This showed the dependence of the total plasma heat load on the fueling rate are qualitatively in good agreement.

Figure 3 shows the plasma heat load to the divertor as function of the outer mid-plane separatrix density *n _{upstream}*. Although JOREK obtains the same power ranges; in a stationary state,

*n*at a given fueling rate is lower in JOREK simulations compared to SOLPS-ITER. This initial drop in

_{upstream}*n*is related to the difficulty of keeping the shape of the initial density profile, which is assumed

_{upstream}*ad hoc*rather than consistent with gas fueling.

The use of diffusive fluid neutrals with a spatially constant *D _{n}* is not sufficient to provide the correct plasma mechanics. However the addition of fluid neutrals allows JOREK simulation to transition from a typically convective (sheath-limited) regime to the conduction limited regime. This relatively simple fluid neutral model captures enhanced energy losses with increased neutral gas puffing, which leads to decreased divertor heat loads. Given the simplicity of the model, we can consider this as reasonable agreement. However, further study for this model is needed to identify the reason for the discrepancies with the SOLPS-ITER heat loads at the low end of upstream densities ( $ n upstream \u2a85 1.5 \xd7 10 19 m \u2212 3$).

### C. Benchmark with kinetic results

The fueling scan using the JOREK model with kinetic neutrals compared to the SOLPS-ITER results is shown in Fig. 4. Only the results without drifts are directly comparable between JOREK and SOLPS-ITER. As opposed to fluid neutrals, kinetic neutrals do not need to assume a diffusivity for the neutral transport. In the kinetic neutral model, the neutrals do have a velocity and can exchange momentum and energy with the plasma via a combination of CX, ionization and recombination.

Figure 4 presents the same scenarios as in Sec. V B. At the lower- to mid-range of the upstream density ( $ n upstream \u2a85 2 \xd7 10 19 m \u2212 3$), there is good agreement in total power of JOREK with SOLPS-ITER. The inner–outer divertor power balance is slightly more symmetric than SOLPS-ITER. Beyond $ n upstream \u2248 2 \xd7 10 19 m \u2212 3$, the SOLPS-ITER plasma transitions from high-recycling to a fully detached divertor plasma. While this happens *n _{upstream}* starts to saturate. In the no-drift JOREK simulation, the transition to detachment occurs at about 10% higher

*n*than SOLPS-ITER. The transition as function of

_{upstream}*n*is less sharp in JOREK. The slope ( $ \u2202 Power / \u2202 n upstream$) is decreasing for higher

_{upstream}*n*This does indicate a similar but less strong form of density saturation from $ n upstream \u2a86 2.4 \xd7 10 19 m \u2212 3$.

_{upstream}Since our kinetic neutral model now only consists of atoms, and neutral–neutral collisions are neglected, one can expect the comparison with SOLPS-ITER to diverge when detaching the plasma from the divertor. We suspect the absence of molecules and neutral–neutral collision together with the exact geometry of the simulation, which is not reproduced in JOREK, can explain the difference in our results from the ones of SOLPS-ITER.

Neutral–neutral collisions combined with the geometry of our simulation grid act to redistribute neutrals in the divertor region. The lack of neutral–neutral collisions allow hot neutrals to be ballistic and reach the wall with a lot of energy. Neutral–neutral collisions would disperse this energy through a lot of neutrals before reaching the wall. The pressure created in the neutrals then drive transport under the divertor dome. This redistribution of neutral particles can effectively change the power distributed between the inner and outer divertor.^{49} In our simulations, there is no grid under the divertor dome.

For plasma detachment, molecules are key around a few eV. They play an important role in getting the detachment behavior right.^{18,50} The molecules create strong (momentum) losses for the plasmas at these low temperature; this creates a stronger ion flux rollover and decreases heat load at the wall. Stronger momentum losses along the field lines in the divertor will alter the transport behavior. This might result in stronger upstream density saturation as fueling from the divertor becomes less efficient.^{18,50} The implementation of molecules in JOREK is left for future work.

A key divertor metric that is important to determine the compatibility of the plasma power fluxes with the plasma facing components at the divertor is the peak heat flux ( $ q p k , \u22a5$). Figure 5 shows $ q p k , \u22a5$ on the inner and outer divertor target. Note that for both SOLPS-ITER and JOREK a toroidally averaged wall is considered (i.e., the tilting and beveling of the actual divertor tiles are not taking into account). To include the shape of the monoblocks, a corrected angle can be taken into account.^{51} Generally, we observe JOREK having a higher $ q p k , \u22a5$. For the same plasma heat load, this indicates that JOREK profiles are narrower than SOLPS-ITER. The slope ( $ \u2202 q p k , \u22a5 / \u2202 n upstream$) is similar between the two codes. However, as with the total power, the reduction of the peak heat flux with an increased fueling rate in JOREK occurs at higher upstream densities compared to SOLPS-ITER. In both simulations at high *n _{upstream}*, the inner and outer target peak heat fluxes become the same.

Generally, the inner target $ q p k , \u22a5$ is lower than the outer target. At the low end of *n _{upstream}*, JOREK has a higher inner $ q p k , \u22a5$. Why there is a strong difference between the behavior between JOREK and SOLPS-ITER in the inner target peak heat flux at low upstream density is not entirely clear. For JOREK, the ion flux is higher on the inner target, and the density is not yet high enough to really decrease the temperature. So the high ion flux with higher temperature combines in a higher heat flux.

Differences in the peak heat fluxes could have several causes. For the same power, a high peak heat flux indicates JOREK has a lower effective cross field heat diffisivity $ \kappa \u22a5 , effective$. There are two notable differences in the simulation that could account for this. For this benchmark, JOREK employs the one-temperature model, while SOLPS-ITER used a two-temperature model, which is important since the same gradients along the field line are implicitly assumed for ion and electrons, which is not the case when two temperatures are considered. Furthermore, molecule–plasma interaction can provide efficient momentum-exchange channels that effectively enhances the cross field particle and energy transport.

### D. Performance of the neutral and kinetic model

We have discussed the physics results of the used fluid neutral and kinetic neutral model. The kinetic neutral model obtains quantitatively better results. The fluid neutral model with a spatially constant diffusion coefficient does not suffice to obtain the desired divertor solutions. Although the fluid model is much simpler to implement, i.e., add one extra equation and variable to the system of equations, it is typically not computational more advantageous. On typically four nodes, with 36 cpus per node, the kinetic neutral model using $ \u223c 10 6$ simulation particles using a time step of $ 2 \xd7 10 \u2212 8$ s, is as fast using the fluid neutral model. This means that in axisymmetric simulations, there is no computational penalty for using the kinetic neutrals over the fluid neutrals. For bigger computational problems (e.g., 3D), the kinetic neutral model scales almost ideally with the number of compute nodes and scales more favorably than the fluid neutral model. To reach steady-state, it takes $ 9000 \u2013 10 \u2009 000$ cpuh, which equates to $ 60 \u2013 70$ h wall-time on said four nodes.

### E. Initial and boundary conditions

The initial conditions of some of the variables are defined by (re-)constructed magnetic equilibrium. This then gives us an initial profiles for *ψ*, *ρ*, and *j*. From experimental profiles (or from other codes in our case), the pressure p is separated into the mass density *ρ* and temperature *T*. If experimentally the parallel velocity is known in the core and edge, we can use it as initial condition for $ v \u2225$. These initial conditions are all functions of the flux, asymmetries will arise naturally. The potential and vorticity u and w are initialized to zero but arise naturally within a few time steps.

For the boundary conditions, in this case, we use a fixed magnetic boundary, *ψ* thus has a Dirichlet boundary condition arising from the (re-)constructed equilibrium. The potential and vorticity u and w are Dirichlet conditions set to the value of 0. *ρ* has a Neuman boundary condition. The boundary conditions for the $ v \u2225$ and the parallel heat flux (i.e., a function T) are described by the usual magnetized sheath boundary conditions with a sheath heat transmission factor.^{52}

### F. Conservation properties with kinetic neutrals

The kinetic particles Boris pusher used for advancing the particles in time preserves the particle energy conservation down to machine precision.^{43} However, for kinetic-fluid interaction not all properties are perfectly numerically conserved. The conservation further depends on the chosen time-stepping scheme. The coupling of moments of the kinetic particle distribution result in explicit source terms in the MHD equations. The time stepping of the JOREK MHD fluid utilizes the fully implicit Gears scheme^{38} (i.e., a BDF2 scheme). The Gears scheme does not conserve energy exactly but is second order accurate in the time step.

The recombination process is calculated from the finite elements of the plasma fluid. This results in an average relative error in momentum $ \u223c 10 \u2212 4$. However, the relative error in energy is $ \u223c 10 \u2212 3$. The ionization, CX, and radiation errors will depend on the error in the projection to the finite elements. The error should scale with the amount of superparticles per finite element. As the neutral–plasma interaction is localized on the centimeter scale, it is non-trivial to assess this *a priori*. However, most neutrals “live” close to the strike points where plasma–neutral interaction is high. To ensure adequate statistical coverage of the plasma–neutral or plasma-wall processes, simulations were performed with $ 10 6 \u2013 10 7$ superparticles with an average weight of $ 5 \xd7 10 13 \u2013 5 \xd7 10 14$ atoms. The weight of each superparticle is self-governed by the physical processes.

The volumetric losses of the plasma and neutrals are the plasma–neutral (pn) interaction and corresponding radiation losses. The surface losses at the boundary are the plasma–wall interaction and neutral–wall interaction. The effective plasma heat load on the wall is the outgoing plasma heat flux minus the energy retained after wall-assisted recombination. The effective neutral heat load is the outgoing kinetic neutral heat flux minus the energy the neutrals retain after colliding with the wall. $ P heating \u2212 d U d t$ is how much energy is supplied to the plasma minus the change in internal energy. In steady-state, the power balance is conserved with a relative error $ \epsilon \u2248 10 \u2212 3 \u2013 10 \u2212 2$. The error depends on the time step size and the amount of simulation particles. For these applications this is a sufficient accuracy.

## VI. THE EFFECTS OF POLOIDAL DRIFTS

### A. Quasi-steady-state plasma heat load and peak heat flux with poloidal drifts

In the benchmark with SOLPS-ITER, V, we focused solely on JOREK results without drifts. By design, the JOREK MHD models include *E* × *B* flows as an essential ingredient. In this section, we compare standard JOREK reduced MHD model (including the evolution of magnetic flux, and vorticity) with JOREK simulations without drifts.

Figure 7 presents the results of the same scenarios as in Secs. V B and V C. At the lower end of *n _{upstream}*, the drifts do not influence the solutions in the divertor. With increasing

*n*, the

_{upstream}*E*×

*B*drifts start redistributing the power balance such that more power goes to the outer divertor target. Simultaneously the inner target heat load decreases faster than in the simulations without drifts. Drifts cause a generally stronger decrease in power with increasing

*n*. Specifically, the

_{upstream}*E*×

*B*flow from the outer target to the inner target through the private region causes the power to drop to the inner target. Around $ n upstream \u2248 2 \xd7 10 19 m \u2212 3$, the plasma is already detached from the inner target. While detachment behavior (roughly $ T e , target < 5 \u2009 eV$) without drifts starts at much higher

*n*.

_{upstream}Around $ n upstream \u2248 2 \u2013 2.2 \xd7 10 19 m \u2212 3$, the upstream density rises with increased fueling while the power to the divertor almost remains constant. This is a phase-like transition of the plasma and coincides with the formation of a so-called high-field-side high-density region (HFSHD-region). The dynamic formation of the HFSHD-region is presented and explained in Sec. VIII.

Figure 8 presents the steady-state peak heat fluxes (*q _{pk}*) as function of

*n*. Although at low

_{upstream}*n*the total power is the same with and without drifts, the heat flux profiles have changed. At the low

_{upstream}*n*, the inner target peak heat flux $ q p k , I T$ is increased while $ q p k , O T$ is decreased. $ q p k , I T$ decreases more quickly due to the

_{upstream}*E*×

*B*drifts. Counter-intuitively with increasing

*n*$ q p k , O T$ increases up until $ n upstream \u2248 1.5 \xd7 10 19 m \u2212 3$. At higher

_{upstream}*n*, $ q p k , O T$ decreases similar to its no-drift counterpart until the phase-like transition of the plasma.

_{upstream}### B. (Quasi-)steady-state ion flux comparison

The ion flux through the boundary in JOREK is $ \Gamma \u221d \rho v \u2225 \xb7 n \u0302$, where the boundary parallel velocity is set to the sound speed because of the Bohm sheath criterion, $ v \u2225 = c s$. This then gives us $ \Gamma \u221d \rho T 1 / 2$. Figure 9 shows the (quasi-) steady-state ion flux profiles over the inner and outer divertor target, with and without drifts, for a range of outer mid-plane plasma densities, *n _{upstream}*. The upper row shows JOREK results (with drifts) and the lower row in the figure shows ion flux results of JOREK without drifts. The results from Fig. 9 correspond to the points in Figs. 7 and 8. As the plasma responds differently to the same amount of gas puffing when the drifts are removed, the comparison cannot easily be done for exactly the same

*n*.

_{upstream}The total plasma heat load to the wall (as shown in Fig. 7) shows that at the low end of *n _{upstream}*, the drifts do not have a big influence. Yet/However with increasing

*n*, the power redistributes differently between the IT and OT. At higher $ n upstream \u2248 1 \xd7 10 19 m \u2212 3$, the location of and height of the peak ion flux is different. With increasing

_{upstream}*n*, we see the qualitative behavior more diverging, i.e., the location and height of the peak, single- vs double-peaked, the rollover and detachment.

_{upstream}In all results (with and without drifts, and at the IT and OT), the ion flux moves upwards away from the strike point (to higher $ s \u2212 s sep$) with increasing *n _{upstream}* as partial detachment sets in. The

*E*×

*B*drift compresses the ion flux on the OT, where the profile remains closer to the separatrix than its “no-drift” counterpart. On the inner target, the

*E*×

*B*drift pushes the profiles upwards, further away from the separatrix.

The no-drift simulations show a clear double-peaked profile. Observing the plasma profiles $ \u2248 1$ cm from the wall, the left peak is due to a combination of peaked *n _{e}* and low $ v \u2225$, and the right peak is at lower

*n*, but $ v \u2225$ sharply increases. This accelerates the plasma differently to the expectations from the simple sheath and we obtain a double-peak ion flux profile. In the steady-state profile with drift, we do not see a strong double-peaked profiles, only more flattened profiles, which is due to the same phenomena, but less strong. Furthermore, the “drift” IT receives significantly lower ion fluxes than the no drift result. The asymmetry between drift OT and IT is due to the poloidal

_{e}*E*×

*B*transport from the OT to the IT and from the IT upwards. The upward

*E*×

*B*flows the IT upwards are an important phenomenon while detaching the plasma, which will be discussed in more depth in Sec. VIII.

The transition from the high-recycling regime toward detachment is typically marked with the “rollover” of the ion flux (i.e., decreasing ion flux with increasing *n _{upstream}*). The “no drift” rollover is fairly symmetrical between the OT and IT. However, as the

*E*×

*B*drift drives asymmetry, it creates a more asymmetrical rollover. On the drift IT rollover begins around $ n upstream \u2248 1.8 \xd7 10 19 m \u2212 3$ while the OT start seeing rollover around $ n upstream \u2248 2.15 \xd7 10 19 m \u2212 3$.

Remarkably, at around $ n upstream \u2248 2.3 \xd7 10 19 m \u2212 3$ in the drifts simulation, there are two different profiles at almost equal *n _{upstream}* ( $ 2.28 \xd7 10 19$ and $ 2.3 \xd7 10 19 m \u2212 3$). This coincides with the formation and stabilization of the HFSHD-region. This non-linear transition has a strong effect on the power distributed to the IT and OT and on the fueling of the plasma from the divertor. This is discussed in more depth in Sec. VIII. During this transition, the ion flux profile first moves downwards (lower $ s \u2212 s sep$) then moves $ \u2248 1$ cm upwards again. This is due to the strong E-field created by the HFSHD-region. This transition only occurs with drifts. Simulation without drifts Fig. 9 show a much smoother behavior on the IT and OT heat fluxes with increasing density.

### C. (Quasi-)steady-state heat flux comparison

The plasma heat flux through the boundary in JOREK, $ q \u221d \gamma sheath \rho T v \u2225$. With $ v \u2225 = c s$, this gives $ q \u221d \gamma sheath \rho T 3 / 2$. Figure 10 shows the (quasi-) steady-state heat flux profiles over the inner and outer divertor targets, with and without drifts, for a range of outer mid-plane plasma densities, *n _{upstream}*. Again, this is only the plasma heat flux. The upper row shows JOREK results (with drifts) and the lower row shows heat flux results of JOREK without

*E*×

*B*drifts. The results from Fig. 10 correspond to the points in Figs. 7–9.

Figures 7 and 8 shows that although the total heat load to the wall is not strongly influenced by drifts at the low end of *n _{upstream}*, the peak heat fluxes are different. As the heat flux $ q \u221d \rho T 3 / 2$ and the ion flux $ \Gamma \u221d \rho T 1 / 2$, it is expected to see stronger changes with heat flux profiles with increasing

*n*.

_{upstream}The no drifts results in Fig. 10 show a similar response on the IT and OT. $ q n o \u2212 drift , I T$ decreases faster with increasing *n _{upstream}* than $ q n o \u2212 drift , O T$; however, qualitatively they behave similarly, and

*q*moves progressively upwards (higher $ s \u2212 s sep$). As with with the ion flux in Fig. 9, the

_{peak}*E*×

*B*drifts compresses and pushes the OT heat flux profile downwards (lower $ s \u2212 s sep$) and pushes the IT heat flux upwards. Due to the $ T 3 / 2$ dependency and higher temperature in the far-SOL, the heat flux peak is always at higher $ s \u2212 s sep$ than the ion flux.

At $ n upstream \u2248 2 \xd7 10 19 m \u2212 3 , \u2009 q I T , n o \u2212 drift$ is still significant, however, $ q I T , drift$ has already completely collapsed and is in a detached state (taking detachment as $ T e , target \u2264 5$ eV). Due to the *E* × *B* drifts, the heat flux profiles on the IT and OT behave strongly asymmetrically. The earlier $ q I T , drift$ decay seems to be a result of two factors, both caused by *E* × *B* drifts. *E* × *B* drifts cause a different distribution of the transport from plasma core through the edge to the wall. Second, *E* × *B* drifts cause transport from the OT to IT via the private flux region. The higher density at the IT lowers the temperature leading to the low heat fluxes on the IT. Double-peaked structures are seen, similarly to those in the ion flux, due to a combination of regions with high *ρ* and low *T* and low *ρ* and high *T*.

The no drift results show a steady progressive decay of the heat flux. There is no sign of phase or regime change into detachment. With drifts, the results are asymmetric and around $ n upstream \u2248 2.3 \xd7 10 19 m \u2212 3$, the heat flux profile moves first downwards (to lower $ s \u2212 s sep$) and upwards. The same as for the ion flux results (Fig. 9).

### D. 2D poloidal flows

So far, 0D and 1D profiles have been presented. To understand how these changes occur due to *E* × *B* flows, let us look at the 2D poloidal flow profiles. 2D flows and divertor geometries have been extensively studied in experiments and simulations to design next generation divertors (such as the ITER divertor).^{53} Even without drifts, neutral recycling associated with the vertical divertor geometry has a deep impact on plasma flows. Herein strong recycling leads to flow reversal away from the target in the separatrix region while it causes strong flows in the outer SOL. Here, we will focus more on the effect of drifts on top of the already existing flows due to the divertor geometry. Figure 11 (left) shows the electron density profile *n _{e}*, with on top of that, the 2D poloidal flows visualized with colored arrows. The color represents the relative importance of the

*E*×

*B*flows compared to the rest of the poloidal flows (the parallel component in the poloidal plane), i.e., $ E \xd7 B factor \u2261 | v E \xd7 B | | v pol \u2212 v E \xd7 B |$. The length of the arrow represents magnitude of $ v pol$.

The quasi-steady-state simulation shown in Fig. 11 is in a high recycling regime ( $ n upstream \u2248 1.95 \xd7 10 19 m \u2212 3$). It is still attached, but clearly there is already some density transport upwards helped by *E* × *B* flows. Figure 11 (right) shows a contour plot of the electric potential; and the *E* × *B* flow direction is along these contours. It is assumed that the wall structures are grounded and we can treat the electric potential with a Dirichlet boundary condition. The plasma tends to form three *E* × *B* vortex regions (i.e., with a closed contour in a small region), all relatively close to the X-point. One on the LFS (low-field-side), one on the HFS (high-field-side), and one in the private region just below the X-point toward the HFS. When an *E* × *B* vortex becomes dominant over non-*E* × *B* flows, it can form a high-density region, in the attached plasma case, these vortices are not dominant.

The *E* × *B* flows are dominant around the X-point, both up- and downstream. The *E* × *B* drift causes the plasma to favor distribution of the upstream plasma toward the OT. The plasma transport from the OT to the IT via the private region is strongly enhanced by *E* × *B* flows. At the IT strike-point, the plasma is partially transported upwards. The low temperature allows more neutrals to travel further upwards before being ionized; allowing for more density above the separatrix strike point, leading to the broadening of the ion flux. The *E* × *B* flows then transport the plasma further along the target upwards. With increased puffing the plasma transports more and more to the HFS *E* × *B* vortex, leading to a transition in the plasma (see Sec. VIII).

From Fig. 11 we can also observe how the core plasma is being fueled from the divertor plasma. In Fig. 11, along the separatrix the plasma is flowing downwards. In the SOL, both on the IT and OT, the plasma moves upwards. However, the flow is not exactly aligned with the flux surfaces; they are slightly pointing inward. The flow inversion along the separatrix occurs around the OMP for this scenario. The flow inversion point varies with different amounts of gas puffing. It typically moves upwards with more puffing.

## VII. DYNAMIC PLASMA RESPONSE WITH KINETIC NEUTRALS

In Secs. V and VI, steady state divertor solutions have been presented. The JOREK steady state solutions are in fact the result of time dependent simulations, with an accurate description of the dynamic evolution of the system. Depending on the timescale of the change in parameters (such as the gas puffing rate), the time dependent solution will not be the same as a series of steady state solutions. In the following, we will discuss some aspects of the time dependent behavior, including the dynamics of detachment and of the formation of the HFSHD region.

### A. Dynamic downstream-upstream plasma response

JOREK simulations (with kinetic neutrals) are typically started from a initial scenario without neutrals and the divertor legs being fully developed with a attached divertor. As we proceed with puffing neutrals, a dynamic plasma response is obtained from an attached to conduction-limited, to high-recycling and eventually detachment, depending on the puffing rate.

As an example, Fig. 12 shows the time evolution of the plasma heat load to the wall, plotted as function on $ n upstream ( t )$ at the OMP. As soon as the neutral gas is introduced (at *t* = 0 ms) in the plasma, we observe a quick decay of $ \u2248 8.3$ MW of power in the first 5 ms. Simultaneously, the upstream density increases from $ \u223c 1$ to $ \u223c 1.5 \xd7 10 19$ m^{−3}. In the following 36.6 ms, the power further decreases by only another 1.3 MW, whereas the upstream density slowly increases, reaching saturation after 30–40 ms. Thus, for this ITER scenario, at least 40 ms (in real time) has to be simulated to obtain a steady-state solution. The figure also includes the power to the divertor from the steady state solutions (red dots). It clearly illustrates the difference between the dynamic and steady state solutions: in the dynamic solution, the power to the divertor decreases about twice as fast (as a function of upstream density) compared to the steady state solution.

The initially rather quick response of the heat load is expected as two simultaneous effects work on decreasing the heat load. First, the plasma-neutral interaction leads to energy losses in the divertor by ionization and radiation. Second, due to the fueling, the plasma becomes more high-recycling, leading to an strong increase in the plasma density at the target. The plasma “tries” to equilibrate the pressure over a magnetic field line, and thus, as the density increases, the plasma temperature decreases.

The time and manner in which this initial power decay occurs is very similar for simulations with different rates of puffing neutrals in the divertor. To which degree the power decreases strongly depends on this rate. From our simulations, typically $ \u223c 50 %$ of the change in *n _{upstream}* occurs in the 5 ms after the neutrals reach the plasma. This behavior might change depending on the waveform used for puffing neutrals and is associated with the propagation of the increased density source upstream toward the mid-plane separatrix whose timescale is dictated by the ion transport timescale along the field line from the divertor to the mid-plane.

If we remove the evolution of the electric potential and vorticity, (i.e., no *E* × *B* drifts), this general pattern remains, although the path $ Power ( n upstream )$ is less concave. Thus, the power decreases slower with an increase in *n _{upstream}* than with

*E*×

*B*drifts included.

### B. Heat flux

As the dynamic progression of the power as a function of *n _{upstream}* and time does not follow the same path as the steady-state path, one might expect the heat and ion fluxes to differ as well. Figure 13 shows the time evolution of the plasma heat flux profile within one simulation along the inner and outer divertor target, with (upper row) and without

*E*×

*B*drifts (lower row), for a range of outer mid-plane plasma densities.

The results with drifts from Fig. 13 are from the same simulation as Fig. 12. Without drifts, the dynamic peak heat flux decreases faster as function of *n _{upstream}* than seen in the steady-state results (i.e., for the same

*n*, the dynamic peak heat flux is lower). The dynamic heat flux profiles also do follow the same shapes as seen in the steady-state results. The height, width, and location of the heat flux peaks are slightly altered.

_{upstream}With drifts, the changes are slightly bigger. The same general behavior where the $ q p k ( t , n upstream )$ decreases faster as function of *n _{upstream}*. In addition the shape of the profile—e.g., singly or doubly peaked profile—is qualitatively different. In the dynamic profiles, much stronger secondary peaks are observed.

In the simulations with and without drifts, the heat flux is decreasing, but without drifts, we do not really see the plasma heat flux or temperature “detaching” from the IT. With drifts, it quickly almost reduces to zero. One could call this a partially detached state, because although the plasma temperature on the IT is below 5 eV, the ion flux just started to decrease.

### C. Ion flux

Figure 14 shows the ion flux profile evolution over the inner and outer divertor target, with and without drifts, for a range of outer mid-plane plasma densities, *n _{upstream}*. Again, this is only the plasma ion flux. The upper row shows JOREK results (with drifts) and the lower row in Fig. 13 show heat flux results of JOREK without

*E*×

*B*drifts.

The results with drifts from Fig. 14 are from the same simulation as Figs. 12 and 13. At low *n _{upstream}*, with and without drifts, the profiles are similar for both the IT (inner target) and OT (outer target). When increasing the neutral fueling; and thus, also

*n*, a strong increase occurs into the ion flux as the plasma starts developing to a high-recycling regime.

_{upstream}In the range below $ n upstream \u2248 2 \xd7 10 19$, for equal *n _{upstream}*, the steady-state peak ion flux is lower than the dynamic peak ion flux. The dynamic ion flux reaches earlier higher recycling stage than would be expected from steady-state solutions due to the dynamic SOL transport, with parallel transport being limited by the parallel sound speed, discussed above. Furthermore, the dynamic ion flux profiles are qualitatively similar to the steady-state solutions.

## VIII. PLASMA DETACHMENT AND THE FORMATION OF THE HFSHD-REGION

With enough gas puffing and a high enough upstream density, the losses in the divertor region become large enough such that the plasma will detach from the wall. In literature, many definitions have been used for describing and declaring the state of plasma detachment (e.g., Ref. 54) In this work, we shall not take a very strict definition of the detachment regime: the plasma is considered detached when the target electron temperature $ T e , target \u2264 5$ eV. In this section, we explore mostly more qualitative effects of detaching the plasma.

In Sec. VI, we have already presented some results of plasmas in detached states (Figs. 7–10). Including drifts, Fig. 7 shows a transition in the power as function of *n _{upstream}*. Figure 9 shows the shift of the ion flux around

*n*indicating a movement of the divertor plasma density in the divertor and at higher

_{upstream}*n*a decay of the ion flux. Figure 10 shows the almost complete decay of the heat flux with high

_{upstream}*n*. In all these figures, we observe a strong dependence on

_{upstream}*E*×

*B*drifts, which eases the detachment on the IT. This section will focus in depth on the JOREK simulations of detached plasmas; and how

*E*×

*B*drifts cause the high-field-side high-density region and a phase-like transition in the plasma when detaching.

### A. Dynamic heat load response toward detachment

During a simulation, the heat load to the wall decreases as expected with increasing *n _{upstream}*. This continues until we reach a critical

*n*, as shown in Fig. 15. The power toward the wall continues decreasing, however

_{upstream}*n*suddenly starts decreasing.

_{upstream}*n*drops slightly and then increases again with almost a constant plasma heat load to the wall. In this particular simulation, the transition lasts for about 6 ms. The stagnation and drop of

_{upstream}*n*indicate a change in the fueling behavior from the divertor. The plasma flows from the divertor going upstream through the sol must have been changed. This then settles, and again start allowing for plasma fueling. By this point, the power on the IT has already been diminished, so this transition is mostly an effect on the OT power, as shown in Fig. 7.

_{upstream}Note, that the stagnation and drop of *n _{upstream}* in this transition is not seen in the steady-state results from Fig. 7. Again, as also shown in Fig. 12, the evolution of the power to the targets as function of

*n*in these dynamic simulations does not follow the same path as the steady-state results. This is due to the short timescale in which the gas puffing is scaled up, and equilibration time between up- and downstream. The neutral gas fueling rate increases more quickly in Fig. 15 than in Fig. 12, this is the reason the power to the divertor is lower for the same

_{upstream}*n*.

_{upstream}This transition (i.e., the stagnation, drop, and re-increase in *n _{upstream}*) occurs over many JOREK time steps and is well resolved in both space and time. If we study 2D results of the simulations, we observe the

*n*back-and-forth transition coincides with the formation and movement of a high-field-side high-density region.

_{upstream}### B. 2D comparison of detachment with and without *E* × *B* drifts

Figure 16 shows a comparison of the electron density in the poloidal plane with drifts (left) and without drifts (right) for comparable *n _{upstream}*. Both simulations are just detached (i.e., $ T e , target \u2264 5$ eV). Including the

*E*×

*B*drifts lead to qualitatively very different results.

*E*×

*B*drifts are essential for transporting plasma across the field and redistributing the plasma flow from the core and divertor. In Fig. 16 (right), without

*E*×

*B*drifts, the density is continuously building up; however, there is no mechanism of transporting it away from the divertor legs. This behavior is quite symmetrical and agrees qualitatively well with SOLPS-ITER without drifts.

^{18}As explained in Sec. VI, the

*E*×

*B*drifts push the OT plasma leg downwards and transport plasma via the private region from the OT to the IT. On the IT, the

*E*×

*B*drifts push the plasma upwards enhancing the cross field transport. With more fueling from the divertor, the HFS

*E*×

*B*vortex becomes dominant and starts trapping increasingly more density forming the so-called HFSHD-region.

The HFSHD-region is not a passive structure, it influences the flows in the divertor region while transitioning to a detached plasma state. This can be observed more clearly by examining the dynamic ion flux profiles toward the wall. We note that SOLPS simulations of ITER divertor plasmas with drifts at high SOL power and with impurities show a similar behavior with the appearance of an HFSHD-region because of similar flow patterns.^{29,31} In our simulations, however, we find that this phenomenology is driven by drifts and main plasma flow dynamics not specifically requiring high power operation and the injection of impurities, as described below.

### C. Dynamic ion flux response toward detachment with and without *E* × *B* drifts

Figure 17 shows the ion flux profiles as function of $ s \u2212 s sep$ progressing over time for comparable simulations with and without drifts. The color indicates *n _{upstream}*. This dynamic response shows how the plasma transition toward detachment. Figure 17 (bottom row) shows the ion flux profiles without drifts. The ion flux first increases with increasing

*n*. The peak at the IT occurs at a slightly lower

_{upstream}*n*than the OT peak. The progression of the ion flux over time and

_{upstream}*n*is qualitatively the same for the IT and OT. Both sides develop double peaked profiles. One could argue the plasma is not fully detached as the ion flux has not fully collapsed; however, for the goal of comparing the physics behavior, this is sufficient.

_{upstream}Figure 17 (top row) shows the ion flux profiles with *E* × *B* drifts. The dynamic response and qualitative results differ strongly from the case without drifts. *E* × *B* drifts favor distribution of ion flux toward the OT creating higher peak ion fluxes. The ion flux rollover point on the OT is around the same *n _{upstream}* for both simulations with and without drifts. The ion flux on the IT with drifts [Fig. 17 (top left)] has lower peak ion fluxes than its no-drift counterpart. With increasing

*n*, the IT ion flux diminishes more quickly and more upwards (to higher $ s \u2212 s sep$). At around $ n upstream \u2248 2.4 \xd7 10 19$ m

_{upstream}^{−3}, the IT ion flux sharply increases around 9 cm above the separatrix strike point. At even higher

*n*, the IT ion flux decays rather quickly again. During this transition in the plasma, the HFSHD-region starts forming and moving and this has a strong impact on the resulting ion flux at the IT.

_{upstream}The IT ion flux disappears completely at the high end of *n _{upstream}*, thus reaching full detachment, while this is not the case for the ion flux without drifts. In this 20 MW heating scenario, full detachment at the IT always goes paired with the formation of the HFSHD-region.

### D. High-field-side high-density region

For the ITER PFPO-1 H-mode 5MA 20 MW scenario, the JOREK simulations with *E* × *B* drifts always develop a HFSHD-region when going from the conduction-limited regime toward detachment. To fully understand the detachment behavior, we now look more closely on the HFSHD-region, its formation and its interaction with divertor and core plasma.

#### 1. Importance of the HFS E × B vortex

Figure 11 has shown that in an attached high-recycling scenario three *E* × *B* vortices form; one on the HFS, one the LFS and one in the private region, all in the vicinity of the X-point. In the lower range of *n _{upstream}*, the

*E*×

*B*drifts due to the vortices are weaker than the other plasma flows associated with a vertical divertor configuration and the associated recycling. With increasing

*n*, especially the HFS

_{upstream}*E*×

*B*vortex builds up in combination with the increased density coming from the IT strike point and this strongly impacts the divertor plasma flow.

Figure 18 shows the detached scenario with a fully formed HFSHD-region. It shows the electron density profile with the poloidal flow vectors colored with the *E* × *B _{factor}* (the ratio of the

*E*×

*B*flow amplitude divided by the poloidal component of the parallel flow $ | v E \xd7 B | | v \u2225 , pol |$) on Fig. 18 (left) and the electric potential contours on Fig. 18 (right). Specifically, the HFS

*E*×

*B*vortex has grown to become dominant as the $ E \xd7 B factor \u226b 1$. The HFS

*E*×

*B*vortex draws in plasma density from the IT, and this becomes trapped in the vortex. Once the HFSHD-region is stabilized, the electric potential forms a plateau where the density is highest. We observe the potential contour patterns from Fig. 11 become more exaggerated toward detachment. The LFS

*E*×

*B*vortex is lower and does not create a vortex flow pattern. The private region

*E*×

*B*vortex contributes to enhanced transport from the X-point toward the OT and decreases transport from the X-point to the IT. The private region

*E*×

*B*vortex is the dominant actor in transporting the plasma from the OT to the IT via the private region.

The transition toward a dominant HFS *E* × *B* vortex is a rather dynamic phenomenon. This transition is shown as a series of five 2D *n _{e}* profiles in Fig. 19. The arrows here indicate a general movement of the density. There is an initial slow buildup phase of density close to the wall. This phase is followed by a density front building up, moving from the wall toward the X-point. The

*E*×

*B*vortex is significant enough to move this high-density front. Once the front is close to the X-point, it quickly builds up more and more density (till $ n e \u2248 10 \xd7 10 20$ m

^{−3}). It then reaches a critical point when the high density region close the X-point suddenly moves away from the X-point and installs itself at the stable location shown in the last profile in Fig. 19 and also in Fig. 18. Specifically, during this fast movement away from the X-point the flow patterns everywhere are altered. This is most notable on the HFS of the divertor but can also be seen on the OT and in the entire SOL.

The buildup and movement of the HFSHD-region naturally influences its surroundings. This is why we see this dynamic phenomenon also in the ion flux profiles (Fig. 17) and even in the upstream density (Fig. 15). Due to the strong potential difference the HFSHD-region pushes and pull the plasma flow in the divertor legs.

The exact equilibrium position on which the HFSHD-region settles is somewhat dependent on the gas fueling rate in the divertor. Altering the gas fueling rate results in horizontal motion. Exactly how this relation works is not yet fully understood in our simulations.

These simulations only consider as source atomic neutrals, and no molecular species and/or impurities. Since molecules can change the momentum loss along the SOL and, together with impurities, increase energy losses, the formation and position of the HFSHD-region could change. However, the *E* × *B* vortex regions always form with plasma density build-up in the divertor and are, thus, almost unavoidable. Second, this simulation does not contain the diamagnetic drift (i.e., $ \u2207 p$ drift). The transport in the SOL will be altered with the addition of diamagnetic drifts. However, as the HFSHD-region forms during detachment in regions of low pressure, the diamagnetic drift is not likely to alter the transport significantly. The HFSHD-region itself, although containing very high densities, is a very low pressure region.

#### 2. Fueling from the divertor

The poloidal velocity profile changes strongly in the presence of the HFSHD-region. There is a clear qualitative difference in the $ v pol$ profiles between Figs. 11 and 18. This means, that with enough divertor gas puff, the flows in the divertor, and from the divertor to the core plasma will become significantly different. The ion mass flux in the poloidal plane $ \Gamma \u2192 pol = \rho v pol$, then tells us in this scenario how the core is fueled from the divertor. This subsection highlights the changes in $ \Gamma \u2192 pol$ after the formation of the HFSHD-region.

Figure 20 shows two schematics of the ion flux directions. Figure 20 (left) represents the attached plasma case, and Fig. 20 (right) the detached plasma after the formation of the HFSHD-region. Although the ion flux profiles are more complicated, this representation illustrates the general behavior. The color of the arrows indicates from where the flows originate, the size of the arrows do not indicate the magnitude; thus, they only trace the flow direction. The flux surface chosen are the separatrix and an approximate flux surface to indicate separation between the SOL and far-SOL regions.

Along the separatrix the ion flow from the core flows through (and slightly around) the X-point toward the strike point. Figure 20 (left) shows that due to *E* × *B* drift the ion flux is mostly distributed to the outer divertor leg. In the private flux region, a small part is redirected to the inner divertor leg. Another part of the private flux region transport comes from recycling on the outer target. The blue arrows indicate flows from the divertor plasma. These arise from plasma recycling and ionization of the gas puff. Figure 20 (left) shows the flux in the SOL comes mostly from the divertor plasma. Higher up along the SOL the flows are almost parallel to the separatrix, but it has a small component crossing the flux surfaces ( $ \Gamma \u2192 pol \xb7 \u2207 \psi \u2260 0$). Depending on the exact gas puffing, in the attached case, the flow inversion point on the separatrix is somewhere between the OMP and the top of the plasma.

After the HFSHD-region has formed, the $ \Gamma \u2192 pol$ profile has qualitative strongly changed. Figure 20 (right) shows in purple the HFSHD-region and with it the HFS *E* × *B* vortex. The recycled flux (cyan) from the IT is interrupted traveling upwards and circles back. The flow from the core (orange) travels from the X-point to halfway along the divertor plasma leg at the OT. Along the leg it gets deflected into the private region toward the IT. Due to multiple flow inversions the result is a flux along the inner divertor plasma leg toward the X-point that then goes into the private region. Generally, the flow patterns become more complicated.

These changes due to the HFSHD-region are not only confined to the divertor region. Due to the interrupted flow from the IT, there is no plasma traveling upstream in the HFS SOL. The flow reversal point along the separatrix is now approximately where the blue arrow crosses the separatrix in Fig. 20 (right). That entails that the flows in the SOL come for the OT and travel all along the SOL to finally reach the X-point from the HFS and the HFSHD-region. It should be noted that strong changes of plasma flows with the onset of detachment have been long seen in divertor tokamaks.^{55,56}

To observe whether the plasma retains this behavior on a longer timescale, some simulations have been extended to $ \u2248 150$ ms. From the typical timescales covered by JOREK, this simulation is completely steady-state. However, as in the simulation the upstream density was ramping up due to increased fueling, it takes of the order of a confinement time before this is fully equilibrated with the core of the plasma at the magnetic axis. The flows and effective fueling will still change on the timescale of a confinement time when considering equilibration of the core plasma density or the applied fueling.

#### 3. Ionization patterns in the presence of HFSHD-region

As a final dive investigating detachment and the HFSHD-region with our kinetic neutral model, we investigate the ionization patterns and the neutral density distribution. In the method, how the kinetic particles are coupled to JOREK, moments of the particle distribution can be projected to the same finite element grid. The ionization rate moment in Fig. 21 (left) is normalized and, for visibility reasons, both plots only show four orders of magnitude. The neutral density is presented as the number density (atoms/m^{3}).

As this simulation only fuels kinetic neutrals from the divertor region, there are only two ways by which neutrals can be found outside the divertor region. One is by recombination of the core plasma, and more importantly by transport across the divertor and SOL plasma. In the attached plasma case, nearly all neutrals are ionized and do not cross the divertor plasma leg. *E* × *B* drifts transport the plasma higher up the inner target where they undergo wall-assisted recombination. This provides the main source of neutrals outside of the private region. Figure 21 (right) show that in the detached case there is a significant neutral density close to the inner target at the height of the X-point. The ionizing region is still close to the wall on IT. *T _{e}* at the IT strike point is approximately 1 eV; this is low enough such that neutrals can transport past this region. In this particular case the outer target strike point temperature is approximately 10 eV, thus, not detached. On the OT, the ionization region is strongest close the wall and, due to the even higher temperature in the OT, SOL transport do not transport neutrals significantly upwards on the OT.

It is expected that when we introduce molecular hydrogen to the kinetic neutral model, we will obtain a more realistic ionization front. Since molecular hydrogen–plasma interaction introduces many more energy loss channels, it is expected that the plasma will detach at lower *n _{upstream}*. Specifically close to the strike points, molecules should lower

*T*for the 20 MW scenario since this is a rather low power level for these ITER plasmas. For the ITER Q = 10 scenarios with 100 MW, solely adding neutrals shall not be sufficient to detach the plasma. For studying detachment in these conditions, the puffing of impurity species is needed such as Neon.

_{e}Although the HFSHD-region is a high density region, the pressure in this region is low. Ionization occurs most significantly on the edge of the HFSHD-region, where the temperature sharply increases again. This allows the HFSHD-region as well to gather neutral density.

Note that if a significant part of the gas fueling would be provided by the top gas valve (as shown in Fig. 1), this would cooldown the (far-) SOL temperature. This could then allow for stronger neutral transport across the OT strike point, although this is not found to have a major impact on divertor behavior according to Ref. 18.

## IX. CONCLUSION

Realistic divertor conditions are important to predict impurity sputtering and transport and to study the consequences due transient MHD activity on the divertor. We have derived a plasma-atomic neutral interaction model based on the effective rate coefficients (From OpenADAS); this model has been implemented in JOREK as both a fluid neutral and a kinetic neutral model. The fluid neutral model evolves only the neutral mass density equation. The kinetic neutral model evolves individual superparticles in 3D geometry and includes ionization, recombination, and charge-exchange. The kinetic model further includes plasma and neutral boundary interactions based on the SDTRIM coefficients.

To benchmark our new models, we have performed a neutral gas puffing scan for an ITER early operations H-mode scenario (ITER PFPO-1, 5MA, 20 MW, H-mode) and compared it to SOLPS-ITER simulations (without poloidal drifts).^{18} To compare relatively similar physics models, we reduced the JOREK physics model by removing the evolution the magnetic fields and by removing *E* × *B* drifts. For all simulations, a divertor gas puff was used.

Using the fluid neutral model, we observe a power decay with increased gas puffing, but quantitatively, it does not match well. A one-variable fluid neutral model with spatially constant diffusivity is not sufficient for realistic divertor plasma modeling.

Using the kinetic neutral model, at low-to-medium upstream densities ( $ n upstream = 1 \xd7 10 19 \u2013 2 \xd7 10 19 m \u2212 3$), the total plasma heat load of the no-drift JOREK models agrees well (difference $ < 10 %$) with SOLPS-ITER simulation.^{18} At higher upstream densities ( $ n upstream > 2 \xd7 10 19 m \u2212 3$), the transition to detachment occurs at a $ \u223c 10 %$ higher upstream density and is less sharp than the detachment transition in SOLPS-ITER. The JOREK kinetic neutral model only includes atomic species. The extra energy loss channels provided by molecular species are not included, and thus, JOREK obtains detachment condition ( $ T target < 5 \u2009 eV$) at higher upstream densities than SOLPS-ITER. The JOREK peak heat fluxes in the conduction-limited regime are approximately twice as high as SOLPS-ITER simulation. The exact underlying reason for this is not fully identified, but it is expected to be due to the use of a single temperature model that does not distinguish electrons and ions for parallel transport in the SOL.

The kinetic neutral model using $ \u223c 10 6$ simulation particles is as fast as the neutral fluid model and scales more favorably with increasing computational nodes.

As a first application of the developed kinetic neutral model, we have looked at the effects of *E* × *B* drifts on the divertor solutions and the plasma heat load. The electric field forms vortices (on either side of the X-point and in the private flux region), leading to predominantly poloidal *E* × *B* flows. In the steady-state, *E* × *B* drifts result in a faster, compared to no-drifts, decrease in power to the wall as a function of the outer mid-plane separatrix density. This comes paired with a more asymmetric distribution of power between the inner and outer targets and allows the inner target to detach more easily with a rapidly diminishing heat flux. At higher upstream density, $ n upstream \u2248 2.2 \xd7 10 19 m \u2212 3$, a phase-like transition in the plasma develops. This coincides with the formation of a high-field-side high density (HFSHD)-region. This phase-like transition only develops with *E* × *B* drifts.

Understanding this transition requires time-dependent simulations since the formation of the HFSHD-region is a transient phenomenon at a timescale of $ 10 \u2212 20$ ms. While the plasma density in the divertor is ramping up, an *E* × *B* vortex close to the X-point on the HFS is becoming dominant. Simultaneously, *E* × *B* flows enhance cross field transport upwards at the IT. Plasma density is collected in the center of the HFS *E* × *B* vortex and the newly formed HFSHD-region moves horizontally toward the X-point and back, reaching a stable position. The HFS *E* × *B* vortex alters the poloidal flows not only in the divertor region, but in the SOL as a whole. This may have consequences for plasma fueling.

Next steps to further improve the divertor solutions in JOREK include the addition of hydrogenic molecules. JOREK models already exist including two temperatures, diamagnetic flows, and radiative kinetic impurities, but have not been used here. The effects of these will be left for further studies. The molecules increase the energy and momentum losses in the divertor changing the precise upstream density at which detachment occurs. Since the HFSHD-region is a low pressure region, the diamagnetic drift is not expected to have a strong influence there. On the contrary, for the H-mode edge plasmas in which large pressure gradients are expected, diamagnetic flows should play an important role.

With the implementation of a kinetic neutral model and, thus, more realistic divertor solutions, we have improved the capability to study the consequences of MHD activity (ELMs and disruptions) on the divertor plasmas and associated power fluxes including the application of 3D fields for ELM control which render the divertor plasma 3D as well.

## ACKNOWLEDGMENTS

This work was part of the research programme (No. 19KE01) of the Foundation for Nederlandse Wetenschappelijk Onderzoek Instituten (NWO-I), which was part of the NWO and has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200 – EUROfusion) and by the ITER Organization under Implementing Agreement No. 43-2174 with the Eindhoven University of Technology (TU/e) and with TU/e work order No. 107640. ITER is the Nuclear Facility, INB No. 174. Computations were performed using computer clusters Marconi-Fusion and the ITER Scientific Data and Computing Center (SDCC). The views and opinions expressed herein do not necessarily reflect those of the European Commission nor the ITER Organization nor NWO-I. Discussions with members of the ITER Organization (R. A. Pitts, S. D. Pinches, and X. Bonnin) are acknowledged.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Sven Quirijn Korving:** Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (lead). **Guido T. A. Huijsmans:** Funding acquisition (equal); Methodology (equal); Resources (equal); Software (equal); Supervision (equal); Writing – review & editing (equal). **Jae-Sun Park:** Data curation (equal); Methodology (equal); Resources (equal). **Alberto Loarte:** Funding acquisition (equal); Methodology (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

Raw data were generated at the ITER Scientific Data and Computing Center (SDCC) and Marconi-Fusion large scale facility. Derived data supporting the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: CONSERVATION EQUATIONS IN JOREK FORM

This appendix contains the JOREK form of the conservation equations used in Sec. III. JOREK uses a different form of the conservation equations than the Eularian form. As we rewrite the terms the JOREK form of the equations are still conservative. The assumptions used in JOREK are not included in this sections. The projections of the momentum equations to obtain the reduced MHD model are not included either.

##### 1. Mass conservation

##### 2. Momentum equation

##### 3. Energy equation

###### a. Ionization

###### b. Recombination

## REFERENCES

*E*×

*B*drifts in the scrape-off layer of a divertor tokamak

*E*×

*B*drift in high recycling divertors on target asymmetries

*E*×

*B*drifts for particle and heat transport in divertors

*Reviews of Plasma Physics*

*The Plasma Boundary of Magnetic Fusion Devices*