Self-Consistent Establishment of the Bohm Criterion in Two-Ion-Species Plasmas Based on the Anisotropic-Ion-Pressure Plasma Fluid Scheme

Abstract

The Bohm criterion for two-ion-species plasmas in a scrape-off layer–divertor system is investigated using the anisotropic-ion-pressure plasma fluid scheme and the virtual divertor model, which does not require the boundary condition of the incident flow speed at a divertor target. The flow velocities of the two ion species at the target obtained from our fluid simulations agree well with those obtained from earlier particle-in-cell simulations for both collisionless and collisional cases. The so-called Bohm criterion can be self-consistently established without direct treatment of the sheath region. Examining the time evolution of the flow velocity at the target, we infer that the Bohm criterion is not related to sheath potential formation but is the stability condition for an equilibrium solution of the quasi-neutral plasma region.

1.  Introduction

Heat removal, helium exhaust, and impurity shielding are essential functions of the divertor in a fusion reactor. Many integrated simulation codes for edge regions have been developed for evaluating divertor performance [1–4]. These integrated codes basically adopt the Braginskii plasma fluid model [5] at least for the main ions (e.g. hydrogen isotopes) and electrons from the viewpoint of computational time. Accurately solving the plasma in the sheath region in front of the divertor plate generally requires a kinetic model; thus, the Braginskii plasma fluid model is used only for the quasi-neutral region, with boundary conditions applied at the sheath entrance.

The Bohm criterion for the incident flow speed during plasma–wall interactions is a boundary condition at the front of the sheath entrance. This criterion is derived by considering the conditions for forming a monotonic potential profile in the sheath region [6, 7]. For singly charged multi-ion-species plasmas, the Bohm criterion is described as follows [8]:

  

$$\sum _{\mathrm{j}}^{}\frac{n_{\mathrm{j}}}{m_{\mathrm{j}}{V}_{\mathrm{j}}^2-{\gamma }_{\mathrm{j}}{T}_{\mathrm{j}\mathrm{\parallel }}}\le \frac{n_{\mathrm{e}}}{T_{\mathrm{e}}}.$$(1)

Here, the summation on the left-hand side (LHS) is performed over the ion species j. The variables for the j-species ion, $n_{\mathrm{j}}$, $m_{\mathrm{j}}$, $V_{\mathrm{j}}$, ${\gamma }_{\mathrm{j}}$, and $T_{\mathrm{j}\mathrm{\parallel }}$ represent the density, mass, parallel flow velocity, polytropic index, and parallel temperature, respectively. The electron density and temperature are denoted by $n_{\mathrm{e}}$ and $T_{\mathrm{e}}$, respectively. For a single-ion-species plasma, Eq. (1) is reduced to a simple boundary condition on the ion flow velocity $V_{\mathrm{i}}$; $V_{\mathrm{i}}\ge C_{\mathrm{s}}=\sqrt{({\gamma }_{\mathrm{i}}{T}_{\mathrm{i}\mathrm{\parallel }}+T_{\mathrm{e}})/m_{\mathrm{i}}}$, where $C_{\mathrm{s}}$ denotes the sound speed. This condition has been extensively adopted in many integrated divertor codes. However, scrape-off layer (SOL)–divertor plasmas generally comprise multiple ion species. For example, a DT burning plasma consists of deuterium, tritium, helium, and intrinsic impurity ions. Extrinsic impurity gases such as argon are typically injected to cool the plasma. In addition, in detached (i.e., recombining) plasmas, hydrogen molecular ions play an important role [9–11]. In such a case, as evidenced in Eq. (1), the Bohm criterion only provides a relation among the incident flow velocities of all ion species and is not sufficient to formulate each of them. An early particle-in-cell (PIC) simulation study of two-ion-species plasmas demonstrated that the incident flow velocities depended not only on the ion mass ratio but also on Coulomb collisionality [8], making it difficult to formulate each incident flow velocity. For example, in the multifluid simulation code LINDA-NU, the ion flow velocities at the sheath entrance are given either by the collisionless or collisional limit condition by switching the two conditions depending on the degree of flow velocity relaxation [11]. Therefore, self-consistent modeling of the Bohm criterion considering the effects of Coulomb collisions is a key issue.

In previous research, we used SOL–divertor simulations for single-ion-species plasmas to demonstrate that the Bohm criterion is automatically satisfied by an anisotropic-ion-pressure (AIP) plasma fluid scheme coupled with the virtual divertor (VD) model [12, 13]. Supersonic flows have also been realized naturally in radiation-cooling divertor plasmas [12] and diverging-magnetic-field SOL plasmas [13]. In the AIP plasma fluid scheme, the parallel ion viscosity term, which is approximated as a second-order spatial differential term in the Braginskii plasma fluid model, is directly described in terms of the parallel and perpendicular ion pressures. Thus, the parallel plasma momentum transport equation becomes a hyperbolic partial differential equation that does not require an explicit boundary condition on $V_{\mathrm{i}}$ at the sheath entrance. Instead of imposing an explicit boundary condition, the VD model is incorporated to reproduce the effects of the sheath region on the quasi-neutral plasma region. In this model, the calculation region is extended beyond the sheath entrance or even the actual wall to add an artificial region (as shown in Fig. 1), and the artificial damping terms of the particle, momentum, and energy are set in this region. The Bohm criterion was satisfied without directly solving the non-neutral sheath region, indicating that it is not related to the condition of monotonic potential formation in the sheath region [14].

Fig. 1.  Steady-state profiles of (a) densities, (b) ion flow velocities, and (c) electrostatic potential in the collisionless case ($G_{12}$ = 0). The vertical dashed line at $s/L$ = 0.5 represents the sheath entrance.

Therefore, the combination of the AIP plasma fluid scheme and the VD model is expected to be effective for multi-ion-species plasmas to reproduce the Bohm criterion. In this study, as a first step, we applied this combination to singly charged two-ion-species plasmas, as reported in an earlier PIC simulation study [8]. The numerical model and calculation conditions are explained in Sec. 2. Section 3 presents the simulation results for the collisionless and collisional cases. In Sec. 4.1, by examining the time evolution of the ion flow velocities at the boundary, we discuss the Bohm criterion from the viewpoint of the stability condition for an equilibrium solution of the quasi-neutral plasma. The mechanism of flattening the sourceless-region profiles within the collisional limit is discussed in Sec. 4.2. Section 4.3 discusses plasma behavior in a collisional situation in which supersonic flow is driven by an electron temperature gradient. Finally, a summary is presented in Sec. 5.

2.  Numerical Model and Calculation Conditions

The AIP plasma fluid equations for a singly charged two-ion-species plasma (such as a DT plasma) are described as follows:

  

$$\frac{\mathit{\partial }{n}_{\mathrm{j}}}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }s}(n_{\mathrm{j}}{V}_{\mathrm{j}})=S_{\mathrm{p},\mathrm{j}},$$(2)

  

$$\frac{\mathit{\partial }}{\mathit{\partial }t}\left(m_{\mathrm{j}}{n}_{\mathrm{j}}{V}_{\mathrm{j}}\right)+\frac{\mathit{\partial }}{\mathit{\partial }s}(m_{\mathrm{j}}{n}_{\mathrm{j}}{V}_{\mathrm{j}}^2+p_{\mathrm{j}\mathrm{\parallel }})=S_{\mathrm{m},\mathrm{j}}+F_{\text{fric},\mathrm{j}-{\mathrm{j}}^{\prime }}-\frac{n_{\mathrm{j}}}{n_{\mathrm{e}}}\frac{\partial p_{\mathrm{e}}}{\partial s}.$$(3)

Here, the subscript j represents an ion species ($\mathrm{j}$ = 1 or 2). The time and parallel-to-B coordinates are denoted by $t$ and $s$, respectively. The unknowns include ion densities $n_1$ and $n_2$ and parallel flow velocities $V_1$ and $V_2$. The terms $S_{\mathrm{p},\mathrm{j}}$ and $S_{\mathrm{m},\mathrm{j}}$ denote the volumetric particle and momentum sources, respectively. The frictional force on the j-species due to collisions with j′-species $F_{\text{fric},\mathrm{j}-{\mathrm{j}}^{\prime }}$ works in proportion to the relative flow velocity $V_{{\mathrm{j}}^{\prime }}-V_{\mathrm{j}}$ (see Eq. (6)). The pressure-gradient force is described in terms of the parallel component $p_{\mathrm{j}\mathrm{\parallel }}$; thus, the parallel plasma momentum transport equation (Eq. (3)) is considered a hyperbolic partial differential equation. The electron density is given by the quasi-neutral condition; $n_{\mathrm{e}}=n_1+n_2$. The pressures satisfy the equations of state, $p_{\mathrm{j}\mathrm{\parallel }}=n_{\mathrm{j}}{T}_{\mathrm{j}\mathrm{\parallel }}$ and $p_{\mathrm{e}}=n_{\mathrm{e}}{T}_{\mathrm{e}}$, where both the parallel ion temperature $T_{\mathrm{j}\mathrm{\parallel }}$ and electron temperature $T_{\mathrm{e}}$ are assumed to be constant (i.e., isothermal) in this study for simplicity.

A symmetric one-dimensional (1D) SOL–divertor system with a homogeneous magnetic field is considered. Divertor targets (sheath entrances) are located at $|s|=L/2$, where $L$ represents the parallel connection length between the targets. In the SOL region inside the two X-points, $|s|\le L_{\mathrm{X}}/2$ ($L_{\mathrm{X}}<L$), a homogeneous particle source $S$ is located, which simply reproduces radial particle diffusion from the core plasma in this 1D system. A VD region of sufficient length is set beyond the sheath entrances to reproduce the sheath boundary condition [12, 13, 15].

The volumetric source terms are given as follows:

  

$$S_{\mathrm{p},\mathrm{j}}=\{\begin{array}{cc}{S}_{\mathrm{j}}\hfill & (|s|\le L_{\mathrm{X}}/2)\hfill \\ 0\hfill & (L_{\mathrm{X}}/2<|s|\le L/2)\hfill \\ -\frac{n_{\mathrm{j}}}{{\tau }_{\mathrm{V}\mathrm{D}}}\hfill & (|s|>L/2)\hfill \end{array},$$(4)

  

$$S_{\mathrm{m},\mathrm{j}}=\{\begin{array}{cc}0\hfill & (|s|\le L_{\mathrm{X}}/2)\hfill \\ 0\hfill & (L_{\mathrm{X}}/2<|s|\le L/2)\hfill \\ -\frac{m_{\mathrm{j}}{n}_{\mathrm{j}}{V}_{\mathrm{j}}}{{\tau }_{\mathrm{V}\mathrm{D}}}+\frac{\mathit{\partial }}{\mathit{\partial }s}\left(m_{\mathrm{j}}{n}_{\mathrm{j}}{D}_{\mathrm{m}}^{\mathrm{V}\mathrm{D}}\frac{\mathit{\partial }{V}_{\mathrm{j}}}{\mathit{\partial }s}\right)\hfill & (|s|>L/2)\hfill \end{array}.$$(5)

Here, ${\tau }_{\mathrm{V}\mathrm{D}}$ denotes an input characteristic time of the artificial damping terms in the VD region [12]. As reported in Ref. [12], both ends of the VD region are connected via periodic boundary conditions. The artificial viscosity $D_{\mathrm{m}}^{\mathrm{V}\mathrm{D}}$ is set to be finite only near the center of the VD region to connect the flow velocity profiles. For the frictional force, we adopt the following model that was used in the SOLPS-ITER code [16] and LINDA-NU code [11]:

  

$$F_{\text{fric},\mathrm{j}-{\mathrm{j}}^{\prime }}=m_{\mathrm{j}}{n}_{\mathrm{j}}{\nu }_{\mathrm{s},\mathrm{j}-{\mathrm{j}}^{\prime }}(V_{{\mathrm{j}}^{\prime }}-V_{\mathrm{j}}),$$(6)

  

$${\nu }_{\mathrm{s},\mathrm{j}-{\mathrm{j}}^{\prime }}=\frac{G_{{\mathrm{j}\mathrm{j}}^{\prime }}}{\frac{3}{4\sqrt{2\mathrm{\pi }}}\frac{\sqrt{m_{\mathrm{j}}}{T}_{\mathrm{j}}^{3/2}}{ln\mathrm{\Lambda }}{\left(\frac{4\mathrm{\pi }{ϵ}_0}{e^2}\right)}^2}{n}_{{\mathrm{j}}^{\prime }}\sqrt{\frac{m_{{\mathrm{j}}^{\prime }}}{m_{\mathrm{j}}+m_{{\mathrm{j}}^{\prime }}}}.$$(7)

The slowing-down frequency of the j-species due to collisions with j′-species is denoted by ${\nu }_{\mathrm{s},\mathrm{j}-{\mathrm{j}}^{\prime }}$. The notations $e$ and ${ϵ}_0$ denote the elementary charge and vacuum permittivity, respectively, and $\ln \mathrm{\Lambda }$ represents the Coulomb logarithm. The ion temperature $T_{\mathrm{j}}$ is assumed to be equal to $T_{\mathrm{j}\mathrm{\parallel }}$ for simplicity. The dimensionless parameter $G_{{\mathrm{j}\mathrm{j}}^{\prime }}$ is uniquely determined by the colliding ion species j and j′ and is of the order of unity in reality [17]. However, in this study, we vary $G_{{\mathrm{j}\mathrm{j}}^{\prime }}$ artificially as a control parameter to change collisionality and thus the dominance of $F_{\text{fric},\mathrm{j}-{\mathrm{j}}^{\prime }}$.

Hereafter, we normalize the physics parameters based on the thermal velocity of the light ion, $v_{1,\mathrm{th}}=\sqrt{T_{1\mathrm{\parallel }}/m_1}$; $V_{\mathrm{j}}$ is normalized by $v_{1,\mathrm{th}}$, $n_{\mathrm{j}}$ is normalized by a characteristic density $n_{\mathrm{c}}=S_1{L}_{\mathrm{X}}/\left(2v_{1,\mathrm{th}}\right)$, $t$ is normalized by a characteristic time ${\tau }_{\mathrm{c}}=L/v_{1,\mathrm{th}}$, and the collisionality is estimated by $L{\nu }_{\mathrm{s},12}/v_{1,\mathrm{th}}$.

Table 1 summarizes the calculation conditions. We select $m_2/m_1$ = 4 as reported in an earlier PIC simulation study [8]. For example, ion species j = 1 and 2 may correspond to H⁺ and He⁺, respectively. However, under the present calculation conditions, the electron mass does not affect the results as long as it is sufficiently smaller than the ion masses. Therefore, the analysis is not limited to the specific case of H⁺ and He⁺. The spatial and time resolutions, $\mathrm{\Delta }s$ and $\mathrm{\Delta }t$, are selected (Table 1) such that $\mathrm{\Delta }s\gg v_{1,\mathrm{th}}\mathrm{\Delta }t$ and $\mathrm{\Delta }s\ll v_{1,\mathrm{th}}{\tau }_{\mathrm{V}\mathrm{D}}$ (in the VD region) are satisfied. The present system is fully symmetric about $s$ = 0; thus, we only present the profiles in the $s$ ≥ 0 region in Sec. 3. However, in an earlier PIC simulation study [8], small asymmetry was observed. This asymmetric problem in the two-ion-species plasma is left for future work. The initial profiles at $t$ = 0 are uniformly given as $n_1/n_{\mathrm{c}}$ = $n_2/n_{\mathrm{c}}$ ≈ 0.536, and $V_1$ = $V_2$ = 0.

Table 1. Calculation conditions.

parameters	values
$L_{\mathrm{X}}/L$	1/2
$S_2/S_1$	1
$T_{1\mathrm{\parallel }}/T_{\mathrm{e}}=T_{2\mathrm{\parallel }}/T_{\mathrm{e}}$	1
$m_2/m_1$	4
${\tau }_{\mathrm{V}\mathrm{D}}/{\tau }_{\mathrm{c}}$	3.35 × 10−3
$\mathrm{\Delta }s/L$	6.25 × 10−5
$\mathrm{\Delta }t/{\tau }_{\mathrm{c}}$	6.70 × 10−6
$G_{12}$	0–100

3.  Results

3.1  Collisionless condition

First, we set $G_{12}$ = 0 to simulate a collisionless condition [15]. Figure 1 shows the steady-state profiles of (a) densities, (b) flow velocities, and (c) electrostatic potential. The densities are hilly in the particle source region ($s/L$ ≤ $L_{\mathrm{X}}/(2L)$ = 0.25) and are flat in the sourceless region ($s/L$ > 0.25). The flow velocities increase in the source region and are constant in the sourceless region. The potential is evaluated by the Boltzmann relation $e\varphi /T_{\mathrm{e}}=\ln (n_{\mathrm{e}}/n_{\mathrm{e}}(s=0))$. The presheath $\varphi$ profile emerges in the source region. This common presheath potential accelerates both ion species; thus, the ratio of flow speed, $V_2/V_1$, becomes approximately $\sqrt{m_1/m_2}$. Because the particle sources are set the same, $S_1=S_2$, the ratio of densities becomes inversely proportional to $V_2/V_1$, i.e., $n_2/n_1=\sqrt{m_2/m_1}$.

The steady-state profiles of $V_1$ and $V_2$ in the sourceless region are almost flat, with values of approximately $1.4v_{1,\mathrm{th}}$ and $0.7v_{1,\mathrm{th}}$, respectively. This is theoretically valid because there is no physical factor that causes gradients in any of the plasma parameters in this case. However, the ion flow velocity very near the sheath entrance $s/L$ = 0.5 is increased by numerical errors due to the finite grid size $\mathrm{\Delta }s$ effect [12]. Therefore, we select the data at $s/L$ = 0.49 instead. These values are plotted as a point on the $V_1$-$V_2$ plane (the cross point in Fig. 2).

Fig. 2.  Steady-state ion flow velocities at $s/L$ = 0.49 near the sheath entrance (cross) in the collisionless case ($G_{12}$ = 0). The solid lines represent the theoretical curves of the Bohm criterion (Eq. (8)). The dashed and dotted lines denote $V_2=\sqrt{m_1/m_2}{V}_1$ and $V_2=V_1$, respectively.

The theoretical curves of the Bohm criterion in two-ion-species plasmas have been reported in Ref. [8]. Equation (1) can be transformed into a handy formula as follows:

  

$$(V_1^2-C_{\mathrm{s},1}^2)(V_2^2-C_{\mathrm{s},2}^2)=\frac{n_1{n}_2{T}_{\mathrm{e}}^2}{n_{\mathrm{e}}^2{m}_1{m}_2},$$(8)

where $C_{\mathrm{s},\mathrm{j}}^2=(n_{\mathrm{j}}{T}_{\mathrm{e}}+n_{\mathrm{e}}{T}_{\mathrm{j}\mathrm{\parallel }})/(n_{\mathrm{e}}{m}_{\mathrm{j}})$ for the isothermal case (i.e., ${\gamma }_{\mathrm{j}}$ = 1). Equation (8) represents a hyperbola with asymptotic lines of $V_{\mathrm{j}}=C_{\mathrm{s},\mathrm{j}}$ as shown in Fig. 2. For convenience, we call one curve in the right upper region an “upper branch” ($V_{\mathrm{j}}>C_{\mathrm{s},\mathrm{j}}$), and the other curve in the left lower region a “lower branch” ($V_{\mathrm{j}}<C_{\mathrm{s},\mathrm{j}}$), as reported in Refs. [14, 15].

Next, we compare the simulation and theoretical curves. We find that the measured point ($V_1$, $V_2$) is located very near the intersection of the “upper branch” of the hyperbolic curves with the linear line of $V_2=\sqrt{m_1/m_2}{V}_1$. This result agrees well with the lowest collisionality case in the earlier PIC simulation study [8]. If the Bohm criterion is governed by the sheath formation condition, the intersection point of the “lower branch” with the $V_2=\sqrt{m_1/m_2}{V}_1$ line would also be observed for a steady-state solution. Contrary to this supposition, the “lower branch” solution was never obtained in either the present fluid simulations or the earlier PIC simulations. The branch selection mechanism is discussed in Sec. 4.1.

3.2  Collisional conditions

Next, we set $G_{12}$ = 1, 10, and 100 to investigate the effects of the frictional force. The collisionality estimated at $s/L$ = 0.375 (sourceless region) is $L{\nu }_{\mathrm{s},12}/v_{1,\mathrm{th}}$ ≈ 3, 30, and 300, respectively. Figures 3, 4, and 5 show the steady-state profiles in the $G_{12}$ = 1, 10, and 100 cases, respectively. As $G_{12}$ and frictional force increase, the difference between $V_1$ and $V_2$ reduces. In the case of $G_{12}$ = 100, $V_1\approx V_2$ holds, indicating that the plasma is sufficiently close to the collisional limit. The presheath $\varphi$ profile emerges even in the sourceless region (i.e., 0.25 ≤ $s/L$ ≤ 0.5). However, its dominance weakens as the plasma approaches the collisional limit. Consequently, in the $G_{12}$ = 100 case, the $\varphi$ profile looks similar to that in Fig. 1(c) (i.e., the collisionless case). The reason for this is discussed in Sec. 4.2.

Fig. 3.  Steady-state profiles of (a) densities, (b) ion flow velocities, and (c) electrostatic potential in the $G_{12}$ = 1 case. The vertical dashed line at $s/L$ = 0.5 represents the sheath entrance.

Fig. 4.  Steady-state profiles of (a) densities, (b) ion flow velocities, and (c) electrostatic potential in the $G_{12}$ = 10 case. The vertical dashed line at $s/L$ = 0.5 represents the sheath entrance.

Fig. 5.  Steady-state profiles of (a) densities, (b) ion flow velocities, and (c) electrostatic potential in the $G_{12}$ = 100 case. The vertical dashed line at $s/L$ = 0.5 represents the sheath entrance.

Figure 6 shows the relationship between the steady-state values of the ion flow velocities and the theoretical curves of the Bohm criterion at $s/L$ = 0.49 for all finite $G_{12}$ settings. As collisionality increases ($G_{12}$ ≫ 1), the point of ($V_1$, $V_2$) moves away from the $V_2=\sqrt{m_1/m_2}{V}_1$ line for the collisionless case and asymptotically approaches the direct proportion line, $V_2=V_1$, due to the strong frictional force. Within the collisional limit, the two-species-ion fluids behave as single-species-ion fluids with an effective temperature, $T_{\mathrm{i}\mathrm{\parallel },\text{eff}}=(n_1{T}_{1\mathrm{\parallel }}+n_2{T}_{2\mathrm{\parallel }})/n_{\mathrm{e}}$, and an effective mass, $m_{\text{eff}}=(n_1{m}_1+n_2{m}_2)/n_{\mathrm{e}}$. Note that the temperatures $T_{1\mathrm{\parallel }}$ and $T_{2\mathrm{\parallel }}$ are relaxed to a unified isotropic ion temperature $T_{\mathrm{i}}$ in the collisional limit. Thus, the plasma tended to satisfy the Bohm criterion for a single-ion species plasma; $V_1=V_2=C_{\mathrm{s},\text{eff}}=\sqrt{(T_{\mathrm{e}}+T_{\mathrm{i}})/m_{\text{eff}}}$. The point ($V_1=C_{\mathrm{s},\text{eff}}$, $V_2=C_{\mathrm{s},\text{eff}}$) is plotted as a triangle in Fig. 6 and moves away from the “upper branch.” This tendency also reproduces the earlier PIC simulation study, except for supersonic plasma flows driven from temperature gradients [8]. The plasma behavior under a finite electron temperature gradient is discussed in Sec. 4.3. The results presented in Secs. 3.1 and 3.2 demonstrate that at least in isothermal cases, the AIP plasma fluid scheme coupled with the VD model succeeded in the self-consistent establishment of the Bohm criterion for two-ion-species plasma in a wide range of collisionality.

Fig. 6.  Steady-state ion flow velocities at $s/L$ = 0.49 near the sheath entrance (crosses) in (a) $G_{12}$ = 1, (b) 10, and (c) 100 cases. The triangles represent the sound speed with the effective mass $C_{\mathrm{s},\text{eff}}=\sqrt{(T_{\mathrm{i}}+T_{\mathrm{e}})/m_{\text{eff}}}$. The solid lines represent the theoretical curves of the Bohm criterion (Eq. (8)). The dashed and dotted lines represent $V_2=\sqrt{m_1/m_2}{V}_1$ and $V_2=V_1$, respectively.

4.  Discussion

4.1  Equilibrium solutions of system equations and their stability

Assuming the steady state ($\mathit{\partial }/\mathit{\partial }t$ = 0) and sourceless and collisionless conditions, the system equations (Eqs. (2) and (3)) are reduced for equilibrium solutions (denoted by the subscript 0) as follows:

  

$$\left(\begin{array}{cccc}\hfill a_{11}\hfill & \hfill a_{12}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill a_{23}\hfill & \hfill a_{24}\hfill \\ \hfill a_{31}\hfill & \hfill a_{32}\hfill & \hfill a_{33}\hfill & \hfill 0\hfill \\ \hfill a_{41}\hfill & \hfill 0\hfill & \hfill a_{43}\hfill & \hfill a_{44}\hfill \end{array}\right)\left(\begin{array}{c}\hfill \mathit{\partial }{n}_{10}/\mathit{\partial }s\hfill \\ \hfill \mathit{\partial }{V}_{10}/\mathit{\partial }s\hfill \\ \hfill \mathit{\partial }{n}_{20}/\mathit{\partial }s\hfill \\ \hfill \mathit{\partial }{V}_{20}/\mathit{\partial }s\hfill \end{array}\right)=\mathbf{0},$$(9)

  

$$a_{11}=V_{10},a_{12}=n_{10},$$(10)

  

$$a_{23}=V_{20},a_{24}=n_{20},$$(11)

  

$$\begin{array}{c}{a}_{31}=n_{\mathrm{e}0}(m_1{V}_{10}^2+T_{1\mathrm{\parallel }})+n_{10}{T}_{\mathrm{e}},\hfill \\ a_{32}=2n_{\mathrm{e}0}{m}_1{n}_{10}{V}_{10},\hfill \\ a_{33}=n_{10}{T}_{\mathrm{e}},\hfill \end{array}$$(12)

  

$$\begin{array}{c}{a}_{41}=n_{20}{T}_{\mathrm{e}},\hfill \\ a_{43}=n_{\mathrm{e}0}(m_2{V}_{20}^2+T_{2\mathrm{\parallel }})+n_{20}{T}_{\mathrm{e}},\hfill \\ a_{44}=2n_{\mathrm{e}0}{m}_2{n}_{20}{V}_{20}.\hfill \end{array}$$(13)

The condition that the determinant of the coefficient matrix becomes zero results in the same equation as Eq. (8), which was derived from the analysis of the sheath region [6–8]. However, a viewpoint of the stability of equilibrium solutions can be introduced, as discussed in Ref. [14]. As shown in Fig. 2, although the intersection of the “lower branch” and the $V_2=\sqrt{m_1/m_2}{V}_1$ line is the other equilibrium, it is never numerically obtained. Figure 7(a) shows the time evolutions of $V_1$ and $V_2$ at $s/L$ = 0.49 in the collisionless case ($G_{12}$ = 0). The ratio of LHS to RHS in Eq. (8) at $s/L$ = 0.49 is also shown in Fig. 7(b). Equation (8) is satisfied instantaneously at $t/{\tau }_{\mathrm{c}}$ ≈ 0.03 and $t/{\tau }_{\mathrm{c}}$ ≈ 0.2 and asymptotically as the plasma approaches the steady state. At $t/{\tau }_{\mathrm{c}}$ ≈ 0.03, $V_1<C_{\mathrm{s},1}$ and $V_2<C_{\mathrm{s},2}$ hold, as illustrated in Fig. 7(a), meaning that the point of ($V_1$, $V_2$) is instantaneously on the “lower branch.” In contrast, at $t/{\tau }_{\mathrm{c}}$ ≈ 0.2, $V_1>C_{\mathrm{s},1}$ and $V_2>C_{\mathrm{s},2}$ hold, meaning that the point of ($V_1$, $V_2$) is instantaneously on the “upper branch.” In addition, the point of ($V_1$, $V_2$) asymptotically returns to the “upper branch” as the plasma approaches the steady state. As shown in Fig. 7, the point of ($V_1$, $V_2$) passes the “lower branch” equilibrium at $t/{\tau }_{\mathrm{c}}$ ≈ 0.03 and does not stay there, indicating that the “lower branch” equilibrium is unstable. Stability analysis of the equilibrium solutions by analyzing the real part of the eigenvalues of the Jacobian of this system will be performed in future studies.

Fig. 7.  Time evolutions of (a) ion flow velocities $V_{\mathrm{j}}$ and asymptotic velocities $C_{\mathrm{s},\mathrm{j}}$, and (b) the ratio of LHS to right-hand side (RHS) of Eq. (8) at $s/L$ = 0.49 near the sheath entrance in the collisionless case ($G_{12}$ = 0). The horizontal dashed line in (b) represents LHS/RHS = 1, and the vertical lines represent the timing at which Eq. (8) holds.

Based on the above discussion and the fact that the combination of the AIP fluid scheme and the VD model self-consistently establishes the Bohm criterion without directly solving the non-neutral sheath region, we can posit that the Bohm criterion is not related to sheath potential formation but is a condition of stable equilibrium solutions in the quasi-neutral plasma region.

4.2  Flattening of sourceless-region profiles within collisional limit

Here, we discuss the mechanism of flattening the sourceless-region profiles within the collisional limit, as shown in Fig. 5. Based on the system equations (Eqs. (2) and (3)), $\mathrm{\delta }V=V_1-V_2$ is roughly governed as follows:

  

$$V_1\frac{\partial \mathrm{\delta }V}{\partial s}\approx -\frac{e}{m_1}\frac{\partial \varphi }{\partial s}-{\nu }_{\mathrm{s},12}\mathrm{\delta }V.$$(14)

Based on the above equation, $\mathrm{\delta }V$ is driven by the presheath electric field in the source region ($s\le L_{\mathrm{X}}/2$). In the collisionless case ($G_{12}$ = 0), $\mathrm{\delta }V\approx 0.67v_{1,\mathrm{th}}$ at the X-point ($s=L_{\mathrm{X}}/2$). However, $\mathrm{\delta }V$ at the X-point reduces as $G_{12}$ and ${\nu }_{\mathrm{s},12}$ increase because of the second term: $\mathrm{\delta }V\approx 0.51v_{1,\mathrm{th}}$ for $G_{12}$ = 1 (Fig. 3) and $\mathrm{\delta }V\approx 0.21v_{1,\mathrm{th}}$ for $G_{12}$ = 10 (Fig. 4). In the sourceless region ($s>L_{\mathrm{X}}/2$), $\mathrm{\delta }V$ decays with the spatial scale of $V_1/{\nu }_{\mathrm{s},12}$. Thus, in the collisional limit, $F_{\text{fric},12}=-F_{\text{fric},21}=-m_1{n}_1{\nu }_{\mathrm{s},12}\mathrm{\delta }V$ converges to zero. Consequently, a collisional two-ion-species plasma can be regarded as a collisionless single-ion-species plasma. This is why the sourceless-region profiles in the collisional limit case ($G_{12}$ = 100) are flattened.

4.3  Supersonic flow driven by electron temperature gradient

Here, we discuss the effects of an electron temperature gradient on the two-ion-species plasma behavior in a collisional case (i.e., $V_1\approx V_2\approx V$). The following $T_{\mathrm{e}}$ profile is introduced only in this subsection:

  

$$T_{\mathrm{e}}\left(s\right)=\{\begin{array}{cc}{T}_{\mathrm{e}0}\hfill & (|s|\le L_{\mathrm{X}}/2)\hfill \\ \frac{T_{\mathrm{e}0}(L-2|s|)}{L-L_{\mathrm{X}}}+\frac{T_{\mathrm{ed}}(2|s|-L_{\mathrm{X}})}{L-L_{\mathrm{X}}}\hfill & (L_{\mathrm{X}}/2<|s|\le L/2)\hfill \\ T_{\mathrm{ed}}\hfill & (|s|>L/2)\hfill \end{array},$$(15)

where $T_{\mathrm{ed}}$ ($<T_{\mathrm{e}0}$) denotes the electron temperature at the sheath entrance ($|s|=L/2$). The electron temperature in the source region $T_{\mathrm{e}0}$ is set in the same way as in the other sections (i.e., $T_{1\mathrm{\parallel }}/T_{\mathrm{e}0}=T_{2\mathrm{\parallel }}/T_{\mathrm{e}0}=1$). By considering the particle and momentum balances and assuming that the Mach number at the X-point is unity (i.e., $V=C_{\mathrm{s},\text{eff}}=\sqrt{(T_{\mathrm{e}0}+T_{\mathrm{i}})/m_{\text{eff}}}$), which is justified later in Fig. 8, we obtain the following formula to estimate the Mach number at the sheath entrance $M_{\mathrm{t}}$ [8]:

Fig. 8.  Steady-state profiles of (a) densities and (b) ion flow velocities with $T_{\mathrm{ed}}=T_{\mathrm{e}0}/4$ and $G_{12}$ = 100. The profile of the sound speed with the effective mass $C_{\mathrm{s},\text{eff}}=\sqrt{(T_{\mathrm{i}}+T_{\mathrm{e}})/m_{\text{eff}}}$ is also shown in (b). The vertical dashed line at $s/L$ = 0.5 represents the sheath entrance.

  

$$M_{\mathrm{t}}=K^{\mathrm{-1}}+\sqrt{K^{\mathrm{-2}}-1},$$(16)

where the cooling factor is given by $K=\sqrt{(T_{\mathrm{ed}}+T_{\mathrm{i}})/(T_{\mathrm{e}0}+T_{\mathrm{i}})}$ (< 1). A numerical simulation is performed by setting $G_{12}$ = 100 and $T_{\mathrm{ed}}=T_{\mathrm{e}0}/4$ (i.e., $K$ ≈ 0.79 and $M_t$ ≈ 2.0). Figure 8 shows the steady-state profiles. The flow velocities become larger than $C_{\mathrm{s},\text{eff}}$ (i.e., supersonic) in the sourceless region (0.25 ≤ $s/L$ ≤ 0.5) due to the electron temperature gradient. In addition, the Mach number is approximately unity at the X-point ($s/L=L_{\mathrm{X}}/(2L)$ = 0.25), validating the assumption of Eq. (16). Figure 9 plots the values of $V_1$ and $V_2$ at $s/L$ = 0.49 by a cross on the $V_1$-$V_2$ plane. The point ($V_1=M_{\mathrm{t}}{C}_{\mathrm{s},\text{eff}}$, $V_2=M_{\mathrm{t}}{C}_{\mathrm{s},\text{eff}}$) is also plotted as an open triangle. As shown in Fig. 9, the cross is located very near the open triangle, demonstrating that the supersonic flow expected in collisional two-ion-species plasmas can be successfully reproduced. Therefore, as demonstrated in single-ion-species plasmas [12, 13], the AIP plasma fluid scheme coupled with the VD model successfully establishes the Bohm criterion for two-ion-species plasmas, including the expected supersonic flows in collisional cases.

Fig. 9.  Steady-state ion flow velocities at $s/L$ = 0.49 near the sheath entrance (cross) with $T_{\mathrm{ed}}=T_{\mathrm{e}0}/4$ and $G_{12}$ = 100. The closed triangle represents the sound speed with the effective mass at the sheath entrance $C_{\mathrm{s},\text{eff}}=\sqrt{(T_{\mathrm{i}}+T_{\mathrm{ed}})/m_{\text{eff}}}$. The open triangle represents the target flow velocity estimated by Eq. (16). The solid lines represent the theoretical curves of the Bohm criterion (Eq. (8)). The dashed and dotted lines denote $V_2=\sqrt{m_1/m_2}{V}_1$ and $V_2=V_1$, respectively.

5.  Summary

Self-consistent modeling of the Bohm criterion in multi-ion-species plasmas, considering the effects of Coulomb collisions, is a key issue improving the accuracy of the plasma fluid model for SOL–divertor integrated codes. In this study, the combination of the AIP plasma fluid scheme and the VD model [12, 13] was applied to two-ion-species plasmas over a wide range of collisionality to investigate the reproducibility of the Bohm criterion.

The AIP plasma fluid equations were applied to isothermal (${\gamma }_{\mathrm{j}}$ = 1) and singly charged two-ion-species plasmas in a 1D SOL–divertor system with a homogeneous magnetic field. The effect of collisionality is reflected by the frictional force term in the parallel plasma momentum transport equation. The mass ratio of the two-ion species was set to 4, as simulated in an earlier PIC simulation study [8].

The theoretical curves of the Bohm criterion for two-ion-species plasmas form a hyperbola with asymptotic lines $V_{\mathrm{j}}=C_{\mathrm{s},\mathrm{j}}=\sqrt{(n_{\mathrm{j}}{T}_{\mathrm{e}}+n_{\mathrm{e}}{T}_{\mathrm{j}\mathrm{\parallel }})/(n_{\mathrm{e}}{m}_{\mathrm{j}})}$ [8]. For convenience, one curve in the upper right region ($V_{\mathrm{j}}>C_{\mathrm{s},\mathrm{j}}$) is referred to as the “upper branch,” and the other curve in the lower left region ($V_{\mathrm{j}}<C_{\mathrm{s},\mathrm{j}}$) is referred to as the “lower branch.” In the collisionless case [15], the point of ($V_1$, $V_2$) appeared very near the intersection of the “upper branch” and the line of $V_2=\sqrt{m_1/m_2}{V}_1$. As the plasma approached the collisional limit, the point of ($V_1$, $V_2$) approached the sound speed $C_{\mathrm{s},\text{eff}}=\sqrt{(T_{\mathrm{i}}+T_{\mathrm{e}})/m_{\text{eff}}}$ with the effective mass $m_{\text{eff}}=(n_1{m}_1+n_2{m}_2)/n_{\mathrm{e}}$. These tendencies reproduced the earlier PIC simulation study [8]; thus, the proposed AIP/VD model self-consistently establishes the Bohm criterion in two-ion-species plasmas without directly solving the sheath region in a wide range of collisionality.

In the collisionless case, the equilibrium solutions of the proposed system equations satisfy the same equation of the Bohm criterion (Eq. (8)). However, it is possible to introduce a viewpoint of the stability of equilibrium solutions, as discussed in Ref. [14]. If the Bohm criterion is governed by the sheath formation condition, the intersection point of the “lower branch” with the $V_2=\sqrt{m_1/m_2}{V}_1$ line would be observed, along with a steady-state solution. However, as shown in the numerical results, this “lower branch” equilibrium was never obtained, as reported in Ref. [8]. Time evolutions of the point of ($V_1$, $V_2$) show that it passes the “lower branch” equilibrium and does not stay there, indicating that the “lower branch” equilibrium is unstable. These results indicate that the Bohm criterion is not related to sheath potential formation but is a condition of stable equilibrium solutions in the quasi-neutral plasma region.

For the collisional limit case, flattening of the sourceless-region profiles was observed. Because the frictional force converges to zero within the collisional limit, the plasma can be regarded as a collisionless single-ion-species plasma. This is the flattening mechanism. In addition, using an electron temperature gradient, supersonic flows can be successfully reproduced by the proposed AIP/VD model.

In future research, we plan to perform stability analysis of the equilibrium solutions by analyzing the real part of the eigenvalues of the Jacobian of the proposed system. In addition, the energy transport equations will be incorporated into the AIP plasma fluid scheme.

 Acknowledgments

This work was conducted with the support and under the auspices of the NIFS Collaboration Research Program (NIFS23KUGM174, NIFS23KIST050). The authors would like to thank Enago (www.enago.jp) for the English language review.

References

• [1]   H.  Kawashima et al., Plasma Fusion Res. 1, 031 (2006).
• [2]   K.  Shimizu et al., Nucl. Fusion 49, 065028 (2009).
• [3]   R.  Schneider et al., Contrib. Plasma Phys. 46, 3 (2006).
• [4]   X.  Bonnin et al., Plasma Fusion Res. 11, 1403102 (2016).
• [5]   S.I.  Braginskii, Rev. Plasma Phys. 1, 205 (1965).
• [6]   D.  Bohm, The Characteristics of Electrical Discharges in Magnetic Fields, edited by A. Guthrie and R.K. Wakerling, (New York, McGraw-Hill, 1949), chapter 3, p. 77.
• [7]   K.-U.  Riemann, IEEE Trans. Plasma Sci. 23, 709 (1995).
• [8]   S.  Azuma et al., Contrib. Plasma Phys. 52, 512 (2012).
• [9]   E.M.  Hollmann et al., Rev. Sci. Instrum. 72, 623 (2001).
• [10]   H.  Yazawa et al., Jpn. J. Appl. Phys. 45, 8208 (2006).
• [11]   K.  Sugiura et al., Contrib. Plasma Phys. 64, e202300150 (2024).
• [12]   S.  Togo et al., J. Comput. Phys. 310, 109 (2016).
• [13]   S.  Togo et al., Nucl. Fusion 59, 076041 (2019).
• [14]   T.  Takizuka et al., Proceedings 13th Burning Plasma Simulation Initiative Meeting (2015) 60. https://www.riam.kyushu-u.ac.jp/sosei/bpsi/2015/BPSI2015_proc.pdf.
• [15]   S.  Togo et al., Proceedings 43rd JSST Annual International Conference on Simulation Technology & the 23rd Asia Simulation Conference (2024) 60.
• [16]   E.  Sytova et al., Contrib. Plasma Phys. 58, 622 (2018).
• [17]   S.O.  Makarov et al., Phys. Plasmas 28, 062308 (2021).