Dependence on ion temperature of shallow-angle magnetic presheaths with adiabatic electrons

Abstract

The magnetic presheath is a boundary layer occurring when magnetized plasma is in contact with a wall and the angle $\unicode[STIX]{x1D6FC}$ between the wall and the magnetic field $\boldsymbol{B}$ is oblique. Here, we consider the fusion-relevant case of a shallow-angle, $\unicode[STIX]{x1D6FC}\ll 1$, electron-repelling sheath, with the electron density given by a Boltzmann distribution, valid for $\unicode[STIX]{x1D6FC}/\sqrt{\unicode[STIX]{x1D70F}+1}\gg \sqrt{m_{\text{e}}/m_{\text{i}}}$, where $m_{\text{e}}$ is the electron mass, $m_{\text{i}}$ is the ion mass, $\unicode[STIX]{x1D70F}=T_{\text{i}}/ZT_{\text{e}}$, $T_{\text{e}}$ is the electron temperature, $T_{\text{i}}$ is the ion temperature and $Z$ is the ionic charge state. The thickness of the magnetic presheath is of the order of a few ion sound Larmor radii $\unicode[STIX]{x1D70C}_{\text{s}}=\sqrt{m_{\text{i}}(ZT_{\text{e}}+T_{\text{i}})}/ZeB$, where e is the proton charge and $B=|\boldsymbol{B}|$ is the magnitude of the magnetic field. We study the dependence on $\unicode[STIX]{x1D70F}$ of the electrostatic potential and ion distribution function in the magnetic presheath by using a set of prescribed ion distribution functions at the magnetic presheath entrance, parameterized by $\unicode[STIX]{x1D70F}$. The kinetic model is shown to be asymptotically equivalent to Chodura’s fluid model at small ion temperature, $\unicode[STIX]{x1D70F}\ll 1$, for $|\text{ln}\,\unicode[STIX]{x1D6FC}|>3|\text{ln}\,\unicode[STIX]{x1D70F}|\gg 1$. In this limit, despite the fact that fluid equations give a reasonable approximation to the potential, ion gyro-orbits acquire a spatial extent that occupies a large portion of the magnetic presheath. At large ion temperature, $\unicode[STIX]{x1D70F}\gg 1$, relevant because $T_{\text{i}}$ is measured to be a few times larger than $T_{\text{e}}$ near divertor targets of fusion devices, ions reach the Debye sheath entrance (and subsequently the wall) at a shallow angle whose size is given by $\sqrt{\unicode[STIX]{x1D6FC}}$ or $1/\sqrt{\unicode[STIX]{x1D70F}}$, depending on which is largest.

1 Introduction

Plasma–wall interaction is important in systems such as plasma discharges (Lieberman & Lichtenberg 2005), fusion devices (Stangeby 2000), magnetic filters (Anders, Anders & Brown 1995), plasma probes (Hutchinson 2002) and thrusters (Martinez-Sanchez & Pollard 1998). In the context of nuclear fusion research, the plasma–wall interaction at the divertor or limiter targets of fusion devices is directly related to the boundary conditions to be imposed (Loizu et al. 2012) on models of plasma in the open-field line region (the scrape-off layer). The heat flux reaching the wall of the device must be minimized and one way to do so is to make the magnetic field lines reach the divertor or limiter target at a shallow angle $\unicode[STIX]{x1D6FC}\ll 1$ ( $\unicode[STIX]{x1D6FC}$ is measured in radians unless otherwise indicated) (Loarte et al. 2007). In typical devices, $\unicode[STIX]{x1D6FC}\sim 0.05{-}0.2~\text{radians}({\sim}3{-}12^{\circ })$ , and in ITER it is expected that $\unicode[STIX]{x1D6FC}\sim 0.04~\text{radians}\sim 2.5^{\circ }$ (Pitts et al. 2009). Hence, it is crucial to understand plasma–wall interaction at such small angles in order to address the problem of exhaust in fusion plasmas.

The magnetic presheath (Chodura 1982) is a boundary layer with a width of a few ion sound Larmor radii, $\unicode[STIX]{x1D70C}_{\text{s}}=\sqrt{m_{\text{i}}(ZT_{\text{e}}+T_{\text{i}})}/ZeB$ , next to the wall, where $T_{\text{i}}$ and $T_{\text{e}}$ are the ion and the electron temperatures respectively, $m_{\text{i}}$ is the ion mass, $Z$ is the ionic charge state, e is the proton charge and $B$ is the magnetic field strength. This region is characterized by a balance between electric and magnetic forces on the ions. Closer to the wall, in steady state, there is a non-neutral layer called the Debye sheath which typically repels electrons. The Debye sheath has a thickness of a few Debye lengths, $\unicode[STIX]{x1D706}_{\text{D}}=\sqrt{\unicode[STIX]{x1D716}_{0}T_{\text{e}}/\text{e}^{2}n_{\text{e}}}$ , where $n_{\text{e}}$ is the electron density and $\unicode[STIX]{x1D716}_{0}$ is the permittivity of free space, and is characterized by the electric forces dominating the ion dynamics. The Debye length is generally much smaller than the ion sound gyroradius, $\unicode[STIX]{x1D706}_{\text{D}}\ll \unicode[STIX]{x1D70C}_{\text{s}}$ , and therefore the magnetic presheath can be solved as a separate quasineutral system. Moreover, we assume that ions collide for the last time when they are a distance $d_{\text{coll}}\gg \unicode[STIX]{x1D70C}_{\text{s}}$ away from the wall, and therefore the magnetic presheath is collisionless. The latter assumption is expected to hold in attached divertor regimes of operation, whereas in detached divertors the temperature is so low that the collisional scale may be small enough to make $d_{\text{coll}}\sim \unicode[STIX]{x1D70C}_{\text{s}}$ (Tskhakaya 2017).

Due to their small mass relative to the ions, most electrons are usually repelled by the sheath electric field, and we thus assume the electrons to be in thermal equilibrium. This assumption becomes less valid when the angle between the magnetic field and the wall is very small and when the ion temperature is sufficiently large compared to the electron temperature (as we will see in § 2). For the ions, many magnetic presheath models use fluid equations, which rely on $T_{\text{i}}=0$ (Chodura 1982; Riemann 1994; Ahedo 1997; Ahedo & Carralero 2009). However, in the vicinity of the divertor target of a typical tokamak plasma, the ion temperature is at least as large as the electron temperature, $T_{\text{i}}\sim T_{\text{e}}$ (Mosetto et al. 2015), making a kinetic treatment of the ions necessary (Siddiqui et al. 2016). In this paper, we study the dependence of the magnetic presheath on the parameter

(1.1) $$\begin{eqnarray}\unicode[STIX]{x1D70F}=\frac{T_{\text{i}}}{ZT_{\text{e}}},\end{eqnarray}$$

which is of fundamental importance in kinetic models of turbulence. For $Z=1$ , $\unicode[STIX]{x1D70F}$ is simply the ratio of ion to electron temperature. Early attempts to solve the magnetic presheath by retaining the ion distribution function made use of analytical solutions of the ion trajectories (Holland, Fried & Morales 1993; Parks & Lippmann 1994; Cohen & Ryutov 1998; Daube & Riemann 1999) with a variety of assumptions, giving valuable insight into the characteristics of the ion motion in the magnetic presheath. Later, there were several particle-in-cell (PIC) studies of the Chodura and Debye sheaths (Tskhakaya & Kuhn 2003, 2004; Khaziev & Curreli 2015), as well as some kinetic simulations using a Eulerian–Vlasov approach (Coulette & Manfredi 2014, 2016). Here, we use analytical solutions of the ion trajectories in a magnetic field whose angle with the wall is small (Holland et al. 1993; Cohen & Ryutov 1998). An asymptotic theory of magnetic presheaths with $\unicode[STIX]{x1D6FC}\ll 1$ , and an associated numerical scheme to obtain self-consistent solutions of the electrostatic potential, was presented in detail in Geraldini, Parra & Militello (2017, 2018). The method also determines the ion distribution function at the Debye sheath entrance. Although only valid for grazing angles, this method has yielded several analytical results, is valid within the current paradigm of plasma exhaust in a fusion device and is computationally fast.

Figure 1. An ion gyro-orbit is shown schematically at a distance of approximately an ion gyroradius $\unicode[STIX]{x1D70C}_{\text{i}}$ from the wall (grey horizontal surface). The magnetic field is constant and the angle between the magnetic field and the wall is small, $\unicode[STIX]{x1D6FC}\ll 1$ (in radians). The electric field is directed towards the wall and is a function of the coordinate $x$ .

This paper is structured as follows. The orderings and geometry of the magnetic presheath are discussed in § 2. We use the shallow-angle ( $\unicode[STIX]{x1D6FC}\ll 1$ ) kinetic model described in Geraldini et al. (2017, 2018) which we briefly review in § 3. In § 4 we discuss the two limits of small $\unicode[STIX]{x1D70F}$ ( $\unicode[STIX]{x1D70F}\ll 1$ ) and large $\unicode[STIX]{x1D70F}$ ( $\unicode[STIX]{x1D70F}\gg 1/\unicode[STIX]{x1D6FC}$ ) analytically using our kinetic model. In particular, we show that our kinetic model is consistent with: the fluid model of Chodura (1982) for $\unicode[STIX]{x1D70F}\ll 1$ ; a kinetic model that assumes a half-Maxwellian ion distribution function, briefly discussed in § 3B of Cohen & Ryutov (1998), for $\unicode[STIX]{x1D70F}\gg 1/\unicode[STIX]{x1D6FC}$ . Using a set of boundary conditions that recovers those used in the small and large temperature limits, numerical results of the shallow-angle kinetic model are obtained for finite values of $\unicode[STIX]{x1D70F}$ . The boundary conditions and numerical results are presented in § 5. We conclude by summarizing and discussing our results in § 6.

In order to help the reader keep track of the several symbols used in this paper (many of which were introduced in references Geraldini et al. (2017, 2018)), we include a glossary in appendix A. For each symbol, the glossary includes a brief verbal definition (or an equation) and a reference to the equation where it first appears in the main text.

2 Orderings

Consider a magnetized plasma in steady state, in the region $x\geqslant 0$ , in contact with a wall, defined as the plane $x=0$ . We use a set of orthogonal axes, depicted in the top-right corner of figure 1, with the $x$ -axis aligned normal to the wall, and the $y$ - and $z$ -axes aligned in the two directions parallel to the wall. The magnetic field is uniform and given by

(2.1) $$\begin{eqnarray}\boldsymbol{B}=B\cos \unicode[STIX]{x1D6FC}\hat{\boldsymbol{z}}-B\sin \unicode[STIX]{x1D6FC}\hat{\boldsymbol{x}}.\end{eqnarray}$$

In (2.1), $\hat{\boldsymbol{x}}$ and $\hat{\boldsymbol{z}}$ denote unit vectors parallel to the $x$ and $z$ -axes and $\unicode[STIX]{x1D6FC}\ll 1$ is the small angle between the magnetic field and the wall. The components of the ion velocity in the three directions are $v_{x}$ , $v_{y}$ and $v_{z}$ . The system is uniform in the plane parallel to the wall, and thus every quantity is independent of the value of $y$ and $z$ . The ion motion can therefore be described using four coordinates: $x$ , $v_{x}$ , $v_{y}$ , and $v_{z}$ .

We consider a plasma with a single ion species and an electron species. An electric field normal to the wall is present to repel the most mobile of the plasma species – the electrons – away from the wall,

(2.2) $$\begin{eqnarray}\boldsymbol{E}=-\unicode[STIX]{x1D719}^{\prime }(x)\hat{\boldsymbol{x}},\end{eqnarray}$$

where $\unicode[STIX]{x1D719}$ is the electrostatic potential and a prime denotes differentiation with respect to $x$ . The electrostatic potential is assumed to monotonically converge to some value at $x\rightarrow \infty$ , and this value is set to be $\unicode[STIX]{x1D719}=0$ . Moreover, it has been shown that $\unicode[STIX]{x1D719}(x)-\unicode[STIX]{x1D719}(0)\propto \sqrt{x}$ at $x\rightarrow 0$ (see (141) and (142) in Geraldini et al. (2018)), so that the magnetic presheath electric field diverges at the Debye sheath entrance. ¹ The coordinate system and the geometry are depicted in figure 1.

Since the electric field is present to repel electrons from the wall, the characteristic size of the electrostatic potential $\unicode[STIX]{x1D719}$ is given by

(2.3) $$\begin{eqnarray}e\unicode[STIX]{x1D719}\sim T_{\text{e}}.\end{eqnarray}$$

Ions gain energies of the order of $Ze\unicode[STIX]{x1D719}\sim ZT_{\text{e}}$ ; at such energies, they have a velocity of the order of the Bohm speed,

(2.4) $$\begin{eqnarray}v_{\text{B}}=\sqrt{\frac{ZT_{\text{e}}}{m_{\text{i}}}}.\end{eqnarray}$$

If the energy gained by the ions during this acceleration is smaller than their thermal energy, $ZT_{\text{e}}\lesssim T_{\text{i}}$ , the typical ion velocity is the ion thermal speed,

(2.5) $$\begin{eqnarray}v_{\text{t},\text{i}}=\sqrt{\frac{2T_{\text{i}}}{m_{\text{i}}}}.\end{eqnarray}$$

From (2.4) and (2.5) it follows that, in general, the ion’s speed has a characteristic size equal to the ion sound speed ²   $c_{\text{s}}$ ,

(2.6) $$\begin{eqnarray}|\boldsymbol{v}|\sim c_{s}=\sqrt{v_{\text{B}}^{2}+\frac{1}{2}v_{\text{t},\text{i}}^{2}}=\sqrt{\frac{ZT_{\text{e}}+T_{\text{i}}}{m_{\text{i}}}}.\end{eqnarray}$$

Note that $c_{\text{s}}=v_{\text{B}}$ when $\unicode[STIX]{x1D70F}=0$ and $c_{\text{s}}=\sqrt{T_{\text{i}}/m_{\text{i}}}=v_{\text{t},\text{i}}/\sqrt{2}$ when $\unicode[STIX]{x1D70F}=\infty$ .

We proceed to argue that the typical size of the magnetic presheath, denoted $d_{\text{mps}}$ , is the ion sound gyroradius (Chodura 1982),

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{\text{s}}=\frac{c_{\text{s}}}{\unicode[STIX]{x1D6FA}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}=ZeB/m_{\text{i}}$ is the typical ion gyrofrequency. We consider the two limits $\unicode[STIX]{x1D70F}\ll 1$ and $\unicode[STIX]{x1D70F}\gg 1$ separately. When the ion temperature is much smaller than the electron temperature, $\unicode[STIX]{x1D70F}\ll 1$ , the only way by which ions can acquire the Bohm velocity $v_{\text{B}}$ in the direction normal to the wall – necessary to satisfy the Bohm condition at the Debye sheath entrance (Riemann 1991) – is if the electric field becomes large enough that it demagnetizes the ion orbits. From the ordering $|\boldsymbol{v}|\sim v_{\text{B}}$ for the ion speed and by balancing the magnetic and electric forces, we obtain $\unicode[STIX]{x1D719}^{\prime }(x)\sim T_{\text{e}}/ed_{\text{mps}}\sim v_{\text{B}}B$ ; hence, $d_{\text{mps}}\sim \unicode[STIX]{x1D70C}_{\text{B}}$ , where

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{\text{B}}=\frac{v_{\text{B}}}{\unicode[STIX]{x1D6FA}}.\end{eqnarray}$$

Since $\unicode[STIX]{x1D70C}_{\text{B}}\simeq \unicode[STIX]{x1D70C}_{\text{s}}$ for $\unicode[STIX]{x1D70F}\ll 1$ , this is consistent with $d_{\text{mps}}\sim \unicode[STIX]{x1D70C}_{\text{s}}$ . When the ion temperature is large, $\unicode[STIX]{x1D70F}\gg 1$ , the radius of gyration of the ions is larger than $\unicode[STIX]{x1D70C}_{\text{B}}$ . The length scale of the magnetic presheath is set by the ion density variation, and therefore must satisfy $d_{\text{mps}}\sim \unicode[STIX]{x1D70C}_{\text{i}}=v_{\text{t},\text{i}}/\unicode[STIX]{x1D6FA}$ , where $\unicode[STIX]{x1D70C}_{\text{i}}$ is the ion gyroradius. Since $\unicode[STIX]{x1D70C}_{\text{i}}\simeq \sqrt{2}\unicode[STIX]{x1D70C}_{\text{s}}$ for $\unicode[STIX]{x1D70F}\gg 1$ , this is consistent with $d_{\text{mps}}\sim \unicode[STIX]{x1D70C}_{\text{s}}$ . When $\unicode[STIX]{x1D70F}\sim 1$ , both arguments are valid, since $\unicode[STIX]{x1D70C}_{\text{i}}\sim \unicode[STIX]{x1D70C}_{\text{B}}\sim \unicode[STIX]{x1D70C}_{\text{s}}$ .

The assumption of an electron-repelling wall is not valid for any value of $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D70F}$ . We proceed to obtain the condition on $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D70F}$ for this assumption to be valid. We expect electrons to travel at characteristic velocities equal to their thermal speed,

(2.9) $$\begin{eqnarray}v_{\text{t},\text{e}}=\sqrt{\frac{2T_{\text{e}}}{m_{\text{e}}}},\end{eqnarray}$$

where $m_{\text{e}}$ is the electron mass. The typical electron velocity is so large, $v_{\text{t},\text{e}}\gg v_{\text{B}}$ , that electrons are virtually unaffected by the electric field, since they are subject to magnetic forces, $ev_{y}B\sim ev_{\text{t},\text{e}}B$ , much larger than electric forces, $e\unicode[STIX]{x1D719}^{\prime }\lesssim ev_{\text{B}}B$ . Moreover, electron gyro-orbits are small, $\unicode[STIX]{x1D70C}_{\text{e}}\ll \unicode[STIX]{x1D70C}_{\text{s}}$ . Hence, averaging over the small-scale gyro-motion, the electrons in the magnetic presheath stream parallel to the magnetic field at a velocity of the order of $v_{\text{t},\text{e}}$ . Conversely, the ion motion close to the wall in the magnetic presheath consists of gyro-orbits distorted by the electric field, and so the ions reach the wall travelling at a velocity of the order of $c_{\text{s}}$ . Considering an ion and an electron initially at a distance ${\sim}\unicode[STIX]{x1D70C}_{\text{s}}$ from the wall, and remembering that the electron motion is constrained to be parallel to the magnetic field, the electron has to travel a longer distance than the ion by a factor of $1/\unicode[STIX]{x1D6FC}$ . However, the electron travels this distance at a speed larger than the ion’s by a factor $v_{\text{t},\text{e}}/c_{\text{s}}=\sqrt{m_{\text{i}}/(m_{\text{e}}(1+\unicode[STIX]{x1D70F}))}$ . Hence, the electron reaches the wall in a shorter time than the ion if

(2.10) $$\begin{eqnarray}\sqrt{\frac{m_{\text{e}}}{m_{\text{i}}}}\sqrt{\unicode[STIX]{x1D70F}+1}\ll \unicode[STIX]{x1D6FC}.\end{eqnarray}$$

If condition (2.10) is satisfied, the wall repels most of the electrons back into the plasma, and the ordering for the magnitude of the ion velocity, equation (2.6), is self-consistent.

For an electron-repelling wall, the electron distribution function is typically considered to be well approximated by a Maxwellian. The reason for this is that the collisional processes outside of the collisionless sheath and presheath drive it to a Maxwellian, and the sheath repels most of the electrons back into the plasma. Hence, the electron density is assumed to be given by a Boltzmann distribution.

3 Kinetic ion model

In this section we briefly review the shallow-angle kinetic model presented in detail in Geraldini et al. (2017, 2018). In § 3.1 we use the asymptotic expansion in $\unicode[STIX]{x1D6FC}\ll 1$ to write the ion velocity in terms of slowly varying orbit parameters, finding that there are approximately periodic solutions to the ion motion. In moving across the magnetic presheath, ions conserve two quantities to lowest order in $\unicode[STIX]{x1D6FC}$ : the total energy $U$ and an adiabatic invariant $\unicode[STIX]{x1D707}$ (Cohen & Ryutov 1998). The adiabatic invariant is directly related to the approximately periodic nature of the ion motion, and coincides with the usual magnetic moment only when the electric field variation over the ion gyroradius scale is small.

When written as a function of $\unicode[STIX]{x1D707}$ and $U$ , the distribution function is constant across the magnetic presheath, to lowest order in $\unicode[STIX]{x1D6FC}$ . This is used, in § 3.2, to write an expression for the ion density. In § 3.3, we write the quasineutrality equation and summarize the main equations of the shallow-angle kinetic model.

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. 2018). 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. (2017)). 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. ³ 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 1968; Taylor & Hastie 1968; Catto 1978; Antonsen & Lane 1980; Frieman & Chen 1982). 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 1998; Geraldini et al. 2017). 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. (2017, 2018). 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. (2017, 2018), 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. (2018).

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

(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. (2018) 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 (1991), it is shown that in most presheath models the Bohm condition is self-consistently satisfied in the equality form, as in (3.42).

4 The limits of small and large ion temperature

In order to study the effect of ion temperature in a kinetic model of the magnetic presheath, it is essential to check that the model is consistent with expected results in appropriate limits of $\unicode[STIX]{x1D70F}$ . Here, we study the limits $\unicode[STIX]{x1D6FC}^{1/3}\ll \unicode[STIX]{x1D70F}\ll 1$ (cold ions) and $\unicode[STIX]{x1D70F}\gg 1/\unicode[STIX]{x1D6FC}\gg 1$ (hot ions) of the kinetic model introduced in § 3. For cold ions, $\unicode[STIX]{x1D70F}\ll 1$ and so the ion distribution function is narrow when compared to the Bohm velocity, as $v_{\text{t},\text{i}}=\sqrt{2\unicode[STIX]{x1D70F}}v_{\text{B}}\ll v_{\text{B}}$ . Since we have argued in § 2 that the typical ion speed in the magnetic presheath is, for $\unicode[STIX]{x1D70F}\ll 1$ , the Bohm velocity, the ion distribution function can be taken to be a delta function to lowest order in $\unicode[STIX]{x1D70F}$ . The result of this approximation is the fluid theory first presented in Chodura (1982). Conversely, for hot ions, $\unicode[STIX]{x1D70F}\gg 1$ and the size of the Bohm velocity is negligible compared to the thermal velocity, $v_{\text{B}}=v_{\text{t},\text{i}}/\sqrt{2\unicode[STIX]{x1D70F}}\ll v_{\text{t},\text{i}}$ . This means that an accurate knowledge of the whole of the ion distribution function is important when studying the limit $\unicode[STIX]{x1D70F}\gg 1$ . For the purpose of this paper, we assume that the ion distribution function is a half-Maxwellian when entering the collisionless magnetic presheath, as was done in Cohen & Ryutov (1998).

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 1982) 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 (1982) (and later in Riemann (1994)), 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. (2017). 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 (1998) 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. (2018), 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}}$ . ⁴ 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 Numerical results

In this section, we study the magnetic presheath at finite values of $\unicode[STIX]{x1D70F}$ using numerical simulations. First, in § 5.1, we parameterize a set of magnetic presheath entrance distribution functions using $\unicode[STIX]{x1D70F}$ in a way that is consistent with the limits of small ( $\unicode[STIX]{x1D70F}\ll 1$ ) and large ( $\unicode[STIX]{x1D70F}\gg 1/\unicode[STIX]{x1D6FC}$ ) ion temperature studied in the previous section. Then, in § 5.2, we present numerical solutions of the electrostatic potential profile and of the ion distribution function at the Debye sheath entrance.

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. (2018) 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. (2018), 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 (1991), 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 1991) 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. 2018). 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$ .

6 Conclusion

In this paper we have studied the dependence of a grazing-angle electron-repelling magnetic presheath on ion temperature using the kinetic model in Geraldini et al. (2017, 2018). The cold ion limit, $\unicode[STIX]{x1D70F}=T_{\text{i}}/ZT_{\text{e}}\ll 1$ , is described by Chodura’s fluid model, giving the solution (4.6) to lowest order in $\unicode[STIX]{x1D6FC}$ . In the limit $3/|\text{ln}\,\unicode[STIX]{x1D6FC}|<1/|\text{ln}\,\unicode[STIX]{x1D70F}|\ll 1$ , we have analytically shown that the solution of the shallow-angle kinetic model is asymptotically equivalent to the fluid solution in (4.6) to lowest order in $\unicode[STIX]{x1D70F}$ and $\unicode[STIX]{x1D6FC}$ . The numerical results for $\unicode[STIX]{x1D70F}=0.2$ , shown in figure 5, confirm that the kinetic solution tends to the fluid solution at small $\unicode[STIX]{x1D70F}$ . We have also shown that, despite the ordering $\unicode[STIX]{x1D70C}_{\text{i}}\ll \unicode[STIX]{x1D70C}_{\text{B}}$ for $\unicode[STIX]{x1D70F}\ll 1$ , the characteristic spatial extent of ion gyromotion in the direction normal to the wall grows to $\unicode[STIX]{x1D70C}_{\text{B}}\sqrt{|\text{ln}\,\unicode[STIX]{x1D6FC}|}$ as the ion approaches the wall, thus becoming comparable to the size of the magnetic presheath. The growth of ion gyro-orbits is accompanied by a decrease in the gyration velocity in order to conserve the adiabatic invariant, as can be seen in figure 3. Hence, if the ion thermal energy is too small, the gyration velocity of ion orbits becomes comparable to the orbit drift, thus invalidating the gyrokinetic assumption underlying our kinetic model. For the largest orbits, our kinetic model breaks down if $\unicode[STIX]{x1D70F}^{3}\lesssim \unicode[STIX]{x1D6FC}$ .

In the hot ion limit, $\unicode[STIX]{x1D70F}\rightarrow \infty$ , our model corresponds to a model briefly studied in Cohen & Ryutov (1998), which we described in § 4.2. From the electrostatic potential results shown in figure 5, the largest values of ion temperature, $\unicode[STIX]{x1D70F}=5$ and $\unicode[STIX]{x1D70F}=10$ , are consistent with the large ion temperature limit. Our results for the distribution function at the Debye sheath entrance (shown in figures 6 and 7, for $\unicode[STIX]{x1D6FC}=0.05$ ) show that the angle between a typical ion trajectory and the wall is smaller at large values of $\unicode[STIX]{x1D70F}$ . Correspondingly, ions that have traversed the magnetic presheath tend to have a smaller spread of the normal component of the velocity, $v_{x}$ . The latter effect, which is also present for $\unicode[STIX]{x1D70F}\sim 1$ and $|\text{ln}\,\unicode[STIX]{x1D6FC}|\gg 1$ (to be treated in a future publication), is particularly prominent for $\unicode[STIX]{x1D70F}\gg 1$ . For $1\ll \unicode[STIX]{x1D70F}\ll 1/\unicode[STIX]{x1D6FC}$ ions reach the wall with a range of velocities that is centred at $v_{x}\approx v_{\text{B}}$ (consistent with the kinetic Bohm condition (3.42)) and whose spread is $v_{\text{B}}\sim v_{\text{t},\text{i}}/\sqrt{\unicode[STIX]{x1D70F}}$ (see, for example, $\unicode[STIX]{x1D6FC}=0.05$ and $\unicode[STIX]{x1D70F}=5$ in figure 6). For $\unicode[STIX]{x1D70F}\gg 1/\unicode[STIX]{x1D6FC}$ , ions reach the wall with a range of velocities that is peaked at $v_{x}\sim v_{\text{B}}/\sqrt{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70F}}\ll v_{\text{B}}$ (essentially $v_{x}\simeq 0$ ), and whose spread is $\unicode[STIX]{x1D6FC}^{1/2}v_{\text{t},\text{i}}$ (see the plot for $\unicode[STIX]{x1D6FC}=0.05$ and $\unicode[STIX]{x1D70F}\rightarrow \infty$ in figure 6).

Chodura’s fluid model of the magnetic presheath can give electrostatic potential profiles that are qualitatively similar to the ones obtained using our kinetic model for $\unicode[STIX]{x1D70F}\lesssim 1$ (see figure 5). At larger values of $\unicode[STIX]{x1D70F}$ , the quantitative difference between the fluid profile and the kinetic profile becomes more evident. For very large values of $\unicode[STIX]{x1D70F}$ , the potential drop normalized to electron temperature is up to 30 % smaller than for $\unicode[STIX]{x1D70F}=0$ . However, at such large values of $\unicode[STIX]{x1D70F}$ the electrons would not be adiabatic, as was assumed here, since the assumption $\unicode[STIX]{x1D6FC}/\sqrt{1+\unicode[STIX]{x1D70F}}\gg \sqrt{m_{\text{e}}/m_{\text{i}}}$ would not be satisfied. In this case, the Debye sheath would not repel most of the electrons back into the magnetic presheath, and a kinetic treatment of both ions and electrons would be necessary. The ordering $\unicode[STIX]{x1D6FC}/\sqrt{1+\unicode[STIX]{x1D70F}}\sim \sqrt{m_{\text{e}}/m_{\text{i}}}$ has mostly been avoided in the literature to date, but is becoming more relevant for fusion devices since $\sqrt{m_{\text{e}}/m_{\text{i}}}\sim 0.02~\text{rad}\approx 1^{\circ }$ for deuterium plasmas, $\unicode[STIX]{x1D70F}\gtrsim 1$ near divertor targets (Mosetto et al. 2015) and $\unicode[STIX]{x1D6FC}\sim 2.5^{\circ }$ is expected in ITER (Pitts et al. 2009).