Effects of recycling neutral on density shoulder formation in tokamak plasmas

Towards a physical understanding of the formation of flattened upstream scrape-off-layer (SOL) density profiles, namely ‘density shoulders’, a self-consistent one-dimensional radial transport model has been developed to estimate the upstream profiles covering both the core plasma and SOL region at the tokamak midplane. For the SOL region, the effective density and temperature profiles for the ionization process are obtained by the weighted averaging of the upstream and downstream profiles, which can distinguish the open-target operation from the closed-target operation by a weighting factor. Compared with enhanced turbulent convective transport, it is complementary for the model to study the competition between the effective source S eff and the parallel particle loss L SOL. It indicates that: (1) an appropriate S eff intensity controlled by the neutral pressure due to divertor or wall recycling and (2) an appropriate S eff peak position in a far SOL region adjusted by the plasma current as well as the weighting factor could offset the damping effect of L SOL on the density profile. Then S eff over L SOL in a far SOL region could be the sole process involved in bringing about SOL density shoulders.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. divergent state, requiring more detailed experimental and modeling studies.
Basically, the SOL density radial profile in steady state would result from a balance among the three main processes (namely the ionization source, parallel loss and perpendicular transport) [26]. The upstream SOL shoulder formation could be simply characterized by a substantial increase of the far SOL density radial e-folding length λ n = |n e /∇ r n e | [9,17,27], however, where the corresponding mechanism is not explicit. According to previous research, the density shoulder with flattened SOL profile is believed to be a consequence of the increased ratio between perpendicular and parallel transport [18,28]. As for the perpendicular transport, the advection rather than diffusion transport is the likely mechanism, in close compliance with relevant experiment results [12,29,30]. The related analyses imply that the advective transport is correlated with large amplitude fluctuations [31,32] manifested as filamentary structures in the SOL [33][34][35]. Meanwhile, it is noteworthy that in particular the intermittent filamentary structures (i.e. plasma blobs) are found to have strong interactions with the neutrals [36][37][38]. Moreover, since the density shoulder usually occurs at high density and low current [8,9,18,25], a common conjecture has been concluded that the filament characteristics could be changed via electrical disconnection from the divertor target sheath corresponding to the high divertor collisionality Λ div [39] for a high plasma density, and ultimately the advective transport could be improved [27,40,41]. Although collisionality seems to be the key factor, there are still some experimental phenomena that deviate from the expected dependence of λ n on Λ div [42], which implies that Λ div might be a necessary but not sufficient condition to obtain a flatter density profile [22,23,42]. As an important complement of the previous studies, which centered on upstream SOL turbulence, more recent reports on JET [9], AUG [24] and EAST [25] plasmas emphasized that divertor recycling or neutral pressure may be the sole process involved in SOL density shoulders, where JET plasmas in particular indicated that the character of shoulder formation has completely changed for different divertor configurations (vertical target (VT) or horizontal target (HT)). So, in short, the three main processes (again, the ionization source, parallel loss and perpendicular transport) are not independent, but affect each other intricately, which brings challenges to revealing the physical mechanism behind shoulder formation.
Inspired by the analyses of neutrals on JET [9], AUG [24] and EAST [25] plasmas, the aim of this paper is also to explore the formation of density shoulders. But for simplicity disentangled from complex convective turbulent transport, the main content is to model a very simple idea as shown in figure 1, which is that, without additional radial convective transport, the density profile can be flattened in a certain SOL region where the effective particle source S eff is dominant over the parallel loss L SOL . To demonstrate this idea, a self-consistent one-dimensional radial transport model has been developed, whose basic features are as follows: (a) Compared with radial convective transport for the SOL region, this is a complementary modeling study that focuses on the competition between the ionization source and parallel loss. (b) The effective particle ionization source is governed by the main atomic and molecular collision processes, where the processes of gas puffing or recycling are reduced as boundary conditions of molecular transport. (c) The model covers both the core plasma and SOL region with the parameterized radial diffusion coefficient for the upstream transport at the tokamak midplane, where the transport coefficients are constrained by the typical experimental scaling of tokamak energy confinement. (d) The downstream plasma profiles near the target are obtained via the basic two-point-model module [26,43], which make it possible for the model to be used to investigate the contribution of upstream and downstream parameters to the effective ionization.
This paper is structured as follows. In section 2, the complete physical model is presented, including the radial particle and heat transport and neutral ionization physics. In section 3, Figure 1. Conceptual diagram for a mechanism of density shoulder caused by the competition between SOL particle loss term L SOL and effective ionization source term S eff . the numerical investigations for the dependence of the density shoulder on the effective particle ionization are demonstrated. Finally, conclusions and remaining issues are summarized in section 4.

The model structure
According to the toroidal symmetry of the tokamak plasma as well as the divertor configuration, the transport inside the separatrix with closed magnetic flux surfaces could be reduced to one-dimensional radial transport [44]. On the other side, the heat and particle fluxes from plasma core across the separatrix will encounter SOL physics with the open field line geometry. In the SOL region, in addition to a certain radial transport, the parallel transport along the magnetic field line and neutral transport cannot be neglectedble [26]. Therefore, these transports should be coupled with each other to make the model self-consistent. Considering that the density shoulder is located at the LFS SOL of the tokamak midplane, the solution domain of the model is limited to the outer SOL region as schematically shown in figure 2, where the VT and HT of the outer divertor correspond to the closed-target operation (CTO) and the open-target operation (OTO), respectively. Specifically, in the model, the downstream profiles contribute more to ionization for CTO and contribute less for OTO. The upstream and downstream are connected by the SOL parallel transport, and have been assumed to be governed by the twopoint model [26,43]. The details of the model structure are presented in the following sections.

The radial transport for the upstream plasma
Assuming that the density shoulder is dominated by bulk plasma species transport, the impurity transport has been neglected in the model for simplicity. Similar to the transport equations in the COREDIV code [45], in this model the radial transport equations are solved for the electron density n e and   (1) and (2): According to the quasi-neutrality condition, the main plasma ion density n i is approximately equal to n e , and the electron temperature is given by T e ≡ p e /n e . There are many methods for plasma heating [2,46,47], but in the model the effective heat source is assumed to be loaded on the electron transport channel by the term of Q eff = Q 0 exp −r 2 / 2σ 2 Q , where the profile normalized with the amplitude Q 0 has been shown in figure 3 for the deposition width σ Q = 0.25a. a is the minor radius.
In the transport equations, Γ n = −D e ∇ r n e = n e v eff r is the radial particle flux with its corresponding effective particle convection velocity v eff r . The radial heat fluxes in the electron and ion channel are given by the formula Γ p e = −n e χ e ∇ r T e + 5 2 p e v eff r . The radial transport is anomalous with the prescribed radial transport coefficients of the order of Bohm diffusion [44,45,47]. In the model, a simple parabolic dependence on the minor radius has been assumed for the radial thermal diffusivity χ e in the core region [48,49], as follows: where ρ = r/a is the normalized radius, and the form factor α χ = 1/3 is adopted in the model. According to the scalings of tokamak energy confinement [44,50] for both L-mode and H-mode discharge, the effective energy confinement time τ E,eff is proportional to the plasma current I p , when other characteristic parameters of the device are fixed. Typically, χ sep is defined as the radial thermal diffusivity at the separatrix with a typical value χ 0 = 1 m 2 s −1 for the reference plasma current I ref,χ . Meanwhile, the ratio between D e and the thermal diffusivity χ e has a simple fixed relation of 4D e ∼ χ e [47,48]. The radial transport coefficients in the SOL region (namely χ SOL or D SOL ) for ρ > 1 are set to be constant as shown in figure 3, which is necessary for the model to eliminate the influence of the transport coefficient on the distribution of the SOL density profile. Generally, reasonable density and pressure profiles can be obtained by the transport of core plasma inside the separatrix, and can be regarded as equivalent boundary conditions for the transport in the SOL region. Q e,SOL and L SOL are the heat and particle losses dominated by the parallel transport in the SOL region. Ifτ nu and τ qu are taken as the effective dwell time for the particles and the heat traveling from the upstream plasma to the downstream target [26], the explicit expressions for Q e,SOL and L SOL in equations (1) and (2) can be derived in the form Q e,SOL = −p e /τ qu and L SOL = −n e /τ nu [17,48,51,52]. In the SOL region, the connection length has a value in the order of L c ∼ πq 95 R 0 , where q 95 is the edge safety factor and R 0 is the major radius. Based on the local plasma parameters, the effective dwell times can be estimated as τ nu = L c / (Mc s ) and τ qu = L 2 c /χ eff e , where M is the effective parallel Mach number [17], c s is the thermal ion sound speed and χ eff e is the effective parallel thermal diffusivity of electrons in the SOL. According to the Spitzer-Harm thermal th,e υ e in the collisional limit and the freestreaming expression χ fs = v th,e L c π for kinetic motions in low collisionality, χ eff e is defined as χ eff e = χ SP e χ fs χ SP e +χ fs [53][54][55], where v th,e is the thermal velocity and υ e is the electron collisional rate. Additionally, it should be noted that the outer side of the SOL is shadowed by the limiter as shown in figure 2. In the limiter shadow, the transport form stays the same, but the connection length is assumed to be shortened to L c /3 for all the cases analyzed by the model. Finally, S eff is the effective electron source, which depends on the plasma parameter profiles and the background neutrals in the SOL. However, differing from the magnetic surface averaged profiles in the core plasma, the profiles in the SOL change gradually from upstream to downstream. Although a complete two-dimensional profile distribution for the whole SOL region cannot be given by the current model, the effective plasma profiles for the neutral ionization can be estimated by the two-point model [26,43], which will be discussed in the next section.

The particle and power balance in the SOL region
In the SOL region, the plasma flow toward the target may experience momentum loss due to the frictional collisions with neutrals, viscous forces and volume recombination. For simplicity, those processes can be included by introducing a momentum loss factor f mom . Then according to the basic two-point model [26,43,[56][57][58][59] without the effects of magnetic flux expansion, the relation between the upstream and downstream total pressures can readily be extended as equation (4) for the pressure balance: where n eu = n e and n d are the upstream and downstream plasma density, T eu = T e and T ed are the upstream and downstream plasma electron temperature, T iu = T i and T id are the upstream and downstream plasma ion temperature, respectively. The momentum loss factor f mom is a function of the target temperature [43,56,60]. Meanwhile, the parallel heat flux will strike on the target through the sheath heat transmission. Then for the electron transport channel, the power balance between the upstream and target can be expressed as equation (5): where q ed is the parallel electron heat flux density entering the sheath to the target and γ s ∼ 7.8 is the sheath heat transmission coefficient. According to the parallel heat transport, both the parallel heat conduction flux q e,cond and the parallel heat convection flux q e,conv contribute to the parallel electron heat flux, q e = q e,cond + q e,conv . Furthermore, for simplicity, the momentum loss is expected to be a function of the target temperature. For a first approximation, f mom could be estimated by an arbitrarily chosen form [43,60]: where A, T * and n are the fitting parameters. In particular, according to the report in references [43,61], (A, T * , n) mom−loss = (1, 0.8 eV, 2.1) is adopted to evaluate f mom for deuterium plasma without impurities. It should be noted that there is little or no momentum loss for T ed > 10 eV. For a pure hydrogenic plasma, the heat conduction coefficients for electrons are κ 0e ∼ 2000 [26,62], for q e W m −2 , T e, i (eV), L c (m), which give the corresponding parallel heat conduction flux for electrons:  By introducing the specific heat at constant pressure C p ≡ 2.5, the corresponding parallel heat convection flux for electrons is given by equation (8): For high collisionality, electrons and ions have the same target temperature T ed = T id , so the local collisional adiabatic sound speed can be expressed as c sd = 2eT ed /m i , which gives the downstream particle flux Γ d = n d c sd .
Since the relationship between the upstream and downstream radial profiles has been established by the two-point model, the weighted average profiles to the upstream and downstream can be obtained as the effective profiles for neutral ionization, which reduces the atomic and molecular processes in the SOL into one-dimensional radial processes. The contribution fraction of the upstream or downstream profile to the effective profile can be represented by introducing a specific Ultimately, the main rate coefficients for the effective source S eff become a function of the weighted density n αW and the weighted temperature T αW . It is necessary to clarify the physical meaning of this weighting factor α W in the model. The neutral particles from the gas puffing and the recycling of the first wall or divertor will be eventually ionized in the SOL region. If the ionization near the divertor target region is dominant for the CTO, then the downstream profiles should have a greater contribution to the effective source S eff associated with a relatively large α W . Otherwise, for an OTO with a relatively small α W , the ionization processes related to upstream parameters cannot be ignored.
Besides the coupling between the upstream and downstream plasma parameters, the transport of neutrals also should be considered in the model for the ionization sources towards self-consistency, which will be discussed in the next section.

The reduced transport of neutrals
Due to the gas puffing and the main-chamber or the divertor recycling [7,12,63], there should be an effective background neutral pressure in the SOL. For the operating conditions of the tokamak, the ionization of neutrals belongs to the paradigm of low-pressure gas discharge in a toroidal vacuum chamber.
To estimate the effective source S eff at the SOL, only the main atomic and molecular collision processes in plasma [26,37] are considered in the model, which yield three dominant neutral species: (1) hydrogenic molecules from the edge gas pressure, (2) cold neutral gas atoms (Frank-Condon atoms) from the dissociation of molecules, (3) hot neutral atoms from the charge exchange of the ion with the Frank-Condon atom and electron-ion recombination. Also for the toroidal symmetry, the reduced neutral radial transports are given by equations (9)- (11), which are necessary to estimate a self-consistent radial profile for the electron source S eff : where n M , n FC , n A are the density of the hydrogenic molecule, Frank-Condon atom and hot neutral atom, respectively. According to the database of the atomic data and analysis structure [64,65], the main rate coefficients of ionization σv ion , recombination σv rec , dissociation σv dis , dissociative ionization σv dion,M and charge exchange σv cex are estimated for the atomic and molecular processes, as shown in figure 4. Then the effective source has an approximate relation S eff = σv ion (n FC + n A ) n e + σv dion,M n M n e − σv rec n i n e , where n i is the ion density. The hydrogenic molecule has a loss term due to the dissociation reactions given by the relation Correspondingly, the source of the Frank-Condon atom and the source of the hot neutral atom are given by S FC = 2 σv dis n M n e + σv dion,M n M n e − σv cex n i n FC − σv ion n e n FC and S A = σv cex n i n FC + σv rec n i n e − σv ion n e n A . Furthermore, the neutral particle fluxes of the hydrogenic molecule, Frank-Condon atom and hot neutral atom are given by respectively. The diffusivities are given by where υ mean = σv cex + σv ion n i and υ dis = σv dis + σv dion,M n e are the mean interaction frequencies for the atom and molecule. v M = eT M /m M is the molecular thermal velocity with m M = 2m i , and [26] for the Frank-Condon atom and the hot neutral atom.
Taking the deuterium discharge as an example, it is assumed that the deuterium molecules have room gas temperature (T M ∼ 300 K ∼ 0.026 eV). The temperature of T FC ∼ 2.5 eV is adopted for Frank-Condon atoms, which is in the production range of 2 to 4 eV by electron-impact disassociation of molecules [26], while the hot atoms have the same temperature as background ions, T A = T i . Apparently, compared with the cold atom or molecule, the hot atom has a longer mean free path, λ A = v A /υ mean . As shown in figure 5, in a wide range of temperature, the values of λ A are in centimeter order when the density is greater than 1e 19 m −3 . In the SOL region, the connection length L c is in the order of tens of meters. Therefore, the condition of λ A L c is generally valid for the atomic and molecular processes.

The boundary condition and the free modeling parameter declaration
For the boundary conditions in the model, according to the solution domain as shown in figure 2, ρ = 0 is the position of the core boundary with zero fluxes, and ρ > 1 is for the SOL region. Due to the limiter shadow, ρ b is set to be the position of the outer boundary, where the range of 1 < ρ < ρ b is much bigger than the SOL width. Then, except the outer boundary condition of the molecule n Mb , the natural boundary conditions are applied for all the equations at the core and the outer boundaries. Due to the gas puffing and the recycling, an effective constant source from the background neutral pressure can be assumed for the hydrogenic molecule at the outer boundary. The modeling parameters in the model are grouped into three types, as shown in table 1. First, if taking some devices as typical analysis objects, the typical geometric parameters, heating power of transport and so on for type-I parameters are fixed by referring to the typical discharges. Then the edge safety factor q 95 can be estimated by the formula q 95 = q 95,ref I p for different devices, which is important for the connection length L c varying with the plasma current. Secondly, based on the profile distribution in figure 3 for simplicity, the model mainly analyzes the variation of the upstream density profile in the SOL with different boundary neutral pressure and plasma current, which requires that the radial diffusion coefficient in the SOL should be consistent with the usual experimental diagnosis. Therefore, the strategy applied in the model is to calibrate the type-II parameters with the experimental scaling law, where the reference current for separatrix diffusivity, I ref,χ , can be calibrated within a certain range of plasma current variation to guarantee the core transport governed by equations (1) and (2) to be broadly in line with the scaling of tokamak energy confinement. Last, the effective parallel Mach number, M, and the weighting factor, α W , are sorted to type-III parameters, which directly affect the parallel loss and effective profiles for the neutral ionization. Due to the fact that the mechanisms for determining their values are out of the scope of this study, the model will incompletely analyze their sensitivity to the upstream density profile by parameter scanning for a typical device discharge.

Calibration and sensitivity of corresponding parameters
In principle, the physical model above can be applied to any tokamak device. However, the specific device parameters in table 2 could help to visualize the modeling results.
Referring to a specific scaling law, the model can calibrate the reference current for separatrix diffusivity, I ref,χ , for different devices as shown in table 2, so that the energy confinement time of the model is consistent with the experimental scaling within the corresponding plasma current variation. The details are described as follows.
The global energy confinement time is defined by τ E,model = 1.5n(T i +T e )dV P net , where P net is the total power for transport. In the model, n i = n e = n i and T i = T e are assumed, which implies that half of the power will be consumed in the electron transport channel. Taking H-mode scaling (i.e. IPB98(y,2) scaling [44]) of deuterium plasma as a reference, it should be noted that the bulk radiation power P rad,bulk is not subtracted in the computations of total loss power P Loss for the confinement time τ E,98y2 ∝ P −0. 69 Loss . If the bulk radiation can be assessed by comparing the radiation normalized to P Loss , i.e. the radiation fraction f rad , then the total loss power can be estimated by P Loss = P net / (1 − f rad ). In reality, the f rad is within an uncertain range [70], but for simplicity, the compromise estimates of f rad are specified as listed in table 2. Finally, when I re f ,χ of  Where JET-like parameters in table 2 are adopted, for the low plasma current I p = 1.7 MA regime. different devices are calibrated and fixed as a reference for the thermal diffusivity in equation (3), the confinement time calculated by the model is consistent with the experimental scaling for different plasma currents and different separatrix densities, as shown in figure 6. Generally, this procedure could also be applied to L-mode scaling with a different reference current of I ref,χ , which does not change the radial transport properties, and the corresponding radial diffusivity χ sep is still in the order of 1 m 2 s −1 . The purpose of the calibration is to obtain more reasonable radial diffusion coefficients. Then, the model could further study the effect of effective neutral ionization on the density profile.
Besides χ sep in the model, the effective parallel Mach number M and weighting factor α W could directly affect the steadystate distribution of the final upstream density profile. Through parameter scanning for a specific device, figures 7 and 8 demonstrate the sensitivity of the density shoulder on these parameters, where for better comparability, all the steadystate density profiles reach the same separatrix density n sep by feedback control of the boundary molecular density n Mb . By definition, the amplitude of characteristic length, Max(λ n ), could indicate the degree of density flattening except for the situation of ∇ r n e 0 for a shoulder rolled over in the range of r sep < r < r lim iter . It is shown that the shoulder formation is sensitive to α W rather than to the increase of M and χ SOL /χ sep . There is a similar result in that if S eff + L SOL are large enough in a certain region of S eff + L SOL > 0, they will lead to a density shoulder or even shoulder roll-over. However, when the ionizations are going to be fully dominated by the downstream profiles for α W → 1, the S eff peak will be too close to the separatrix, or when the upstream profiles are going to be fully dominant for α W → 0, the S eff peak will be too close to the limiter, which both lead to an insufficient intensity of S eff + L SOL to form the density shoulder. Additionally, as shown in figure 8, the density shoulder is not sensitive to the change of transport coefficient ratio χ SOL /χ sep . But if the transport coefficient ratio χ SOL /χ sep is too small, it will weaken the formation of the density shoulder, which implies that the H-mode density shoulder would be much weaker than in L-mode if the transport coefficient χ SOL of the H-mode is smaller than that of the L-mode.
In the following sections, as a compromise, M = 0.2 is adopted as the effective parallel Mach number [17], and χ SOL /χ sep = 1 is specified in the model, but they do not strictly distinguish the scalings of tokamak energy confinement [44,50] for L-mode or H-mode discharge, where in fact the real H-mode physics have exceeded the scope of this study. Meanwhile, in order to illustrate the difference of the closure degree of the divertor, α W = 0.8 and α W = 0.5 are specified typically for the CTO and for the OTO, respectively.

Typical case of density shoulder formation in the upstream SOL
In the model, one way to control the plasma density regime in steady state is to change the outer boundary value of molecular density, n Mb , which is associated with the background neutral pressure due to the gas puffing or the recycling. For each increased value of n Mb , the model reaches a steady state. The upstream density could evolve into a flattened profile with the increasing separatrix density, which is more likely to occur in relatively low current conditions. In addition, taking the JET-like parameters in table 2 as an example, for the low plasma current I p = 1.7 MA regime with the weighting factor α W = 0.5, there is a shoulder-like profile formed by a flattened upstream density profile as shown in figures 9(a) and (b), where the near (r − r sep = 0-1.5 cm) and far (r − r sep = 1.5-3 cm) SOL are specified. The radial distributions of the density e-folding length, λ n = |n e /∇ r n e |, are shown in figure 9(c), whose amplitudes indicate the degree of the flattening. According to the definition, the 'near' SOL represents the region of short e-folding length near the separatrix and the 'far' SOL refers to the region relatively far from the separatrix with the longer e-folding length [9,17]. Furthermore, due to the downstream profiles of n d and T ed , the effective collisionality Λ div for the divertor region parameters could be estimated by the definition of Λ div = L c υ e Ω i C sd Ω e = L c υ e m e C sd m i [9,27], as shown in figure 9(d). Due to the step-up of n Mb , the increase in n d and the decrease oinT ed will increase Λ div dramatically. By averaging over the near and far SOL regions, the regional averaged collisionality Λ div can be obtained for each interval region, respectively. As shown in figure 10(a), for the low plasma current I p = 1.7 MA regime, there is a sharp increase in value of the λ n amplitude, i.e. Max(λ n ), at a critical point of Λ div for the OTO-like regime with α W = 0.5, which indicates that the flattening of the density profile has been enhanced. However, for the CTO-like regime with α W = 0.8 or for the high plasma current I p = 3.5 MA regime with α W = 0.5, the values of Max(λ n ) remain nearly constant when Λ div increases, which is consistent with the modeling result of the single exponential profile without flattening. In particular, as shown schematically in the topologies of HT and VT in figure 2, this result also agrees with the statement from the JET experiment that 'there is no shoulder formation for the vertical (closed) target operation, while there are shoulders formed for the horizontal (open) target operation' [9].
According to the definition, the sheath-limited (SL) regime is characterized by the isothermality along each individual flux tube, and the conduction-limited (CL) regime is characterized by significant parallel temperature drops due to the finite value of heat conductivity [26]. According to the temperature ratios between upstream and downstream at separatrix versus Λ div are shown in figure 10(b), there are the divertor plasma transitions from SL to CL conditions during the increase of divertor collisionality. But this transition trend is an inevitable consequence as the divertor target n d increases, and T ed reduces strongly (leading to Λ div increasing), which neither depends on the closure degree of the divertor nor on the shape of the upstream density profile. Another phenomenon is the saturation of core density. Compared to the CTO-like regime, the Greenwald density fraction of core line-averaged density f GW will saturate at a lower level for the OTO-like regime (for a fixed plasma current), as shown in figure 10(c). This is because the flattened SOL density profile will weaken the depth of the ionization source entering the core plasma. There may be more phenomena associated with density flattening, but for the current model, the causal mechanism for SOL density flattening must be included in the transport equation.

Understanding of the modeling of the density shoulder formation in the upstream SOL
In the model, the density profile is completely determined by the local transport equation. As a complementary modeling study on recycling, the model does not focus on turbulent evolution. Accordingly, the convection transport mechanism is not included in the model. It should be noted that the radial thermal diffusivity is of the order of 1 m 2 s −1 by the energy confinement calibration, which is inversely proportional to the plasma current but uniform in the SOL region. It can be considered that the uniformly distributed diffusion coefficient would not be the reason for the local flattening of the density profile. Therefore, according to equation (1), the only mechanism to flatten the SOL density profile in the model is the balance between the effective source S eff and the local particle losses L SOL .
From the same cases as shown in figure 10, the typical radial profiles of S eff and L SOL are shown in figure 11. For the purposes of comparison, three density regimes with different separatrix values of upstream density (i.e. low, medium and high n sep ) are selected for each case, where the same density regime has similar n sep between each case. Again, as shown in figures 11(a1)-(a3), the upstream density flattening in the far SOL only occurred at high density and low current for the OTO-like regime with α W = 0.5, while there is no significant evidence that the density profile is flattening in the other two cases.
The mechanism leading to the formation of the flattened density profile is primarily revealed in figure 11(b1), where a certain region in which the algebraic sum of S eff + L SOL should be significantly greater than a critical value |C| cri is formed when the separatrix density increases to a certain level, and this region is also associated with the far SOL. In contrast, the profiles of S eff + L SOL are generally below zero in the other two cases shown in figures 11(b2) and (b3), which means that the loss terms L SOL are dominant and lead to almost monotonic attenuation of the density profiles. In other words, the profile of S eff or L SOL itself alone is not sufficient for the upstream density shoulder formation.
More specifically, first, the shape of the S eff profile in the model is determined by the weighted density n αW and the weighted temperature T αW , which could be broadened and be peaked further away from the separatrix, as shown in figures 11(c1)-(c3). Second, the absolute value of loss term |L SOL | ∝ n e √ T e L −1 c is dominated by the temperature profile with shorter decay length, which always peaks at the separatrix and has a step increase in the limiter shadow region due to the shorter connection length, as shown in figures 11(d1)-(d3). Apparently, the profile of S eff + L SOL is dominated by L SOL near the separatrix, which corresponds to the near SOL region where the upstream density has a short e-folding length. In the limiter shadow region, the L SOL is usually dominant, and could block any mechanism to flatten the density profile. Therefore, between the near SOL and the limiter shadow region, S eff can become dominant due to the possibility that the S eff profile could peak away from the separatrix. It should be noted that T eu is generally higher and has a longer decay length than T ed . Compared to the CTO-like regime, the weighted pressure p αW = n αW T αW with a broader profile is consistent with a greater contribution of T eu to T αW , as shown in figures 11(e1) and (e2), which leads to a positive feedback mechanism that makes S eff + L SOL |C| cri strong enough to flatten the density profile in the far SOL region.
It should be noted that, in the model, the connection length L c and the radial separatrix diffusivity χ sep are both inversely proportional to the plasma current I p . Due to the parallel heat loss term Q e,SOL in SOL region, the decay length of T eu will decrease with the decrease in connection length for the high plasma current regime, where the corresponding p αW profile has a narrower distribution, as shown in figures 11(e1) and (e3). So relatively, the S eff profile will peak close to the separatrix for the high plasma current regime. Therefore, besides the CTO-like regime, the high plasma current regime with a shorter connection length could also inhibit the positive feedback mechanism for the density profile flattening.
In fact, it is difficult to manifest the positive feedback mechanism of shoulder formation and to quantify the trigger threshold |C| cri , especially for the previous cases with the discrete step-up of n Mb . In order to promote the understanding of the positive feedback mechanism, figure 12 illustrates the spatiotemporal evolution of the upstream density profile associated with other profiles by slowly and continuously increasing the boundary molecular density, n Mb , in time. The same result is found for the CTO-like regime with α W = 0.8; there is no sign of shoulder formation yet with the same ramp rate of n Mb . This implies that the threshold condition |C| cri or threshold pressure of neutral to achieve a density shoulder in the CTO-like regime will be much greater than that in the OTO-like regime. In particular, for the OTO-like regime with α W = 0.5, when the boundary value n Mb increases to a certain value, the effective ionization source increases significantly, and the corresponding density profile gradually broadens to form a shoulder and extends towards the limiter. Limited by the current understanding, the threshold of |C| cri in the model could be a research topic in future work.

The modeling extrapolation for ITER-like parameters
The mechanism that makes S eff + L SOL strong enough to flatten the upstream density profile is applicable to all device parameters shown in table 2. However, when the model is extrapolated to an ITER-like large-size device, the density flattening behaves differently from the current medium-size devices. Besides the main parameters listed in table 2, the limiter radius r − r sep ∼ 10 cm [16,71] is specified for the ITERlike plasma in the model. Basically, although there is a large gap between the limiter and the separatrix, the density profile flattening can expand towards the limiter, which could occur in relatively low density or high current regimes. The details of typical cases are as follows. Compared with the plasma current I p = 15 MA, I p = 9 MA could be a low current regime for ITER-like parameters. The common phenomena are similar to those for JET-like parameters, where the mechanism for density flattening is completely suppressed in the CTO-like regime with α W = 0.8, as shown in figures 13(a2) and 14(a2). Also, there are near and far SOL regions with longer widths for the flattened density profiles in the OTO-like regime with α W = 0.5, as shown in figures 13(a1) and 14(a1). But for the low current regime with α W = 0.5, the flattened density profile almost always exists or occurs in the initial low-density condition, which implies that high density may be a necessary but not sufficient condition to obtain a flatter density profile for a large-size device. The reason is explicit in the model that the flattening is only determined by the region where the S eff + L SOL > 0 is strong enough to shield the effect of particle loss, as compared between figures 13(b1) and (b2). In the same way, by the comparison between figures 14(b1) and (b2), it can be understood that the profile flattening could occur with increasing density to enhance the term of S eff for the high current regime, which means that the inhibition mechanism of high current would be ineffective when the ionization region is far away from the separatrix.
Furthermore, for the device with ITER-like parameters, the separatrix temperature could reach hundreds of eV under the normal conditions of gas puffing or recycling, which means that the loss term |L SOL | ∝ √ T e will play a leading role in a wider area to form a wider near-SOL region. According to the weighted temperature and pressure as shown in figures 13 and 14, the effective ionization temperature could be maintained about tens of eV in the far SOL region (r − r sep = 2-10 cm). Therefore, compared with the CTO-like regime (i.e. the loss term always dominates), the effective source term S eff would certainly be dominant in a region far away from the separatrix in the OTO-like regime, which finally leads to the flattening of the density profile.
It should be noted that the saturation of core density also occurs in the case with ITER-like parameters. Compared to medium-size devices (e.g. the cases with JET-like parameters in table 2), the core density for cases with ITER-like parameters saturated at a very low separatrix divertor collisionality Λ div,sep corresponding to low separatrix density and high separatrix temperature, as the Greenwald density fraction f GW shown in figure 15(a). Before the saturation of f GW , although the ratios of T eu /T ed sep keep increasing with the increase of Λ div,sep , the ratios are still around 1 for all the cases shown in figure 15(b), which implies that all cases operate in the SL regime, regardless of whether the density profile is flattened or not.
Finally, it should not be ignored that there is a side effect of a significant drop of f GW after the saturation, especially for the case with the flattened density profile, as shown in figures 15(a) and 10(c). The reason for the f GW drop is that almost all of the neutral particles are ionized in the far-SOL region, resulting in insufficient effective ionization into the core plasma. For the case with ITER-like parameters with high separatrix temperature, if the effective source S eff in the far SOL is too strong, the upstream density profile could roll over, which would more strongly prevent the penetration of neutral  particles into the core plasma. Therefore, with the normal gas puffing or recycling, the model cannot be used to study the high core density regime of f GW ∼ 1. However, the positive feedback mechanism to flatten the upstream density profile is mainly determined by the SOL physics, so the results above could provide some reference for predicting the shape of the upstream density profile for a large-size device.

Conclusions and discussions
To investigate the recycling effects on the upstream SOL density profile for tokamak divertor plasmas, the main atomic and molecular collision processes have been considered as the particle ionization sources in the model, where gas puffing or recycling from the wall or divertor has been reduced as an effective molecular pressure boundary condition. Coupled with neutral processes, a self-consistent one-dimensional radial transport model has been developed to estimate the upstream profiles covering both the core plasma and SOL region at the tokamak midplane. Meanwhile, the basic two-point model has been also coupled to the model for parallel transport in the SOL region to determine the downstream parameters near the target. The main conclusions and discussions in this paper are briefly summarized as follows: (a) Compared with turbulent convective transport, it is complementary for the model to study the competition between the effective source S eff and the parallel particle losses L SOL . The model has verified the mechanism shown in figure 1, which is that the algebraic sum of S eff + L SOL should be strong enough in a certain SOL region to form a flattened density profile. Although the model claims that a larger value of α W means a CTO-like regime, on the other hand, as shown in the figure 16, changing the value of α W actually changes the peak location of S eff . Since the absolute value of loss term |L SOL | ∝ n e √ T e L −1 c always peaks at the separatrix, it could more effectively offset the effect of the loss term on the density profile when the effective source peak of S eff is located in the far SOL region. (b) For the effect of plasma current in the model, besides the loss term |L SOL | directly increasing with the shorter connection length L c for the high plasma current regime, the characteristic decay length of the upstream temperature profile, λ T eu ∝ L c χ sep /χ eff e , will also decrease with the increase in plasma current for the shorter connection length L c . Therefore, the changing of the plasma current I p as the changing of the weighting factor α W is to change the peak location of S eff . Finally, an appropriate ionization source intensity associated with neutral pressure plus an appropriate peak position eventually leads to a density shoulder.
(c) The model could be extrapolated to an ITER-like largesize device without any additional assumptions. For a device with ITER-like parameters, the separatrix temperature could reach hundreds of eV corresponding to insufficient fueling into the core plasma by gas puffing or recycling. Although the mechanism for density flattening is suppressed in the CTO-like regime, compared to medium-size devices in the OTO-like regime, the effective ionization could be dominant further away from the separatrix, where the flattened density profile almost always exists in the initial low-density low-current condition or occurs in the high-current regime. Besides the divertor collisionality, it seems that high density and low plasma current also become necessary but not sufficient conditions to obtain a flatter density profile in the OTO-like regime for large-size devices. (d) It should be pointed out that the physical quantities mentioned in the model (such as the effective divertor collisionality Λ div , Greenwald density fraction f GW , and the transition from SL to CL) are not directly related to the density shoulder mechanism shown in figure 1. They may affect the density shoulder for other mechanisms, which is not within the analytical capability of this model.
It is necessary to clarify the limitations related to the recycling and the radial diffusion coefficient. First, since this is a one-dimensional model, it cannot distinguish the real divertor geometry, nor can it identify the specific location of neutral recycling. As a strategy in the model, the gas puffing or recycling from the wall or divertor have been reduced as an effective molecular pressure boundary condition, namely the boundary molecular density. The contribution of upstream and downstream profiles to ionization is distinguished by the weighting factor. Therefore, there is no real recycling process in the model, only the reduced effective boundary molecular density. Second, although the model does not consider radial convective transport, the model still needs a reasonable radial diffusion coefficient. Therefore, the model uses experimental scaling to calibrate the transport coefficient. The model is not currently able to strictly distinguish between L-mode or H-mode discharge, but can only give the calibration that the coefficients are in the order of 1 m 2 s −1 . Towards promoting an understanding of recycling mechanisms, the current model results could provide a reference for the 2D model with real divertor geometry in future.
Finally, it should be noted that it is almost impossible for a simple model to cover all the energy loss processes. Then, there are still several issues remaining for further improving the model. For the two-point model, the complex transport processes in the target region are simplified. Without the crossfield in two-dimensional transport, the conduction and convection of heat flux along the magnetic field line are reduced by associating with the connection length and the sheath heat transmission coefficient. As a compromise, the factor f mom is introduced to identify other synthesis processes of momentum loss, but the effects of impurities on those factors could give another multi-species correction. Therefore, the effects of impurity and detachment are beyond the capability of the model. Also, other specific core plasma heating methods might further distinguish the electron and ion transport channels, which could give a different ratio between the electron and ion temperature. Since the model analysis has included the main energy transport and loss mechanism, the results could be qualitatively valuable for the design reference of the device.
Grant No. BJPY2019A01 and the K C Wong Education Foundation, the Science Foundation of Institute of Plasma Physics, Chinese Academy of Sciences under Grant No. DSJJ-18-02.