3.1 Ion trajectories in terms of slowly changing orbit parameters

The equations of motion of an ion in the magnetic presheath are

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{v}_{x}=-\frac{\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}^{\prime }}{B}+\unicode[STIX]{x1D6FA}v_{y}\cos \unicode[STIX]{x1D6FC}, & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{v}_{y}=-\unicode[STIX]{x1D6FA}v_{x}\cos \unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FA}v_{z}\sin \unicode[STIX]{x1D6FC}, & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{v}_{z}=\unicode[STIX]{x1D6FA}v_{y}\sin \unicode[STIX]{x1D6FC}. & \displaystyle\end{eqnarray}$$

Expanding (3.1)–(3.3) in $\unicode[STIX]{x1D6FC}\ll 1$ and neglecting second-order terms, we obtain

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{v}_{x}\simeq -\frac{\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}^{\prime }}{B}+\unicode[STIX]{x1D6FA}v_{y}, & \displaystyle\end{eqnarray}$$

(3.5) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{v}_{y}\simeq -\unicode[STIX]{x1D6FA}v_{x}-\unicode[STIX]{x1D6FA}v_{z}\unicode[STIX]{x1D6FC}, & \displaystyle\end{eqnarray}$$

(3.6) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{v}_{z}\simeq \unicode[STIX]{x1D6FA}v_{y}\unicode[STIX]{x1D6FC}. & \displaystyle\end{eqnarray}$$

We introduce three orbit parameters: the orbit position

(3.7) $$\begin{eqnarray}\bar{x}=x+\frac{1}{\unicode[STIX]{x1D6FA}}v_{y},\end{eqnarray}$$

the perpendicular energy

(3.8) $$\begin{eqnarray}U_{\bot }=\frac{1}{2}v_{x}^{2}+\frac{1}{2}v_{y}^{2}+\frac{\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x)}{B},\end{eqnarray}$$

and the total energy

(3.9) $$\begin{eqnarray}U=\frac{1}{2}v_{x}^{2}+\frac{1}{2}v_{y}^{2}+\frac{1}{2}v_{z}^{2}+\frac{\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x)}{B}.\end{eqnarray}$$

The orbit parameters vary over a time scale which is longer by a factor of $1/\unicode[STIX]{x1D6FC}$ than the time scale $1/\unicode[STIX]{x1D6FA}$ over which $x$ , $v_{x}$ and $v_{y}$ vary, $\bar{x}/\dot{\bar{x}}\sim U_{\bot }/\dot{U}_{\bot }\sim 1/\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FA}\gg |\boldsymbol{v}|/|\dot{\boldsymbol{v}}|\sim 1/\unicode[STIX]{x1D6FA}$ . The total energy $U$ is exactly constant, $\dot{U}=0$ . The instantaneous particle velocities can be expressed in terms of the instantaneous position $x$ and the orbit parameters,

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle v_{x}=\pm V_{x}(x,\bar{x},U_{\bot })=\pm \sqrt{2(U_{\bot }-\unicode[STIX]{x1D712}(x,\bar{x}))}, & \displaystyle\end{eqnarray}$$

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle v_{y}=\unicode[STIX]{x1D6FA}(\bar{x}-x), & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle v_{z}=V_{\Vert }(U_{\bot },U)=\sqrt{2(U-U_{\bot })}, & \displaystyle\end{eqnarray}$$

where

(3.13) $$\begin{eqnarray}\unicode[STIX]{x1D712}(x,\bar{x})=\frac{1}{2}\unicode[STIX]{x1D6FA}^{2}(x-\bar{x})^{2}+\frac{\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x)}{B}\end{eqnarray}$$

is an effective potential function. In (3.12) we assumed $v_{z}>0$ because all ions enter the magnetic presheath with $v_{z}>0$ , are accelerated to larger values of $v_{z}$ , reach the Debye sheath and are then absorbed by the wall (Geraldini et al. Reference Geraldini, Parra and Militello2018). For convenience, in (3.10) we introduced the symbol $V_{x}$ to denote the absolute value of $v_{x}$ as a function of $x$ , $\bar{x}$ and $U_{\bot }$ , and in (3.12) we introduced the symbol $V_{\Vert }$ to denote $v_{z}$ as a function of $U_{\bot }$ and $U$ .

For times comparable to the typical ion gyroperiod, $2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FA}$ , the orbit parameters are constant to lowest order in $\unicode[STIX]{x1D6FC}$ and (3.10)–(3.13) can be used to infer the approximate particle trajectory. From (3.10), the ion motion is periodic to lowest order in $\unicode[STIX]{x1D6FC}$ if, for some $\bar{x}$ and $U_{\bot }$ , turning points $x_{\text{b}}$ (bottom) and $x_{\text{t}}$ (top) exist such that: (i)  $V_{x}(x_{\text{b}},\bar{x},U_{\bot })=V_{x}(x_{\text{t}},\bar{x},U_{\bot })=0$ and (ii)  $\unicode[STIX]{x1D712}(x,\bar{x})\leqslant U_{\bot }$ in the interval $x_{\text{b}}\leqslant x\leqslant x_{\text{t}}$ . Then, the ion will move back and forth between $x_{\text{b}}$ and $x_{\text{t}}$ with period ${\sim}2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FA}$ . In order to satisfy (ii), the turning points must lie on either side of an effective potential minimum $x_{\text{m}}$ which, by definition, satisfies

(3.14) $$\begin{eqnarray}\unicode[STIX]{x1D712}^{\prime }(x_{\text{m}},\bar{x})=\unicode[STIX]{x1D6FA}^{2}(x_{\text{m}}-\bar{x})+\frac{\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}^{\prime }(x_{\text{m}})}{B}=0\end{eqnarray}$$

and

(3.15) $$\begin{eqnarray}\unicode[STIX]{x1D712}^{\prime \prime }(x_{\text{m}})=\unicode[STIX]{x1D6FA}^{2}+\frac{\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}^{\prime \prime }(x_{\text{m}})}{B}>0.\end{eqnarray}$$

The value of $\unicode[STIX]{x1D712}$ evaluated at the effective potential minimum is, using (3.13) and (3.14),

(3.16) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{\text{m}}(\bar{x})=\unicode[STIX]{x1D712}(x_{\text{m}},\bar{x})=\frac{1}{2}\left(\frac{\unicode[STIX]{x1D719}^{\prime }(x_{\text{m}})}{B}\right)^{2}+\frac{\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x_{\text{m}})}{B}.\end{eqnarray}$$

The ion motion in the $x$ direction (normal to the wall) is exactly periodic for $\unicode[STIX]{x1D6FC}=0$ , with constant orbit parameters (in the $y$ direction, the motion is a sum of an exactly periodic motion and a constant $\boldsymbol{E}\times \boldsymbol{B}$ drift, as explained in Geraldini et al. (Reference Geraldini, Parra and Militello2017)). The small angle $\unicode[STIX]{x1D6FC}$ perturbs the periodic motion by a small amount, since the orbit parameters become slowly changing in time. Thus, the ion velocity can be approximately decomposed into a periodic piece, with period ${\sim}2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D6FA}$ , and a piece that is approximately constant over the time scale of the periodic motion. Under such circumstances, there is a quantity related to the underlying periodic motion, called an adiabatic invariant, which is a constant of the overall quasi-periodic motion to lowest order in the perturbation parameter.Footnote ³ The adiabatic invariant in this system is given by

(3.17) $$\begin{eqnarray}\unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{\text{gk}}(\bar{x},U_{\bot })=\frac{1}{\unicode[STIX]{x03C0}}\int _{x_{\text{b}}}^{x_{\text{t}}}V_{x}(x,\bar{x},U_{\bot })\,\text{d}x\sim \frac{v_{\text{t},\text{i}}^{2}}{\unicode[STIX]{x1D6FA}},\end{eqnarray}$$

and is constant to lowest order in $\unicode[STIX]{x1D6FC}$ . The ordering $\unicode[STIX]{x1D707}\sim v_{\text{t},\text{i}}^{2}/\unicode[STIX]{x1D6FA}$ on the far right is obtained in the following way. We define the quantities

(3.18) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{x}=x-x_{\text{m}}\sim x_{\text{t}}-x_{\text{b}},\end{eqnarray}$$

and

(3.19) $$\begin{eqnarray}w_{x}=\dot{\unicode[STIX]{x1D70C}}_{x}\sim \sqrt{2(U_{\bot }-\unicode[STIX]{x1D712}_{\text{m}})}.\end{eqnarray}$$

Note that the size of $w_{x}$ is the characteristic orbital velocity, and the size of $\unicode[STIX]{x1D70C}_{x}$ is the characteristic spatial extent of the orbit in the $x$ direction (normal to the wall). From (3.17) we estimate $\unicode[STIX]{x1D707}\sim w_{x}\unicode[STIX]{x1D70C}_{x}$ . At the magnetic presheath entrance, where the electric field is very small, the ion gyro-orbit is circular to a good approximation. Hence, the orbital velocity is of the order of the ion thermal velocity, $w_{x}\sim v_{\text{t},\text{i}}$ , and the ion orbit size is of the order of the ion thermal gyroradius, $\unicode[STIX]{x1D70C}_{x}\sim \unicode[STIX]{x1D70C}_{\text{i}}=v_{\text{t},\text{i}}/\unicode[STIX]{x1D6FA}$ . The ordering in (3.17) follows because $\unicode[STIX]{x1D707}$ is an adiabatic invariant, and so $\unicode[STIX]{x1D707}$ is conserved to lowest order in $\unicode[STIX]{x1D6FC}$ as the ion moves across the magnetic presheath. Recall that the ion motion must retain an approximate periodicity for $\unicode[STIX]{x1D707}$ to be an adiabatic invariant.

3.2 Ion density

Treating the ion motion as periodic to lowest order in some expansion parameter is akin to conventional gyrokinetics (Rutherford & Frieman Reference Rutherford and Frieman1968; Taylor & Hastie Reference Taylor and Hastie1968; Catto Reference Catto1978; Antonsen & Lane Reference Antonsen and Lane1980; Frieman & Chen Reference Frieman and Chen1982). At every point, the ion’s trajectory can be approximated to lowest order by a periodic orbit whose period is faster than any other time scale of interest. As in gyrokinetic theory, the ion distribution function can be shown to be independent of the fast time scale to lowest order in $\unicode[STIX]{x1D6FC}$ . Moreover, since $\unicode[STIX]{x1D707}$ and $U$ are both constants of the perturbed motion (at least to lowest order in $\unicode[STIX]{x1D6FC}$ ), the distribution function written in terms of the variables $\unicode[STIX]{x1D707}$ and $U$ , $F(\unicode[STIX]{x1D707},U)$ , can be shown to be constant across the magnetic presheath (Cohen & Ryutov Reference Cohen and Ryutov1998; Geraldini et al. Reference Geraldini, Parra and Militello2017). Therefore, the function $F(\unicode[STIX]{x1D707},U)$ is completely determined by ions entering the magnetic presheath at $x\rightarrow \infty$ . In order to write $F(\unicode[STIX]{x1D707},U)$ from the distribution function at $x\rightarrow \infty$ expressed in terms of $\boldsymbol{v}$ , denoted $f_{\infty }(\boldsymbol{v})$ , we use the equations

(3.20) $$\begin{eqnarray}\unicode[STIX]{x1D707}=\frac{v_{x}^{2}+v_{y}^{2}}{2\unicode[STIX]{x1D6FA}}\end{eqnarray}$$

and

(3.21) $$\begin{eqnarray}U=\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D707}+{\textstyle \frac{1}{2}}v_{z}^{2}.\end{eqnarray}$$

These equations are obtained by setting $\unicode[STIX]{x1D719}=0$ in (3.9) and (3.17), and are thus valid at $x\rightarrow \infty$ . Note that the self-consistent form of $f_{\infty }(\boldsymbol{v})$ should be independent of the gyrophase angle, which at $x\rightarrow \infty$ is $\tan ^{-1}(v_{x}/v_{y})$ .

The ion density, $n_{\text{i}}$ , can be obtained by taking an integral in the velocity space variables $\bar{x}$ , $U_{\bot }$ and $U$ , as explained in Geraldini et al. (Reference Geraldini, Parra and Militello2017, Reference Geraldini, Parra and Militello2018). There are two distinct contributions to the ion density: one due to ions in quasiperiodic orbits

(3.22) $$\begin{eqnarray}n_{\text{i},\text{cl}}[\unicode[STIX]{x1D719}](x)=\int _{\bar{x}_{\text{m}}(x)}^{\infty }\unicode[STIX]{x1D6FA}\,\text{d}\bar{x}\int _{\unicode[STIX]{x1D712}(x,\bar{x})}^{\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})}\frac{2\,\text{d}U_{\bot }}{V_{x}(x,\bar{x},U_{\bot })}\int _{U_{\bot }}^{\infty }\frac{F(\unicode[STIX]{x1D707}_{\text{gk}}(\bar{x},U_{\bot }),U)\,\text{d}U}{V_{\Vert }(U_{\bot },U)},\end{eqnarray}$$

and another due ions that are about to intersect the wall,

(3.23) $$\begin{eqnarray}\displaystyle n_{\text{i},\text{op}}[\unicode[STIX]{x1D719}](x) & = & \displaystyle \int _{\bar{x}_{\text{m},\text{o}}(x)}^{\infty }\unicode[STIX]{x1D6FA}\,\text{d}\bar{x}\int _{\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})}^{\infty }\frac{F(\unicode[STIX]{x1D707}_{\text{gk}}(\bar{x},\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})),U)\,\text{d}U}{V_{\Vert }(\unicode[STIX]{x1D712}_{\text{M}}(\bar{x}),U)}\nonumber\\ \displaystyle & & \displaystyle \times \,[V_{x}(x,\bar{x},\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})+\unicode[STIX]{x1D6E5}_{\text{M}}(\bar{x},U))-V_{x}(x,\bar{x},\unicode[STIX]{x1D712}_{\text{M}}(\bar{x}))].\end{eqnarray}$$

The notation $f[\unicode[STIX]{x1D719}](x)$ represents a functional $f$ that depends on the whole function $\unicode[STIX]{x1D719}$ , and not just on its value at a particular position $x$ . In (3.22), the subscript ‘cl’ stands for ‘closed’ and in (3.23) the subscript ‘op’ stands for ‘open’, corresponding to ions whose trajectory can be approximated by a closed orbit (i.e. periodic) and an open orbit (i.e. terminating at the wall). The total ion density is the sum of the closed and open orbit densities of equations (3.22) and (3.23), respectively,

(3.24) $$\begin{eqnarray}n_{\text{i}}(x)=n_{\text{i},\text{cl}}[\unicode[STIX]{x1D719}](x)+n_{\text{i},\text{op}}[\unicode[STIX]{x1D719}](x).\end{eqnarray}$$

In (3.22) and (3.23), we have introduced several quantities which are derived and explained in detail in Geraldini et al. (Reference Geraldini, Parra and Militello2017, Reference Geraldini, Parra and Militello2018), and we have assumed that $\unicode[STIX]{x1D719}(x)$ , $\unicode[STIX]{x1D719}^{\prime }(x)$ and $\unicode[STIX]{x1D719}^{\prime \prime }(x)$ are all monotonic functions of $x$ . The minimum allowed orbit position $\bar{x}_{\text{m}}$ for an ion at position $x$ to be in an orbit that is periodic to lowest order in $\unicode[STIX]{x1D6FC}$ is

(3.25) $$\begin{eqnarray}\bar{x}_{\text{m}}(x)=\min _{s\in [0,x)}\frac{1}{2}(x+s)+\frac{\unicode[STIX]{x1D719}(x)-\unicode[STIX]{x1D719}(s)}{\unicode[STIX]{x1D6FA}B(x-s)}.\end{eqnarray}$$

The minimum allowed orbit position $\bar{x}_{\text{m},\text{o}}$ for an ion at position $x$ to be in an orbit that is not periodic to lowest order in $\unicode[STIX]{x1D6FC}$ is

(3.26) $$\begin{eqnarray}\bar{x}_{\text{m},\text{o}}(x)=\left\{\begin{array}{@{}ll@{}}\bar{x}_{\text{c}}\quad & \text{for }x<x_{\text{c}},\\ \bar{x}_{\text{m}}(x)\quad & \text{for }x\geqslant x_{\text{c}}.\end{array}\right.\end{eqnarray}$$

In (3.26) we have introduced the two quantities $\bar{x}_{\text{c}}$ and $x_{\text{c}}$ , defined via

(3.27) $$\begin{eqnarray}\bar{x}_{\text{c}}=\min _{x\in [0,\infty ]}\left(x+\frac{\unicode[STIX]{x1D719}^{\prime }(x)}{\unicode[STIX]{x1D6FA}B}\right)=x_{\text{c}}+\frac{\unicode[STIX]{x1D719}^{\prime }(x_{\text{c}})}{\unicode[STIX]{x1D6FA}B}.\end{eqnarray}$$

The effective potential maximum $\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})$ is the largest value of $\unicode[STIX]{x1D712}(s,\bar{x})$ for a given value of $\bar{x}$ and for values of $s$ smaller than the position of the effective potential minimum $x_{\text{m}}$ ,

(3.28) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})=\unicode[STIX]{x1D712}(x_{\text{M}},\bar{x})=\max _{s\in [0,x_{\text{m}}]}\unicode[STIX]{x1D712}(s,\bar{x}).\end{eqnarray}$$

The quantity $x_{\text{M}}$ is the position of the effective potential maximum at a given value of $\bar{x}$ . For $\bar{x}=\bar{x}_{\text{c}}$ , the values of $\unicode[STIX]{x1D712}_{\text{M}}$ and $\unicode[STIX]{x1D712}_{\text{m}}$ coincide with

(3.29) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{\text{c}}\equiv \unicode[STIX]{x1D712}(x_{\text{c}},\bar{x}_{\text{c}}).\end{eqnarray}$$

For $\bar{x}\geqslant \bar{x}_{\text{m}}(x)$ , we are guaranteed to find $\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})\geqslant \unicode[STIX]{x1D712}(x,\bar{x})$ and $x\geqslant x_{\text{M}}$ , so that there are closed orbit solutions (to lowest order in $\unicode[STIX]{x1D6FC}$ ) passing through $x$ . For $\bar{x}\geqslant \bar{x}_{\text{m},\text{o}}(x)$ , we are guaranteed to find $\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})\geqslant \unicode[STIX]{x1D712}(x,\bar{x})$ , so that there are open orbit solutions with $U_{\bot }\simeq \unicode[STIX]{x1D712}_{\text{M}}(\bar{x})$ passing through $x$ . Finally, the quantity $\unicode[STIX]{x1D6E5}_{\text{M}}$ is the range of possible values of $v_{x}^{2}/2$ that an ion in an open orbit can have at a given value of $\bar{x}$ and $U$ , and is given by

(3.30) $$\begin{eqnarray}\unicode[STIX]{x1D6E5}_{\text{M}}(\bar{x},U)=2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FC}V_{\Vert }(\unicode[STIX]{x1D712}_{\text{M}}(\bar{x}),U)\left.\frac{\text{d}\unicode[STIX]{x1D707}}{\text{d}\bar{x}}\right|_{\text{open}},\end{eqnarray}$$

where

(3.31) $$\begin{eqnarray}\left.\frac{\text{d}\unicode[STIX]{x1D707}}{\text{d}\bar{x}}\right|_{\text{open}}=\frac{\text{d}}{\text{d}\bar{x}}[\unicode[STIX]{x1D707}_{\text{gk}}(\bar{x},\unicode[STIX]{x1D712}_{\text{M}}(\bar{x}))].\end{eqnarray}$$

Equation (3.30) is derived in appendix B from the expression for $\unicode[STIX]{x1D6E5}_{\text{M}}$ given in Geraldini et al. (Reference Geraldini, Parra and Militello2018).

3.3 Quasineutrality and summary of equations

The magnetic presheath is quasineutral: the ion charge density is equal and opposite to the electron charge density, and so

(3.32) $$\begin{eqnarray}Zn_{\text{i}}(x)=n_{\text{e}}(x).\end{eqnarray}$$

Since the electrons are assumed to be in thermal equilibrium, the electron number density is

(3.33) $$\begin{eqnarray}n_{\text{e}}(x)=Zn_{\infty }\exp \left(\frac{e\unicode[STIX]{x1D719}(x)}{T_{\text{e}}}\right),\end{eqnarray}$$

where $n_{\infty }$ is the ion density for $x\rightarrow \infty$ . Using (3.24) and (3.33), the quasineutrality equation of our kinetic model can be written as

(3.34) $$\begin{eqnarray}n_{\text{i},\text{cl}}[\unicode[STIX]{x1D719}](x)+n_{\text{i},\text{op}}[\unicode[STIX]{x1D719}](x)=n_{\infty }\exp \left(\frac{e\unicode[STIX]{x1D719}(x)}{T_{\text{e}}}\right).\end{eqnarray}$$

Equation (3.34) is used to determine the self-consistent electrostatic potential $\unicode[STIX]{x1D719}(x)$ across the magnetic presheath. A condition that must be satisfied in order for (3.34) to have a solution is (Geraldini et al. Reference Geraldini, Parra and Militello2018)

(3.35) $$\begin{eqnarray}\int \frac{f_{\infty }(\boldsymbol{v})}{v_{z}^{2}}\,\text{d}^{3}v\leqslant \frac{n_{\infty }}{v_{\text{B}}^{2}},\end{eqnarray}$$

which we refer to as the kinetic Chodura condition.

Once $\unicode[STIX]{x1D719}(x)$ is calculated, we can obtain several interesting quantities. The component $u_{x}$ of the ion fluid velocity in the direction normal to the wall is obtained by use of the steady-state ion continuity equation $\text{d}/\text{d}x(n_{\text{i}}u_{x})=0$ . The quasineutrality equation (3.32) and the expression for the electron density (3.33) lead to $n_{\text{i}}=n_{\infty }\exp (e\unicode[STIX]{x1D719}/T_{\text{e}})$ . Hence, using the boundary conditions $n_{\text{i}}(\infty )=n_{\infty }$ and $u_{x}(\infty )=u_{x\infty }$ , we obtain

(3.36) $$\begin{eqnarray}u_{x}(x)=u_{x\infty }\exp \left(-\frac{e\unicode[STIX]{x1D719}(x)}{T_{\text{e}}}\right).\end{eqnarray}$$

The value of $u_{x\infty }$ is obtained from the flow velocity in the direction parallel to the magnetic field at $x\rightarrow \infty$ , projected in the direction normal to the wall. Since to lowest order in $\unicode[STIX]{x1D6FC}$ the velocity component $u_{z\infty }$ is equal to the component of the velocity parallel to the wall, we have

(3.37) $$\begin{eqnarray}u_{x\infty }=-\unicode[STIX]{x1D6FC}u_{z\infty }=-\frac{\unicode[STIX]{x1D6FC}}{n_{\infty }}\int v_{z}f_{\infty }(\boldsymbol{v})\,\text{d}^{3}v.\end{eqnarray}$$

The ion distribution function at the Debye sheath entrance, $x=0$ , is given by

(3.38) $$\begin{eqnarray}f_{0}(\boldsymbol{v})=F(\unicode[STIX]{x1D707},U)\hat{\unicode[STIX]{x1D6F1}}(v_{x},-V_{x}(0,\bar{x},\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})+\unicode[STIX]{x1D6E5}_{M}(\bar{x},U)),-V_{x}(0,\bar{x},\unicode[STIX]{x1D712}_{\text{M}}(\bar{x}))),\end{eqnarray}$$

where $\bar{x}=v_{y}/\unicode[STIX]{x1D6FA}$ at $x=0$ , and $\hat{\unicode[STIX]{x1D6F1}}$ is the top-hat function defined by

(3.39) $$\begin{eqnarray}\hat{\unicode[STIX]{x1D6F1}}\left(y,h_{1},h_{2}\right)=\left\{\begin{array}{@{}ll@{}}1\quad & \text{if }h_{1}<y\leqslant h_{2},\\ 0\quad & \text{else.}\end{array}\right.\end{eqnarray}$$

In order to study the three-dimensional distribution function at the Debye sheath entrance of (3.38), we define the distribution of the velocity component normal to the wall,

(3.40) $$\begin{eqnarray}\displaystyle f_{0x}(v_{x}) & = & \displaystyle \int _{0}^{\infty }\text{d}v_{y}\int _{0}^{\infty }f_{0}(\boldsymbol{v})\,\text{d}v_{z}=\int _{\bar{x}_{\text{c}}}^{\infty }\unicode[STIX]{x1D6FA}\,\text{d}\bar{x}\int _{\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})}^{\infty }\frac{F(\unicode[STIX]{x1D707}(\bar{x},\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})),U)\,\text{d}U}{V_{\Vert }(\unicode[STIX]{x1D712}_{\text{M}}(\bar{x}),U)}\nonumber\\ \displaystyle & & \displaystyle \times \,\hat{\unicode[STIX]{x1D6F1}}[v_{x},-V_{x}(0,\bar{x},\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})+\unicode[STIX]{x1D6E5}_{M}(\bar{x},U)),-V_{x}(0,\bar{x},\unicode[STIX]{x1D712}_{\text{M}}(\bar{x}))]\,\text{d}U,\end{eqnarray}$$

and the two-dimensional distribution of the velocity components tangential to the wall,

(3.41) $$\begin{eqnarray}\displaystyle f_{0yz}(v_{y},v_{z}) & = & \displaystyle \int _{-\infty }^{0}f_{0}(\boldsymbol{v})\,\text{d}v_{x}\nonumber\\ \displaystyle & = & \displaystyle F(\unicode[STIX]{x1D707}(\bar{x},\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})),U)[V_{x}(x,\bar{x},\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})+\unicode[STIX]{x1D6E5}_{\text{M}}(\bar{x},U))-V_{x}(x,\bar{x},\unicode[STIX]{x1D712}_{\text{M}})].\qquad\end{eqnarray}$$

Equation (3.40) is obtained by integrating (3.38) over $\bar{x}$ and $U$ , without integrating over $v_{x}$ . Equation (3.41) is obtained by integrating (3.38) over $v_{x}$ and re-expressing the distribution as a function of $v_{y}$ and $v_{z}$ , using $v_{y}=\unicode[STIX]{x1D6FA}\bar{x}$ (valid at $x=0$ ) and $v_{z}=\sqrt{2(U-\unicode[STIX]{x1D712}_{\text{M}}(\bar{x}))}$ . In Geraldini et al. (Reference Geraldini, Parra and Militello2018) it was shown that the equation

(3.42) $$\begin{eqnarray}\int \frac{f_{0}(\boldsymbol{v})}{v_{x}^{2}}\,\text{d}^{3}v\equiv \int _{-\infty }^{0}\frac{f_{0x}(v_{x})}{v_{x}^{2}}\,\text{d}v_{x}=\frac{n_{\text{i}}(0)}{v_{\text{B}}^{2}},\end{eqnarray}$$

which corresponds to the equality form of the well-known kinetic Bohm condition, is satisfied self-consistently by the magnetic presheath solution. In the review paper Riemann (Reference Riemann1991), it is shown that in most presheath models the Bohm condition is self-consistently satisfied in the equality form, as in (3.42).

4.1 Cold ions ( $\unicode[STIX]{x1D70F}\ll 1$ )

In this subsection, we argue that our kinetic model is equivalent to Chodura’s fluid model, which is valid for $\unicode[STIX]{x1D70F}=0$ , in an appropriate limit for $\unicode[STIX]{x1D70F}\ll 1$ .

In order to compare the fluid and kinetic models with each other, we first briefly recap the fluid analysis. We start by generalizing to arbitrary values of $\unicode[STIX]{x1D6FC}$ . All ions are assumed to have the same velocity at a given position $x$ , such that $\boldsymbol{v}=\boldsymbol{u}(x)$ , where $\boldsymbol{u}$ is the fluid velocity vector (a function of position only). The fluid velocity at $x\rightarrow \infty$ is chosen to be

(4.1a-c ) $$\begin{eqnarray}u_{x\infty }=-v_{\text{B}}\sin \unicode[STIX]{x1D6FC},\quad u_{y\infty }=0,\quad \text{and}\quad u_{z\infty }=v_{\text{B}}\cos \unicode[STIX]{x1D6FC}.\end{eqnarray}$$

The choice (4.1) corresponds to flow parallel to the magnetic field satisfying the Chodura (or Bohm–Chodura) condition (Chodura Reference Chodura1982) with the equality sign. Using (3.36) and (4.1a ), the ion fluid velocity at every position can be written in terms of the electrostatic potential at that position,

(4.2) $$\begin{eqnarray}u_{x}=-v_{\text{B}}\exp \left(-\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right)\sin \unicode[STIX]{x1D6FC}.\end{eqnarray}$$

As shown in § C.1, from the momentum equations and equation (4.2), one obtains a first-order differential equation for the electrostatic potential,

(4.3) $$\begin{eqnarray}\displaystyle & & \displaystyle \left(\sin ^{2}\unicode[STIX]{x1D6FC}\exp \left(-\frac{2e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right)-1\right)^{2}\frac{v_{\text{B}}^{2}}{\unicode[STIX]{x1D6FA}^{2}\cos ^{2}\unicode[STIX]{x1D6FC}}\left(\frac{e\unicode[STIX]{x1D719}^{\prime }}{T_{\text{e}}}\right)^{2}=1-\sin ^{2}\unicode[STIX]{x1D6FC}\exp \left(-\frac{2e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right)-\frac{2e\unicode[STIX]{x1D719}}{T_{\text{e}}}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\frac{1}{\cos ^{2}\unicode[STIX]{x1D6FC}}\left[2-\exp \left(\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right)-\exp \left(-\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right)\sin ^{2}\unicode[STIX]{x1D6FC}\right]^{2}.\end{eqnarray}$$

Equation (4.3) was originally derived in Chodura (Reference Chodura1982) (and later in Riemann (Reference Riemann1994)), in terms of $u_{x}$ instead of $\unicode[STIX]{x1D719}$ , and is valid for all values of $\unicode[STIX]{x1D6FC}$ (provided that $\unicode[STIX]{x1D6FC}\gg \sqrt{m_{\text{e}}/m_{\text{i}}}$ as discussed in § 2).

For $\unicode[STIX]{x1D6FC}\ll 1$ , the relationship between electrostatic potential and fluid velocity, equation (4.2), simplifies to

(4.4) $$\begin{eqnarray}u_{x}\simeq -\unicode[STIX]{x1D6FC}v_{\text{B}}\exp \left(-\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right).\end{eqnarray}$$

In the fluid model, the equality form of the Bohm condition is $u_{x}(0)=-v_{\text{B}}$ , which, from (4.4), leads to equation

(4.5) $$\begin{eqnarray}\frac{e\unicode[STIX]{x1D719}(0)}{T_{\text{e}}}\simeq \ln \unicode[STIX]{x1D6FC}.\end{eqnarray}$$

In § C.2, we expand (4.3) for $\unicode[STIX]{x1D6FC}\ll 1$ , thus obtaining an equation for the electrostatic potential for $\unicode[STIX]{x1D70F}=0$ and $\unicode[STIX]{x1D6FC}\ll 1$ ,

(4.6) $$\begin{eqnarray}x\simeq \unicode[STIX]{x1D70C}_{\text{B}}\int _{\ln \unicode[STIX]{x1D6FC}+1/2-(1/2)\unicode[STIX]{x1D6FC}^{2}}^{e\unicode[STIX]{x1D719}/T_{\text{e}}+(1/2)\unicode[STIX]{x1D6FC}^{2}[\exp (-2e\unicode[STIX]{x1D719}/T_{\text{e}})-1]}\frac{\text{d}p}{\sqrt{-3-2p+4e^{p}-e^{2p}}}.\end{eqnarray}$$

In order to obtain a correspondence between our kinetic model and Chodura’s fluid model, we define a new expansion parameter, $\unicode[STIX]{x1D716}\equiv 1/|\text{ln}\,\unicode[STIX]{x1D6FC}|$ , and take the ordering

(4.7) $$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D716}}\equiv |\text{ln}\,\unicode[STIX]{x1D6FC}|\sim |\text{ln}\,\unicode[STIX]{x1D70F}|\sim |\text{ln}\,\unicode[STIX]{x1D6FC}|-3|\text{ln}\,\unicode[STIX]{x1D70F}|\gg |\ln \unicode[STIX]{x1D716}|\sim 1.\end{eqnarray}$$

The ordering (4.7) restricts $\unicode[STIX]{x1D6FC}$ (and $\unicode[STIX]{x1D70F}$ ) to be exponentially small, $\unicode[STIX]{x1D6FC}=\exp (-1/\unicode[STIX]{x1D716})$ . However, with the ordering (4.7), the kinetic model is asymptotically equivalent to the fluid model with an exponentially small error in $\unicode[STIX]{x1D716}$ (or, equivalently, a small error in $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D70F}$ ). Hence, since the error is so small, in practice $\unicode[STIX]{x1D6FC}$ need not be excessively small (we require $\unicode[STIX]{x1D6FC}<0.1$ ).

Before analysing the kinetic model for $\unicode[STIX]{x1D70F}\ll 1$ using the ordering (4.7), it is instructive to solve (4.6) explicitly for $\unicode[STIX]{x1D716}\ll 1$ . Using the boundary condition (4.5), we order $e\unicode[STIX]{x1D719}/T_{\text{e}}\sim 1/\unicode[STIX]{x1D716}$ in the magnetic presheath. Then, we expand (4.6) in $\unicode[STIX]{x1D716}\ll 1$ to obtain

(4.8) $$\begin{eqnarray}x\simeq \unicode[STIX]{x1D70C}_{\text{B}}\int _{\ln \unicode[STIX]{x1D6FC}}^{e\unicode[STIX]{x1D719}/T_{\text{e}}}\frac{\text{d}p}{\sqrt{-2p}}.\end{eqnarray}$$

Carrying out the integral in (4.8), the electrostatic potential in the magnetic presheath is

(4.9) $$\begin{eqnarray}\frac{e\unicode[STIX]{x1D719}(x)}{T_{\text{e}}}\simeq \left\{\begin{array}{@{}ll@{}}-{\textstyle \frac{1}{2}}(x/\unicode[STIX]{x1D70C}_{\text{B}}-\sqrt{2/\unicode[STIX]{x1D716}})^{2}\quad & \text{for }x/\unicode[STIX]{x1D70C}_{\text{B}}<\sqrt{2/\unicode[STIX]{x1D716}},\\ 0\quad & \text{for }x/\unicode[STIX]{x1D70C}_{\text{B}}\geqslant \sqrt{2/\unicode[STIX]{x1D716}}.\end{array}\right.\end{eqnarray}$$

From (4.9), the length scale of the magnetic presheath is ${\sim}\unicode[STIX]{x1D70C}_{\text{B}}/\sqrt{\unicode[STIX]{x1D716}}$ .

We will find that, in the ordering (4.7), the magnetic presheath can be divided into three regions where different types of ion trajectories are dominant:

1. (i) a region far from the wall,

   (4.10) $$\begin{eqnarray}\frac{x}{\unicode[STIX]{x1D70C}_{\text{B}}}>\sqrt{\frac{2}{\unicode[STIX]{x1D716}}}-\sqrt{4|\text{ln}\,\unicode[STIX]{x1D70F}|},\end{eqnarray}$$
   where all ions are in small approximately periodic orbits (closed orbits);

2. (ii) a region close to the wall,

   (4.11) $$\begin{eqnarray}\frac{x}{\unicode[STIX]{x1D70C}_{\text{B}}}<\sqrt{\frac{2}{\unicode[STIX]{x1D716}}}-\sqrt{\frac{2}{\unicode[STIX]{x1D716}}-2|\text{ln}\,\unicode[STIX]{x1D70F}|},\end{eqnarray}$$
   where all ions are in open orbits;

3. (iii) an intermediate region,

   (4.12) $$\begin{eqnarray}\sqrt{\unicode[STIX]{x1D716}}\ll \frac{x}{\unicode[STIX]{x1D70C}_{\text{B}}}<\sqrt{\frac{2}{\unicode[STIX]{x1D716}}},\end{eqnarray}$$
   where ions moving towards the wall transition from small closed orbits to larger, distorted closed orbits, and finally to open orbits.

In §§ 4.1.1–4.1.3, we study the three regions in the order listed above. Instead of taking the limit $\unicode[STIX]{x1D70F}\ll 1$ of (3.22) and (3.23) directly, which we leave to appendix E, in §§ 4.1.1 and 4.1.2 we derive the flow velocity of ions in closed and open orbits. For ions in closed orbits, the flow velocity is much smaller than the particle velocity, as most of the particle velocity is periodic in time and gives no contribution to the flow (because the periodic motion is averaged over); hence, the flow velocity is equal to the drift velocity of the ion gyro-orbits. However, for ions in open orbits sufficiently close to the wall the motion has no periodic piece, and so the flow velocity is equal to the individual particle velocity. From the flow velocity $u_{x}$ and (4.4), we derive equations for the electrostatic potential $\unicode[STIX]{x1D719}$ in the regions (4.10) and (4.11). To lowest order in $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D70F}$ , the solution for the electrostatic potential in the part of these two regions that overlaps with the intermediate region (4.12) is a parabola. Therefore, we assume that the lowest-order solution for $\unicode[STIX]{x1D719}(x)$ in the whole intermediate region (4.12) is a parabola, and use this to write an approximate kinetic quasineutrality equation for the region (4.12). Finally, in § 4.1.4 we write an approximate differential equation for the electrostatic potential, whose solution is (4.6), and show that it is equivalent to the equations describing the electrostatic potential in the three regions.

Figure 2. Effective potential curves $\unicode[STIX]{x1D712}(x,\bar{x})$ (solid lines), corresponding to the electrostatic potential profile $\unicode[STIX]{x1D719}(x)$ (dashed line) given by the approximation (4.6) (valid for $\unicode[STIX]{x1D70F}=0$ ) with $\unicode[STIX]{x1D6FC}=0.05$ , shown for five different values of $\bar{x}$ . From the equation $\unicode[STIX]{x1D712}(\bar{x},\bar{x})=\unicode[STIX]{x1D719}(\bar{x})$ , the values of $\bar{x}$ are where the dashed line intersects the solid lines. For the different values of $\bar{x}$ , the values of $U_{\bot }$ (horizontal dotted lines) of an ion with $\unicode[STIX]{x1D707}\ll v_{\text{B}}^{2}/\unicode[STIX]{x1D6FA}$ are $U_{\bot }\simeq \unicode[STIX]{x1D712}_{\text{m}}(\bar{x})$ . When the difference between $\unicode[STIX]{x1D712}_{\text{m}}(\bar{x})$ (local minimum) and $\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})$ (local maximum) becomes so small that $U_{\bot }\simeq \unicode[STIX]{x1D712}_{\text{M}}(\bar{x})$ (shaded region around the solid vertical line, $x=x_{\text{c}}$ ), the ion gyro-orbit is distorted and enlarged.

4.1.1 Far from the wall

From (4.9), the characteristic size of the magnetic presheath is $\unicode[STIX]{x1D70C}_{\text{B}}/\sqrt{\unicode[STIX]{x1D716}}$ . Hence, sufficiently far away from the wall, all ions are in closed orbits with a radius of gyration, $\unicode[STIX]{x1D70C}_{\text{i}}$ , that is small compared with the size of the magnetic presheath, $\unicode[STIX]{x1D70C}_{\text{B}}/\sqrt{\unicode[STIX]{x1D716}}$ . The motion of the ions is thus drift kinetic. As shown in figure 2, for $\unicode[STIX]{x1D712}^{\prime \prime }(x)\neq 0$ the effective potential $\unicode[STIX]{x1D712}$ looks like a parabola locally near the minimum,

(4.13) $$\begin{eqnarray}\unicode[STIX]{x1D712}(x,\bar{x})-\unicode[STIX]{x1D712}_{\text{m}}(\bar{x})=\frac{1}{2}\unicode[STIX]{x1D712}^{\prime \prime }(x_{\text{m}})(x-x_{\text{m}})^{2}\left(1+O\left(\frac{\unicode[STIX]{x1D70C}_{x}}{l}\right)\right),\end{eqnarray}$$

where, in the error term, we have introduced the characteristic length scale over which the second derivative of the effective potential, $\unicode[STIX]{x1D712}^{\prime \prime }$ , changes,

(4.14) $$\begin{eqnarray}l=\left|\frac{\unicode[STIX]{x1D712}^{\prime \prime }(x)}{\unicode[STIX]{x1D712}^{\prime \prime \prime }(x)}\right|.\end{eqnarray}$$

Consider an ion moving in an effective potential given by (4.13). The turning points $x_{\text{b}}$ and $x_{\text{t}}$ are solutions of the equation $U_{\bot }=\unicode[STIX]{x1D712}(x,\bar{x})$ , and so

(4.15) $$\begin{eqnarray}U_{\bot }-\unicode[STIX]{x1D712}_{\text{m}}(\bar{x})=\frac{1}{8}\unicode[STIX]{x1D712}^{\prime \prime }(x_{\text{m}})(x_{\text{t}}-x_{\text{b}})^{2}\left(1+O\left(\frac{\unicode[STIX]{x1D70C}_{x}}{l}\right)\right).\end{eqnarray}$$

Recalling the definitions and orderings in (3.18)–(3.19), equation (4.15) corresponds to $w_{x}^{2}\sim \unicode[STIX]{x1D712}^{\prime \prime }(x)\unicode[STIX]{x1D70C}_{x}^{2}$ . Thus, we obtain the ordering $\unicode[STIX]{x1D70C}_{x}\sim w_{x}/\sqrt{\unicode[STIX]{x1D712}^{\prime \prime }(x)}$ relating the typical spatial extent in the $x$ direction (normal to the wall) of the ion orbit, $\unicode[STIX]{x1D70C}_{x}$ , to the typical orbital velocity component in the same direction, $w_{x}$ . Note that $w_{x}/\unicode[STIX]{x1D70C}_{x}\sim \sqrt{\unicode[STIX]{x1D712}^{\prime \prime }(x_{\text{m}})}$ is the characteristic gyrofrequency of the approximately periodic motion of the ion. This is consistent with the elliptical gyro-orbits studied in the Appendix of Geraldini et al. (Reference Geraldini, Parra and Militello2017). Moreover, from (3.17) we have the relationship $\unicode[STIX]{x1D707}\sim w_{x}\unicode[STIX]{x1D70C}_{x}\sim v_{\text{t},\text{i}}^{2}/\unicode[STIX]{x1D6FA}\sim \unicode[STIX]{x1D70F}v_{\text{B}}^{2}/\unicode[STIX]{x1D6FA}$ , from which we obtain the estimates

(4.16) $$\begin{eqnarray}w_{x}\sim \left(\frac{\unicode[STIX]{x1D712}^{\prime \prime }(x_{\text{m}})}{\unicode[STIX]{x1D6FA}^{2}}\right)^{1/4}\sqrt{\unicode[STIX]{x1D70F}}v_{\text{B}}\end{eqnarray}$$

and

(4.17) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{x}\sim \left(\frac{\unicode[STIX]{x1D712}^{\prime \prime }(x_{\text{m}})}{\unicode[STIX]{x1D6FA}^{2}}\right)^{-1/4}\sqrt{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D70C}_{\text{B}}.\end{eqnarray}$$

The electrostatic potential $\unicode[STIX]{x1D719}$ given in (4.9) has a discontinuous second derivative: for $x/\unicode[STIX]{x1D70C}_{\text{B}}\geqslant \sqrt{2/\unicode[STIX]{x1D716}}$ , we have $\unicode[STIX]{x1D719}^{\prime \prime }(x)\simeq 0$ and $\unicode[STIX]{x1D712}^{\prime \prime }(x)\simeq \unicode[STIX]{x1D6FA}^{2}$ , while for $x/\unicode[STIX]{x1D70C}_{\text{B}}<\sqrt{2/\unicode[STIX]{x1D716}}$ we have $\unicode[STIX]{x1D719}^{\prime \prime }(x)\simeq -\unicode[STIX]{x1D6FA}^{2}$ and $\unicode[STIX]{x1D712}^{\prime \prime }(x)\simeq 0$ . Hence, to lowest order in $\unicode[STIX]{x1D716}$ , the second derivative of the electrostatic potential is not determined (note that (4.3) does not specify $\unicode[STIX]{x1D719}^{\prime \prime }(x)$ ). However, the abrupt jump in the value of $\unicode[STIX]{x1D719}^{\prime \prime }(x)$ and $\unicode[STIX]{x1D712}^{\prime \prime }(x)$ occurring at $x/\unicode[STIX]{x1D70C}_{\text{B}}=\sqrt{2/\unicode[STIX]{x1D716}}$ is a reflection of a decrease of $\unicode[STIX]{x1D719}^{\prime \prime }(x)$ and $\unicode[STIX]{x1D712}^{\prime \prime }(x)$ in going from $x\rightarrow \infty$ to $x/\unicode[STIX]{x1D70C}_{\text{B}}<\sqrt{2/\unicode[STIX]{x1D716}}$ . From (4.17), the size of ion orbits is $\unicode[STIX]{x1D70C}_{x}\sim \sqrt{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D70C}_{\text{B}}\sim \unicode[STIX]{x1D70C}_{\text{i}}$ when $\unicode[STIX]{x1D712}^{\prime \prime }(x_{\text{m}})\simeq \unicode[STIX]{x1D6FA}^{2}$ . Conversely, when $\unicode[STIX]{x1D712}^{\prime \prime }(x_{\text{m}})\ll \unicode[STIX]{x1D6FA}^{2}$ , the spatial extent of the ion orbits is larger, $\unicode[STIX]{x1D70C}_{x}\gg \unicode[STIX]{x1D70C}_{\text{i}}$ . The growth of the ion orbit as it approaches the wall in the magnetic presheath is shown in figure 3. Note that as $\unicode[STIX]{x1D70C}_{x}$ becomes larger, the typical orbital velocity $w_{x}\simeq v_{x}$ becomes smaller (see (4.16)). When $\unicode[STIX]{x1D70C}_{x}$ grows so large that $\unicode[STIX]{x1D70C}_{x}\sim l$ , equations (4.13) and (4.15)–(4.17) cease to be valid as the effective potential can no longer be Taylor expanded near its minimum. This happens when the ion reaches the shaded region in figure 2.

Figure 3. An example of an ion orbit shown at two different positions: far from the wall (green), and in the intermediate region (blue). (a) The approximate trajectory is shown in the coordinates $({\tilde{y}},x)$ , where ${\tilde{y}}$ is a $y$ -coordinate in a frame of reference that is moving with the average $v_{y}$ of the ion. (b) The trajectory is shown in phase space coordinates $(v_{x},x)$ . The invariance of $\unicode[STIX]{x1D707}$ ensures that the area of the closed orbits in (b) is constant.

We proceed to solve for the ion motion by assuming that $\unicode[STIX]{x1D70C}_{x}$ is small,

(4.18) $$\begin{eqnarray}x=x_{\text{m}}+\unicode[STIX]{x1D70C}_{x}\simeq x_{\text{m}}.\end{eqnarray}$$

From (3.14) and (4.18), we obtain the value of $v_{y}=\unicode[STIX]{x1D6FA}(\bar{x}-x)$ ,

(4.19) $$\begin{eqnarray}v_{y}=\frac{\unicode[STIX]{x1D719}^{\prime }(x_{\text{m}})}{B}-\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D70C}_{x}\simeq \frac{\unicode[STIX]{x1D719}^{\prime }(x_{\text{m}})}{B}.\end{eqnarray}$$

Indeed, since the orbital velocity is small, by (4.16), and the angle between the magnetic field and the wall is shallow, the motion of the ion is approximately parallel to the wall. Thus, the magnetic force away from the wall, $ZeBv_{y}$ , is approximately equal to the electric force towards the wall, $Ze\unicode[STIX]{x1D719}^{\prime }(x)\simeq Ze\unicode[STIX]{x1D719}^{\prime }(x_{\text{m}})$ . Using (3.8) and (4.19), the perpendicular energy of an ion at a position $x\simeq x_{\text{m}}$ is given by

(4.20) $$\begin{eqnarray}U_{\bot }=\frac{1}{2}\left(\frac{\unicode[STIX]{x1D719}^{\prime }(x_{\text{m}})}{B}\right)^{2}+\frac{\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x_{\text{m}})}{B}+O(\unicode[STIX]{x1D712}^{\prime \prime }(x_{\text{m}})\unicode[STIX]{x1D70C}_{x}^{2},w_{x}^{2}).\end{eqnarray}$$

The first error in (4.20) is a combination of the orbital component of $(1/2)v_{y}^{2}$ , $\unicode[STIX]{x1D6FA}^{2}\unicode[STIX]{x1D70C}_{x}^{2}/2$ , and the quadratic term, $\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}^{\prime \prime }(x_{\text{m}})\unicode[STIX]{x1D70C}_{x}^{2}/2B$ , of the Taylor expansion of $\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x)/B$ near $x=x_{\text{m}}$ . Note that the term $O(\unicode[STIX]{x1D719}^{\prime }(x_{\text{m}})\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D70C}_{x}/B)$ arising from taking the square of (4.19) has cancelled with the linear term of the Taylor expansion of $\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x)/B$ . The second error in (4.20) comes from neglecting $(1/2)v_{x}^{2}$ . From (4.16) and (4.17), the two errors have the same size. From (4.1), the total energy of an ion at $x\rightarrow \infty$ is given by $U=v_{\text{B}}^{2}/2$ . Hence, the $z$ -component of the ion velocity is, using (3.12) and (4.20) with $U=v_{\text{B}}^{2}/2$ ,

(4.21) $$\begin{eqnarray}v_{z}=\sqrt{v_{\text{B}}^{2}-\left(\frac{\unicode[STIX]{x1D719}^{\prime }(x_{\text{m}})}{B}\right)^{2}-\frac{2\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x_{\text{m}})}{B}}+O\left(\frac{w_{x}^{2}}{v_{\text{B}}}\right).\end{eqnarray}$$

In order to obtain the ion fluid velocity $u_{x}$ , we do not require the exact velocity of an ion, $v_{x}\simeq w_{x}$ , as most of this velocity gives a quasi-periodic motion at the small length scale $\unicode[STIX]{x1D70C}_{x}\ll \unicode[STIX]{x1D70C}_{\text{B}}$ . Instead, we require a drift velocity, denoted $v_{\text{d}}$ , defined as

(4.22) $$\begin{eqnarray}v_{\text{d}}={\dot{x}}_{\text{m}}.\end{eqnarray}$$

Using this definition for the drift velocity and (3.19) and (4.18), the ion velocity ${\dot{x}}=v_{x}$ can be split into two distinct contributions,

(4.23) $$\begin{eqnarray}v_{x}=w_{x}+v_{\text{d}}.\end{eqnarray}$$

As mentioned at the end of § 3.1, it is important that the ion motion be approximately periodic for our shallow-angle kinetic model to be valid. The ion motion can be considered as an approximately periodic orbit only if the characteristic period ${\sim}\unicode[STIX]{x1D70C}_{x}/w_{x}$ is much smaller than the characteristic time it takes for the ion orbit to drift (at speed  $v_{\text{d}}$ ) by a distance $l$ such that the electrostatic potential has changed by a significant amount, ${\sim}l/v_{\text{d}}$ . Thus, the criterion for approximate periodicity is

(4.24) $$\begin{eqnarray}\frac{v_{\text{d}}}{l}\ll \frac{w_{x}}{\unicode[STIX]{x1D70C}_{x}}.\end{eqnarray}$$

We proceed to obtain an expression for $v_{\text{d}}$ . From (3.5), we obtain

(4.25) $$\begin{eqnarray}v_{x}=-\unicode[STIX]{x1D6FC}v_{z}-\frac{\dot{v}_{y}}{\unicode[STIX]{x1D6FA}}.\end{eqnarray}$$

Taking the derivative of (4.19), we obtain

(4.26) $$\begin{eqnarray}\dot{v}_{y}=v_{\text{d}}\frac{\unicode[STIX]{x1D719}^{\prime \prime }(x_{\text{m}})}{B}-\unicode[STIX]{x1D6FA}w_{x}.\end{eqnarray}$$

Inserting (4.21), (4.23) and (4.26) into (4.25), the terms proportional to $w_{x}$ on the left- and right-hand sides cancel and we obtain an implicit equation for $v_{\text{d}}$ ,

(4.27) $$\begin{eqnarray}v_{\text{d}}=-\unicode[STIX]{x1D6FC}\sqrt{v_{\text{B}}^{2}-\left(\frac{\unicode[STIX]{x1D719}^{\prime }(x_{\text{m}})}{B}\right)^{2}-\frac{2\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x_{\text{m}})}{B}}-\frac{v_{\text{d}}\unicode[STIX]{x1D719}^{\prime \prime }(x_{\text{m}})}{\unicode[STIX]{x1D6FA}B}+O\left(\unicode[STIX]{x1D6FC}\frac{w_{x}^{2}}{v_{\text{B}}}\right).\end{eqnarray}$$

The right-hand side of (4.27) consists of the small component of parallel streaming in the $x$ direction, approximately given by $-\unicode[STIX]{x1D6FC}v_{z}$ , a polarization drift, approximately given by $-v_{\text{d}}\unicode[STIX]{x1D719}^{\prime \prime }/\unicode[STIX]{x1D6FA}B$ , and the error term coming from the error in $v_{z}$ . By manipulating (4.27), we obtain

(4.28) $$\begin{eqnarray}v_{\text{d}}(x_{\text{m}})=\frac{-\unicode[STIX]{x1D6FC}\sqrt{v_{\text{B}}^{2}-(\unicode[STIX]{x1D719}^{\prime }(x_{\text{m}})/B)^{2}-2\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x_{\text{m}})/B}}{1+\unicode[STIX]{x1D719}^{\prime \prime }(x_{\text{m}})/\unicode[STIX]{x1D6FA}B}\left(1+O\left(\frac{w_{x}^{2}}{v_{\text{B}}^{2}}\right)\right).\end{eqnarray}$$

An alternative procedure to derive (4.28) is to obtain the time derivative of $x_{\text{m}}$ by using the chain rule, $v_{\text{d}}={\dot{x}}_{\text{m}}=\dot{\bar{x}}\,\text{d}x_{\text{m}}/\text{d}\bar{x}$ , as shown in appendix D. Equation (4.28) is divergent for $\unicode[STIX]{x1D719}^{\prime \prime }(x_{\text{m}})=-\unicode[STIX]{x1D6FA}B$ , but approximating the ion motion as a closed orbit becomes invalid close to the divergence, as it requires $v_{\text{d}}$ to be small by (4.24).

The ion fluid velocity $u_{x}$ is the average value of $v_{x}$ at a fixed position $x$ , not at a fixed guiding centre position $x_{\text{m}}$ . The orbital velocity $w_{x}$ averages to zero provided that the motion is approximately periodic (condition (4.24)). Moreover, writing $v_{\text{d}}(x_{\text{m}})\simeq v_{\text{d}}(x)-v_{\text{d}}^{\prime }(x)\unicode[STIX]{x1D70C}_{x}+O(v_{\text{d}}\unicode[STIX]{x1D70C}_{x}^{2}/l^{2})$ and using the fact that the linear piece in $\unicode[STIX]{x1D70C}_{x}$ averages to zero for approximately periodic motion (condition (4.24)), we obtain

(4.29) $$\begin{eqnarray}u_{x}(x)=\frac{-\unicode[STIX]{x1D6FC}\sqrt{v_{\text{B}}^{2}-(\unicode[STIX]{x1D719}^{\prime }(x)/B)^{2}+2\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x)/B}}{1+\unicode[STIX]{x1D719}^{\prime \prime }(x)/\unicode[STIX]{x1D6FA}B}\left(1+O\left(\frac{\unicode[STIX]{x1D70C}_{x}^{2}}{l^{2}}\right)\right).\end{eqnarray}$$

The $O(w_{x}^{2}/v_{\text{B}}^{2})$ error in (4.28) is neglected in (4.29) as it is smaller than the $O(\unicode[STIX]{x1D70C}_{x}^{2}/l^{2})$ error.

The assumption that $w_{x}$ averages to zero also implies that we have neglected the contribution from the open orbits, $n_{\text{i},\text{op}}(x)\simeq 0$ , so that $n_{\text{i}}(x)\simeq n_{\text{i},\text{cl}}(x)$ . The closed orbit density can then be obtained from (4.29) and the fact that $n_{\text{i}}u_{x}=-\unicode[STIX]{x1D6FC}n_{\infty }v_{\text{B}}$ ,

(4.30) $$\begin{eqnarray}n_{\text{i},\text{cl}}(x)=\frac{n_{\infty }v_{\text{B}}(1+\unicode[STIX]{x1D719}^{\prime \prime }(x)/\unicode[STIX]{x1D6FA}B)}{\sqrt{v_{\text{B}}^{2}-(\unicode[STIX]{x1D719}^{\prime }(x)/B)^{2}-2\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x)/B}}\left(1+O\left(\frac{\unicode[STIX]{x1D70F}\unicode[STIX]{x1D70C}_{\text{B}}^{2}}{l^{2}}\sqrt{\frac{\unicode[STIX]{x1D6FA}^{2}}{|\unicode[STIX]{x1D712}^{\prime \prime }(x)|}}\right)\right),\end{eqnarray}$$

where the error has been rewritten using the ordering (4.17). This result can also be derived by taking the limit $\unicode[STIX]{x1D70F}\ll 1$ in (3.22), which is a more direct though perhaps less intuitive approach (see appendix E). We can substitute either of (4.29) or (4.30) into (4.4) or (3.34), respectively, to obtain a differential equation for the electrostatic potential,

(4.31) $$\begin{eqnarray}\frac{v_{\text{B}}(1+\unicode[STIX]{x1D719}^{\prime \prime }/\unicode[STIX]{x1D6FA}B)}{\sqrt{v_{\text{B}}^{2}-(\unicode[STIX]{x1D719}^{\prime }/B)^{2}-2\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}/B}}\left(1+O\left(\frac{\unicode[STIX]{x1D70C}_{\text{B}}^{2}}{l^{2}}\left(\frac{|\unicode[STIX]{x1D712}^{\prime \prime }|}{\unicode[STIX]{x1D6FA}^{2}}\right)^{-1/2}\unicode[STIX]{x1D70F}\right)\right)=\exp \left(\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right).\end{eqnarray}$$

Recalling that $1+\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}^{\prime \prime }(x)/B=\unicode[STIX]{x1D712}^{\prime \prime }(x)/\unicode[STIX]{x1D6FA}^{2}$ , the ordering that results from (4.31) is $\unicode[STIX]{x1D712}^{\prime \prime }(x)/\unicode[STIX]{x1D6FA}^{2}\sim \exp (e\unicode[STIX]{x1D719}/T_{\text{e}})$ , which also leads to the ordering $1/l\sim \unicode[STIX]{x1D712}^{\prime \prime \prime }/\unicode[STIX]{x1D712}^{\prime \prime }\sim e\unicode[STIX]{x1D719}^{\prime }/T_{\text{e}}$ . Moreover, using the fact that balancing the terms in the denominator of (4.31) gives $\unicode[STIX]{x1D70C}_{\text{B}}^{2}(e\unicode[STIX]{x1D719}^{\prime }/T_{\text{e}})^{2}\sim e\unicode[STIX]{x1D719}/T_{\text{e}}$ , we obtain the ordering

(4.32) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D70C}_{\text{B}}^{2}}{l^{2}}\sim \frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}.\end{eqnarray}$$

Then, upon rearranging (4.31) and re-expressing the error, we obtain

(4.33) $$\begin{eqnarray}\frac{v_{\text{B}}(1+\unicode[STIX]{x1D719}^{\prime \prime }/\unicode[STIX]{x1D6FA}B)}{\sqrt{v_{\text{B}}^{2}-(\unicode[STIX]{x1D719}^{\prime }/B)^{2}-2\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}/B}}=\exp \left(\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right)+O\left(\unicode[STIX]{x1D70F}\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\exp \left(\frac{e\unicode[STIX]{x1D719}}{2T_{\text{e}}}\right)\right).\end{eqnarray}$$

Multiplying (4.33) by $e\unicode[STIX]{x1D719}^{\prime }/T_{\text{e}}$ , integrating once and using the boundary condition $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}^{\prime }=0$ at $x\rightarrow \infty$ , gives

(4.34) $$\begin{eqnarray}2-\sqrt{1-\unicode[STIX]{x1D70C}_{\text{B}}^{2}\left(\frac{e\unicode[STIX]{x1D719}^{\prime }}{T_{\text{e}}}\right)^{2}-\frac{2e\unicode[STIX]{x1D719}}{T_{\text{e}}}}=\exp \left(\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right)+O\left(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D70F}\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\exp \left(\frac{e\unicode[STIX]{x1D719}}{2T_{\text{e}}}\right)\right).\end{eqnarray}$$

Upon integrating the error on the right-hand side of (4.33), we obtain two distinct contributions to the error in (4.34): one is $O(\unicode[STIX]{x1D70F})$ and the other is $O(\unicode[STIX]{x1D70F}(e\unicode[STIX]{x1D719}/T_{\text{e}})\exp (e\unicode[STIX]{x1D719}/2T_{\text{e}}))$ . Note that both these error terms are exactly equal to zero for $\unicode[STIX]{x1D719}=0$ and they are comparable in size for $-e\unicode[STIX]{x1D719}/T_{\text{e}}\lesssim 1$ , while for $-e\unicode[STIX]{x1D719}/T_{\text{e}}\sim 1/\unicode[STIX]{x1D716}\gg 1$ the term $O(\unicode[STIX]{x1D70F})$ is larger. However, this larger term tends to a constant, which as we will see does not affect the functional form of the solution $\unicode[STIX]{x1D719}(x)$ , but only shifts the value of the constant of integration by a small amount. Equation (4.34) can be rearranged to obtain

(4.35) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{\text{B}}^{2}\left(\frac{e\unicode[STIX]{x1D719}^{\prime }}{T_{\text{e}}}\right)^{2}+3+\frac{2e\unicode[STIX]{x1D719}}{T_{\text{e}}}=4\exp \left(\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right)-\exp \left(\frac{2e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right)+O\left(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D70F}\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\exp \left(\frac{e\unicode[STIX]{x1D719}}{2T_{\text{e}}}\right)\right).\end{eqnarray}$$

Finally, equation (4.35) can be integrated to obtain the electrostatic potential far away from the wall, although a boundary condition in the intermediate region, which we have not yet specified, is required to carry out the integration.

For $-e\unicode[STIX]{x1D719}/T_{\text{e}}\gg 1$ , all the terms on the right-hand side of (4.35) become small except for the $O(\unicode[STIX]{x1D70F})$ term, which approaches a constant, and the solution approaches the parabola

(4.36) $$\begin{eqnarray}\frac{e\unicode[STIX]{x1D719}(x)}{T_{\text{e}}}\simeq \frac{e\unicode[STIX]{x1D719}_{\text{p}}(x)}{T_{\text{e}}}=-\frac{3}{2}+\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}-\frac{1}{2}\frac{(x-C)^{2}}{\unicode[STIX]{x1D70C}_{\text{B}}^{2}}.\end{eqnarray}$$

Here, $C$ is a constant determined by boundary conditions, and we denoted the constant $O(\unicode[STIX]{x1D70F})$ error coming from the right-hand side of (4.35) as $\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}$ , where $\unicode[STIX]{x1D705}$ is an unknown constant of order unity. The electrostatic potential at the wall is large, $-e\unicode[STIX]{x1D719}(0)/T_{\text{e}}=|\text{ln}\,\unicode[STIX]{x1D6FC}|=1/\unicode[STIX]{x1D716}\gg 1$ , and so we expect (4.36) to become valid closer to the wall. If we assume that (4.36) is valid at $x=0$ to lowest order in $\unicode[STIX]{x1D716}$ and impose $-e\unicode[STIX]{x1D719}(0)/T_{\text{e}}=|\text{ln}\,\unicode[STIX]{x1D6FC}|=1/\unicode[STIX]{x1D716}$ , we obtain

(4.37) $$\begin{eqnarray}C\simeq \unicode[STIX]{x1D70C}_{\text{B}}\sqrt{\frac{2}{\unicode[STIX]{x1D716}}}.\end{eqnarray}$$

To lowest order in $\unicode[STIX]{x1D716}$ , equation (4.36), with $C$ given by (4.37), is equivalent to (4.9), which was obtained from the fluid model. However, note that the non-constant piece of the error in (4.35) becomes comparable to the first term on the right-hand side when $\exp (e\unicode[STIX]{x1D719}/T_{\text{e}})\sim \unicode[STIX]{x1D70F}(e\unicode[STIX]{x1D719}/T_{\text{e}})\exp (e\unicode[STIX]{x1D719}/2T_{\text{e}})$ . Hence, equation (4.35) fails to correctly determine the potential when $\exp (e\unicode[STIX]{x1D719}/T_{\text{e}})\sim \unicode[STIX]{x1D70F}^{2}|\text{ln}\,\unicode[STIX]{x1D70F}|^{2}\sim \unicode[STIX]{x1D70F}^{2}/\unicode[STIX]{x1D716}^{2}$ . From (4.17), with the ordering $\unicode[STIX]{x1D712}^{\prime \prime }/\unicode[STIX]{x1D6FA}^{2}\sim \exp (e\unicode[STIX]{x1D719}/T_{\text{e}})$ , and (4.32), this value of $\exp (e\unicode[STIX]{x1D719}/T_{\text{e}})$ corresponds to $\unicode[STIX]{x1D70C}_{x}\sim l\sim \sqrt{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D70C}_{\text{B}}$ , which is the point at which the approximation in (4.13) ceases to be valid. The validity of (4.35) is thus restricted to

(4.38) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D70F}^{2}}{\unicode[STIX]{x1D716}^{2}}\ll \exp \left(\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right).\end{eqnarray}$$

Note that, from (4.36), (4.37) and (4.38), the validity region is given by (4.10) to lowest order in $\unicode[STIX]{x1D716}$ . Ion gyro-orbits grow in size as they approach the wall, as shown in figure 3, making the treatment of this section invalid for $x/\unicode[STIX]{x1D70C}_{\text{B}}\leqslant \sqrt{2/\unicode[STIX]{x1D716}}-\sqrt{4|\text{ln}\,\unicode[STIX]{x1D70F}|}$ , where $\unicode[STIX]{x1D70C}_{x}$ is no longer small.

4.1.2 Near the wall

When $U_{\bot }\simeq \unicode[STIX]{x1D712}_{\text{M}}$ , ions transition to open orbits and thereafter reach the wall in a time scale of the order of a gyroperiod, $\unicode[STIX]{x1D70C}_{x}/w_{x}\sim 1/\sqrt{\unicode[STIX]{x1D712}^{\prime \prime }(x_{\text{m}})}$ . Previously, we saw that the spatial extent of a closed ion orbit enlarges as the ion approaches the wall; for the moment, we take $\unicode[STIX]{x1D70C}_{x}\sim \unicode[STIX]{x1D70C}_{\text{B}}$ for ions transitioning from closed to open orbits, ignoring any potential scaling with $\unicode[STIX]{x1D716}$ . For such transitioning ions, we expect that $\unicode[STIX]{x1D712}_{\text{M}}-\unicode[STIX]{x1D712}_{\text{m}}\sim w_{x}^{2}\sim \unicode[STIX]{x1D70F}^{2}v_{\text{B}}^{2}/\unicode[STIX]{x1D70C}_{x}^{2}$ , since $\unicode[STIX]{x1D707}\sim \unicode[STIX]{x1D70C}_{x}w_{x}\sim \unicode[STIX]{x1D70F}v_{\text{B}}^{2}$ . Thus, $\unicode[STIX]{x1D712}_{\text{M}}-\unicode[STIX]{x1D712}_{\text{m}}$ is small in $\unicode[STIX]{x1D70F}$ . Recall, from (3.29), that $\unicode[STIX]{x1D712}_{\text{c}}$ is defined to be the value of the effective potential at $\bar{x}=\bar{x}_{\text{c}}$ such that $\unicode[STIX]{x1D712}_{\text{c}}=\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})=\unicode[STIX]{x1D712}_{\text{m}}(\bar{x})$ . Hence, it follows that $\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})\simeq \unicode[STIX]{x1D712}_{\text{c}}$ and $\bar{x}\simeq \bar{x}_{\text{c}}$ for all ions in open orbits, as can be seen in figure 2. The error in approximating $\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})-\unicode[STIX]{x1D712}(x,\bar{x})\simeq \unicode[STIX]{x1D712}_{\text{c}}-\unicode[STIX]{x1D712}(x,\bar{x}_{\text{c}})$ can be obtained by calculating

(4.39) $$\begin{eqnarray}\left.\frac{\text{d}}{\text{d}\unicode[STIX]{x1D707}}\right|_{\text{open}}(\unicode[STIX]{x1D712}_{\text{M}}-\unicode[STIX]{x1D712}(x,\bar{x}))=\left(\frac{\text{d}\unicode[STIX]{x1D707}}{\text{d}\bar{x}}\right)_{\text{open}}^{-1}\unicode[STIX]{x1D6FA}^{2}(x-x_{\text{M}})\sim \frac{\unicode[STIX]{x1D70F}\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D70C}_{\text{B}}^{2}}{\unicode[STIX]{x1D70C}_{x}^{2}},\end{eqnarray}$$

where we used $\text{d}(\unicode[STIX]{x1D712}_{\text{M}}-\unicode[STIX]{x1D712}(x,\bar{x}))/\text{d}\bar{x}=\unicode[STIX]{x1D6FA}^{2}(x-x_{\text{M}})\sim \unicode[STIX]{x1D6FA}^{2}\unicode[STIX]{x1D70C}_{x}$ and estimated $(\text{d}\unicode[STIX]{x1D707}/\text{d}\bar{x})_{\text{open}}\sim \unicode[STIX]{x1D6FA}^{2}\unicode[STIX]{x1D70C}_{x}^{2}/w_{x}\sim \unicode[STIX]{x1D6FA}\unicode[STIX]{x1D70C}_{x}^{3}/\unicode[STIX]{x1D70F}\unicode[STIX]{x1D70C}_{\text{B}}^{2}$ from (B 7). Since typical ion orbits have values of $\unicode[STIX]{x1D707}$ differing by $O(\unicode[STIX]{x1D70F}v_{\text{B}}^{2}/\unicode[STIX]{x1D6FA})$ , the values of $\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})-\unicode[STIX]{x1D712}(x,\bar{x})$ of such orbits change by $O(\unicode[STIX]{x1D70F}^{2}v_{\text{B}}^{2}\unicode[STIX]{x1D70C}_{\text{B}}^{2}/\unicode[STIX]{x1D70C}_{x}^{2})$ . Recall that an open orbit has $U_{\bot }-\unicode[STIX]{x1D712}(x,\bar{x})=\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})-\unicode[STIX]{x1D712}(x,\bar{x})+O(\unicode[STIX]{x1D6E5}_{\text{M}})$ : from (3.30) and the previous estimate for $(\text{d}\unicode[STIX]{x1D707}/\text{d}\bar{x})_{\text{open}}$ , we obtain the scaling $\unicode[STIX]{x1D6E5}_{\text{M}}\sim \unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FA}v_{\text{B}}\unicode[STIX]{x1D70C}_{x}^{3}/\unicode[STIX]{x1D70F}\unicode[STIX]{x1D70C}_{\text{B}}^{2}$ . If the condition

(4.40) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D70F}^{3}}\sim \unicode[STIX]{x1D6FC}^{A}=\exp \left(-\frac{A}{\unicode[STIX]{x1D716}}\right)\ll 1\end{eqnarray}$$

is satisfied, with $A$ a positive constant, the $O(\unicode[STIX]{x1D6E5}_{\text{M}})$ error term is small in $\unicode[STIX]{x1D6FC}$ compared to the $O(\unicode[STIX]{x1D70F}^{2}v_{\text{B}}^{2}\unicode[STIX]{x1D70C}_{\text{B}}^{2}/\unicode[STIX]{x1D70C}_{x}^{2})$ error. Then, using (3.10) with $U_{\bot }-\unicode[STIX]{x1D712}(x,\bar{x})\simeq \unicode[STIX]{x1D712}_{\text{c}}-\unicode[STIX]{x1D712}(x,\bar{x}_{\text{c}})+O(\unicode[STIX]{x1D70F}^{2}v_{\text{B}}^{2}\unicode[STIX]{x1D70C}_{\text{B}}^{2}/\unicode[STIX]{x1D70C}_{x}^{2})$ , the velocity of an ion in an open orbit near the wall is

(4.41) $$\begin{eqnarray}v_{x}=-\sqrt{2(\unicode[STIX]{x1D712}_{\text{c}}-\unicode[STIX]{x1D712}(x,\bar{x}_{\text{c}}))+O(\unicode[STIX]{x1D70F}^{2}v_{\text{B}}^{2}\unicode[STIX]{x1D70C}_{\text{B}}^{2}/\unicode[STIX]{x1D70C}_{x}^{2})}.\end{eqnarray}$$

As we will see at the end of this section, condition (4.40) must be satisfied for our kinetic model to be valid; otherwise, the velocity of all ions transitioning from closed to open orbits is not known to lowest order in $\unicode[STIX]{x1D6FC}$ , making the ion density incorrect in a large region. The ordering (4.7) includes the validity condition (4.40).

Assuming that $x$ is sufficiently close to the wall that most ions are in open orbits, $n_{\text{i}}(x)\simeq n_{\text{i},\text{op}}(x)$ , the ion fluid velocity is

(4.42) $$\begin{eqnarray}u_{x}(x)=-\sqrt{2(\unicode[STIX]{x1D712}_{\text{c}}-\unicode[STIX]{x1D712}(x,\bar{x}_{\text{c}}))+O(\unicode[STIX]{x1D70F}^{2}v_{\text{B}}^{2}\unicode[STIX]{x1D70C}_{\text{B}}^{2}/\unicode[STIX]{x1D70C}_{x}^{2})}.\end{eqnarray}$$

Then, from (4.42) and the continuity equation $n_{\text{i},\text{op}}(x)u_{x}(x)=-\unicode[STIX]{x1D6FC}n_{\infty }v_{\text{B}}$ , we obtain an expression for the open orbit density,

(4.43) $$\begin{eqnarray}n_{\text{i},\text{op}}(x)=\frac{\unicode[STIX]{x1D6FC}n_{\infty }v_{\text{B}}}{\sqrt{2(\unicode[STIX]{x1D712}_{\text{c}}-\frac{1}{2}\unicode[STIX]{x1D6FA}^{2}(x-\bar{x}_{\text{c}})^{2}-\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x)/B)+O(\unicode[STIX]{x1D70F}^{2}v_{\text{B}}^{2}\unicode[STIX]{x1D70C}_{\text{B}}^{2}/\unicode[STIX]{x1D70C}_{x}^{2})}}.\end{eqnarray}$$

In § E.2, we derive (4.43) by expanding $n_{\text{i},\text{op}}(x)$ in (3.23) to lowest order in $\unicode[STIX]{x1D70F}\ll 1$ . In order for (4.43) to be valid, we require the error term to be small, implying $n_{\text{i},\text{op}}\ll (\unicode[STIX]{x1D6FC}n_{\infty }/\unicode[STIX]{x1D70F})(\unicode[STIX]{x1D70C}_{x}/\unicode[STIX]{x1D70C}_{\text{B}})$ . Moreover, since $\unicode[STIX]{x1D70C}_{x}$ here quantifies the characteristic size of closed orbits while the ion is transitioning from a closed to an open orbit, the ion density changes from $n_{\text{i},\text{cl}}\gg n_{\infty }\unicode[STIX]{x1D70F}^{2}/\unicode[STIX]{x1D716}^{2}$ (recall (4.38) and the fact that $n_{\text{i},\text{cl}}\sim n_{\infty }\exp (e\unicode[STIX]{x1D719}/T_{\text{e}})$ ) to $n_{\text{i},\text{op}}\ll (\unicode[STIX]{x1D6FC}n_{\infty }/\unicode[STIX]{x1D70F})(\unicode[STIX]{x1D70C}_{x}/\unicode[STIX]{x1D70C}_{\text{B}})$ over a length scale of $\unicode[STIX]{x1D70C}_{x}$ . If the validity condition (4.40) is satisfied, this drop in density corresponds to a decrease of $\ln (\unicode[STIX]{x1D70F}^{3}/\unicode[STIX]{x1D6FC})\sim 1/\unicode[STIX]{x1D716}$ in the normalized electrostatic potential $e\unicode[STIX]{x1D719}/T_{\text{e}}$ . In order to be consistent with the lowest-order electric field obtained from (4.36) and (4.37), $e\unicode[STIX]{x1D719}^{\prime }/T_{\text{e}}=-(x-C)/\unicode[STIX]{x1D70C}_{\text{B}}^{2}\sim 1/\sqrt{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D70C}_{\text{B}}$ , transitioning ions must have $\unicode[STIX]{x1D70C}_{x}\sim \unicode[STIX]{x1D70C}_{\text{B}}/\sqrt{\unicode[STIX]{x1D716}}$ .

Inserting (4.43) into (3.34) with $n_{\text{i},\text{cl}}(x)=0$ , or inserting (4.42) into (4.4), we obtain

(4.44) $$\begin{eqnarray}\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}+\frac{1}{2}\unicode[STIX]{x1D6FC}^{2}\exp \left(-\frac{2e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right)=K-\frac{3}{2}-\frac{(x-\bar{x}_{\text{c}})^{2}}{2\unicode[STIX]{x1D70C}_{\text{B}}^{2}}+O(\unicode[STIX]{x1D716}\unicode[STIX]{x1D70F}^{2}).\end{eqnarray}$$

The constants

(4.45) $$\begin{eqnarray}K=\frac{\unicode[STIX]{x1D712}_{\text{c}}}{v_{\text{B}}^{2}}+\frac{3}{2}\end{eqnarray}$$

and $\bar{x}_{\text{c}}$ are to be determined; they are related by the boundary condition $e\unicode[STIX]{x1D719}(0)/T_{\text{e}}=\ln \unicode[STIX]{x1D6FC}$ , giving

(4.46) $$\begin{eqnarray}\bar{x}_{\text{c}}\simeq \unicode[STIX]{x1D70C}_{\text{B}}\sqrt{-2\ln \unicode[STIX]{x1D6FC}-4+2K}\sim \frac{\unicode[STIX]{x1D70C}_{\text{B}}}{\sqrt{\unicode[STIX]{x1D716}}}.\end{eqnarray}$$

Note that far from the wall, equation (4.44) gives

(4.47) $$\begin{eqnarray}\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\simeq K-\frac{3}{2}-\frac{(x-\bar{x}_{\text{c}})^{2}}{2\unicode[STIX]{x1D70C}_{\text{B}}^{2}}.\end{eqnarray}$$

For (4.47) to be valid, we require $\unicode[STIX]{x1D6FC}^{2}\exp (-2e\unicode[STIX]{x1D719}/T_{\text{e}})\sim \unicode[STIX]{x1D6FC}^{2}\exp ((x-\bar{x}_{\text{c}})^{2}/\unicode[STIX]{x1D70C}_{\text{B}}^{2})\sim \exp (-x\bar{x}_{\text{c}}/\unicode[STIX]{x1D70C}_{\text{B}}^{2}+x^{2}/2\unicode[STIX]{x1D70C}_{\text{B}}^{2})\ll 1$ , where we have used $\bar{x}_{\text{c}}\simeq \unicode[STIX]{x1D70C}_{\text{B}}\sqrt{2/\unicode[STIX]{x1D716}}$ ; hence, the electrostatic potential becomes well approximated by (4.47) for $x\gg \unicode[STIX]{x1D70C}_{\text{B}}\sqrt{\unicode[STIX]{x1D716}}$ . The derivation of (4.44) fails when $n_{\text{i},\text{op}}\sim \unicode[STIX]{x1D6FC}n_{\infty }/\sqrt{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D70F}$ (recall (4.43) with $\unicode[STIX]{x1D70C}_{x}\sim \unicode[STIX]{x1D70C}_{\text{B}}/\sqrt{\unicode[STIX]{x1D716}}$ ), and so the validity of (4.44) is restricted to

(4.48) $$\begin{eqnarray}\exp \left(\frac{e\unicode[STIX]{x1D719}}{T_{\text{e}}}\right)\ll \frac{\unicode[STIX]{x1D6FC}}{\sqrt{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D70F}}.\end{eqnarray}$$

From (4.47) and (4.48), we obtain the estimate (4.11) for the region where (4.44) is valid. Outside of the validity region (4.11), the velocity of a typical ion is of the order of the gyration velocity, $\sqrt{2(\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})-\unicode[STIX]{x1D712}(x,\bar{x}))}\sim w_{x}$ , and so the assumption that all ions are in open orbits is invalid.

4.1.3 Intermediate region

With the ordering (4.40), there is a finite region where (4.35) and (4.44) are not valid: from (4.10) and (4.11), this region is

(4.49) $$\begin{eqnarray}\sqrt{\frac{2}{\unicode[STIX]{x1D716}}}-\sqrt{\frac{2}{\unicode[STIX]{x1D716}}-2|\text{ln}\,\unicode[STIX]{x1D70F}|}\leqslant \frac{x}{\unicode[STIX]{x1D70C}_{\text{B}}}\leqslant \sqrt{\frac{2}{\unicode[STIX]{x1D716}}}-\sqrt{4|\text{ln}\,\unicode[STIX]{x1D70F}|}.\end{eqnarray}$$

However, the solution of (4.35) tends to (4.36) for $\sqrt{2/\unicode[STIX]{x1D716}}-x/\unicode[STIX]{x1D70C}_{\text{B}}\gg 1$ and (4.44) tends to (4.47) for $x/\unicode[STIX]{x1D70C}_{\text{B}}\gg \sqrt{\unicode[STIX]{x1D716}}$ . Hence, we proceed by assuming that in the intermediate region (4.12), which includes the region (4.49), the electrostatic potential is simultaneously given by the parabolas in (4.36) and (4.47) to lowest order in $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D70F}$ . This provides the value of $K$ , $K=\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}\ll 1$ , and the equality $C=\bar{x}_{\text{c}}$ . Using (4.46) with $K\simeq 0$ , the value of $C$ and $\bar{x}_{\text{c}}$ is

(4.50) $$\begin{eqnarray}C=\bar{x}_{\text{c}}\simeq \unicode[STIX]{x1D70C}_{\text{B}}\sqrt{-2\ln \unicode[STIX]{x1D6FC}-4}.\end{eqnarray}$$

The neglected term $\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}$ causes a small constant correction to the value of $C$ , as we had claimed in the discussion following (4.34).

From (3.13), the effective potential curves associated with the parabolic electrostatic potential of (4.36) are a set of straight lines,

(4.51) $$\begin{eqnarray}\unicode[STIX]{x1D712}(x;\bar{x})\simeq -{\textstyle \frac{3}{2}}v_{\text{B}}^{2}+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D6FA}^{2}(\bar{x}^{2}-C^{2})-\unicode[STIX]{x1D6FA}^{2}(\bar{x}-C)x.\end{eqnarray}$$

In figure 2, a family of effective potential curves $\unicode[STIX]{x1D712}(x;\bar{x})$ are plotted for different values of the orbit position $\bar{x}$ for $\unicode[STIX]{x1D6FC}=0.05$ : the curves shown are indeed close to straight lines in the shaded region, as (4.51) suggests. Since straight lines do not have a local minimum – which is necessary to approximate the ion motion as a periodic orbit – the small non-parabolic piece of the electrostatic potential,

(4.52) $$\begin{eqnarray}\unicode[STIX]{x1D719}_{\text{np}}(x)=\unicode[STIX]{x1D719}(x)-\unicode[STIX]{x1D719}_{\text{p}}(x),\end{eqnarray}$$

must be retained. In (4.52), $\unicode[STIX]{x1D719}(x)$ is the solution to the quasineutrality equation (3.34) for a given value of $\unicode[STIX]{x1D70F}$ and $\unicode[STIX]{x1D6FC}$ . In the intermediate region we take $\unicode[STIX]{x1D719}(x)\simeq \unicode[STIX]{x1D719}_{\text{p}}(x)$ and calculate $\unicode[STIX]{x1D719}_{\text{np}}(x)$ as a higher-order asymptotic correction from the following equation:

(4.53) $$\begin{eqnarray}n_{\text{i},\text{cl}}[\unicode[STIX]{x1D719}_{\text{p}}+\unicode[STIX]{x1D719}_{\text{np}}](x)+n_{\text{i},\text{op}}[\unicode[STIX]{x1D719}_{\text{p}}+\unicode[STIX]{x1D719}_{\text{np}}](x)=n_{\infty }\exp \left(\frac{e\unicode[STIX]{x1D719}_{\text{p}}(x)}{T_{\text{e}}}\right)\left(1+O\left(\frac{e\unicode[STIX]{x1D719}_{\text{np}}}{T_{\text{e}}}\right)\right).\end{eqnarray}$$

On the right-hand side of (4.53), we neglected terms small in $e\unicode[STIX]{x1D719}_{\text{np}}/T_{\text{e}}\ll 1$ to simplify the expression for the electron density. On the left-hand side, we included the non-parabolic piece $\unicode[STIX]{x1D719}_{\text{np}}$ because no effective potential minima exist – and so no closed or open ion orbits can be solved for – when $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{\text{p}}$ .

At the beginning of this section, we noted that (4.35) and (4.44) do not have a common region of validity. Therefore, it is crucial that (4.35) and (4.53) be simultaneously valid in some overlap region of finite size; the same has to be true for (4.44) and (4.53). Equation (4.35) is valid in region (4.10), and equation (4.53) is valid in the region (4.12). Hence, the overlap region in which both equations are valid is

(4.54) $$\begin{eqnarray}1\ll \sqrt{\frac{2}{\unicode[STIX]{x1D716}}}-\frac{x}{\unicode[STIX]{x1D70C}_{\text{B}}}<\sqrt{4|\text{ln}\,\unicode[STIX]{x1D70F}|},\end{eqnarray}$$

where we re-expressed the lowest-order inequality $x/\unicode[STIX]{x1D70C}_{\text{B}}<\sqrt{2/\unicode[STIX]{x1D716}}$ to the more precise form $1\ll \sqrt{2/\unicode[STIX]{x1D716}}-x/\unicode[STIX]{x1D70C}_{\text{B}}$ in order to emphasize the necessity of the ordering $|\text{ln}\,\unicode[STIX]{x1D70F}|\sim 1/\unicode[STIX]{x1D716}\gg 1$ . We proceed to calculate $\unicode[STIX]{x1D719}_{\text{np}}(x)$ in this region. Inserting $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{\text{p}}+\unicode[STIX]{x1D719}_{\text{np}}$ in (4.31) and rearranging the error term, we obtain

(4.55) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D719}_{\text{np}}^{\prime \prime }(x)}{2\unicode[STIX]{x1D6FA}B}=\exp \left(-\frac{3}{2}-\frac{(x-C)^{2}}{2\unicode[STIX]{x1D70C}_{\text{B}}^{2}}\right)+O\left(\frac{\unicode[STIX]{x1D70F}\unicode[STIX]{x1D70C}_{\text{B}}}{l^{2}}\sqrt{\frac{e\unicode[STIX]{x1D719}_{\text{np}}^{\prime \prime }}{T_{\text{e}}}}\right),\end{eqnarray}$$

where we have used $\unicode[STIX]{x1D719}\simeq \unicode[STIX]{x1D719}_{\text{p}}$ in the denominator to get $\sqrt{v_{\text{B}}^{2}-(\unicode[STIX]{x1D719}^{\prime }(x)/B)^{2}-2\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}(x)/B}\simeq 2v_{\text{B}}$ . From the definition of $l$ in (4.14) and using $\unicode[STIX]{x1D712}^{\prime \prime }=\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D719}_{\text{np}}^{\prime \prime }/B$ with equation (4.55), we obtain $l\sim \unicode[STIX]{x1D70C}_{\text{B}}^{2}/(C-x)\sim \unicode[STIX]{x1D70C}_{\text{B}}\sqrt{\unicode[STIX]{x1D716}}$ , consistent with our previous estimate for $l$ in this region (before equation (4.38)). Integrating (4.55) twice and imposing $\unicode[STIX]{x1D719}_{\text{np}}^{\prime }(x)=\unicode[STIX]{x1D719}_{\text{np}}(x)=0$ at $(C-x)/\unicode[STIX]{x1D70C}_{\text{B}}\rightarrow \infty$ (where the electrostatic potential becomes more parabolic) gives

(4.56) $$\begin{eqnarray}\displaystyle \frac{e\unicode[STIX]{x1D719}_{\text{np}}(x)}{T_{\text{e}}} & = & \displaystyle 2\exp \left(-\frac{3}{2}-\frac{(x-C)^{2}}{2\unicode[STIX]{x1D70C}_{\text{B}}^{2}}\right)-\sqrt{2\unicode[STIX]{x03C0}}(C-x)\exp \left(-\frac{3}{2}\right)\left(1-\text{erf}\left(\frac{C-x}{\sqrt{2}\unicode[STIX]{x1D70C}_{\text{B}}}\right)\right)\nonumber\\ \displaystyle & & \displaystyle +\,O\left(\frac{\unicode[STIX]{x1D70F}}{\sqrt{\unicode[STIX]{x1D716}}}\sqrt{\frac{e\unicode[STIX]{x1D719}_{\text{np}}}{T_{\text{e}}}}\right),\end{eqnarray}$$

where we have used that the double integral of the term $O((\unicode[STIX]{x1D70F}\unicode[STIX]{x1D70C}_{\text{B}}/l^{2})\sqrt{e\unicode[STIX]{x1D719}_{\text{np}}^{\prime \prime }/T_{\text{e}}})$ is $O((\unicode[STIX]{x1D70F}\unicode[STIX]{x1D70C}_{\text{B}}/l)\sqrt{e\unicode[STIX]{x1D719}_{\text{np}}/T_{\text{e}}})$ . We proceed to consider the part of the intermediate region that is closest to the wall. Equation (4.44) is valid in the region (4.11) near the wall and equation (4.53) is valid in the intermediate region (4.12). Hence, the region in which these two equations are both valid is the overlap of (4.11) and (4.12),

(4.57) $$\begin{eqnarray}\sqrt{\unicode[STIX]{x1D716}}\ll \frac{x}{\unicode[STIX]{x1D70C}_{\text{B}}}<\sqrt{\frac{2}{\unicode[STIX]{x1D716}}}-\sqrt{\frac{2}{\unicode[STIX]{x1D716}}-2|\text{ln}\,\unicode[STIX]{x1D70F}|}.\end{eqnarray}$$

From (4.44) and using $e\unicode[STIX]{x1D719}_{\text{np}}/T_{\text{e}}\ll 1$ (the assumption behind equation (4.53)), we extract

(4.58) $$\begin{eqnarray}\frac{e\unicode[STIX]{x1D719}_{\text{np}}(x)}{T_{\text{e}}}=-\frac{1}{2}\unicode[STIX]{x1D6FC}^{2}\exp \left(3+\frac{(x-C)^{2}}{\unicode[STIX]{x1D70C}_{\text{B}}^{2}}\right)+O(\unicode[STIX]{x1D70F}^{2}\unicode[STIX]{x1D716}).\end{eqnarray}$$

In the region (4.49), equation (4.53) cannot be simplified further. Ion orbits are large, $\unicode[STIX]{x1D70C}_{\text{B}}\sqrt{\unicode[STIX]{x1D716}}\lesssim \unicode[STIX]{x1D70C}_{x}\lesssim \unicode[STIX]{x1D70C}_{\text{B}}/\sqrt{\unicode[STIX]{x1D716}}$ , and the non-parabolic piece of the electrostatic potential is small, $-\unicode[STIX]{x1D70F}^{2}\unicode[STIX]{x1D716}\lesssim e\unicode[STIX]{x1D719}_{\text{np}}(x)/T_{\text{e}}\lesssim \unicode[STIX]{x1D70F}^{2}/\unicode[STIX]{x1D716}$ (consistent with the errors in (4.56) and (4.58), noting that $\unicode[STIX]{x1D719}_{\text{np}}$ is positive in (4.56) and negative in (4.58)). From the previous paragraph we deduce that $l\sim \sqrt{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D70C}_{\text{B}}$ , and so the double derivative of the effective potential is small, $-\unicode[STIX]{x1D70F}^{2}\lesssim \unicode[STIX]{x1D712}^{\prime \prime }(x)/\unicode[STIX]{x1D6FA}^{2}=\unicode[STIX]{x1D70C}_{\text{B}}^{2}e\unicode[STIX]{x1D719}_{\text{np}}^{\prime \prime }/T_{\text{e}}\lesssim \unicode[STIX]{x1D70F}^{2}/\unicode[STIX]{x1D716}^{2}$ . The characteristic size of the periodic piece of the ion velocity is obtained from the relation $\unicode[STIX]{x1D707}\sim \unicode[STIX]{x1D70C}_{x}w_{x}\sim \unicode[STIX]{x1D70F}v_{\text{B}}^{2}/\unicode[STIX]{x1D6FA}$ (which holds provided the adiabatic invariant is still conserved), giving $\unicode[STIX]{x1D70F}\sqrt{\unicode[STIX]{x1D716}}\lesssim w_{x}/v_{\text{B}}\lesssim \unicode[STIX]{x1D70F}/\sqrt{\unicode[STIX]{x1D716}}$ . The characteristic size of the drift of closed ion orbits is obtained from (4.28), giving $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D716}^{2}/\unicode[STIX]{x1D70F}^{2}\lesssim v_{\text{d}}/v_{\text{B}}\lesssim \unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D70F}^{2}$ . Ignoring the dependences on $\unicode[STIX]{x1D716}$ and focusing only on how quantities scale with $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D70F}$ (which, from (4.7), are both exponentially small in $\unicode[STIX]{x1D716}$ ), we obtain $\unicode[STIX]{x1D70C}_{x}\sim l\sim \unicode[STIX]{x1D70C}_{\text{B}}$ , $w_{x}\sim \unicode[STIX]{x1D70F}v_{\text{B}}$ and $v_{\text{d}}\sim \unicode[STIX]{x1D6FC}v_{\text{B}}/\unicode[STIX]{x1D70F}^{2}$ . Hence, the condition $v_{\text{d}}/l\ll w_{x}/\unicode[STIX]{x1D70C}_{x}$ (4.24), which is necessary for the ion motion to be approximately periodic, implies that $v_{\text{d}}/w_{x}\sim \unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D70F}^{3}\ll 1$ . In (4.40) $\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D70F}^{3}$ is required to be small in $\unicode[STIX]{x1D6FC}$ to ensure that the motion remains periodic to lowest order in $\unicode[STIX]{x1D6FC}$ . The scaling with $\unicode[STIX]{x1D70F}^{3}$ in (4.40) implies that the kinetic model is not valid for relatively large values of $\unicode[STIX]{x1D70F}\ll 1$ . This unfortunate scaling arises because of the growth of the ion orbits: if small ion orbits reached $x=0$ , we would expect $l\sim \unicode[STIX]{x1D70C}_{x}\sim \unicode[STIX]{x1D70C}_{\text{i}}$ , $\unicode[STIX]{x1D712}^{\prime \prime }(x)\sim \unicode[STIX]{x1D6FA}^{2}$ , $w_{x}\sim \sqrt{\unicode[STIX]{x1D70F}}v_{\text{B}}$ , $v_{\text{d}}\sim \unicode[STIX]{x1D6FC}v_{\text{B}}$ , and so the condition (4.24) would give the weaker requirement $\unicode[STIX]{x1D6FC}\ll \sqrt{\unicode[STIX]{x1D70F}}$ for the model to be valid near $x=0$ . Hence, the orbit growth and the associated large polarization drift have a strong negative effect on the condition for validity of the model, multiplying the power by which $\unicode[STIX]{x1D70F}$ is raised by a factor of six. It is for this reason that, as we will see in § 5, we do not obtain numerical solutions of (3.34) for values of $\unicode[STIX]{x1D70F}$ lower than $\unicode[STIX]{x1D70F}=0.2$ .

4.1.4 Uniformly valid solution

We proceed to obtain an expression for $\unicode[STIX]{x1D719}(x)$ , equation (4.6), that is uniformly valid across the whole magnetic presheath to lowest order in $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D70F}$ . In order to do this, we first make a change of variables: guided by the form of (4.44), we introduce the function

(4.59) $$\begin{eqnarray}\unicode[STIX]{x1D713}(x)=\frac{e\unicode[STIX]{x1D719}(x)}{T_{\text{e}}}+\frac{1}{2}\unicode[STIX]{x1D6FC}^{2}\left(\exp \left(-\frac{2e\unicode[STIX]{x1D719}(x)}{T_{\text{e}}}\right)-1\right).\end{eqnarray}$$

The term $-\unicode[STIX]{x1D6FC}^{2}/2$ is small but is included in the definition (4.59) in order to have the desirable exact property that $\unicode[STIX]{x1D719}=0$ when $\unicode[STIX]{x1D713}=0$ : then, far from the wall, where $-e\unicode[STIX]{x1D719}/T_{\text{e}}\ll 1/\unicode[STIX]{x1D716}$ , the relation $\unicode[STIX]{x1D713}=(e\unicode[STIX]{x1D719}/T_{\text{e}})(1+O(\unicode[STIX]{x1D6FC}^{2}))$ is satisfied.

We proceed to show that the equation

(4.60) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{\text{B}}^{2}\unicode[STIX]{x1D713}^{\prime 2}+3+2\unicode[STIX]{x1D713}=4\exp (\unicode[STIX]{x1D713})-\exp (2\unicode[STIX]{x1D713})+O(\unicode[STIX]{x1D70F},\unicode[STIX]{x1D70F}\unicode[STIX]{x1D713}\exp ({\textstyle \frac{1}{2}}\unicode[STIX]{x1D713}),\unicode[STIX]{x1D70F}^{2})\end{eqnarray}$$

is equivalent to the equations describing the electrostatic potential in the three regions of the magnetic presheath. All the errors on the right-hand side of (4.60) are exactly equal to zero at $\unicode[STIX]{x1D713}=0$ . Moreover, the $O(\unicode[STIX]{x1D70F})$ error tends to a constant for $-\unicode[STIX]{x1D713}\gg 1$ , and this constant does not have an effect on the functional form of $\unicode[STIX]{x1D713}(x)$ . Hence, we retain the smaller errors $O(\unicode[STIX]{x1D70F}\unicode[STIX]{x1D713}\exp (\unicode[STIX]{x1D713}/2),\unicode[STIX]{x1D70F}^{2})$ . First, we compare equation (4.60) with equation (4.35), valid in the region (4.10). Since, from (4.38), $\unicode[STIX]{x1D6FC}^{2}\exp (-2e\unicode[STIX]{x1D719}/T_{\text{e}})\ll \unicode[STIX]{x1D6FC}^{2}/\unicode[STIX]{x1D70F}^{4}$ in this region, it follows that $e\unicode[STIX]{x1D719}/T_{\text{e}}=\unicode[STIX]{x1D713}+O(\unicode[STIX]{x1D6FC}^{2}/\unicode[STIX]{x1D70F}^{4})$ . Hence, equation (4.60) directly follows from (4.35), after noting that the $O(\unicode[STIX]{x1D6FC}^{2}/\unicode[STIX]{x1D70F}^{4})$ error term is smaller than the $O(\unicode[STIX]{x1D70F}^{2})$ error term (the smallest in (4.60)) because of the validity condition (4.40). Next, we compare the solution to equation (4.60) with equation (4.44) (recall that $K=\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}$ , where $\unicode[STIX]{x1D705}$ is an unknown constant of order unity, and $\bar{x}_{\text{c}}=C$ ), valid in the region (4.11) close to the wall. From (4.59), in this region $-e\unicode[STIX]{x1D719}/T_{\text{e}}=-\unicode[STIX]{x1D713}+O(1)\sim 1/\unicode[STIX]{x1D716}$ and so, from (4.48), $\exp (\unicode[STIX]{x1D713}/2)\ll \sqrt{\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D70F}}\ll \unicode[STIX]{x1D70F}$ . Therefore, all terms on the right-hand side of (4.60), except for the constant $O(\unicode[STIX]{x1D70F})$ term, become smaller than the $O(\unicode[STIX]{x1D70F}^{2})$ term. Re-expressing the constant $O(\unicode[STIX]{x1D70F})$ error term as $\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}$ and integrating gives

(4.61) $$\begin{eqnarray}\unicode[STIX]{x1D713}=-\frac{3}{2}+\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}-\frac{(x-C)^{2}}{2\unicode[STIX]{x1D70C}_{\text{B}}^{2}}+O\left(\unicode[STIX]{x1D716}\unicode[STIX]{x1D70F}\unicode[STIX]{x1D713}\exp \left(\frac{1}{2}\unicode[STIX]{x1D713}\right),\unicode[STIX]{x1D716}\unicode[STIX]{x1D70F}^{2}\right).\end{eqnarray}$$

The size of the error terms in (4.61) is obtained as follows. From $\unicode[STIX]{x1D713}\sim 1/\unicode[STIX]{x1D716}$ and $l\sim \unicode[STIX]{x1D713}^{\prime \prime }/\unicode[STIX]{x1D713}^{\prime \prime \prime }\sim \unicode[STIX]{x1D70C}_{\text{B}}^{2}/(x-C)\sim \sqrt{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D70C}_{\text{B}}$ , the error terms in the expression for $\unicode[STIX]{x1D713}^{\prime }$ are smaller by a factor of $\sqrt{\unicode[STIX]{x1D716}}$ compared with the error terms in the expression for $\unicode[STIX]{x1D713}^{\prime 2}$ , equation (4.60). Upon integrating the expression for $\unicode[STIX]{x1D713}^{\prime }$ to obtain $\unicode[STIX]{x1D713}$ , we multiply these error terms by a factor of $l\sim \sqrt{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D70C}_{\text{B}}$ . Using (4.59), and remembering that in the region (4.11) the size of the largest error term in (4.61) is $O(\unicode[STIX]{x1D716}\unicode[STIX]{x1D70F}^{2})$ , observe that equation (4.61) is equivalent to (4.44).

It only remains to be shown that (4.60) is valid in the part of the intermediate region (4.49) where neither (4.35) (the equation determining $\unicode[STIX]{x1D719}$ far from the wall) nor (4.44) (the equation determining $\unicode[STIX]{x1D719}$ close to the wall) are valid. From the equation $\unicode[STIX]{x1D719}=\unicode[STIX]{x1D719}_{\text{p}}+\unicode[STIX]{x1D719}_{\text{np}}$ , equation (4.36) for $\unicode[STIX]{x1D719}_{\text{p}}$ , and the fact that $e|\unicode[STIX]{x1D719}_{\text{np}}|/T_{\text{e}}\lesssim \unicode[STIX]{x1D70F}^{2}/\unicode[STIX]{x1D716}$ , the electrostatic potential in this region is given by $e\unicode[STIX]{x1D719}/T_{\text{e}}=-3/2+\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}-(x-C)^{2}/2\unicode[STIX]{x1D70C}_{\text{B}}^{2}+O(\unicode[STIX]{x1D70F}^{2}/\unicode[STIX]{x1D716})$ . From (4.38) for the size of the electrostatic potential in the region far from the wall, and equation (4.48) for the size of the electrostatic potential in the region close to the wall, we obtain the ordering $\unicode[STIX]{x1D6FC}/\sqrt{\unicode[STIX]{x1D716}}\unicode[STIX]{x1D70F}\lesssim \exp (e\unicode[STIX]{x1D719}/T_{\text{e}})\lesssim \unicode[STIX]{x1D70F}^{2}/\unicode[STIX]{x1D716}^{2}$ for the size of the electrostatic potential in this region. Hence, the first error term in (4.61) becomes $O(\unicode[STIX]{x1D70F}^{2}/\unicode[STIX]{x1D716})$ . Moreover, from $\unicode[STIX]{x1D716}^{4}\unicode[STIX]{x1D6FC}^{2}/\unicode[STIX]{x1D70F}^{4}\lesssim \unicode[STIX]{x1D6FC}^{2}\exp (-2e\unicode[STIX]{x1D719}/T_{\text{e}})\lesssim \unicode[STIX]{x1D716}\unicode[STIX]{x1D70F}^{2}$ we obtain $e\unicode[STIX]{x1D719}/T_{\text{e}}=\unicode[STIX]{x1D713}+O(\unicode[STIX]{x1D716}\unicode[STIX]{x1D70F}^{2})$ . Hence, from equation (4.61) and the associated $O(\unicode[STIX]{x1D70F}^{2}/\unicode[STIX]{x1D716})$ error in this region, we obtain $e\unicode[STIX]{x1D719}/T_{\text{e}}=-3/2+\unicode[STIX]{x1D705}\unicode[STIX]{x1D70F}-(x-C)^{2}/2\unicode[STIX]{x1D70C}_{\text{B}}^{2}+O(\unicode[STIX]{x1D70F}^{2}/\unicode[STIX]{x1D716})$ . Equation (4.60) is thus a good approximation also in the region (4.49).

Using the definition (4.59), and (4.5), the boundary condition at the wall is $\unicode[STIX]{x1D713}(0)=\ln \unicode[STIX]{x1D6FC}+1/2-\unicode[STIX]{x1D6FC}^{2}/2$ . This can be used to integrate equation (4.60) (neglecting the error terms) and obtain the approximate electrostatic potential solution, as in (4.6).

4.2 Hot ions ( $\unicode[STIX]{x1D70F}\gg 1/\unicode[STIX]{x1D6FC}$ )

In the limit of very hot ions, $\unicode[STIX]{x1D70F}\gg 1/\unicode[STIX]{x1D6FC}$ , we assume that the ion distribution function is a half-Maxwellian at the magnetic presheath entrance,

(4.62) $$\begin{eqnarray}f_{\infty }(\boldsymbol{v})=\frac{2n_{\infty }}{\unicode[STIX]{x03C0}^{3/2}v_{\text{t},\text{i}}^{3}}\exp \left(-\frac{|\boldsymbol{v}|^{2}}{v_{\text{t},\text{i}}^{2}}\right)\unicode[STIX]{x1D6E9}(v_{z}),\end{eqnarray}$$

where we introduced the Heaviside step function

(4.63) $$\begin{eqnarray}\unicode[STIX]{x1D6E9}(s)=\left\{\begin{array}{@{}ll@{}}1\quad & \text{for }s\geqslant 0,\\ 0\quad & \text{for }s<0.\end{array}\right.\end{eqnarray}$$

Since $\unicode[STIX]{x1D719}(\infty )=0$ , $U=|\boldsymbol{v}|^{2}/2$ and we re-express (4.62) to

(4.64) $$\begin{eqnarray}F=2n_{\infty }\left(\frac{m_{\text{i}}}{2\unicode[STIX]{x03C0}T_{\text{i}}}\right)^{3/2}\exp \left(-\frac{m_{\text{i}}U}{T_{\text{i}}}\right).\end{eqnarray}$$

Equation (4.62) is one of many choices that could be made. The reason we choose this distribution function is that it was also used in Cohen & Ryutov (Reference Cohen and Ryutov1998) in the equivalent limit of small electron temperature. We consider the limit $\unicode[STIX]{x1D70F}\ll \unicode[STIX]{x1D6FC}^{2}m_{\text{i}}/m_{\text{e}}$ in order to be consistent with condition (2.10) for an electron-repelling sheath.

For $\unicode[STIX]{x1D70F}\rightarrow \infty$ , ion orbits are undistorted by the presheath potential drop necessary to repel the electrons. We expect $e\unicode[STIX]{x1D719}(x)/T_{\text{e}}\sim 1$ , and therefore the ion flow and density can be computed using $Ze\unicode[STIX]{x1D719}(x)/T_{\text{i}}=(1/\unicode[STIX]{x1D70F})e\unicode[STIX]{x1D719}(x)/T_{\text{e}}\simeq 0$ across the magnetic presheath. The effective potential is a parabola with its minimum at $x_{\text{m}}=\bar{x}$ ,

(4.65) $$\begin{eqnarray}\unicode[STIX]{x1D712}(x,\bar{x})={\textstyle \frac{1}{2}}\unicode[STIX]{x1D6FA}^{2}(x-\bar{x})^{2}.\end{eqnarray}$$

This is an effective potential whose maximum for $x<x_{\text{m}}$ is given by

(4.66) $$\begin{eqnarray}\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})=\unicode[STIX]{x1D712}(0,\bar{x})={\textstyle \frac{1}{2}}\unicode[STIX]{x1D6FA}^{2}\bar{x}^{2}.\end{eqnarray}$$

The minimum value of $\bar{x}$ necessary for an ion at position $x$ to be in a closed orbit or an open orbit is, using (3.25) and (3.26) with $\unicode[STIX]{x1D719}(x)=0$ ,

(4.67) $$\begin{eqnarray}\bar{x}_{\text{m},\text{o}}(x)=\bar{x}_{\text{m}}(x)={\textstyle \frac{1}{2}}x.\end{eqnarray}$$

Moreover, the adiabatic invariant is $\unicode[STIX]{x1D707}=U_{\bot }/\unicode[STIX]{x1D6FA}$ .

Inserting the distribution function (4.64) into (3.22), the closed orbit density is

(4.68) $$\begin{eqnarray}\displaystyle & & \displaystyle n_{\text{i},\text{cl}}(x)=2n_{\infty }\left(\frac{m_{\text{i}}}{2\unicode[STIX]{x03C0}T_{\text{i}}}\right)^{3/2}\nonumber\\ \displaystyle & & \displaystyle \quad \times \int _{x/2}^{\infty }\unicode[STIX]{x1D6FA}\,\text{d}\bar{x}\int _{(1/2)\unicode[STIX]{x1D6FA}^{2}(x-\bar{x})^{2}}^{(1/2)\unicode[STIX]{x1D6FA}^{2}\bar{x}^{2}}\frac{2\,\text{d}U_{\bot }}{\sqrt{2(U_{\bot }-\unicode[STIX]{x1D712}(x,\bar{x}))}}\int _{U_{\bot }}^{\infty }\frac{\exp (-m_{\text{i}}U/T_{\text{i}})\,\text{d}U}{\sqrt{2(U-U_{\bot })}}.\hspace{12.0pt}\end{eqnarray}$$

Changing variables to $\tilde{v}_{y}=(\bar{x}-x)/\unicode[STIX]{x1D70C}_{\text{i}}$ , $\tilde{U} _{\bot }=m_{\text{i}}(U_{\bot }-(1/2)\unicode[STIX]{x1D6FA}^{2}(x-\bar{x})^{2})/T_{\text{i}}$ and $\tilde{U} =m_{\text{i}}(U-U_{\bot })/T_{\text{i}}$ gives

(4.69) $$\begin{eqnarray}\displaystyle n_{\text{i},\text{cl}}(x) & = & \displaystyle \frac{n_{\infty }}{\unicode[STIX]{x03C0}^{3/2}}\int _{-x/2\unicode[STIX]{x1D70C}_{\text{i}}}^{\infty }\,\text{d}\tilde{v}_{y}\exp (-\tilde{v}_{y}^{2})\int _{0}^{(x/\unicode[STIX]{x1D70C}_{\text{i}})(2\tilde{v}_{y}+x/\unicode[STIX]{x1D70C}_{\text{i}})}\tilde{U} _{\bot }^{-1/2}\exp (-\tilde{U} _{\bot })\,\text{d}\tilde{U} _{\bot }\nonumber\\ \displaystyle & & \displaystyle \times \int _{0}^{\infty }\tilde{U} ^{-1/2}\exp (-\tilde{U} )\,\text{d}\tilde{U} .\end{eqnarray}$$

Evaluating the integral over $\tilde{U}$ and the integral over $\tilde{U} _{\bot }$ leads to

(4.70) $$\begin{eqnarray}n_{\text{i},\text{cl}}(x)=\frac{n_{\infty }}{\sqrt{\unicode[STIX]{x03C0}}}\int _{-x/2\unicode[STIX]{x1D70C}_{\text{i}}}^{\infty }\exp (-\tilde{v}_{y}^{2})\text{erf}\left(\sqrt{\frac{x}{\unicode[STIX]{x1D70C}_{\text{i}}}\left(2\tilde{v}_{y}+\frac{x}{\unicode[STIX]{x1D70C}_{\text{i}}}\right)}\right)\,\text{d}\tilde{v}_{y},\end{eqnarray}$$

where we introduced the error function,

(4.71) $$\begin{eqnarray}\text{erf}(s)=\frac{2}{\sqrt{\unicode[STIX]{x03C0}}}\int _{0}^{s}\exp (-s^{\prime 2})\,\text{d}s^{\prime }.\end{eqnarray}$$

For $x\ll \unicode[STIX]{x1D70C}_{\text{i}}$ , the integral in (4.70) simplifies in the following ways: (i) the lower limit of integration can be set to $\tilde{v}_{y}=0$ , since the contribution to the integral from the integration range $[0,\infty ]$ is dominant; (ii) the factor $\sqrt{(x/\unicode[STIX]{x1D70C}_{\text{i}})(2\tilde{v}_{y}+x/\unicode[STIX]{x1D70C}_{\text{i}})}$ in the argument of the error function can be replaced by $\sqrt{2\tilde{v}_{y}x/\unicode[STIX]{x1D70C}_{\text{i}}}$ , since this replacement is accurate in most of the integration range except where $\tilde{v}_{y}\sim x/\unicode[STIX]{x1D70C}_{\text{i}}\ll 1$ ; (iii) the error function can be approximated by $\text{erf}(\sqrt{2\tilde{v}_{y}x/\unicode[STIX]{x1D70C}_{\text{i}}})\simeq 2\sqrt{2\tilde{v}_{y}x/\unicode[STIX]{x1D70C}_{\text{i}}}/\sqrt{\unicode[STIX]{x03C0}}$ for $x/\unicode[STIX]{x1D70C}_{\text{i}}\ll 1/\tilde{v}_{y}$ , which holds everywhere except in the region $\tilde{v}_{y}\gtrsim \unicode[STIX]{x1D70C}_{\text{i}}/x\gg 1$ , where the integrand is exponentially small. Hence, equation (4.70) becomes, for $x\ll \unicode[STIX]{x1D70C}_{\text{i}}$ ,

(4.72) $$\begin{eqnarray}n_{\text{i},\text{cl}}(x)\simeq \frac{2n_{\infty }}{\unicode[STIX]{x03C0}}\sqrt{\frac{2x}{\unicode[STIX]{x1D70C}_{\text{i}}}}\int _{0}^{\infty }\sqrt{\tilde{v}_{y}}\exp (-\tilde{v}_{y}^{2})\,\text{d}\tilde{v}_{y}.\end{eqnarray}$$

Re-expressing the integral over $\tilde{v}_{y}$ in terms of the standard gamma function $\unicode[STIX]{x1D6E4}$ ,

(4.73) $$\begin{eqnarray}\int _{0}^{\infty }\sqrt{\tilde{v}_{y}}\exp (-\tilde{v}_{y}^{2})\,\text{d}\tilde{v}_{y}=\frac{1}{2}\int _{0}^{\infty }\frac{\exp (-\unicode[STIX]{x1D709})}{\unicode[STIX]{x1D709}^{1/4}}\,\text{d}\unicode[STIX]{x1D709}=\frac{1}{2}\unicode[STIX]{x1D6E4}\left(\frac{3}{4}\right),\end{eqnarray}$$

we simplify (4.72) to

(4.74) $$\begin{eqnarray}n_{\text{i},\text{cl}}(x)\simeq \frac{\sqrt{2}}{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D6E4}\left(\frac{3}{4}\right)\sqrt{\frac{x}{\unicode[STIX]{x1D70C}_{\text{i}}}}n_{\infty }.\end{eqnarray}$$

The density of open orbits is given by

(4.75) $$\begin{eqnarray}\displaystyle & & \displaystyle n_{\text{i},\text{op}}(x)=\int _{(1/2)x}^{\infty }\unicode[STIX]{x1D6FA}\,\text{d}\bar{x}\int _{(1/2)\unicode[STIX]{x1D6FA}^{2}\bar{x}^{2}}^{\infty }\frac{F(\unicode[STIX]{x1D6FA}^{2}\bar{x}^{2}/2,U)}{\sqrt{2(U-\unicode[STIX]{x1D712}_{\text{M}}(\bar{x}))}}\nonumber\\ \displaystyle & & \displaystyle \quad \times \,(\sqrt{2(\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})-\unicode[STIX]{x1D712}(x,\bar{x})+\unicode[STIX]{x1D6E5}_{\text{M}}(\bar{x},U))}-\sqrt{2(\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})-\unicode[STIX]{x1D712}(x,\bar{x}))})\,\text{d}U.\end{eqnarray}$$

Note that, in (4.75), we have used

(4.76) $$\begin{eqnarray}\unicode[STIX]{x1D707}={\textstyle \frac{1}{2}}\unicode[STIX]{x1D6FA}\bar{x}^{2}\end{eqnarray}$$

for the adiabatic invariant of ions with $U_{\bot }=\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})=\unicode[STIX]{x1D712}(0,\bar{x})=\unicode[STIX]{x1D6FA}^{2}\bar{x}^{2}/2$ . Using (3.30), we obtain

(4.77) $$\begin{eqnarray}\unicode[STIX]{x1D6E5}_{\text{M}}=2\unicode[STIX]{x1D6FC}\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FA}\bar{x}\sqrt{2(U-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D6FA}^{2}\bar{x}^{2})}.\end{eqnarray}$$

Then, using the dimensionless integration variables $\tilde{v}_{z}=\sqrt{m_{\text{i}}(U-\unicode[STIX]{x1D6FA}^{2}\bar{x}^{2}/2)/T_{\text{i}}}$ and $\tilde{\bar{x}}=\bar{x}/\unicode[STIX]{x1D70C}_{\text{i}}$ , equation (4.75) reduces to

(4.78) $$\begin{eqnarray}\displaystyle n_{\text{i},\text{op}}(x) & = & \displaystyle \frac{2n_{\infty }}{\unicode[STIX]{x03C0}^{3/2}}\int _{x/2\unicode[STIX]{x1D70C}_{\text{i}}}^{\infty }\,\text{d}\tilde{\bar{x}}\exp (-\tilde{\bar{x}}^{2})\int _{0}^{\infty }\exp (-\tilde{v}_{z}^{2})\nonumber\\ \displaystyle & & \displaystyle \times \left(\sqrt{\frac{x}{\unicode[STIX]{x1D70C}_{\text{i}}}\left(2\tilde{\bar{x}}-\frac{x}{\unicode[STIX]{x1D70C}_{\text{i}}}\right)+4\unicode[STIX]{x1D6FC}\unicode[STIX]{x03C0}\tilde{\bar{x}}\tilde{v}_{z}}-\sqrt{\frac{x}{\unicode[STIX]{x1D70C}_{\text{i}}}\left(2\tilde{\bar{x}}-\frac{x}{\unicode[STIX]{x1D70C}_{\text{i}}}\right)}\right)\,\text{d}\tilde{v}_{z}.\end{eqnarray}$$

Equation (4.78) does not simplify further for general values of $x$ , but can be simplified for $x\ll \unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{\text{i}}$ . Evaluating the ion density at $x=0$ using (4.78), we obtain

(4.79) $$\begin{eqnarray}n_{\text{i},\text{op}}(0)=\frac{1}{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D6E4}^{2}\left(\frac{3}{4}\right)\sqrt{\unicode[STIX]{x1D6FC}}n_{\infty }.\end{eqnarray}$$

We then proceed to evaluate $n_{\text{i},\text{op}}(x)-n_{\text{i},\text{op}}(0)$ for $x\ll \unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{\text{i}}$ . In this ordering, $(x/\unicode[STIX]{x1D70C}_{\text{i}})(2\tilde{\bar{x}}-x/\unicode[STIX]{x1D70C}_{\text{i}})\ll 4\unicode[STIX]{x1D6FC}\unicode[STIX]{x03C0}\tilde{\bar{x}}\tilde{v}_{z}$ and so

(4.80) $$\begin{eqnarray}n_{\text{i},\text{op}}(x)-n_{\text{i},\text{op}}(0)\simeq -\frac{2n_{\infty }}{\unicode[STIX]{x03C0}^{3/2}}\int _{x/2\unicode[STIX]{x1D70C}_{\text{i}}}^{\infty }\,\text{d}\tilde{\bar{x}}\exp (-\tilde{\bar{x}}^{2})\sqrt{\frac{x}{\unicode[STIX]{x1D70C}_{\text{i}}}\left(2\tilde{\bar{x}}-\frac{x}{\unicode[STIX]{x1D70C}_{\text{i}}}\right)}\int _{0}^{\infty }\exp (-\tilde{v}_{z}^{2})\,\text{d}\tilde{v}_{z}.\end{eqnarray}$$

Note that there is a small integration region, $\tilde{v}_{z}\lesssim x/\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{\text{i}}$ , where $(x/\unicode[STIX]{x1D70C}_{\text{i}})(2\tilde{\bar{x}}-x/\unicode[STIX]{x1D70C}_{\text{i}})\gtrsim 4\unicode[STIX]{x1D6FC}\unicode[STIX]{x03C0}\tilde{\bar{x}}\tilde{v}_{z}$ , but the contribution to the integral from this region is higher order in $x/\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{\text{i}}\ll 1$ . In (4.80) we take $\sqrt{(x/\unicode[STIX]{x1D70C}_{\text{i}})(2\tilde{\bar{x}}+x/\unicode[STIX]{x1D70C}_{\text{i}})}\simeq \sqrt{2\tilde{\bar{x}}x/\unicode[STIX]{x1D70C}_{\text{i}}}$ , since this is accurate everywhere except where $\tilde{\bar{x}}\sim x/\unicode[STIX]{x1D70C}_{\text{i}}\ll \unicode[STIX]{x1D6FC}$ , and evaluate the integrals over $\tilde{v}_{z}$ and $\tilde{\bar{x}}$ to obtain

(4.81) $$\begin{eqnarray}n_{\text{i},\text{op}}(x)-n_{\text{i},\text{op}}(0)\simeq -\frac{\sqrt{2}}{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D6E4}\left(\frac{3}{4}\right)\sqrt{\frac{x}{\unicode[STIX]{x1D70C}_{\text{i}}}}n_{\infty }.\end{eqnarray}$$

From (4.81), the open orbit density near $x=0$ decreases proportionally to $\sqrt{x}$ , but the increase of the closed orbit density in (4.74) is faster by a factor of $2$ , leading to the total ion density increasing proportionally to $\sqrt{x}$ ,

(4.82) $$\begin{eqnarray}n_{\text{i}}(x)-n_{\text{i}}(0)\simeq \frac{\sqrt{2}}{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D6E4}\left(\frac{3}{4}\right)\sqrt{\frac{x}{\unicode[STIX]{x1D70C}_{\text{i}}}}n_{\infty },\end{eqnarray}$$

for $x\ll \unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}_{\text{i}}$ .

The ion density profile for $\unicode[STIX]{x1D70F}\rightarrow \infty$ is, according to (3.24), the sum of (4.70) and (4.78). The potential profile is obtained by imposing quasineutrality and inverting the Boltzmann relation for the electron density, to find

(4.83) $$\begin{eqnarray}\frac{e\unicode[STIX]{x1D719}(x)}{T_{\text{e}}}=\ln \left(\frac{n_{\text{i}}(x)}{n_{\infty }}\right).\end{eqnarray}$$

The potential drop across the magnetic presheath can be calculated by using $n_{\text{i},\text{cl}}(0)=0$ (from (4.70)) and (4.79),

(4.84) $$\begin{eqnarray}\frac{e\unicode[STIX]{x1D719}(0)}{T_{\text{e}}}=\ln \left(\frac{\unicode[STIX]{x1D6E4}^{2}(3/4)}{\unicode[STIX]{x03C0}}\sqrt{\unicode[STIX]{x1D6FC}}\right)\simeq \ln (0.48\sqrt{\unicode[STIX]{x1D6FC}}).\end{eqnarray}$$

Inserting the distribution function (4.64) and the value of $\bar{x}_{\text{m},\text{o}}$ in (4.67) into (3.40), the distribution of the ion velocity component perpendicular to the wall at $x=0$ is

(4.85) $$\begin{eqnarray}f_{0x}(v_{x})=\frac{n_{\infty }}{v_{\text{t},\text{i}}\unicode[STIX]{x03C0}}\unicode[STIX]{x1D6E9}(-v_{x})\int _{0}^{\infty }\exp (-\tilde{\bar{x}}^{2})\left[1-\text{erf}\left(\frac{v_{x}^{2}}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FC}\tilde{\bar{x}}v_{\text{t},\text{i}}^{2}}\right)\right]\,\text{d}\tilde{\bar{x}}.\end{eqnarray}$$

Inserting the distribution function (4.64) into (3.41), the distribution of the ion velocity components parallel to the wall at $x=0$ is

(4.86) $$\begin{eqnarray}f_{0yz}(v_{y},v_{z})=\frac{4\sqrt{\unicode[STIX]{x1D6FC}}n_{\infty }}{\unicode[STIX]{x03C0}}\frac{\sqrt{v_{y}v_{z}}}{v_{\text{t},\text{i}}^{3}}\exp \left(-\frac{v_{y}^{2}+v_{z}^{2}}{v_{\text{t},\text{i}}^{2}}\right)\unicode[STIX]{x1D6E9}(v_{y})\unicode[STIX]{x1D6E9}(v_{z}).\end{eqnarray}$$

To conclude, we briefly point out and resolve an apparent contradiction in the validity of our kinetic model when $\unicode[STIX]{x1D70F}\gg 1/\unicode[STIX]{x1D6FC}$ . In Geraldini et al. (Reference Geraldini, Parra and Militello2018), we found that the self-consistent electrostatic potential prohibits the presence of ions entering the Debye sheath with zero velocity normal to the wall. This is in apparent contradiction with the situation described in this section: when undistorted circular orbits reach the wall, there are ion trajectories tangential to the wall and thus there is a finite number of ions which have a normal component of the velocity equal to zero. This is reflected in the fact that, from (4.85), $f_{0x}(0)\neq 0$ . In reality, there is a small region near $x=0$ in which the electric field distorts ion orbits just before they reach the wall, so that $\unicode[STIX]{x1D712}_{\text{M}}(\bar{x})=\unicode[STIX]{x1D712}(x_{\text{M}},\bar{x})$ with $x_{\text{M}}\ll \unicode[STIX]{x1D70C}_{\text{i}}$ . The quasi-tangential ions (with $v_{x}\simeq 0$ ) must be accelerated to values of $v_{x}$ such that the Bohm condition (3.42) is satisfied with the equality sign. If these very slow ions do not accelerate to large enough values of $v_{x}$ , the integral on the left-hand side of (3.42) becomes too large and the Bohm condition cannot be satisfied. Conversely, if these ions are accelerated too much towards the wall, the Bohm condition cannot be satisfied with the equality sign, as in (3.42), which is in contradiction with our theory. Thus, one can think of the real distribution function as the distribution function in (4.85) (which is plotted as a dashed line in figure 6 c), but shifted in such a way that the peak of the distribution function is at $v_{x}=-\bar{v}$ instead of $v_{x}=0$ , and the distribution function is effectively equal to zero for $|v_{x}|<\bar{v}\ll \sqrt{\unicode[STIX]{x1D6FC}}v_{\text{t},\text{i}}$ .Footnote ⁴ Since the width of the distribution function, $\sqrt{\unicode[STIX]{x1D6FC}}v_{\text{t},\text{i}}$ , is much larger than $v_{\text{B}}$ for $\unicode[STIX]{x1D70F}\gg 1/\unicode[STIX]{x1D6FC}$ , the Bohm integral on the right-hand side evaluates approximately to $f_{0x}(\bar{v})/\bar{v}$ for the real distribution function. Then, approximating $f_{0x}(\bar{v})\sim n_{\text{i}}(0)/\sqrt{\unicode[STIX]{x1D6FC}}v_{\text{t},\text{i}}$ , we obtain the estimate $\bar{v}\sim v_{\text{B}}/\sqrt{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}}$ to satisfy the Bohm condition (3.42). Hence, the final piece of the electrostatic potential drop, which is responsible for distorting the ion orbits enough to satisfy the kinetic Bohm condition, is smaller than the total electrostatic potential drop by a factor of $m_{\text{i}}\bar{v}^{2}/T_{\text{e}}\sim 1/\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}\ll 1$ . Note that the pair of conditions (2.10) and $1/\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}\ll 1$ require $m_{\text{i}}/m_{\text{e}}\gg \unicode[STIX]{x1D70F}^{3}$ to be satisfied. The size $h$ of the region near $x=0$ where this final potential drop occurs is obtained by balancing the electric force, $Ze\unicode[STIX]{x1D719}^{\prime }\sim ZT_{\text{e}}/h\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}$ , with the magnetic force $Zev_{y}B\sim m_{\text{i}}\unicode[STIX]{x1D6FA}v_{\text{t},\text{i}}$ , giving $h/\unicode[STIX]{x1D70C}_{\text{i}}\sim 1/\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}^{2}\ll 1$ . The spatial resolution necessary to resolve this region can be prohibitively high even for $1/\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}\sim 1$ , since $\unicode[STIX]{x1D70F}\gg 1$ , and it is for this reason that, as we will see in § 5, we do not obtain numerical solutions for values of $\unicode[STIX]{x1D70F}$ larger than $\unicode[STIX]{x1D70F}=10$ .

5.1 Boundary conditions

The ion distribution function, $f_{\infty }(\boldsymbol{v})$ , that enters the magnetic presheath is determined by a kinetic solution of the bulk plasma or of the collisional presheath. Without such a solution, there is an infinite possible number of distribution functions we could choose as boundary conditions. We proceed to parameterize a set of such distribution functions using $\unicode[STIX]{x1D70F}=T_{\text{i}}/ZT_{\text{e}}$ . We design them to recover the two limits studied in § 4.

We proceed to make a number of observations about the properties that an appropriate set of distribution functions must satisfy. Considering the strong resemblance of the kinetic Chodura condition (3.35) with the kinetic Bohm condition, whose equality form is (3.42), we choose that (3.35) be satisfied with the equality sign,

(5.1) $$\begin{eqnarray}\int \frac{f_{\infty }(\boldsymbol{v})}{v_{z}^{2}}\,\text{d}^{3}v=\frac{n_{\infty }}{v_{\text{B}}^{2}}.\end{eqnarray}$$

The assumption behind equation (5.1) is that, just as the magnetic presheath solution self-consistently satisfies the kinetic Bohm condition with the equality sign, the collisional presheath will self-consistently satisfy the kinetic Chodura condition with the equality sign. In order to be consistent with the models in § 4 in the limits $\unicode[STIX]{x1D70F}\rightarrow 0$ and $\unicode[STIX]{x1D70F}\rightarrow \infty$ , we also choose a set of distribution functions that,

1. (i) for $\unicode[STIX]{x1D70F}\rightarrow 0$ is a Maxwellian that peaks at $v_{z}=v_{\text{B}}$ ;

2. (ii) for $\unicode[STIX]{x1D70F}\rightarrow \infty$ is a half-Maxwellian that peaks at $v_{z}=0$ .

A set of distribution functions that has all the above properties is

(5.2) $$\begin{eqnarray}f_{\infty }(\boldsymbol{v})=\left\{\begin{array}{@{}ll@{}}\displaystyle {\mathcal{N}}n_{\infty }\frac{4v_{z}^{2}}{\unicode[STIX]{x03C0}^{3/2}v_{\text{t},\text{i}}^{5}}\exp \left(-\frac{|\boldsymbol{v}-uv_{\text{t},\text{i}}\hat{\boldsymbol{z}}|^{2}}{v_{\text{t},\text{i}}^{2}}\right)\unicode[STIX]{x1D6E9}(v_{z})\quad & \text{for }\unicode[STIX]{x1D70F}\leqslant 1,\\ \displaystyle {\mathcal{N}}n_{\infty }\frac{4v_{z}^{2}}{\unicode[STIX]{x03C0}^{3/2}v_{\text{t},\text{i}}^{3}(v_{\text{t},\text{i}}^{2}+rv_{z}^{2})}\exp \left(-\frac{|\boldsymbol{v}|^{2}}{v_{\text{t},\text{i}}^{2}}\right)\unicode[STIX]{x1D6E9}(v_{z})\quad & \text{for }\unicode[STIX]{x1D70F}>1,\end{array}\right.\end{eqnarray}$$

where $\unicode[STIX]{x1D6E9}$ is the Heaviside step function defined in (4.63). The values of $u$ and $r$ in (5.2) are chosen such that condition (5.1) is satisfied. For $\unicode[STIX]{x1D70F}\leqslant 1$ , decreasing $\unicode[STIX]{x1D70F}$ increases the parameter $u$ , which increases the flow velocity of the distribution function. For $\unicode[STIX]{x1D70F}\ll 1$ and $u\gg 1$ , the distribution function tends to a shifted Maxwellian with flow velocity given by $uv_{\text{t},\text{i}}$ . For $\unicode[STIX]{x1D70F}>1$ , the parameter $r$ increases from $0$ to $\infty$ for increasing $\unicode[STIX]{x1D70F}$ . For $r>1$ , the distribution function becomes small for values of $v_{z}$ smaller than $v_{\text{t},\text{i}}/\sqrt{r}$ , and thus the parameter $r$ determines the region of velocity space around $v_{z}=0$ where there are almost no particles. For $\unicode[STIX]{x1D70F}\gg 1$ and $r\gg 1$ , the distribution function at the entrance of the magnetic presheath is a half-Maxwellian with a very narrow region around $v_{z}=0$ where the distribution function vanishes. The quantity ${\mathcal{N}}$ is a normalization constant that ensures that

(5.3) $$\begin{eqnarray}n_{\infty }=\int f_{\infty }(\boldsymbol{v})\,\text{d}^{3}v.\end{eqnarray}$$

Note that, from (3.20), (3.21) and (5.2), we can write the distribution function in the form $F(\unicode[STIX]{x1D707},U)$ ,

(5.4) $$\begin{eqnarray}F(\unicode[STIX]{x1D707},U)=\left\{\begin{array}{@{}ll@{}}\displaystyle {\mathcal{N}}n_{\infty }\frac{8(U-\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D707})}{\unicode[STIX]{x03C0}^{3/2}v_{\text{t},\text{i}}^{5}}\exp \left[-\frac{2}{v_{\text{t},\text{i}}^{2}}(\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D707}+(\sqrt{2(U-\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D707})}-uv_{\text{t},\text{i}})^{2})\right]\quad & \text{for }\unicode[STIX]{x1D70F}\leqslant 1,\\ \displaystyle {\mathcal{N}}n_{\infty }\frac{8(U-\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D707})}{\unicode[STIX]{x03C0}^{3/2}v_{\text{t},\text{i}}^{3}(v_{\text{t},\text{i}}^{2}+2r(U-\unicode[STIX]{x1D6FA}\unicode[STIX]{x1D707}))}\exp \left(-\frac{2U}{v_{\text{t},\text{i}}^{2}}\right)\quad & \text{for }\unicode[STIX]{x1D70F}>1.\end{array}\right.\end{eqnarray}$$

The value of the normalization constant ${\mathcal{N}}$ is, from (5.3),

(5.5) $$\begin{eqnarray}{\mathcal{N}}=\left\{\begin{array}{@{}ll@{}}\displaystyle \left[(1+2u^{2})(1+\text{erf}(u))+\frac{2u}{\sqrt{\unicode[STIX]{x03C0}}}\exp (-u^{2})\right]^{-1}\quad & \text{for }\unicode[STIX]{x1D70F}\leqslant 1,\\ \displaystyle r^{3/2}\left[2\sqrt{r}-2\sqrt{\unicode[STIX]{x03C0}}\exp \left(\frac{1}{r}\right)\left(1-\text{erf}\left(\frac{1}{\sqrt{r}}\right)\right)\right]^{-1}\quad & \text{for }\unicode[STIX]{x1D70F}>1.\end{array}\right.\end{eqnarray}$$

The values of $u$ and $r$ are, from (5.1), given by

(5.6) $$\begin{eqnarray}1+\text{erf}(u)=\unicode[STIX]{x1D70F}\left[(1+2u^{2})(1+\text{erf}(u))+\frac{2u}{\sqrt{\unicode[STIX]{x03C0}}}\exp (-u)\right],\end{eqnarray}$$

(5.7) $$\begin{eqnarray}\displaystyle & & \displaystyle r\sqrt{\unicode[STIX]{x03C0}}\exp \left(\frac{1}{r}\right)\left(1-\text{erf}\left(\frac{1}{\sqrt{r}}\right)\right)\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D70F}\left[2\sqrt{r}-2\sqrt{\unicode[STIX]{x03C0}}\exp \left(\frac{1}{r}\right)\left(1-\text{erf}\left(\frac{1}{\sqrt{r}}\right)\right)\right],\end{eqnarray}$$

and are plotted as functions of $\unicode[STIX]{x1D70F}$ in figure 4. The fluid velocity in the $z$ direction at the magnetic presheath entrance, $u_{z\infty }$ , is given by the equations

(5.8) $$\begin{eqnarray}\frac{u_{z\infty }}{v_{\text{t},\text{i}}}=\frac{\displaystyle u(3+2u^{2})(1+\text{erf}(u))+\frac{2}{\sqrt{\unicode[STIX]{x03C0}}}\exp (-u^{2})(1+u^{2})}{\displaystyle (1+2u^{2})(1+\text{erf}(u))+\frac{2}{\sqrt{\unicode[STIX]{x03C0}}}u\exp (-u^{2})}\quad \text{for }\unicode[STIX]{x1D70F}\leqslant 1,\end{eqnarray}$$

and

(5.9) $$\begin{eqnarray}\frac{u_{z\infty }}{v_{\text{t},\text{i}}}=\frac{2}{\sqrt{\unicode[STIX]{x03C0}r}}\frac{\displaystyle r-\exp \left(\frac{1}{r}\right)E_{1}\left(\frac{1}{r}\right)}{\displaystyle 2\sqrt{r}-2\sqrt{\unicode[STIX]{x03C0}}\exp \left(\frac{1}{r}\right)\left(1-\text{erf}\left(\frac{1}{\sqrt{r}}\right)\right)}\quad \text{for }\unicode[STIX]{x1D70F}>1.\end{eqnarray}$$

In (5.9), we have introduced the exponential integral,

(5.10) $$\begin{eqnarray}E_{1}(\unicode[STIX]{x1D709})=\int _{\unicode[STIX]{x1D709}}^{\infty }\frac{\exp (-\unicode[STIX]{x1D702})}{\unicode[STIX]{x1D702}}\,\text{d}\unicode[STIX]{x1D702}.\end{eqnarray}$$

Using (5.6)–(5.9), in figure 4 we plot the value of $u_{z\infty }$ as a function of $\unicode[STIX]{x1D70F}$ . Equations (5.5)–(5.9) are derived in appendix F.

Figure 4. (a) The numbers $\ln r$ and $\ln u$ as a function of the parameter $\ln \unicode[STIX]{x1D70F}$ . (b) The flow velocity at the magnetic presheath entrance, $u_{z\infty }$ , as a function of the parameter $\unicode[STIX]{x1D70F}$ . The dashed line corresponds to $\unicode[STIX]{x1D70F}>1$ , where $r$ (instead of $u$ ) is used to parameterize the distribution functions in (5.2). Note that $u_{z\infty }/v_{\text{t},\text{i}}\rightarrow 1/\sqrt{2\unicode[STIX]{x1D70F}}$ for $\unicode[STIX]{x1D70F}\rightarrow 0$ , $u_{z\infty }/v_{\text{t},\text{i}}=2/\sqrt{\unicode[STIX]{x03C0}}\approx 1.13$ for $\unicode[STIX]{x1D70F}=1$ and $u_{z\infty }/v_{\text{t},\text{i}}\rightarrow 1/\sqrt{\unicode[STIX]{x03C0}}\approx 0.56$ for $\unicode[STIX]{x1D70F}\rightarrow \infty$ .

To conclude this subsection, we verify that the distribution functions have the required properties at $\unicode[STIX]{x1D70F}\rightarrow 0$ and $\unicode[STIX]{x1D70F}\rightarrow \infty$ . From (5.6), note that taking the limit $\unicode[STIX]{x1D70F}\rightarrow 0$ leads to $u\simeq \sqrt{1/2\unicode[STIX]{x1D70F}}\gg 1$ , so that the ion distribution function $f_{\infty }$ in (5.2) is indeed a Maxwellian that peaks at $v_{z}=v_{\text{t},\text{i}}/\sqrt{2\unicode[STIX]{x1D70F}}=v_{\text{B}}$ . Moreover, note that taking the limit $\unicode[STIX]{x1D70F}\rightarrow \infty$ in (5.7) leads to $r\simeq (2\unicode[STIX]{x1D70F})^{2}/\unicode[STIX]{x03C0}\gg 1$ , so that $f_{\infty }$ is a half Maxwellian that peaks at $v_{z}=0$ . In the next subsection, we present the numerical results obtained for finite values of $\unicode[STIX]{x1D70F}$ .

5.2 Numerical solutions

The numerical scheme presented in Geraldini et al. (Reference Geraldini, Parra and Militello2018) is used to obtain numerical solutions to the quasineutrality equation (3.34) for values of $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D70F}$ in the range $0.01\leqslant \unicode[STIX]{x1D6FC}\leqslant 0.2$ (approximately corresponding to $0.57^{\circ }\leqslant \unicode[STIX]{x1D6FC}\leqslant 11^{\circ }$ ) and $0.2\leqslant \unicode[STIX]{x1D70F}\leqslant 10$ . We define a quantity

(5.11) $$\begin{eqnarray}{\tilde{n}}(x)=1-\frac{Zn_{\text{i}}(x)}{n_{\text{e}}(x)}.\end{eqnarray}$$

In the numerical scheme, all quantities are discretized and so ${\tilde{n}}_{\unicode[STIX]{x1D708}}={\tilde{n}}(x_{\unicode[STIX]{x1D708}})$ is a set of values defined on a grid of values of $x_{\unicode[STIX]{x1D708}}$ , where $\unicode[STIX]{x1D708}$ is an index running from $0$ to some value $\unicode[STIX]{x1D702}$ . The exact solution to equation (3.34) has ${\tilde{n}}(x)=0$ everywhere, but numerically ${\tilde{n}}_{\unicode[STIX]{x1D708}}$ cannot be made to be arbitrarily small at all grid points. Hence, we use the following convergence criterion to define what constitutes a valid numerical solution to equation (3.34),

(5.12) $$\begin{eqnarray}\left(\frac{1}{\unicode[STIX]{x1D702}+1}\mathop{\sum }_{\unicode[STIX]{x1D708}=0}^{\unicode[STIX]{x1D702}}{\tilde{n}}_{\unicode[STIX]{x1D708}}^{2}\right)^{1/2}<E,\end{eqnarray}$$

where $E$ is a small number. An iteration scheme, outlined in Geraldini et al. (Reference Geraldini, Parra and Militello2018), is performed to find the numerical electrostatic potential solution $\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D708}}=\unicode[STIX]{x1D719}(x_{\unicode[STIX]{x1D708}})$ for a given value of $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D70F}$ . The solution numerically satisfies the quasineutrality equation with an error $E=0.7\,\%$ for all values of $\unicode[STIX]{x1D70F}$ except for $\unicode[STIX]{x1D70F}=0.2$ , where $E=1.2\,\%$ .

The electrostatic potential drop across the magnetic presheath is shown in figure 5(a) as a function of $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D70F}$ . The numerical results approaching $\unicode[STIX]{x1D70F}=0.2$ and $\unicode[STIX]{x1D70F}=10$ are consistent with the results obtained using (4.5) (valid for small $\unicode[STIX]{x1D70F}$ , $3/|\text{ln}\,\unicode[STIX]{x1D6FC}|<1/|\text{ln}\,\unicode[STIX]{x1D70F}|\ll 1$ ) and using (4.84) (valid for $\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}\gg 1$ ), shown with dashed lines. The shaded region is where we expect the assumption of an electron-repelling wall not to be suitable for deuterium ions, $\unicode[STIX]{x1D6FC}\lesssim \sqrt{1+\unicode[STIX]{x1D70F}}\sqrt{m_{\text{e}}/m_{\text{i}}}\sim 0.02\sqrt{1+\unicode[STIX]{x1D70F}}$ . Considering the unshaded region in figure 5(a), the potential drop with finite ion temperature is up to 10 %–15 % smaller than the cold ion ( $\unicode[STIX]{x1D70F}=0$ ) potential drop. For a fixed angle, $\unicode[STIX]{x1D6FC}=0.05~\text{rad}\approx 3^{\circ }$ , the electrostatic potential profiles for different values of $\unicode[STIX]{x1D70F}$ are shown in figure 5(b). The blue dashed curve labelled ‘0’ in figure 5(b) is obtained from (4.6), while the red dashed curve marked ‘ $\infty$ ’ is obtained from (4.83). The numerical profiles are consistent with the limits $\unicode[STIX]{x1D70F}=0$ and $\unicode[STIX]{x1D70F}=\infty$ .

Figure 5. (a) The electrostatic potential drop across the magnetic presheath $\unicode[STIX]{x1D719}(0)$ is shown as a function of the angle $\unicode[STIX]{x1D6FC}$ and the parameter $\unicode[STIX]{x1D70F}$ . The region where $\unicode[STIX]{x1D6FC}\lesssim \sqrt{1+\unicode[STIX]{x1D70F}}\sqrt{m_{\text{e}}/m_{\text{i}}}$ , and therefore the ordering (2.10) breaks down, is shaded. (b) Electrostatic potential profiles for $\unicode[STIX]{x1D6FC}=0.05$ at different values of $\unicode[STIX]{x1D70F}$ , marked on the curves.

Figure 6. The distributions of the component $v_{z}$ of the ion velocity at the magnetic presheath entrance $x\rightarrow \infty$  (a,b) and the component $v_{x}$ of the velocity at the Debye sheath entrance $x=0$  (c,d) are shown for $\unicode[STIX]{x1D6FC}=0.05$ for three different values of the parameter $\unicode[STIX]{x1D70F}$ , labelled next to the corresponding curve. The velocities are normalized to $v_{\text{t},\text{i}}$ in (a,c) and to $v_{\text{B}}$ in (b,d). Magnetized ions at the magnetic presheath entrance move parallel to the magnetic field. Hence, $v_{z}$ is responsible for the flow of ions to the wall. At the Debye sheath entrance, the ion flow towards the wall is determined by $|v_{x}|$ . The red dashed lines in (a,c) are the distribution functions in the limit $\unicode[STIX]{x1D70F}\rightarrow \infty$ . The blue vertical dashed lines in (b,d) are the cold ion distribution functions, $\unicode[STIX]{x1D70F}=0$ .

Figure 7. The ion distribution functions $f_{\infty yz}(v_{y},v_{z})$ (a) and $f_{0yz}(v_{y},v_{z})$ (b) for $\unicode[STIX]{x1D6FC}=0.05$ and, from top to bottom, for $\unicode[STIX]{x1D70F}=0.2$ , $\unicode[STIX]{x1D70F}=1$ , $\unicode[STIX]{x1D70F}=5$ and $\unicode[STIX]{x1D70F}=\infty$ (see § 4.2). The Bohm speed $v_{\text{B}}/v_{\text{t},\text{i}}=1/\sqrt{2\unicode[STIX]{x1D70F}}$ is marked as a horizontal line in all panels, and also as a vertical line on the right panels.

While the solution to a fluid model can give a good estimate of the electrostatic potential profile in the magnetic presheath at some range of finite temperatures, it provides no information on the velocity distribution of the ions. The ions hitting the wall can cause sputtered neutral impurities to be thrown back into the plasma, and the sputtering yield is sensitively dependent on the kinetic energy and angle of incidence of the ion on the target. Hence, it is important to predict the ion distribution function at the wall. Since in the Debye sheath ions only undergo an acceleration towards the wall, see e.g. Riemann (Reference Riemann1991), the distribution function of ions at the Debye sheath entrance is expected to be similar in shape to the distribution function at the wall. For different values of $\unicode[STIX]{x1D70F}$ , in figure 6 we plot the distribution function $f_{0x}(v_{x})$ (defined in (3.40)) and compare it with the boundary condition $f_{\infty z}(v_{z})=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }f_{\infty }(\boldsymbol{v})\,\text{d}v_{y}\,\text{d}v_{x}$ . Equation (4.85) is the dashed curve in figure 6(c). The equality form of the kinetic Bohm condition (3.42) (Riemann Reference Riemann1991) is approximately numerically satisfied for all distribution functions in the parameter range of the presented simulations; recall that (3.42) is an analytical property of the self-consistent solution of (3.34) (Geraldini et al. Reference Geraldini, Parra and Militello2018). Note that at values of $\unicode[STIX]{x1D70F}$ larger than $\unicode[STIX]{x1D70F}=10$ , it becomes computationally expensive to resolve the sharp gradient of the distribution function near $v_{x}=0$ , as discussed at the end of § 4.2. In all of our simulations, the distribution $f_{0x}(v_{x})$ is found to be both narrower and more centred around $v_{\text{B}}$ than $f_{\infty z}(v_{z})$ . In figure 7, we plot the functions $f_{\infty yz}(v_{y},v_{z})=\int _{-\infty }^{\infty }f_{\infty }(\boldsymbol{v})\,\text{d}v_{x}$ and $f_{0yz}(v_{y},v_{z})$ . Equation (4.86) is the bottom right panel in figure 7. For $\unicode[STIX]{x1D70F}\lesssim 1$ , the ions have very large tangential velocities at $x=0$ (compared with $x=\infty$ ) due to the large increase in the $y$ -component of the velocity, related to the $\boldsymbol{E}\times \boldsymbol{B}$ drift acquired by the ion orbit in the magnetic presheath.

We can summarize the numerical results for the distribution function as follows:

1. (i) for $\unicode[STIX]{x1D70F}\gg 1$ , the velocity components tangential to the wall, $v_{y}$ and $v_{z}$ , remain unaffected while the velocity component normal to the wall, $v_{x}$ , becomes of the order of whichever is largest between $\sqrt{\unicode[STIX]{x1D6FC}}v_{\text{t},\text{i}}$ and $v_{\text{B}}\sim v_{\text{t},\text{i}}/\sqrt{\unicode[STIX]{x1D70F}}$ ;

2. (ii) for $\unicode[STIX]{x1D70F}\lesssim 1$ , all velocity components are affected by the magnetic presheath electric field and become of order $v_{\text{B}}$ (ignoring factors of $|\text{ln}\,\unicode[STIX]{x1D6FC}|$ ).

For large ion temperatures, $\unicode[STIX]{x1D70F}\gtrsim 5$ , the velocity component normal to the wall at the Debye sheath entrance is small because the electrostatic potential necessary to repel electrons barely affects the ions. In this case, there are two regimes of interest. Firstly, if $1\ll \unicode[STIX]{x1D70F}\ll 1/\unicode[STIX]{x1D6FC}$ , most ions are accelerated to $|v_{x}|\simeq v_{\text{B}}\sim v_{\text{t},\text{i}}/\sqrt{\unicode[STIX]{x1D70F}}\ll v_{\text{t},\text{i}}$ , as expected if the Bohm condition (3.42) is to be satisfied, and the spread of the ion distribution function in the $x$ direction, $f_{0x}(v_{x})$ , is $v_{\text{B}}$ . The numerical solution for $\unicode[STIX]{x1D70F}=5$ and $\unicode[STIX]{x1D6FC}=0.05$ , where $f_{0x}(v_{x})$ is shown in figure 6(c,d), is adequately described by this regime. Secondly, if $\unicode[STIX]{x1D70F}$ is such that $\unicode[STIX]{x1D70F}\gg 1/\unicode[STIX]{x1D6FC}$ , the velocity spread of the distribution function is $\sqrt{\unicode[STIX]{x1D6FC}}v_{\text{t},\text{i}}$ , satisfying $v_{\text{B}}\ll \sqrt{\unicode[STIX]{x1D6FC}}v_{\text{t},\text{i}}\ll v_{\text{t},\text{i}}$ ; this regime corresponds to the limit taken in § 4.2, where $f_{0x}(v_{x})$ is given in (4.85) and plotted in figure 6(c) as a red dashed line. For $\unicode[STIX]{x1D6FC}\sim 1/\unicode[STIX]{x1D70F}$ , the velocity spread is $\sqrt{\unicode[STIX]{x1D6FC}}v_{\text{t},\text{i}}\sim v_{\text{B}}\ll v_{\text{t},\text{i}}$ , as both of the estimates above are valid. The tangential velocity of a typical ion with $\unicode[STIX]{x1D70F}\gtrsim 5$ remains roughly of the same size, $v_{y}\sim v_{z}\sim v_{\text{t},\text{i}}$ , and therefore the angle between the ion trajectory and the wall is shallow at the Debye sheath entrance. For $\unicode[STIX]{x1D70F}\lesssim 1$ , the typical size of all the velocity components is $v_{\text{B}}$ and thus the angle between the ion trajectory and the wall is of order unity. Hence, an ion is expected to impinge on the wall at an angle whose size is small when $\unicode[STIX]{x1D70F}\gg 1$ and order unity when $\unicode[STIX]{x1D70F}\lesssim 1$ .