Development of the neutral model in the nonlinear MHD code JOREK: Application to E × B drifts in ITER PFPO-1 plasmas

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⁸ 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¹² and ELM control by RMPs⁵ and by pellet triggering.¹³ 

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¹⁴ 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.¹⁵ Kinetic particles have also been used to simulate the fast particle destabilization of TAE modes.¹⁶ 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.¹⁷ 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²³ 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.³⁶ 

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

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

The physics variables are Discretized using finite elements based on bicubic Bézier surfaces¹⁰ in the poloidal plane and by a Fourier series in the toroidal direction. These higher order elements provide $G 1 −$ continuity, i.e., continuity of the variables and their gradients in the R,Z plane. Even higher order continuity is available.³⁷ 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³⁸ (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.

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.¹⁷ 

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.

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 − + i + → 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.³⁹ 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.

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 $ρ ≡ m i n i + m e n e$⁠, and n_e = n_i in the absence of impurities. We then assume the limit of $m e → 0$⁠. Then, $ρ = 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.

The kinetic neutral model is part of the kinetic particle framework developed in JOREK.¹⁴ 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.

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.

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_e), it is possible to improve the order of accuracy by making the coupling terms (semi-) implicit in future work.

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 (⁠ $Δ 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.

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, $ϕ = ϕ plane + 1 2 U rng [ 0 , 1 ] Δ ϕ$ to place it in the sector associated with the poloidal plane with coordinate $ϕ plane$ and distance between planes $Δ ϕ$⁠.

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.¹⁸) 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 , I p = 5 M A$⁠, P_heating = 20 MW, $D ⊥ = 0.3$ m²s⁻¹, $χ i , ⊥ = χ e , ⊥ = 1.0$ m²s⁻¹. 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 ≈ 2 × 10 2$ m²s⁻¹. For the kinetic neutrals, such an assumption is not needed. To reach a quasi-steady-state, JOREK simulation were run between $11 and 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.¹⁸ 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 $ψ , j , u , w , ρ , T = T i + T e , v ∥$⁠, and ρ_n, the neutral density for the addition of fluid neutrals.

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, $ϕ$⁠) and local finite element coordinates.

As the fueling rate scan in the study by Park et al.¹⁸ 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 $ψ , 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.¹⁸ 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.

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_upstream at a given fueling rate is lower in JOREK simulations compared to SOLPS-ITER. This initial drop in n_upstream is related to the difficulty of keeping the shape of the initial density profile, which is assumed 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 ⪅ 1.5 × 10 19 m − 3$⁠).

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 ⪅ 2 × 10 19 m − 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 ≈ 2 × 10 19 m − 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_upstream than SOLPS-ITER. The transition as function of n_upstream is less sharp in JOREK. The slope (⁠ $∂ Power / ∂ n upstream$⁠) is decreasing for higher n_upstream This does indicate a similar but less strong form of density saturation from $n upstream ⪆ 2.4 × 10 19 m − 3$⁠.

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.⁴⁹ 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 , ⊥$⁠). Figure 5 shows $q p k , ⊥$ 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.⁵¹ Generally, we observe JOREK having a higher $q p k , ⊥$⁠. For the same plasma heat load, this indicates that JOREK profiles are narrower than SOLPS-ITER. The slope (⁠ $∂ q p k , ⊥ / ∂ 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 , ⊥$ is lower than the outer target. At the low end of n_upstream, JOREK has a higher inner $q p k , ⊥$⁠. 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 $κ ⊥ , 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.

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 $∼ 10 6$ simulation particles using a time step of $2 × 10 − 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 – 10 000$ cpuh, which equates to $60 – 70$ h wall-time on said four nodes.

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 ∥$⁠. 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 ∥$ 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.⁵² 

The kinetic particles Boris pusher used for advancing the particles in time preserves the particle energy conservation down to machine precision.⁴³ 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³⁸ (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 $∼ 10 − 4$⁠. However, the relative error in energy is $∼ 10 − 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 – 10 7$ superparticles with an average weight of $5 × 10 13 – 5 × 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 − 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 $ε ≈ 10 − 3 – 10 − 2$⁠. The error depends on the time step size and the amount of simulation particles. For these applications this is a sufficient accuracy.

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_upstream, the 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_upstream. Specifically, the 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 ≈ 2 × 10 19 m − 3$⁠, the plasma is already detached from the inner target. While detachment behavior (roughly $T e , target < 5 eV$⁠) without drifts starts at much higher n_upstream.

Around $n upstream ≈ 2 – 2.2 × 10 19 m − 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_upstream. Although at low n_upstream the total power is the same with and without drifts, the heat flux profiles have changed. At the low n_upstream, 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 E × B drifts. Counter-intuitively with increasing n_upstream $q p k , O T$ increases up until $n upstream ≈ 1.5 × 10 19 m − 3$⁠. At higher n_upstream, $q p k , O T$ decreases similar to its no-drift counterpart until the phase-like transition of the plasma.

The ion flux through the boundary in JOREK is $Γ ∝ ρ v ∥ · n ̂$⁠, where the boundary parallel velocity is set to the sound speed because of the Bohm sheath criterion, $v ∥ = c s$⁠. This then gives us $Γ ∝ ρ 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_upstream, the power redistributes differently between the IT and OT. At higher $n upstream ≈ 1 × 10 19 m − 3$⁠, the location of and height of the peak ion flux is different. With increasing n_upstream, we see the qualitative behavior more diverging, i.e., the location and height of the peak, single- vs double-peaked, the rollover and detachment.

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 − 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 $≈ 1$ cm from the wall, the left peak is due to a combination of peaked n_e and low $v ∥$⁠, and the right peak is at lower n_e, but $v ∥$ 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 × 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 ≈ 1.8 × 10 19 m − 3$ while the OT start seeing rollover around $n upstream ≈ 2.15 × 10 19 m − 3$⁠.

Remarkably, at around $n upstream ≈ 2.3 × 10 19 m − 3$ in the drifts simulation, there are two different profiles at almost equal n_upstream (⁠ $2.28 × 10 19$ and $2.3 × 10 19 m − 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 − s sep$⁠) then moves $≈ 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.

The plasma heat flux through the boundary in JOREK, $q ∝ γ sheath ρ T v ∥$⁠. With $v ∥ = c s$⁠, this gives $q ∝ γ sheath ρ 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 ∝ ρ T 3 / 2$ and the ion flux $Γ ∝ ρ 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 − drift , I T$ decreases faster with increasing n_upstream than $q n o − drift , O T$⁠; however, qualitatively they behave similarly, and q_peak moves progressively upwards (higher $s − s sep$⁠). As with with the ion flux in Fig. 9, the E × B drifts compresses and pushes the OT heat flux profile downwards (lower $s − 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 − s sep$ than the ion flux.

At $n upstream ≈ 2 × 10 19 m − 3 , q I T , n o − 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 ≤ 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 ≈ 2.3 × 10 19 m − 3$⁠, the heat flux profile moves first downwards (to lower $s − s sep$⁠) and upwards. The same as for the ion flux results (Fig. 9).

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).⁵³ 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 × B factor ≡ | v E × B | | v pol − v E × 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 ≈ 1.95 × 10 19 m − 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.

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.

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 $≈ 8.3$ MW of power in the first 5 ms. Simultaneously, the upstream density increases from $∼ 1$ to $∼ 1.5 × 10 19$ m⁻³. 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 $∼ 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.

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

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.

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_upstream, a strong increase occurs into the ion flux as the plasma starts developing to a high-recycling regime.

In the range below $n upstream ≈ 2 × 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.

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 ≤ 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_upstream indicating a movement of the divertor plasma density in the divertor and at higher n_upstream a decay of the ion flux. Figure 10 shows the almost complete decay of the heat flux with high n_upstream. In all these figures, we observe a strong dependence on 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.

During a simulation, the heat load to the wall decreases as expected with increasing n_upstream. This continues until we reach a critical n_upstream, as shown in Fig. 15. The power toward the wall continues decreasing, however n_upstream suddenly starts decreasing. n_upstream 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 n_upstream 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.

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_upstream 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 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_upstream back-and-forth transition coincides with the formation and movement of a high-field-side high-density region.

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 ≤ 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.¹⁸ 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.

Figure 17 shows the ion flux profiles as function of $s − 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_upstream. The peak at the IT occurs at a slightly lower n_upstream than the OT peak. The progression of the ion flux over time and n_upstream 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.

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_upstream, the IT ion flux diminishes more quickly and more upwards (to higher $s − s sep$⁠). At around $n upstream ≈ 2.4 × 10 19$ m⁻³, the IT ion flux sharply increases around 9 cm above the separatrix strike point. At even higher n_upstream, 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.

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.

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.

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_upstream, especially the HFS 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 × B | | v ∥ , 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 × B factor ≫ 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 ≈ 10 × 10 20$ m⁻³). 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., $∇ 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.

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 $Γ → pol = ρ v pol$⁠, then tells us in this scenario how the core is fueled from the divertor. This subsection highlights the changes in $Γ → 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 (⁠ $Γ → pol · ∇ ψ ≠ 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 $Γ → 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 $≈ 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.

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³).

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

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.

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).¹⁸ 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 × 10 19 – 2 × 10 19 m − 3$⁠), the total plasma heat load of the no-drift JOREK models agrees well (difference $< 10 %$⁠) with SOLPS-ITER simulation.¹⁸ At higher upstream densities (⁠ $n upstream > 2 × 10 19 m − 3$⁠), the transition to detachment occurs at a $∼ 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 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 $∼ 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 ≈ 2.2 × 10 19 m − 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 − 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.

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.

The authors have no conflicts to disclose.

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

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.

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.