Ion temperature effects on plasma flow in the magnetic mirror configuration

Plasma flow in the magnetic mirror configuration (magnetic nozzle) plays an important role in a number of devices, such as magnetic mirrors used in fusion research and devices for electric propulsion in space. In plasma propulsion, the magnetic nozzle configuration is employed to convert the plasma thermal energy into the ion kinetic energy, thus generating thrust KaganovichPoP2020 ; LongmierPSST2011 . In open mirror fusion devices, the expanding magnetic field of the divertor (expander region) is used to spread the energy over the larger area to reduce the wall heat loads onofri2017magnetohydrodynamic . Plasma flow in the edge region of the divertor tokamak also experiences acceleration to supersonic velocities due to the combined effects of plasma pressure and inhomogeneous magnetic field BufferandPPCF2014 ; KirkPPCF2003 ; GhendrihPPCF2011 . Various aspects of plasma flow and acceleration in the magnetic mirror configurations have been studied. General framework of plasma flow and acceleration in the MHD approximation was formulated in Ref. morozov1980steady, , and subsequently used to study plasma detachment HooperAIAA1993 , acceleration mechanisms and propulsion efficiency AhedoPoP2010 ; MerinoPoP2016 . The role of the magnetic nozzle in the conversion of plasma thermal energy to supersonic flow was demonstrated experimentally InutakePST2004 ; InutakePPCF2007 . The MHD models were used to study the supersonic acceleration in the scrape-off-layer (SOL) of tokamak edge in Refs. GoswamiPoP2014, ; TogoCPP2018, ; TogoNF2019, . Ion pressure anisotropy and related mirror force effects on the plasma flow in the SOL were considered and the predictions of the isotropic and anisotropic pressure models were compared for specific conditions of the advanced divertors in fusion systems in Refs. TogoCPP2018, ; TogoNF2019, . The emphasis of this paper is on the role of the sonic point transition in the formation of the accelerating potential in plasma with anisotropic ion pressure and including atomic processes such as ionization and charge exchange.

It is well known that in the quasineutral approximation the ion inertia and finite temperature will result in the appearance of the sonic point singularity at the point where the local ion velocity is equal to the ion sound speed $c_{s}$. It has recently been emphasized SmolyakovPoP2021 that the conditions at the sonic point in the nozzle region where the magnetic field has the maximum are critical for the existence of the smooth accelerating solution such that the resulting plasma flow is uniquely defined in the whole converging-diverging region, i.e. in the whole range from sub-sonic, $V<c_{s}$, to super-sonic, $V>c_{s}$, velocities. Such smooth accelerating solution is formed by the ambipolar potential supported by the electron pressure gradient. Here we study how the addition of the ion pressure, its anisotropy in particular, the effects of the mirror force, and some dissipative processes modify the conditions for the formation of a smooth accelerating potential. These results are of interest for fusion devices application TogoNF2019 ; onofri2017magnetohydrodynamic as well as for propulsion devices where a large ion temperature is expected LongmierPSST2011 .

The basic model equations are presented in Section II. The Section III discusses general features of the sonic point singularity and ion acceleration with finite ion pressure. Plasma acceleration for cold ions is reviewed in Section IV. The Section V presents the results of the solution of the stationary and dynamic (time-dependent) equations for the base case with anisotropic pressure. The Section VI analyzes modifications due to the ionization and charge-exchange collisions. The summary and discussions are presented in Section VII.

In our model, electrons are considered in the massless isothermal $T_{e}=const$ approximation,

The ions are described by the two-pressure Chew-Golberger-Low (CGL) model chew1956boltzmann in the form

where $p_{\perp}$ and $p_{\|}$ are the perpendicular and parallel ion pressure, respectively. The standard CGL model is generalized here to include the sink terms $S_{\perp}$, $S_{\|}$ due to dissipation related to ionization and charge-exchange. We consider a problem in the paraxial approximation so that the total fluid time derivative for ions is $d/dt=\partial/\partial t+V_{\parallel}{\nabla}_{\parallel}$, where $V_{\parallel}=\mathbf{V}\cdot\mathbf{B}/B$ is the ion parallel velocity, ${\nabla}_{\parallel}=\mathbf{b}\cdot\nabla=\mathbf{B}/B\cdot\nabla$, and $\mathbf{b}=\mathbf{B}/B$ represents a unit vector in the direction of the magnetic field.

In the paraxial approximation, the ion continuity and momentum equations take the form

The effects of anisotropic pressure are included in the momentum equation with anisotropic pressure tensor ${\bm{\pi}}=(p_{\|}-p_{\perp})({\mathbf{b}}{\mathbf{b}}-{\mathbf{I}}/3)$, $p=(2p_{\perp}+p_{\|})/3$,. and where $S_{n},S_{V}$ describe the ionization source/sink effects.

In the context of the plasma flow and acceleration along the open magnetic field lines, the anisotropic MHD equations were considered in Refs. zawaideh1986generalized, ; GuoPoP2014, . Similar equations also follow from the drift-kinetic SmolyakovPoP2010 or gyro-fluid equations SnyderPoP1997 ; RobertsonPoP2016 . In simulations of plasma flow in the SOL tokamak regions, the fluid models with additional source/sink terms in the density, momentum, and energy equations as well as model approximations for the heat fluxes and relaxation terms were used TogoNF2019 . In this paper, our emphasis is on the ion pressure effects, thus we use the standard adiabatic CGL equations chew1956boltzmann modified with the model sink terms $S_{\parallel}$ and $S_{\perp}$ to model the effects of ionization and charge-exchange on the ion pressure evolution, cf. Eqs. (2) and (3).

In general, ionization contributes to the density, ion momentum and pressure evolution, while the charge exchange affects the ion momentum and pressure evolution. In the context of fusion and electric propulsion, the ionization by electron impact is most relevant. The ionization coefficients depend on the neutral species, neutral density, electron temperature and electron density, while charge-exchange coefficients depend on the ions/neutral density and their energy. Thus the ionization and charge-exchange effects may have complex profile dependencies depending on the particular situation and plasma parameters. In this study we are interested in the main parametric trends due to the ionization and charge exchange. Thus, we use a simplified model HelanderPoP1994 ; NgPoP2007 where the effects of ionization and charge exchange are represented by two constant coefficients $\nu_{1}$ and $\nu_{2}$. The $\nu_{1}$ coefficient corresponds to the ionization, while the $\nu_{2}$ coefficient in the ion momentum equation describes the total effects of ionization and charge-exchange. In this study we neglect the heat flux effects in the energy (pressure evolution) as well as any possible ionization heating terms LehnertNF1971 and only the pressure ”decay” terms are included with the same coefficient $\nu_{2}$ for $p_{\|}$ and $p_{\perp}$. In the expanded form, the full system of equations including ionization and charge-exchange effects is then written in the form similarly to Ref. onofri2017magnetohydrodynamic, :

We have to note that here we do not consider the plasma source region but focus on plasma acceleration in the magnetic mirror region, $0<z<L$. Thus, effects of ionization and charge exchange, represented by the $\nu_{1}$ and $\nu_{2}$ coefficients in the continuity, momentum, and pressure equations, are the model approximations for the ionization and charge-exchange effects in the magnetic mirror, while the plasma source is assumed for $z<0$.

Before we proceed with the general case, we consider general conditions for the existence of accelerating solution and describe plasma acceleration for cold and warm ions in the absence of ionization and charge-exchange effects.

In this section we consider a general condition for the existence of the global stationary accelerating solution in the absence of ionization and charge-exchange but taking into account anisotropic ion pressure. Setting $\nu_{1}=\nu_{2}=0$, for the stationary case one can obtain from (1 - 9) the following equation for the ion velocity

This equation has to be solved simultaneously with equations (8) and (9) for $p_{\perp}$ and $p_{\parallel}$. In equation (10), $M=V_{\parallel}/c_{s}$ is the plasma velocity $V_{\parallel}$ normalized to the speed of sound with respect to the electron velocity $c_{s}=\sqrt{(T_{e}/m_{i})}$. Such normalization is convenient because $T_{e}$ is constant and uniform. It is important to note that for warm anisotropic ions, the actual sound velocity includes the parallel ion temperature, $v_{cs}=\sqrt{(T_{e}+3T_{i\|})/m_{i}}$. This is reflected in the modification of the sonic point singularity in equation (10), which occurs at a point where the ion velocity is equal to the local ion sound velocity, $M=M_{cr}=\sqrt{1+3p_{\parallel}/nT_{e}}$. We note that the expression $v_{cs}=\sqrt{(T_{e}+T_{i\|})/m_{i}}$ used for sound speed in some papers TogoNME2019 ; TogoNF2019 is incorrect and thus misrepresents the location of the transition point.

Equation (10) describes two mechanisms of the ion acceleration: the electric field and the mirror force due to the perpendicular ion pressure. The electric field is supported by the electron pressure gradient in Eq. (1), so this mechanism comes from the electron thermal energy. The approximation of isothermal electrons used here assumes an infinite source of electron energy. In practice, electrons can experience cooling along the flow LittlePRL2016 ; MerinoPSST2020 and thus their temperature is not constant. Semi-empirically, this can be described by a general polytropic equation of state for electrons, $p_{e}\sim n^{\gamma}$. The ion acceleration due to the electron pressure with $\gamma\neq 1$ was considered previously SmolyakovPoP2021 (not included here for simplicity).

An additional contribution to the ion acceleration comes from the ion perpendicular pressure (energy), manifested by the $p_{\perp}$ term on the right hand side of (10), and is the result of the mirror force. As it will be discussed below, the effect of the ion parallel pressure modifies the location of the sonic point in Eq. 10 but the impact on the ion velocity is not significant due to the fast decrease of the ion pressure related to the decrease of plasma density along the nozzle and strong density dependence $\simeq n^{3}$ in the CGL equation for $p_{\|}$.

In general, acceleration in the magnetic mirror configuration is similar to the Laval nozzle acceleration. The ion acceleration by the electric field is supported by the density drop. In subsonic region, for $M<M_{cr}$, the ion inertia can be neglected and the acceleration occurs kinematically due to the effective area cross-section decrease as a result of the magnetic field rise, $\partial\ln B/\partial z>0$. Effectively, in this regime for $M\ll 1$, the density gradient is small and the acceleration simply follows from the flow conservation: $V^{-1}\partial V/\partial z\simeq B^{-1}\partial B/\partial z$, under constant density.

Equation (10) exhibits the singularity at $M=M_{cr}$, the point where the ion flow resonates with the ion sound mode. A global smooth accelerating solution can be obtained by regularizing the sonic point at the point $z=z_{m}$ where $\partial\ln B/\partial z=0$, which therefore requires the condition

Expanding the left and right sides of (10) near the singular point and using equation (11) one obtains the expression for the $\partial M/\partial z$ derivative near $z=z_{m}:$

This expression illustrates that the condition $\partial^{2}\ln B/\partial z^{2}<0$, i.e. maximum of the magnetic field, for $z=z_{m}$, is required for the existence of a smooth (regular) solution. Note that condition (12), for $\partial^{2}\ln B/\partial z^{2}<0$, also allows the decelerating solution with $\partial M/\partial z<0$.

Equation (12) shows that there are no regular solutions in case of the minimum of the magnetic field suggesting the possibility of shock solutions. Various dissipative effects however may affect the singularity and allow the flow with $\partial^{2}\ln B/\partial z^{2}>0$ at the sonic point. For example, ionization and charge-exchange modify the regularity condition as shown below in Section VI, Eq. (20). More generally, the singularities may be resolved by collisional MirnovNF1972 ; SkovorodinPPR2012 and finite Debye length PekkerSovJPP1984 effects allowing the flows in complex non-monotonic magnetic field geometries. The finite Debye length effects may also result in formation of the double layers at the singular points.

In this section we overview the case of cold ions providing a simple illustration of the global properties of the accelerating potential formed by the converging-diverging magnetic field configuration. For $p_{\perp}=p_{\parallel}=0$, Eq.(10) can be integrated ManheimerIEEE2001 ; FruchtmanPoP2012 resulting in the implicit equation for the ion velocity in the form

Here $B_{m}=B(z)$ at $z=z_{m}$ and the integration constant was chosen so the solution becomes regular at $z=z_{m}$ where $M=1$.

Generally, acceleration of the magnetized plasma in the magnetic nozzle is analogous to the gas flow in Laval nozzle. It is also similar to the problem of the solar wind acceleration when the effective nozzle is created by the spherical expansion and the gravity force ParkerAJ1958 . The exact solution for equation (13) can be written SmolyakovPoP2021 in terms of the Lambert function which appears in various fields including numerous plasma physics applications DubinovJPP2005 ; RaimbaultPoP2007 ; DubinovPoP2022 . The solution for spherically expanding solar wind can also be written with the Lambert function CranmerAJP2004 .

In the limit of $B(z)/B_{0}\rightarrow 0$, from Eq.(13) one can obtain two asymptotic solutions: $M\simeq B(z)/B_{0}$ and $M\simeq(-2\ln(B(z)/B_{0}))^{1/2}$. These asymptotics correspond to two branches of the regular solution, respectively, in the converging and diverging parts of the nozzle. The singularity at $M=1$, where the ion flow velocity is equal to $c_{s}$, can be removed if $\partial\ln B/\partial z=0$. This condition corresponds to the smooth matching of two branches in Eq.(13) fixing the value of the velocity derivative at the sonic point. Expanding Eq. (13) near $M=1$ gives the expression for the velocity derivative $\partial M/\partial z$, cf. with Eq. (12)

Here, for illustration, we consider the magnetic field mirror with Gaussian profile, $B\left(z\right)=(B_{m}-B_{0})\exp(-\left(z-z_{0}\right)^{2}/(\delta^{2}L^{2}))+B_{0}$, giving a mirror ratio $R=B_{m}/B(z)=4$ at both ends, $z=0$ and $z=L$. Several cases with different widths were considered by changing $\delta$, as shown in Fig. 1a.

As it is shown in Figs. 1b-1c, changing the profile of the magnetic field under the constant mirror ratio modifies the velocity, density and potential profiles but the values at $z=0$ and $z=L$ remain the same. The plasma density is normalized to the value $n_{0}$ at $z=0$: Equation (13) shows that the local value of the plasma velocity is defined by the local value of the magnetic field. It is important to note that this solution has a global character defined by the regularization condition at the sonic point, $V_{\parallel}=c_{s}$, which for cold plasma occurs at the point of the maximum magnetic field, where ${\partial\ln B}/{\partial z}=0$. The condition for the regular (smooth) solution at the sonic point defines the velocity derivative at this point, therefore fixing the velocity profile globally. Changing the values of the mirror ratio $R$ modifies the velocities at $z=0$ and $z=L$ but the value at the maximum $z=z_{m}$ remains equal to $c_{s}$, as shown in Fig. 2. Note that different values of the mirror ratio were considered by changing the value of both $\delta$ and $B_{0}$. The values of $\delta$ were adjusted as to ensure that the magnetic field with different R values would overlap at $z/L=0.5$, i.e. to have the same derivative $\partial V_{\|}/\partial z$ at $z=z_{m}$. Increase of the mirror ratio leads to the decrease of the initial velocity $M_{0}\equiv M(0)<1$, and increase of the exit velocity at $z=L$, $M_{L}\equiv M(L)>1$. The plasma velocity at the nozzle exit only depends on the mirror ratio R (regardless of the details of the magnetic field profile), and has weak logarithmic divergence for $R>>1$.

In this Section, we present the global regular solutions for plasma flow with finite (and anisotropic) ion pressure for the magnetic field profile similar to the C-2U device onofri2017magnetohydrodynamic ; BinderbauerAIP2016 shown in Fig. 3. The global solutions are considered in the region $0<z<L$. The mirror ratios were $B_{m}/B_{0}=8.0$ and $B_{m}/B_{r}=20.0$ at the left, $z=0$, and the right, $z=L$, ends of the nozzle, respectively. The magnetic field setup is described in more detail in the Appendix. The left side of the simulation region represents the transition to the plasma source at $z<0$.

The magnetic geometry used in this paper considers only one side of a typical magnetic mirror system like C-2U. This region has to be matched with the central source (plasma production and heating) region. Similar configurations are used in propulsion applications with the mirror only on one side. It is assumed in our study that the source region is at $z<0$. An important result of our study is the observation that plasma velocity is fixed globally by the condition at the singular point, so that the value of the velocity at $z=0$ cannot be arbitrary. The matching of the accelerating region (as studied here) with the source region is outside of the scope of this paper. However, some comments are provided in Section VII.

The numerical solution of stationary equations (10) and (1 - 9) are obtained by the integration from the the proximity of the sonic point, $z_{m}/L=0.5$, in both $z<z_{m}$ and $z>z_{m}$ directions. The initial conditions for the parameters $M$, $n$, $p_{\parallel}$, and $p_{\perp}$ are obtained from the Taylor series expansion near $z=z_{m}$. Across the entire nozzle electrons were assumed to be isothermal with $T_{e}=200\ eV$. At the left end of the nozzle it was assumed that $T_{i\parallel_{0}}=T_{i\perp_{0}}=200\ eV$ with the density $n_{0}=1.0\times 10^{19}\ m^{-3}$, however all results can be represented in dimensionless units by being normalized to their respective values at $z=0$.

The shooting method was used to achieve the final solution with isotropic ion pressure, $T_{i}{}_{\parallel}{}_{0}=T_{i}{}_{\perp}{}_{0}=200\ eV$ at the left end of the nozzle, which corresponds to the beginning of Region A in Fig. 3.

Similarly to the case of cold ions, plasma acceleration occurs in the regions with the finite gradient of the magnetic field. The perpendicular ion pressure enhances the ion acceleration according to equation (12). The parallel ion pressure shifts the sonic point which is now defined by Eq. (10). However, the parallel pressure decreases fast with distance, so the effects of the finite parallel pressure on the location of the singular point is weak for our parameters.

An interesting observation that can be seen from Fig. 4b is that for a finite ion temperature there is a region where the density is increasing with distance, contrary to the case of cold ions, where the density is always monotonically decreasing in regions of non-constant magnetic field. A similar behavior is observed for the potential, with the small increase with distance being displayed in Fig. 4f.
 The perpendicular plasma pressure follows the increase of the magnetic field $B$ related to the conservation of the adiabatic invariant as $T_{i\perp}(z)/B(z)=const$. The dependence of $T_{i\perp},T_{i\|}$ on $B(z)$ are shown in Fig. 4e. The decrease of $T_{i\|}$ in the outer region is related to strong plasma density dependence $~{}n^{3}$ in the CGL equation of state for $p_{\|}$ and density decrease due to plasma acceleration.

The perpendicular ion temperatures $T_{i\perp 0}=200\ eV$ were changed to study the effect of perpendicular temperature at the nozzle entrance $z=0$. Three additional cases were studies with $T_{i\perp}=400\ eV$, $T_{i\perp}=100\ eV$, and $T_{i\perp}=50\ eV$, while $T_{i\parallel 0}=200\ eV$ and $T_{e}=200\ eV$ remained the same. For the initial anisotropic ion pressure state at $z=0$, the finite perpendicular ion temperature increases the velocity of the accelerated ions as shown in Fig. 5 and displayed in Table I.

Plasma temperatures can be lower in some laboratory devices MakrinichPoP2009 , e.g. $T_{e}=10\ to\ 30\ eV$ and $T_{i}=0.5\ to\ 3\ eV$. The ion temperatures with $T_{i}\ll T_{e}$ will resemble the profiles of the cold ion case. Since, the speed of sound $c_{s}$ has been normalized as $c_{s}=\sqrt{T_{e}/m_{i}}$, the lower values of $T_{e}$ will result in a lower absolute values of the velocities. The ion temperatures profiles will qualitatively remain similar, e.g. as in Fig. 4.

Steady-state solutions of stationary equations obtained by the shooting method were compared with the time-dependent solutions obtained as a solution of the initial value problem for equations (6 - 9), presented in the Appendix with additional details. We have verified that the solution of the time-dependent equations converge well to the stationary solution obtained by the shooting method, as shown in figures 4a, 4b, 4c and 4d. The characteristic time was $t_{c}=L/c_{s}=4.0\times 10^{-5}s$ and each time interval was given as a multiple of $t_{c}$. Each time step in the simulation was equivalent to $1.0\times 10^{-5}t_{c}$. The evolution of the time-dependent problem towards the stationary solution is shown in Fig. 6 for the successive moments in time $t_{1}<t_{2}...<t_{5}$ with the lowest $t_{1}=6.0\times 10^{-2}t_{c}$ and the stationary solution value being reached at $t=2.4t_{c}$. This exercise demonstrates that the stationary accelerating solutions are stable and therefore dynamically accessible in the initial value problem.

In this section we consider the ion flow taking into account ionization and charge-exchange effects. These effects are important in fusion applications due to presence of neutrals in the mirror region. Effects of neutrals were also considered in propulsion applications MakrinichPoP2009 ; FruchtmanJPhysD2017 ; FruchtmanIEEE2011 . Depending on the application, neutral density and therefore ionization and charge-exchange coefficients may vary in a rather wide range lieberman2005principles ; onofri2017magnetohydrodynamic ; BinderbauerAIP2016 ; MakrinichPoP2009 . To estimate possible effects at various conditions, here we use a generic form of the sink terms $S_{\parallel}$ and $S_{\perp}$ with various values of $\nu_{1}$ and $\nu_{2}$ in normalized form $\nu_{2}^{\prime}=\nu_{2}L/c_{s}$, $\nu_{1}^{\prime}=\nu_{1}L/c_{s}$ from low to large values. We have used the values of $L=4\ m$ as the length of the nozzle and $c_{s}=9.78\times 10^{4}\ m/s$ as the speed of sound BinderbauerAIP2016 giving for $\nu^{\prime}_{1}\in[10^{-2}..0.5]$ and $\nu^{\prime}_{2}\in[0.1..2]$ with the dimensional values of $\nu_{1}\in[245\ s^{-1}..1.22\times 10^{4}\ s^{-1}]$ and $\nu_{2}\in[2.45\times 10^{3}\ s^{-1}..4.9\times 10^{4}\ s^{-1}]$, respectively. The ionization and charge-exchange coefficients are $\nu_{1}=\beta_{iz}N_{g}$ and $\nu_{2}=\beta_{cx}N_{g}$, where $N_{g}$ is the density of neutral atoms, and typical rates janev1987collision $\beta_{iz}=3.1\ \times 10^{-14}\ m^{3}\ s^{-1}$ and $\beta_{cx}=6.18\ \times 10^{-}{14}\ m^{3}\ s^{-1}$. For our parameters above, it gives for the values in the range $\nu_{1}\in[245\ s^{-1}..1.22\times 10^{4}\ s^{-1}]$ the neutral density $N_{g}\in[7.8\times 10^{15}\ m^{-3}..3.8\times 10^{17}\ m^{-3}]$, and for $\nu_{2}\in[2.45\times 10^{3}\ s^{-1}...4.9\times 10^{4}\ s^{-1}]$ – the values $N_{g}\in[4.0\times 10^{16}\ m^{-3}..7.9\times 10^{17}\ m^{-3}]$. Note that at lower boundaries of these values the collisional effects are not significant. Therefore, our parameters cover a wide range of typical densities observed in experiments BinderbauerAIP2016 and also extend to the low collisionalities typical for space propulsion applications. It should be noted however that generally in applications the density of neutral atoms is not constant while we assumed the constant and uniform profile. Therefore, our results should be understood as a parametric study towards revealing the main trends.

For the applications in Ref. BinderbauerAIP2016, the lowest value of the ion cyclotron frequency is $\omega_{ci}=7.3\times 10^{6}\ rad\ s^{-1}$ remaining larger than the collisional frequencies, $\omega_{ci}>>(\nu_{1},\nu_{2}),$ therefore justifying the use of the MHD approximation and assumption of magnetized ions. The problem of detachment HooperAIAA1993 ; BreizmanPoP2008 that occurs for space propulsion applications with lower values of the magnetic field is not considered in our paper.

In dimensionless form the stationary equations describing ion dynamics have the form

The plasma parameters here are normalized to their respective values at the left end of the nozzle such that $n^{\prime}=n/n_{0}$, $p_{\parallel}^{\prime}=p_{\parallel}/p_{\parallel_{0}}$, $p_{\perp}^{\prime}=p_{\perp}/p_{\perp_{0}}$, $T^{\prime}_{i_{\parallel}}=T_{i_{\parallel}}/T_{e}$, $T^{\prime}_{i_{\perp}}=T_{i_{\perp}}/T_{e}$, $z^{\prime}=z/L$ and $t^{\prime}=c_{s}t/L$. For the sake of convenience all the primes on parameters will be dropped and it will be assumed that n, $p_{\parallel}$, $p_{\perp}$, $T_{i}{}_{\parallel}$, $T_{i}{}_{\perp}$, $\nu_{1}$, $\nu_{2}$, $z$ and $t$ represent normalized quantities.

Additional terms in Eq. (18) due to the ionization and charge-exchange effects modify the regularization condition at the sonic point. The sonic point is still defined by condition 11. However the location of the sonic point, where the right hand side of equation (18) is zero, is shifted from the position of the maximum magnetic field at $z=z_{m}.$ Expanding near the point, where both sides of Eq. (18) equal zero, we obtain for $\partial M/\partial z$ the following equation

where

Here $T_{i\parallel S}$ and $T_{i\perp S}$ are the normalized values of the parallel and perpendicular temperatures at the sonic point, which no longer occurs at $z=z_{m}$ but is shifted to the right. The shift in the location of the sonic is directly affected by the value of $p_{\parallel}$ that is in turn affected by the values of $\nu_{1}$ and $\nu_{2}$. Thus, higher values of $\nu_{1}$ and $\nu_{2}$ will have a more profound effect on $p_{\parallel}$ and will indirectly shift the sonic point further away from its initial position at $z=z_{m}$. Similar to Section V, we obtain stationary solutions for equations (15) - (19) using the shooting method as described in Section III and compare these solutions with the solution of the initial value problem given by equations (1) - (9). Again, it has been confirmed that time-dependent solutions converge well to the stationary solutions after some relaxation.

There are several features introduced by ionization and charge-exchange effects. One is a non-monotonous behavior of plasma density with the increase due to the ionization, which is especially noticeable in the region of flat magnetic field before any substantial acceleration occurs. This behavior is more pronounced for higher ionization values, as indicated in Fig. 8d. Related to the density behavior, the potential also shows a non-monotonous increase in the region to the left of the maximum of the magnetic field as shown in Fig. 7e.

Another effect introduced by the dissipative terms, is the modification of the perpendicular pressure profile so that it does not follow the increase in the magnetic field. The latter effect is much reduced by dissipation so the perpendicular pressure profile may become similar to the parallel pressure with almost monotonous decrease throughout the whole region, and both the parallel and perpendicular pressure having very similar profiles with $T_{i_{\parallel}}$ and $T_{i_{\perp}}$ as displayed in Figures 7c and 7d. With the above noted modifications, the resulting velocity has a significantly larger value at the left boundary, $M=0.127$, compared with the case in the absence of charge-exchange and ionization, and a lower final acceleration value on the exit side of the nozzle. The reduction in the final value of acceleration is greater for higher charge-exchange values as shown in Fig. 7a but also in the case of ionization alone as shown in Fig. 8a. One has to note also that the presence of the dissipative terms due to the ionization and charge-exchange results in the shift of the position of the sonic point so it is no longer at the magnetic field maximum, as described by equation (20).

The effect of charge-exchange has been further studied for various values of $\nu_{2}$ as shown in Table II and Figures 7a, 7b, 7c, 7d and 7e. In Table II, the values of $M_{0}$ and $M_{L}$ represent the value of the Mach number at the $z/L=0$ and $z/L=1$ ends of the nozzle respectively. Further increase in $\nu_{2}$ results in the non-monotonous behavior of the plasma velocity in the outer region, namely Region C, where the magnetic field is almost flat, between $z/L=0.8$ and $z/L=1$, see Fig. 7a. The width of the region with the density increase on the left side of the mirror is also increasing with $\nu_{2}$ as shown in Fig. 7b. The characteristic increase in the perpendicular pressure so well pronounced for $\nu_{2}=\nu_{1}=0$, is smoothed out by finite values of $\nu_{2}$ and eventually disappears for large values of $\nu_{2}$ as shown in Fig. 7d. A similar behavior is observed for $p_{\parallel}$ and $p_{\perp}$ when only ionization is present, with the decrease being more pronounced for higher values of $\nu_{1}$ as shown in Fig. 8b and Fig. 8c respectively. An interesting observation is that the density curve corresponding to $\nu_{2}=0.35$, is below the curve corresponding to $\nu_{2}=0.1$ in the Region A where the magnetic field is constant. A larger charge-exchange value would result in a greater drag force on the ions, reducing the Mach number and resulting in a higher plasma density. It could be possible that when the rate of ionization is equal to the rate of charge-exchange, as is the case for $\nu_{1}=\nu_{2}=0.1$, the neutrals that formed in the region of constant magnetic field as a result of charge-exchange are themselves ionized, thus enhancing ionization and resulting in a greater increase in plasma density. The same effect is seen for the parallel pressure in Fig. 7c and this effect is likely a result of the behavior of the plasma density. The reduction of the velocity at the exit side, corresponding to the end of Region C, is consistent with the lower overall drop in the electrostatic potential as shown in Fig. 7e and displayed in Table II.

The flow of plasma in a magnetic nozzle was studied in the paraxial approximation using fluid MHD model with fully magnetized ions and anisotropic pressure. Stiffness of the global acceleration profile, which is fully defined by the regularization condition at the sonic point, is one of the most important results of this manuscript. The unique accelerating solution starts with a finite (and fixed) value $V_{\|}<c_{s}$ at the entrance of the magnetic mirror. This solution continues through the sonic point $V_{\|}=c_{s}$ to the diverging region where the ion velocity becomes supersonic. Plasma acceleration occurs as a result of an ambipolar electrostatic potential drop formed in the magnetic mirror, converting finite plasma pressure into the ion kinetic energy. We have shown that the finite perpendicular ion temperature increases the finite ion velocity in the exit region due to the effect of the mirror force. The finite (accelerated) value of the plasma velocity is defined by the electron and ion temperatures and the magnetic mirror aspect ratio. It is important to note that CGL theory is collisionless and the additional ion acceleration due to the perpendicular pressure is a result of the mirror force.

In many applications plasmas are collisionless but coupling of plasma flow with neutrals was noted as an important factor for some fusion onofri2017magnetohydrodynamic ; NgPoP2007 and plasma thruster applications FruchtmanJPhysD2017 ; WachsPSST2020 ; TakahashiAPL2016 . We have studied how the ionization and charge-exchange processes reform the accelerating potential impeding the plasma flow. We have shown that the ionization and charge-exchange shift the position of the sonic point from the maximum magnetic field, reduce plasma acceleration, and may result in non-monotonous behavior of plasma density and electrostatic potential.

We have obtained solutions of stationary equations with the shooting method. We have also performed simulations of time-dependent equations as an initial value problem therefore proving the stability of the obtained stationary solutions.

Similarly to Refs. MerinoPSST2018, ; AhedoPSST2020, , our model is a paraxial approximation which considers the radial variations only in the main order. Though such one-dimensional models are often a good approximation for a full two-dimensional dynamics FruchtmanPoP2012 , two-dimensional effects are expected to be important AhedoPoP2010 for large diameter systems.

In general ion kinetic effects are important for low density collisionless plasma as in Refs. BinderbauerAIP2016, ; onofri2017magnetohydrodynamic, . One of the related limitations of the current study is the use of the standard collisionless two-pressure CGL model which neglects the ion heat fluxes. Further work is required to include heat fluxes either with extended two-pressure models and/or with full kinetic theory AhedoPSST2020 ; MartinezPoP2011 . It is expected that the addition of heat fluxes will result in flattening of the parallel and perpendicular temperature profiles along the entire magnetic nozzle, thus resulting in higher Mach number at the exit point. These effects require further studies. It would be interesting to implement kinetic closures SnyderPoP1997 for the heat flux as proposed in Refs. GuoPoP2014, . The kinetic ion model would also be required to describe the possible ion trapping and demagnetization effects for very low values of the magnetic field in the expander region. We have assumed fully magnetized ions and did not consider the detachment problem which may not be critical for fusion application, but would need to be analyzed for space propulsion applications BreizmanPoP2008 .

It is expected that electron cooling and trapped electrons in the region of the diverging magnetic field (expander) are important MerinoPSST2018 ; WethertonPoP2021 ; SkovorodinPPR2019 . Here, we have focused on ion dynamics and for simplicity assumed isothermal Boltzmann electrons. Effects of general polytropic equation of state for electrons with $\gamma\neq 1$ can be included similarly as in Ref. SmolyakovPoP2021, .

The obtained global solutions do not allow for arbitrary values of plasma velocity at the entrance point of the mirror region and also fully determine the plasma velocity at the mirror exit (fixed by the mirror ratio). The uniqueness of the accelerating solutions will therefore provide the constraints on the matter and energy outflow through the magnetic nozzle (mirror). These effects were not previously considered either in the context of open mirror or propulsion applications. These results raise interesting questions of how the plasma solutions in the source region can be matched to the mirror. Such a solution will require the full analysis of the matter and energy balance, in particular, electron and energy fluxes. Neither of fluid model approximations used here, such as isothermal electrons with infinite electron heat conductivity, nor the ion CGL model which neglects the ion heat fluxes are fully adequate for this purpose. Collisionless regimes of interest for most applications also require kinetic analysis. The adjustment (matching) of plasma source with the magnetic nozzle constraints will occur via partial particle reflections and particle mixing in the source (that have to be described kinetically), possibly enhanced, in part, by non-stationary fluctuations that will occur in regimes when stationary fluid solutions do not exist. It is expected that in general the magnetic nozzle constraints discussed in this paper will reduce plasma losses through the mirror compared to the standard estimates based on the collisional transitions between trapped and passing (loss cone) regions MirnovNF1972 . The particle reflections and mixing will likely smooth out the non-monotonous features in the potential that exist in fluid models in some regimes thus increasing the plasma density in the confinement region. The kinetic analysis of such effects is outside of the scope of the present paper and is left for future studies. It is also important to note that additional drag due to interactions with neutrals (such as charge-exchange interactions) can increase the total thrust and thus be beneficial for the electric propulsion applications as shown experimentally and theoretically MakrinichPoP2009 ; FruchtmanIEEE2011 .

Data generated in this study is available from the authors upon reasonable request.