1. Introduction
When plasma is in contact with a solid surface, such as in fusion experiments (Stangeby Reference Stangeby2000), Hall thrusters (Boeuf Reference Boeuf2017), plasma probes (Hutchinson Reference Hutchinson2002), magnetic filters (Anders, Anders & Brown Reference Anders, Anders and Brown1995) and orbiting spacecraft (Hastings Reference Hastings1995), the resulting interaction affects both the plasma and the surface. Among the many plasma–surface interaction processes, one that is of particular concern is sputtering, where an ion from the plasma reaches the surface material and knocks an atom off the surface. Ionization of sputtered atoms in the plasma produces impurities, thus altering the plasma. Moreover, in the long run sputtering causes erosion of the surface material. The amount of sputtering depends on a wide variety of factors, including surface material, surface roughness, plasma conditions and velocity distributions of particles striking the target (Cohen & Ryutov Reference Cohen and Ryutov1998b; Khaziev & Curreli Reference Khaziev and Curreli2015; Siddiqui et al. Reference Siddiqui, Thompson, Jackson, Kim, Hershkowitz and Scime2016; Drobny et al. Reference Drobny, Hayes, Curreli and Ruzic2017; Krasheninnikov & Kukushkin Reference Krasheninnikov and Kukushkin2017; Lasa et al. Reference Lasa, Canik, Blondel, Younkin, Curreli, Drobny, Roth, Cianciosa, Elwasif and Green2020).
In this paper, we focus on the calculation of the distribution function of plasma ions striking the solid surface. We consider the target surface, or wall, to be smooth, planar and absorbing all incident particles. We consider a plasma magnetized by a uniform magnetic field 
$\boldsymbol {B}$, with one ion species. The angle between the magnetic field and the wall is taken to be small, 
$\alpha \ll 1$ (measured in radians unless otherwise indicated). This situation is particularly relevant in fusion plasmas, where divertors are designed so that the angle between incident magnetic field lines and the target surface is as small as possible. We define a set of right-handed Cartesian axes 
$(x,y,z)$ where 
$x$ measures the distance from the wall, 
$z$ measure displacements in the direction tangential to the wall, such that the magnetic field is in the 
$x$–
$z$ plane, and 
$y$ measures displacements in the remaining direction. The axes are shown on the top-right of figure 1. For simplicity, we assume no gradients tangential to the wall. Thus, the only gradients are in the 
$x$ direction.
Ion gyro-orbits, whose gyro-radius is 
$\rho _{i}$, reaching the target when the angle between the magnetic field 
$\boldsymbol {B}$ and the target is small, 
$\alpha \ll 1$. The axes 
$(x,y,z)$ are labelled. (a) With no normal electric field, the circular orbit moves closer to the target by 
$\alpha \rho _{i}$ after a gyro-period and thus the normal velocity of an ion at the target is 
$v_x \sim \sqrt {\alpha } v_{t,i}$. (b,c) With the magnetic presheath and Debye sheath electric field 
$\boldsymbol {E}$, ions are accelerated to 
$v_x \sim \sqrt {\alpha v_{t,i}^{2} + v_{B}^{2} }$.
The standard picture of the plasma–wall boundary is as follows. Close to the wall, there is a thin positively charged layer called the Debye sheath, with a characteristic size of a few Debye lengths 
$\lambda _{{D}} = \sqrt {\epsilon _0 T_{{e}} / e^{2} n_{{e}} }$, where a strong electric field 
$\boldsymbol {E} = - \boldsymbol {\nabla } \phi$ directed towards the target is present to repel electrons (Riemann Reference Riemann1991; Hershkowitz Reference Hershkowitz2005; Baalrud et al. Reference Baalrud, Scheiner, Yee, Hopkins and Barnat2019). Here, 
$e$ is the proton charge, 
$n_{{e}}$ is the number density of the electrons, 
$\epsilon _0$ is the permittivity of free space, 
$T_{{e}}$ is the temperature of the electrons and 
$\phi (x)$ is the electrostatic potential as a function of the distance from the wall. The purpose of the electric field is to achieve a steady state with comparable (or, in ambipolar conditions, equal) fluxes of ions and electrons to the wall. The size of the electrostatic potential drop necessary to repel electrons is 
$| \phi |\sim T_{e} / e$. The kinetic energy gained by an ion of charge 
$Ze$ in such a potential is 
$Z e |\phi | \sim ZT_{e}$. Hence, the parameter
\begin{equation} \tau = \frac{T_{i}}{ZT_{e}}, \end{equation}
where 
$T_{i}$ is the ion temperature, is a measure of the ratio of ion thermal energy divided by ion kinetic energy gained from the electric field. At the edge of a fusion device one often finds 
$\tau \gtrsim 1$ (Mosetto et al. Reference Mosetto, Halpern, Jolliet, Loizu and Ricci2015). Poisson's equation,
\begin{equation} \varepsilon_0 \phi'' (x) = Zen_{i}(x) - en_{e}(x) \rm , \end{equation}
relates the charge separation to the electrostatic potential in the Debye sheath, where 
$x \sim \lambda _{D}$. Here a prime denotes differentiation with respect to the argument, in this case 
$x$, of the function. At distances from the wall comparable to the ion sound gyro-radius, 
$\rho _{{s}}$, the ion population is depleted owing to a combination of ion gyro-orbit losses and acceleration of ions by the electric field, as schematically shown in figure 1. Here, 
$\rho _{s} = c_{s} / \varOmega$, where 
$c_{s} = \sqrt { ( ZT_{e} + T_{i} ) / m_{i}}$ is the ion sound speed, 
$\varOmega = ZeB / m_{i}$ is the ion gyro-frequency, 
$B = |\boldsymbol {B}|$ and 
$m_{i}$ is the ion mass. As typically 
$\lambda _{D} \ll \rho _{s}$, the region 
$x \sim \rho _{s}$ can be assumed to be quasi-neutral,
\begin{equation} Zn_{i}(x) \simeq n_{e} (x), \end{equation}
and is referred to as the magnetic presheath (and sometimes as the Chodura sheath). A substantial fraction of the electrostatic potential drop between the plasma and the wall must occur in the magnetic presheath, as an electric field is necessary to adjust the electron and ion densities such that (1.3) is preserved. At typically even larger distances from the target, 
$d_{c} \gg \rho _{s}$, ions tend to collide with neutrals or other ions before reaching the target. Thus, the magnetic presheath and Debye sheath can be assumed to be collisionless. In this paper, the form of the ion distribution function in the region 
$\rho _{s} \ll x \ll d_{c}$ is assumed. This region is known as the magnetic presheath entrance.
Several distinct approaches may be used to calculate the velocity distributions of ions reaching the target. An approach that describes all the phenomena at play close to the wall, including the effect of the collisional layer, is to numerically solve the kinetic Vlasov equation for the ions and electrons self-consistently with the Poisson equation for the electrostatic potential (Coulette & Manfredi Reference Coulette and Manfredi2016). An alternative, equally complete, approach is the particle-in-cell (PIC) method (Tskhakaya & Kuhn Reference Tskhakaya and Kuhn2003; Khaziev & Curreli Reference Khaziev and Curreli2015). Both the Vlasov and the PIC approaches offer the most complete description of the plasma, but can be computationally expensive. Simplifying models can offer more immediate calculations. For example, taking into account gyro-orbit losses at the wall, but ignoring the electric field, one can solve for distribution functions at the wall analytically, assuming an incoming Maxwellian (Parks & Lippmann Reference Parks and Lippmann1994) or more refined boundary conditions (Gunn et al. Reference Gunn, Carpentier-Chouchana, Dejarnac, Escourbiac, Hirai, Komm, Kukushkin, Panayotis and Pitts2017). However, in neglecting the electric field this model assumes that some ions can reach the target travelling tangentially,Footnote 1 as the left ion in figure 1(a) does. By introducing an ad hoc analytical electrostatic potential function close to the wall to model the effect of gyro-orbit distortion, Borodkina et al. (Reference Borodkina, Borodin, Kirschner, Tsvetkov, Kurnaev, Komm, Dejarnac and Contributors2016) numerically solved for ion trajectories near the target. The authors found a substantial effect on erosion coefficients, as was also suggested by Siddiqui et al. (Reference Siddiqui, Thompson, Jackson, Kim, Hershkowitz and Scime2016). Daube & Riemann (Reference Daube and Riemann1999) obtained self-consistent solutions of the electrostatic potential and ion distribution function in a magnetic presheath by considering charge exchange collisions with cold neutrals. They calculated the ion density as an integral over characteristics originating at the last collision event. The resulting ion distribution functions exhibit an interesting and involved structure with singularities, which are expected to be smeared out by unstable ion cyclotron modes (Daube, Riemann & Schmitz Reference Daube, Riemann and Schmitz1998) and finite neutral temperature. Tskhakaya Sr & Kos (Reference Tskhakaya and Kos2014) analysed the plasma–wall boundary layers using an asymptotic scale separation and an asymptotic expansion in 
$\alpha \ll 1$. They considered the ion gyro-orbits to have zero spatial extent, but retained all other kinetic effects. In Geraldini, Parra & Militello (Reference Geraldini, Parra and Militello2017), the full approximately periodic ion trajectories in the collisionless magnetic presheath were solved using an expansion in 
$\alpha \ll 1$. This expansion leads to the presence of an adiabatic invariant, as first described by Cohen & Ryutov (Reference Cohen and Ryutov1998a). A numerical scheme to efficiently calculate the self-consistent electrostatic potential was developed by Geraldini, Parra & Militello (Reference Geraldini, Parra and Militello2018). The final open piece of the ion trajectory near the wall was included in the ion density calculation. Velocity distributions of ions reaching the Debye sheath, consistent with a quasi-neutral magnetic presheath, were thus obtained. Although this treatment applies only to grazing angles, it provides an efficient way to solve self-consistently for the effect of the electric field on ion trajectories in the collisionless magnetic presheath.
In this paper, a large gyro-orbit model for the ion distribution function at the target is developed. The full solution of the self-consistent electrostatic potential is bypassed. Instead, the electrostatic potential is assumed to distort ion gyro-orbits only just before ions reach the Debye sheath. This assumption is expected to be more accurate for large gyro-orbits, 
$\tau \gg 1$. The model results are compared with distribution functions obtained using the full self-consistent electrostatic potential solution in the magnetic presheath, with good qualitative agreement for 
$\tau \gtrsim 1$. The agreement between the two methods is better at larger values of 
$\tau$, as expected.
The rest of the paper is structured as follows. In § 2, the orderings assumed in this work are presented and discussed. In § 3, the electron model is introduced. In § 4 ion trajectories in the collisionless magnetic presheath and Debye sheath regions are analysed. Expressions for the velocity distributions of ions reaching the Debye sheath and of ions striking the target are obtained in § 5. These expressions depend on the full electrostatic potential solution in the magnetic presheath, 
$\phi (x)$. The trajectories of ions in large gyro-orbits, for 
$\tau \gg 1$, are analysed in § 6. From this analysis, a model for the ion velocity distribution at the target is developed. In § 7 ion distribution functions obtained from the large gyro-orbit model are compared with those obtained from the full self-consistent electrostatic potential solution 
$\phi (x)$ in the magnetic presheath. Finally, in § 8, the results of the paper are summarized.
2. Orderings
As mentioned in the introduction, the typical electrostatic potential variation across the magnetic presheath and Debye sheath is ordered as 
$|\phi | \sim T_{e} / e$. Hence, the kinetic energy transferred by the electric field to an ion of charge 
$Ze$ is 
$Z e \phi \sim ZT_{e}$ and the characteristic speed of an ion owing to the energy gained from the electric field is the Bohm velocity 
$v_{B} = \sqrt {ZT_{e} / m_{i}}$. The thermal energy of an ion is 
$T_{i}$ and the thermal speed of an ion is 
$v_{t,i} = \sqrt {2T_{i} / m_{i}}$. Adding together the contributions to the energy, the typical kinetic energy of an ion is 
$ZT_{e} + T_{i}$. The ion velocity, denoted by 
$\boldsymbol {v} = (v_x, v_y, v_z)$ where 
$v_k$ is the velocity component in the 
$k$th direction, is therefore ordered such that 
$|\boldsymbol {v}| \sim \sqrt {(ZT_{e} + T_{i} ) / m_{i}} = c_{s}$.
The presence of ion gyro-orbits and the grazing angle of the magnetic field with the target modify the ordering for 
$v_x$ at the target as follows. Consider a circular ion gyro-orbit with no electric field, as shown in figure 1(a). The component of the velocity parallel to the magnetic field is denoted by 
$v_{\parallel }$ and the magnitude of the gyrating component of the velocity is denoted by 
$v_{\perp }$. The gyro-phase angle of the ion is denoted by 
$\varphi$. In the small-angle approximation, 
$\sin \alpha \simeq \alpha$, 
$\cos \alpha \simeq 1$ and the component of the velocity normal to the wall is given by 
$v_x \simeq v_{\perp } \sin \varphi - \alpha v_{\parallel }$. If the gyro-orbit almost touches the wall (
$x \rightarrow 0$) tangentially at a time 
$t=0$, the distance from the wall at a later time 
$t$ is 
$x \simeq (v_{\perp } / \varOmega ) ( 1 - \cos \varphi ) - \alpha v_{\parallel } t$. After a full gyro-period 
$2{\rm \pi} / \varOmega$, the orbit has drifted a little closer to the wall. Therefore, the gyro-phase angle corresponding to 
$x=0$ is no longer 
$\varphi = 0$, yet it has only changed by a small amount. Solving for 
$x=0$ at 
$t = 2{\rm \pi} / \varOmega$ with 
$1 - \cos \varphi \simeq \varphi ^{2} / 2$ gives 
$\varphi \simeq \sqrt { 4 {\rm \pi}\alpha v_{\parallel } / v_{\perp } }$, and thus 
$v_x \simeq - \sqrt {4 {\rm \pi}\alpha v_{\parallel } v_{\perp } }$ (Cohen & Ryutov Reference Cohen and Ryutov1998a). The piece of 
$v_x$ equal to 
$-\alpha v_{\parallel }$ is smaller by a factor of 
$\sqrt { \alpha v_{\parallel } / ( 4{\rm \pi} v_{\perp } ) }$, and can be neglected. Thus, the gyro-phase dependence of ions reaching the target gives rise to an interval in allowed values of normal kinetic energy, 
$0 \leqslant v_x^{2} / 2 < 2{\rm \pi} \alpha v_{\parallel } v_{\perp }$. The electric field, however, can still accelerate the ions by transferring an energy 
$\sim ZT_{e}$ to the normal component of the velocity, as depicted schematically in figure 1(b,c). Note that this additional acceleration towards the target is not obvious. It only happens because, as we will see, the electric field close to the target is sufficiently inhomogeneous (
$|\phi ''(x)|$ is sufficiently large) that it overcomes the magnetic force pulling the ion back away from the target. Combining these two contributions to the normal kinetic energy gives 
$v_x^{2} / 2 \sim ZT_{e} + \alpha T_{i}$. The velocity of the ion at the target therefore satisfies 
$v_x \sim v_{B} \sqrt { 1+\alpha \tau }$ and 
$v_y \sim v_z \sim c_{s}$.
As was discussed in the introduction, the Debye sheath, the magnetic presheath and the collisional region are assumed to satisfy the scale separation 
$\lambda _{D} \ll \rho _{s} \ll d_{c}$. At distances 
$x \sim d_{c} \gg \rho _{s}$, the ion motion is restricted along a field line. Therefore, the size of the collisional region can be expressed as 
$d_{c} \sim \alpha \lambda _{\textrm {mfp}}$, where 
$\lambda _{\textrm {mfp}}$ is the mean free path of an ion near the target. It follows that the angle 
$\alpha$ must satisfy 
$\alpha \gg \rho _{s} / \lambda _{\textrm {mfp}}$ in order for 
$\rho _{s} \ll d_{c}$ to be valid.
In order to simplify the treatment of the electrons, the electron gyro-radius 
$\rho _{e} = \sqrt {2 m_{e} T_{e}} / (eB)$ is assumed to be much smaller than the Debye length, such that 
$\rho _{e} \ll \lambda _{D}$ (Loizu et al. Reference Loizu, Ricci, Halpern and Jolliet2012; Stangeby Reference Stangeby2012). Being tightly bound to the magnetic field lines, electrons have to travel along the magnetic field in order to reach the wall. The typical speed of an electron is the electron thermal speed, 
$v_{t,e} = \sqrt {2T_{e}/m_{e}}$. Conversely, the typical ion velocity close to the wall is 
${\sim }v_{B} \sqrt { 1+\alpha \tau }$ towards the wall. When unopposed by an electric field, the electrons reach the wall much more quickly than the ions provided that 
$\alpha v_{t,e} \gg v_{B} \sqrt { 1+\alpha \tau }$, or 
$\sqrt { (1+\alpha \tau ) Zm_{e} / m_{i} } \ll \alpha$. For 
$\alpha \tau \lesssim 1$, the ordering 
$\alpha \gg \sqrt { Z m_{e} / m_{i} }$ emerges. For 
$\alpha \tau \gg 1$, the ordering 
$\alpha \gg m_{e} \tau Z / m_{i}$ emerges instead. Putting these last two orderings together gives 
$1/ \alpha \ll \tau \ll \alpha m_{i} / m_{e} Z$, which can only be satisfied if, again, 
$\alpha \gg \sqrt { Z m_{e} / m_{i} }$. To summarize, for 
$\alpha \gg \sqrt { Z m_{e} / m_{i} }$ the electrons reach the target much more quickly than the ions. An electric field must therefore set up to repel most of the electrons from the target.
Summarizing the orderings of this work, the physical length scales satisfy
\begin{equation} \rho_{e} \ll \lambda_{D} \ll \rho_{s} \ll d_{c}. \end{equation}The angle and mass ratio satisfy
\begin{equation} \sqrt{\frac{ Zm_{e}}{m_{i}} } \ll \alpha \ll 1. \end{equation}
The validity of these orderings is examined for a current fusion experiment such as JET. In a deuterium plasma, the angle obtained from the square root of mass ratio is 
$\sqrt {Zm_e/m_i} \approx 0.02\ \textrm {rad} \sim 1^{\circ }$. From Militello & Fundamenski (Reference Militello and Fundamenski2011), we estimate for JET: 
$B \sim 2\ \textrm {T}$, 
$T_{e} \sim T_{i} \sim 30\ \textrm {eV}$, 
$n_{e} \sim n_{i} \sim 10^{19} \ \textrm {m}^{-3}$, giving 
$\rho _{s} \sim 1\ \textrm {mm}$, 
$\lambda _{D} \sim \rho _{e} \sim 0.01\ \textrm {mm}$ and 
$\alpha \approx 0.07\ \text {rad} \approx 4^{\circ }$. As, of all the orderings in this paper, 
$\sqrt {Zm_{e} / m_{i} } \ll \alpha$ and 
$\rho _{e} \ll \lambda _{D}$ are the least well-satisfied in fusion devices, it will be necessary to study in more detail the effect of electron inertia and gyro-radius.
3. Electron model
In this work, Maxwellian electrons are assumed to enter the magnetic presheath. We proceed to obtain the relationship between the electron current to the wall and the electrostatic potential at the wall. We also derive, using the ordering (2.2), the Boltzmann expression for the electron density in the magnetic presheath.
According to (2.1), the electron gyro-radius is so small that electrons are essentially tied to the magnetic field line, as shown in figure 2. The electrons stream parallel to the magnetic field with a velocity given by 
$w_{\parallel }$. At the very small length scale 
$\rho _{e} \ll \lambda _{D}$, the electron gyro-motion is unaffected. The electron distribution function entering (that is, for 
$w_{\parallel } > 0$) the magnetic presheath is assumed to be a half-Maxwellian,
\begin{equation} g_{\textrm{MPE}} ( w_{{\parallel}} ) = Z\bar{n}_{\textrm{MPE}} \left( \frac{ m_{e}}{2{\rm \pi} T_e} \right)^{1/2} \exp \left(- \frac{m_e w_{{\parallel}}^{2}}{2T_e} \right),\quad \text{for}\ w_{{\parallel}} > 0 , \end{equation}
with density denoted as 
$Zn_{\textrm {MPE}}$,
\begin{equation} Zn_{\textrm{MPE}}= \int_{-\infty}^{\infty} g_{\textrm{MPE}} ( w_{{\parallel}} )\,\textrm{d}w_{{\parallel}}. \end{equation}
We set the zero of the electrostatic potential to be at the magnetic presheath, 
$\phi _{\textrm {MPE}} = 0$. Assuming the electrostatic potential to be a monotonically increasing function of 
$x$, the number of electrons that enter the magnetic presheath and come back out of it depends on the electrostatic potential at the wall relative to the magnetic presheath entrance, denoted by 
$\phi _{W} = \phi (0) < 0$. Therefore, the constant 
$\bar {n}_{\textrm {MPE}}$ depends on 
$n_{\textrm {MPE}}$ and 
$\phi _{W}$.
An electron gyro-orbit, whose gyro-radius is 
$\rho _{e}$, streaming towards the wall along the magnetic field 
$\boldsymbol {B}$ with velocity 
$w_{\parallel }$.
In the magnetic presheath and Debye sheath, the component of the electron velocity parallel to the magnetic field as a function of 
$x$ is obtained by energy conservation,
\begin{equation} w_{{\parallel}} = \sigma \sqrt{w_{{\parallel} \textrm{MPE}}^{2} + \frac{2e\phi(x)}{m_{e}}} . \end{equation}
Here, 
$w_{\parallel \textrm {MPE}}$ is the electron velocity at the magnetic presheath entrance. The 
$\boldsymbol {E}\times \boldsymbol {B}$ and gyration velocities of an electron remain unaffected by electrostatic potential variations as these have a much longer scale length than the electron gyro-radius, 
$\lambda _{D} \gg \rho _{e}$. In (3.3), 
$\sigma = {\pm }1$ for those electrons reflected before reaching the wall and 
$\sigma = 1$ for those electrons that are not reflected. At 
$x = 0$ the electron velocity is zero if 
$w_{\parallel \textrm {MPE}}^{2} = - 2e\phi _{W} / m_{e}$. Hence, reflected electrons satisfy
\begin{equation} w_{{\parallel} \textrm{MPE}}^{2} <{-} \frac{2e \phi_{W}}{m_{e}} , \end{equation}
as they cannot reach 
$x=0$. Therefore, the full electron distribution function at the magnetic presheath entrance is
\begin{equation} g_{\textrm{MPE}} ( w_{{\parallel}} ) = Z\bar{n}_{\textrm{MPE}} \left( \frac{m_{e}}{2{\rm \pi} T_e} \right)^{1/2} \exp \left(- \frac{m_e w_{{\parallel}}^{2}}{2T_e} \right) \varTheta \left( w_{{\parallel}} + \sqrt{-\frac{2e \phi_{W}}{m_{e}}} \right) , \end{equation}
where 
$\varTheta$ is the Heaviside step function,
\begin{equation} \varTheta ( \xi ) = \begin{cases} 1, & \text{for}\ \xi \geqslant 0 , \\ 0, & \text{for}\ \xi < 0 . \end{cases} \end{equation}
Assuming 
$\text {erf}(\sqrt {-e\phi _{W}/T_{e}}) \simeq 1$, which will be justified in the next paragraph, we obtain
\begin{equation} \bar{n}_{\textrm{MPE}} = \frac{2n_{\textrm{MPE}} }{\left( 1+\text{erf}\left( \sqrt{- e\phi_{W} / T_{e} } \right) \right)} \simeq n_{\textrm{MPE}} . \end{equation} The electron current 
$j_{e\parallel }$ is obtained from the first moment of the distribution function (3.5) (the flux of electrons) multiplied by the electron charge, 
$-e$. The current directed towards the wall is the geometric projection of the parallel current, 
$j_{e, x} = - j_{e\parallel } \sin \alpha \simeq - \alpha j_{e\parallel }$,
\begin{equation} j_{e, x} \simeq \alpha Z e n_{\textrm{MPE}} \left( \frac{T_e}{2{\rm \pi} m_e} \right)^{1/2} \exp \left( \frac{e\phi_{W}}{T_{e}} \right) . \end{equation}
As the electron charge is negative and the electron flow is directed towards the wall (negative), the electron current is directed away from the wall (positive). The electron and ion current are assumed to be similar in size. To be consistent with the Chodura condition (Chodura Reference Chodura1982) at the magnetic presheath entrance, the ion current is assumed to be of the order of the sound speed, giving 
$j_{ e, x} \sim \alpha Z en_{\textrm {MPE}} c_{s}$. Hence, the electrostatic potential at the wall is
\begin{equation} \frac{e\phi_{W} }{ T_{e} } \sim \ln \left( \sqrt{ \frac{2{\rm \pi} m_{e} (1+\tau) }{ m_{i} } } \right), \end{equation}
where 
$\sqrt {2{\rm \pi} m_{e} (1+\tau ) / m_{i} } \ll 1$, justifying 
$\text {erf}(\sqrt {-e\phi _{W}/T_{e}}) \simeq 1$.
The electron distribution function at any point in the magnetic presheath and Debye sheath is (Stangeby Reference Stangeby2012)
\begin{align} g ( x, w_{{\parallel}} ) \simeq Z n_{\textrm{MPE}} \left( \frac{m_{e}}{2{\rm \pi} T_e} \right)^{1/2} \exp \left(\frac{e \phi (x) }{T_{e}} - \frac{ m_e w_{{\parallel}}^{2} }{2T_e} \right) \varTheta \left( w_{{\parallel}} + \sqrt{\frac{2e ( \phi (x) - \phi_{W} ) }{m_{e}}} \right) . \end{align}Hence, the electron density is
\begin{equation} n_{e} ( x ) \simeq \frac{1}{2} \left( 1 + \text{erf} \left( \sqrt{\frac{e ( \phi (x) - \phi_{W} ) }{T_{e}}} \right) \right) Z n_{\textrm{MPE}} \exp \left(\frac{e \phi (x) }{T_{e}} \right) . \end{equation}
In the magnetic presheath the electrostatic potential is at its smallest at the Debye sheath entrance, 
$\lambda _{D} \ll x \ll \rho _{s}$, where 
$\phi (x) \simeq \phi _{\textrm {DSE}}$. Thus, provided 
$\text {erf} ( \sqrt {e ( \phi _{\textrm {DSE}} - \phi _{W} ) / T_{e}} ) \simeq 1$, the electron density in the magnetic presheath is given by the Boltzmann distribution
\begin{equation} n_{e} ( x ) \simeq Z n_{\textrm{MPE}} \exp \left(\frac{e \phi (x) }{T_{e}} \right) . \end{equation} We proceed to justify (3.12). The ion flow speed parallel to the magnetic field at the magnetic presheath entrance, 
$\rho _{s} \ll x \ll d_{c}$, is of the order of the sound speed 
$\sim c_{s}$. Projecting this parallel flow in the direction normal to the target gives 
$\alpha c_{s} \sim \alpha \sqrt {1+\tau } v_{B}$. The ion velocity component perpendicular to the magnetic field averages to zero at the magnetic presheath entrance, as the electric field is small and the target is too far away to capture ions during their gyro-motion. Conversely, at the Debye sheath entrance the size of the ion flow is determined by the ordering for the velocity component normal to the target, 
$v_x \sim \sqrt { 1 + \alpha \tau } v_{B}$. As the number of ions in the magnetic presheath is conserved in steady state, the ion flux into the magnetic presheath, 
$\alpha n_{\textrm {MPE}} \sqrt {1+\tau } v_{B}$, and the ion flux out of the magnetic presheath, 
$n_{\textrm {DSE}} \sqrt {1+\alpha \tau } v_{B}$, are equal. The ion density at the Debye sheath entrance is thus 
$n_{\textrm {DSE}} \sim \alpha n_{\textrm {MPE}} \sqrt {1+\tau } / \sqrt {1+\alpha \tau }$. Hence, we find
\begin{equation} \frac{ e \phi_{\textrm{DSE}} }{ T_{e} } \sim \ln \left( \frac{ \alpha \sqrt{1+\tau} }{ \sqrt{1+\alpha \tau } } \right) \end{equation}and
\begin{equation} \frac{ e ( \phi_{W} - \phi_{\textrm{DSE}} ) }{ T_{e} } \sim \ln \left( \frac{1}{\alpha} \sqrt{\frac{2{\rm \pi} m_{e} (1+\alpha \tau)}{m_{i}} } \right). \end{equation}
Upon neglecting the factors of 
$\alpha \tau$, the estimates in (3.9), (3.13) and (3.14) are consistent with those in Stangeby (Reference Stangeby2012). Equation (3.12) follows from expanding (3.11), with 
$\phi (x) \geqslant \phi _{\textrm {DSE}}$, using the orderings (3.14) and 
$\sqrt {Zm_{e} / m_{i}} \ll \alpha$. Note that (3.9), (3.13) and (3.14) are all negative, with the arguments of the logarithm smaller than unity.
The results of this section that will be used in the rest of the paper are (3.12) for the electron density in the magnetic presheath and (3.8) for the relationship between electron current and wall potential.
4. Ion trajectories
In this section the trajectories of ions in the magnetic presheath and the Debye sheath are analysed in detail. The goal of this section is to relate the velocity of an ion at the target to the energy and magnetic moment of its circular gyro-orbit at the magnetic presheath entrance 
$\rho _{s} \ll x \ll d_{c}$. We analyse the ion trajectories first in the magnetic presheath, § 4.1, and then in the Debye sheath, § 4.2.
4.1. In the magnetic presheath
We proceed to focus on the magnetic presheath, where 
$x\sim \rho _{s}$. Ions move under the influence of a wall-normal electrostatic electric field and a magnetic field at an angle 
$\alpha$ with the wall. The ion equations of motion are
\begin{gather} \dot{v}_x ={-} \frac{\varOmega \phi'(x)}{B} + \varOmega v_y \cos\alpha , \end{gather}
\begin{gather} \dot{v}_y ={-} \varOmega v_x \cos\alpha - \varOmega v_z \sin \alpha, \end{gather}
\begin{gather} \dot{v}_z = \varOmega v_y \sin \alpha. \end{gather}
For grazing angles, 
$\alpha \ll 1$, the equations simplify to
\begin{gather} \dot{v}_x \simeq{-} \frac{\varOmega \phi'(x)}{B} + \varOmega v_y, \end{gather}
\begin{gather} \dot{v}_y \simeq{-} \varOmega v_x - \alpha \varOmega v_z , \end{gather}
\begin{gather} \dot{v}_z \simeq \alpha \varOmega v_y, \end{gather}
where only small terms linear in 
$\alpha$ were retained. It will be useful to introduce two orbit parameters,
\begin{gather} \bar{x} = x + \frac{v_y}{\varOmega}, \end{gather}
\begin{gather} U_{{\perp}} = \frac{1}{2} v_x^{2} + \frac{1}{2} v_y^{2} + \frac{\varOmega \phi(x)}{B}, \end{gather}
whose time derivatives satisfy 
$\dot {\bar {x}} \simeq - \alpha v_z$ and 
$\dot {U}_{\perp } \simeq - \alpha \varOmega v_y v_z$. The third orbit parameter,
\begin{equation} U = \frac{1}{2} v_x^{2} + \frac{1}{2} v_y^{2} + \frac{1}{2} v_z^{2} + \frac{\varOmega \phi(x)}{B} , \end{equation}
is just the total energy of an ion and is exactly conserved, 
$\dot {U} = 0$. From the definitions (4.7)–(4.9), we obtain
\begin{gather} v_z = \sqrt{2\left( U - U_{{\perp}} \right) } , \end{gather}
\begin{gather} v_y = \varOmega (\bar{x} - x) , \end{gather}and
\begin{equation} v_x = {\pm} \sqrt{2\left( U_{{\perp}} - \chi (x, \bar{x}) \right) } . \end{equation}In (4.12) an effective potential function,
\begin{equation} \chi(x, \bar{x}) = \frac{1}{2} \varOmega^{2} \left( x - \bar{x} \right)^{2} + \frac{\varOmega \phi(x)}{B} , \end{equation}
was introduced. Note that, to lowest order in 
$\alpha \ll 1$, 
$v_z$ is equivalent to the velocity component parallel to the magnetic field. The electric field slowly (owing to the grazing angle) pushes ions in the direction parallel to the magnetic field towards larger 
$v_z$ (Geraldini et al. Reference Geraldini, Parra and Militello2017). All ions enter the magnetic presheath with a parallel velocity directed towards the target, and so they have 
$v_z \geqslant 0$ to lowest order in 
$\alpha$. As the parallel velocity towards the wall increases in the magnetic presheath, ions with 
$v_z < 0$ are not present. Therefore, in (4.10) we have set 
$v_z \geqslant 0$.
The orbit parameter 
$\bar {x}$ is referred to as the orbit position, and 
$U_{\perp }$ as the perpendicular energy (perpendicular to the magnetic field). As 
$\dot {\bar {x}}/ \rho _{{s}} \sim \dot {U}_{\perp } / c_{{s}}^{2} \sim \alpha \varOmega \ll \varOmega$, the orbit position and perpendicular energy only change by a very small amount during the timescale 
${\sim } 1/\varOmega$. Neglecting the small change in the orbit parameters (which is a good approximation for a time 
$\ll 1/(\alpha \varOmega )$), particle orbits are solved for as follows. Consider a stationary point of the effective potential, 
$\chi _{\text {st}}(\bar {x}) = \chi (x_{\text {st}}, \bar {x})$, such that 
$\chi '(x_{\text {st}}, \bar {x}) = \varOmega ^{2} (x_{\textrm {st}} - \bar {x} ) + \varOmega \phi ' (x_{\textrm {st}}) / B=0$. Here, it is understood that 
$\chi '(x, \bar {x}) = \partial \chi (x, \bar {x}) / \partial x$. Rearranging this equation gives the orbit parameter as a function of the position of a stationary point,
\begin{equation} \bar{x} = x_{\text{st}} + \frac{ \phi'(x_{\textrm{st}})}{\varOmega B} . \end{equation}
A stationary point is a minimum, 
$x_{\textrm {st}} = x_{m}$, if 
$\chi ''(x_{{m}}, \bar {x}) = \varOmega ^{2} + \varOmega \phi '' (x_{m}) / B > 0$, leading to
\begin{equation} \phi''(x_{m}) >{-}\varOmega B . \end{equation}
At the magnetic presheath entrance, the electrostatic potential is assumed to monotonically converge to the value 
$\phi _{\textrm {MPE}} = 0$. We further assume that 
$\phi ''(x)$ is negative (the magnitude of the electric field, 
$\phi '(x)$, decreases away from the wall) and monotonically converges to zero at the magnetic presheath entrance. Hence, the stationary point is a minimum for 
$x_{\textrm {st}} > x_{c}$, where 
$\phi ''(x_{c}) = -\varOmega B$ if 
$\phi ''_{\textrm {DSE}} \leqslant - \varOmega B$ or 
$\lambda _{D} \ll x_{c} \ll \rho _{s}$ if 
$\phi ''_{\textrm {DSE}} > - \varOmega B$. Here 
$\phi ''_{\textrm {DSE}}$ denotes 
$\phi ''(x)$ at the Debye sheath entrance, 
$\lambda _{D} \ll x \ll \rho _{s}$, and 
$x_{c}$ is a critical point corresponding to the inflection point of 
$\chi$, if it exists, or the Debye sheath entrance 
$\lambda _{D} \ll x_{c} \ll \rho _{s}$. There are either two or one solutions for stationary points of the effective potential according to (4.14), depending on whether the function 
$x + \phi '(x) / (\varOmega B)$ has a stationary point or not. This leads to the distinction between two orbit types in the magnetic presheath. Type I orbits occur when the effective potential 
$\chi (x, \bar {x})$ has only one stationary point: a minimum 
$x_{m}$. Type II orbits occur when 
$\chi (x, \bar {x})$ has two stationary points: a minimum 
$x_{m}$ and a maximum 
$x_{M} < x_{m}$. For 
$\bar {x} > \phi '_{\textrm {DSE}} / (\varOmega B)$, where 
$\phi '_{\textrm {DSE}}$ denotes 
$\phi '(x)$ at the Debye sheath entrance, there is only one solution to (4.14) in the magnetic presheath and therefore there are only type I ion orbits. For both type I and type II orbits, the motion is periodic in the neighbourhood of the minimum. The turning points 
$x_{b}$ (for ‘bottom’) and 
$x_{t}$ (for ‘top’) of the periodic motion satisfy 
$x_{M} \leqslant x_{b} < x_{m}$ and 
$x_{t} > x_{m}$. They are obtained by solving for the positions at which 
$v_x = 0$, that is, 
$U_{\perp } = \chi ( x_{b, t} , \bar {x} )$.
The slow change in 
$\bar {x}$ and 
$U_{\perp }$ cannot be neglected entirely, as it leads to ions eventually reaching the wall. Ion trajectories are approximately periodic over a short timescale, 
${\sim }1/\varOmega$. Over a long enough timescale, 
${\sim }1/(\alpha \varOmega )$, the effect of the slow variation in 
$\bar {x}$ and 
$U_{\perp }$ becomes significant. Nonetheless, the quasi-periodic motion of the ion has an adiabatic invariant
\begin{equation} \mu = \frac{1}{\rm \pi} \int_{x_{b}}^{x_{t}} \sqrt{2\left(U_{{\perp}} - \chi(x, \bar{x}) \right)} \,{\textrm{d}x}, \end{equation}
which is conserved to lowest order in 
$\alpha \ll 1$ during the entire ion trajectory in the magnetic presheath (Cohen & Ryutov Reference Cohen and Ryutov1998a; Geraldini et al. Reference Geraldini, Parra and Militello2017). At the magnetic presheath entrance, 
$\rho _{s} \ll x \ll d_{c}$, 
$\phi (x) = 0$ and so the adiabatic invariant of (4.16) is given by 
$\mu = (1 / {\rm \pi}) \int _{x_{b}}^{x_{t}}\,\textrm {d}s \sqrt {2 U_{\perp } - \varOmega ^{2} (s - \bar {x} )^{2} }$ with 
$x_{b} = \bar {x} - \sqrt {2U_{\perp } }/ \varOmega$ and 
$x_{t} = \bar {x} + \sqrt {2U_{\perp } }/ \varOmega$. Upon changing variables to 
$\varphi$ using 
$s = \bar {x} - ( \sqrt {2U_{\perp }} / \varOmega ) \cos \varphi$, the adiabatic invariant becomes 
$\mu = ( 2U_{\perp } / ({\rm \pi} \varOmega ) ) \int _{0}^{{\rm \pi} } \,\textrm {d}\varphi \sin ^{2} \varphi = U_{\perp } / \varOmega$. Using this result and (4.8) for 
$U_{\perp }$, with 
$\phi (x) = 0$, we obtain 
$\mu = ( v_x^{2} + v_y^{2} )/ ( 2 \varOmega )$. This is equivalent to the magnetic moment to lowest order in 
$\alpha \ll 1$; the small difference is geometric and arises because 
$v_x$ is not exactly perpendicular to the magnetic field.
The ion motion can be described as approximately periodic only insofar as it is not about to be interrupted by the absorbing wall. If the perpendicular energy becomes larger than a threshold value, the ion gyro-orbit becomes sufficiently large that the bottom bounce point disappears. The threshold value of 
$U_{\perp }$ is the maximum value of the effective potential function between the position of the minimum, 
$x=x_{m}$, and the wall, 
$x=0$,
\begin{equation} \chi_{{M}}(\bar{x}) \equiv \chi (x_{{M}}, \bar{x}) = \max_{x\in[0,x_{{m}}]} \chi (x, \bar{x}) . \end{equation}
For type I orbits, the effective potential maximum lies at the Debye sheath entrance 
$\lambda _{D} \ll x_{M} \ll \rho _{s}$, such that 
$\chi _{M} (\bar {x}) \simeq \varOmega ^{2} \bar {x}^{2} / 2 + \varOmega \phi _{\textrm {DSE}} / B$. For type II orbits, the effective potential maximum lies in the magnetic presheath 
$x_{M} \sim \rho _{s}$, such that 
$\chi _{M} (\bar {x}) = \varOmega ^{2} (x_{M}-\bar {x})^{2} / 2 + \varOmega \phi (x_{M}) / B$. In this case, 
$x_{M}$ is a stationary point. As the variation of 
$U_{\perp }$ and 
$\bar {x}$ is slow compared with the timescale of ion motion, ions quickly reach the wall once 
$U_{\perp } > \chi _{M} (\bar {x})$, and therefore these ions have 
$U_{\perp } \simeq \chi _{M} (\bar {x})$. Any ion reaching the wall must, because it comes from an approximately periodic orbit, have a value of orbit position such that an effective potential minimum exists. From (4.14), the smallest value of orbit position, denoted by 
$\bar {x}_{c}$, for ions in the magnetic presheath is
\begin{equation} \bar{x}_{c} = \min \left( x + \frac{\phi'(x)}{\varOmega B} \right) = x_{c} + \frac{ \phi'(x_{c})}{ \varOmega B } . \end{equation}
Note that the second equality defines the value of 
$x_{c}$, which is consistent with the discussion after (4.15) where 
$x_{c}$ is first introduced.
4.2. In the Debye sheath
Here, we focus on ions in the Debye sheath, 
$x \sim \lambda _{D} \ll \rho _{s}$. Considering 
$\bar {x} \sim \rho _{s}$ and neglecting 
$x \ll \rho _{s}$ in (4.11) gives
\begin{equation} v_y \simeq \varOmega \bar{x}. \end{equation}
For every ion in the Debye sheath, we can trace back its trajectory to a quasi-periodic orbit. The associated value of 
$\mu$ is a function of 
$\bar {x} ( \simeq v_y / \varOmega )$ only, because 
$U_{\perp } \simeq \chi _{M}(\bar {x})$ for ions reaching the target,
\begin{equation} \mu_{\textrm{op}} (\bar{x}) = \frac{1}{\rm \pi} \int_{x_{{M}}}^{x_{{t}}} \sqrt{2\left(\chi_{{M}} (\bar{x}) - \chi(x, \bar{x}) \right)} \,{\textrm{d}x}. \end{equation}
Here we have used 
$x_{b} = x_{M}$ for 
$U_{\perp } = \chi _{M} (\bar {x})$. The value of 
$v_z$ is determined by the total energy 
$U$,
\begin{equation} v_z \simeq \sqrt{2\left( U - \chi_{M} (\bar{x} )\right)}. \end{equation} In order to calculate 
$v_x$ in the Debye sheath, the final piece of the ion trajectory in the magnetic presheath must be considered. This is a transition from a quasi-periodic orbit, with at least one turning point in its future trajectory, to an open orbit, with no turning points in its future trajectory. The small change of 
$\bar {x}$ and 
$U_{\perp }$ causes the value of 
$U_{\perp } - \chi _{{M}} (\bar {x})$ to increase until 
$U_{\perp } > \chi _{M} (\bar {x})$. The increase is slow and so the change in 
$U_{\perp } - \chi _{M} (\bar {x})$ incurred by an ion transitioning from 
$U_{\perp } < \chi _{M} (\bar {x})$ to 
$U_{\perp } > \chi _{M} (\bar {x})$ can be calculated approximately by assuming a periodic orbit with fixed 
$U_{\perp } = \chi _{M} (\bar {x})$, as shown in appendix A. Such an orbit is fictitious: it has a bottom turning point coinciding with the position of the effective potential maximum, 
$x_{M}$, and for a type II orbit it takes an infinite time to turn around at 
$x_{M}$. The true orbit turns at 
$x_{b} > x_{M}$ (with 
$U_{\perp } < \chi _{M}$), then once more at 
$x_{t}$ and then passes 
$x_{M}$ (with 
$U_{\perp } > \chi _{M}$) in a finite time 
$\sim \ln (1/\alpha ) / \varOmega$ moving towards the wall. Yet, despite the approximate orbit being qualitatively different from the true orbit, the change in 
$U_{\perp } - \chi _{M} (\bar {x})$ is accurate to lowest order in 
$\alpha$ when calculated from the approximate orbit. This is because the long time spent near 
$x_{M}$ does not contribute to a significant change in 
$U_{\perp } - \chi _{M} (\bar {x})$, as the time derivatives of 
$U_{\perp }$ and of 
$\chi _{M} (\bar {x})$ coincide at 
$x=x_{M}$. The overall change in the quantity 
$U_{\perp } - \chi _{M}(\bar {x})$ during the last gyro-orbit is
\begin{equation} \varDelta_{M}(\bar{x}, U) = 2{\rm \pi} \alpha V_{{\parallel}} \left( \chi_{{M}}(\bar{x}), U \right) \mu_{\textrm{op}}'(\bar{x}) , \end{equation}
where 
$\mu _{\textrm {op}}'(\bar {x})= \textrm {d}\mu _{\textrm {op}}(\bar {x})/\textrm {d}\bar {x}$. Equation (4.22) is derived in appendix A.
The implication of this discussion for ion trajectories in the Debye sheath is that there is a band of possible values of 
$v_x$ for a given value of 
$\bar {x}$ (or 
$\mu$) and 
$U$. Considering 
$v_x^{2} \simeq 2U_{\perp } - \varOmega ^{2} \bar {x}^{2} - 2\varOmega \phi (x) / B$, which follows from (4.12), (4.13) and 
$x\sim \lambda _{D} \ll \rho _{s}$, we obtain the range
\begin{equation} \chi_{M}(\bar{x}) - \frac{1}{2} \varOmega^2 \bar{x}^{2} - \frac{\varOmega \phi (x) }{B} \leqslant \frac{v_x^{2}}{2} < \chi_{M}(\bar{x}) + \varDelta_{M}(\bar{x},U) - \frac{1}{2} \varOmega^2 \bar{x}^{2} - \frac{\varOmega \phi (x) }{B} . \end{equation}
Equation (4.23) is valid at any point in the Debye sheath, including the Debye sheath entrance and the target. For 
$\sqrt {Zm_{e} / m_{i}} \ll 1$ the Debye sheath repels most electrons from the wall and attracts all ions to the wall, so ions in the Debye sheath must have 
$v_x < 0$.
5. Ion velocity distribution
The ion distribution function at the magnetic presheath entrance, 
$\rho _{s} \ll x \ll d_{c}$, is denoted by 
$f_{\textrm {MPE}}(v_x, v_y, v_z )$. The exact distribution function in this region includes a small number of ions with 
$v_z<0$, which are travelling out of the magnetic presheath towards the collisional presheath. However, to lowest order in 
$\rho _{s} \ll d_{c}$ there are no such ions,
\begin{equation} f_{\textrm{MPE}} (v_z < 0) = 0 . \end{equation}
It can be shown that the distribution function is independent of the gyro-phase angle (Cohen & Ryutov Reference Cohen and Ryutov1998a; Geraldini et al. Reference Geraldini, Parra and Militello2017) and therefore can be expressed in the form 
$F(\mu , U )$. The relationship between 
$f_{\textrm {MPE}}$ and 
$F$ is obtained by recalling that 
$\mu = ( v_x^{2} + v_y^{2} )/ ( 2 \varOmega )$ at the magnetic presheath entrance,
\begin{equation} f_{\textrm{MPE}} (v_x, v_y, v_z ) = F\left( \frac{ v_x^{2} + v_y^{2} }{ 2 \varOmega } , \frac{ v_x^{2} + v_y^{2} + v_z^{2}}{ 2 } \right). \end{equation}
The function 
$F(\mu , U )$ is conserved across the magnetic presheath to lowest order in 
$\alpha \ll 1$, because 
$\mu$ and 
$U$ are conserved.
The ion density at the magnetic presheath entrance, denoted by 
$n_{\textrm {MPE}}$, is
\begin{equation} n_{\textrm{MPE}} = 2{\rm \pi} \int_0^{\infty} \varOmega\,\textrm{d}\mu \int_{\varOmega \mu}^{\infty} \frac{F(\mu, U)\, \textrm{d}U}{\sqrt{2\left( U - \varOmega \mu \right) }} = \int f_{\textrm{MPE}} (v_x, v_y, v_z )\,\textrm{d}^{3} v . \end{equation}
The ion current towards the wall, 
$j_{{i},x}$, is obtained from the projection of the flow in the direction parallel to the magnetic field. For 
$\alpha \ll 1$, this is approximately equal to
\begin{equation} \frac{j_{{i},x}}{Ze} \simeq{-} 2{\rm \pi} \alpha \int_0^{\infty} \varOmega \,\textrm{d}\mu \int_{\varOmega \mu}^{\infty} F(\mu, U)\,\textrm{d}U ={-} \alpha \int f_{\textrm{MPE}} (v_x, v_y, v_z ) v_z \,\textrm{d}^{3} v. \end{equation}We define the total current normal to the wall as
\begin{equation} j_{x}= j_{{e},x} + j_{{i}, x}. \end{equation}From (3.8) and (5.5), the electrostatic potential at the wall is
\begin{equation} \exp \left( \frac{e\phi_{W}}{T_{e}} \right) \simeq \frac{ - j_{{i}, x} + j_{x} }{ \alpha Z e n_{\textrm{MPE}} } \sqrt{ \frac{2{\rm \pi} m_e}{T_e} }. \end{equation}The ion current is determined by (5.4), which leads to
\begin{equation} \frac{e\phi_{W} }{ T_{e} } \simeq \ln \left[ \sqrt{\frac{2{\rm \pi} m_e}{T_e}} \left( \frac{1}{n_{\textrm{MPE}}} 2{\rm \pi} \int_0^{\infty} \varOmega\,\textrm{d}\mu \int_{\varOmega\mu}^{\infty}\,\textrm{d}U F \left(\mu, U\right) + \frac{j_x }{\alpha Z en_{\textrm{MPE}}} \right) \right]. \end{equation}
The numerical results of this paper, presented in § 7, are obtained assuming ambipolarity, 
$j_x = 0$.
As was shown in § 4, every value of 
$\mu$ and 
$U$, originally associated with a circular gyro-orbit entering the magnetic presheath, is associated with a specific value of 
$v_y \simeq \varOmega \bar {x}$ and 
$v_z \simeq \sqrt {2( U - \chi _{M} (\bar {x}) ) }$ at the Debye sheath entrance, where 
$\mu = \mu _{\textrm {op}}(\bar {x})$. Here, 
$v_x$ is given by (4.23) with 
$\phi (x) = \phi _{\textrm {DSE}}$. Conservation of the phase space distribution function 
$F(\mu , U)$ leads to the following velocity distribution (Geraldini et al. Reference Geraldini, Parra and Militello2018),
\begin{align} f_{\textrm{DSE}} (v_x, v_y, v_z) &\simeq F \left( \mu_{\textrm{op}} ( \bar{x}) , U \right) \varTheta \left( \bar{x} - \bar{x}_{{c}} \right) \varTheta \left({-}v_x \right) \nonumber\\ &\quad \times \hat{\varPi} \left( \frac{1}{2}v_x^{2} - \chi_{M} ( \bar{x} ) + \frac{1}{2} \varOmega^{2} \bar{x}^{2} + \frac{\varOmega \phi_{\textrm{DSE}}}{B} , 0 , \varDelta_{M}(\bar{x}, U ) \right) . \end{align}Here, we have defined the top-hat function
\begin{equation} \hat{\varPi} ( \xi, \xi_1, \xi_2 ) = \begin{cases} 1, & \text{for } \xi_1 \leqslant \xi < \xi_2 ,\\ 0, & \text{otherwise}. \end{cases} \end{equation}
In appendix B it is shown that the ion current normal to the wall calculated from (5.8) is equal to (5.4), and thus (5.8) satisfies ion conservation. At the wall, where 
$x=0$, the range of possible values of 
$v_x$ associated with each value of 
$\bar {x}$ and 
$U$ is given by (4.23) with 
$\phi (0) = \phi _{W}$,
\begin{align} f_{W} (v_x, v_y, v_z) &\simeq F \left( \mu_{\textrm{op}} ( \bar{x}) , U \right) \varTheta \left( \bar{x} - \bar{x}_{{c}} \right) \varTheta \left({-}v_x \right) \nonumber\\ &\quad \times \hat{\varPi} \left( \frac{1}{2}v_x^{2} - \chi_{M} ( \bar{x} ) + \frac{1}{2} \varOmega^{2} \bar{x}^{2} + \frac{\varOmega \phi_{W}}{B} , 0 , \varDelta_{M}(\bar{x}, U ) \right) . \end{align} In order to obtain 
$f_{\textrm {DSE}}$, and consequently 
$f_{W}$, it is necessary to determine the constants 
$\bar {x}_{c}$ and 
$\phi _{\textrm {DSE}}$, and the functions 
$\chi _{M}(\bar {x})$ and 
$\mu _{\textrm {op}}(\bar {x})$. Recall that, by (4.22), 
$\chi _{M}(\bar {x})$ and 
$\mu _{\textrm {op}}(\bar {x})$ also determine 
$\varDelta _{M}( \bar {x}, U )$. These quantities are specified by the electrostatic potential profile 
$\phi (x)$, which is obtained by solving the quasi-neutrality (1.3). Thus, (5.8) does not per se fully specify 
$f_{\textrm {DSE}} (v_x, v_y, v_z)$. Geraldini et al. (Reference Geraldini, Parra and Militello2018) derived an expression for the ion density 
$n_{i} (x)$ for 
$\alpha \ll 1$, as a functional of the electrostatic potential 
$\phi (x)$. Using this expression, an iterative scheme to obtain the numerical solution 
$\phi (x)$ of the quasi-neutrality (1.3) was presented. In the next section, a model for 
$f_{\textrm {DSE}} (v_x, v_y, v_z)$ is presented, which allows one to bypass obtaining a numerical solution of 
$\phi (x)$ across the whole magnetic presheath.
6. Large ion gyro-orbit model
In this section, we derive a closed set of equations for the quantities 
$\bar {x}_{c}$, 
$\phi _{\textrm {DSE}}$, 
$\chi _{M}(\bar {x})$ and 
$\mu _{\textrm {op}}(\bar {x})$ appearing in (5.8) for 
$f_{\textrm {DSE}}$ and (5.10) for 
$f_{W}$. The derivation assumes 
$\tau \gg 1$ and exploits the approximately undistorted nature of ion gyro-orbits in this limit. In § 6.1, the quasi-neutrality equation is expanded in the magnetic presheath close to the Debye sheath entrance, 
$\lambda _{D} \ll x \ll \rho _{s}$, to obtain a relationship between the distribution function and electric field. Then, in § 6.2, the expression for the electric field is used to derive expressions for the functions 
$\chi _{M}(\bar {x})$ and 
$\mu _{\textrm {op}}(\bar {x})$. This procedure is strictly not self-consistent, as the expression for the electric field derived in the previous subsection is valid closer to the wall than where it is used. To determine the large gyro-orbit distribution function, only the two parameters 
$\bar {x}_{c}$ and 
$\phi _{\textrm {DSE}}$ remain to be specified. In § 6.3, a method to solve for the two parameters is presented.
6.1. Quasi-neutrality at the Debye sheath entrance
In general, solving (1.3) in the magnetic presheath is a numerical task. However, near the Debye sheath entrance the quasi-neutrality equation can be expanded to obtain analytical expressions relating the electric field to the distribution function in this region. This analysis is valid for 
$\sqrt {Zm_{e} / m_{i} } \ll \alpha$, as it assumes (3.12) for the electron density.
The variation in density in the magnetic presheath, close to the Debye sheath entrance, for both ions and electrons is related to the variation in the electrostatic potential, 
$\delta \phi (x) = \phi (x) - \phi _{\textrm {DSE}}$. The Boltzmann distribution (3.12) is expanded near the Debye sheath entrance to obtain
\begin{equation} n_{e} (x) \simeq Z n_{\textrm{MPE}} \exp \left( \frac{e \phi_{\textrm{DSE}} }{T_{e}} \right) \left( 1 + \frac{e\delta \phi}{T_{e}} + \left( \frac{e\delta \phi}{T_{e}} \right)^{2} \right) . \end{equation}
The form of the expansion of the ion density in 
$\delta \phi (x)$ depends on whether ions with 
$v_x=0$ are present or not at the Debye sheath entrance, that is, whether 
$f_{\textrm {DSE}} (v_x= 0) = 0$ or not. If 
$f_{\textrm {DSE}} (v_x= 0) \neq 0$, (5.8) requires that 
$\chi _{M} (\bar {x}) = \varOmega ^{2} \bar {x}^{2} / 2 + \varOmega \phi _{\textrm {DSE}} / B$ for at least some values of 
$\bar {x}$, that is, type I ion orbits must be present. Thus, there are ions whose bottom turning point lies very close to the Debye sheath entrance at 
$x_{b} \leqslant x$. Such ions have a velocity range between 
$|v_x| = 0$ (
$x_{b} = x$) and 
$|v_x| = \sqrt {2( \chi _{M } (\bar {x}) - \chi (x, \bar {x} ) )} \simeq \sqrt {2( \varOmega ^{2} \bar {x} x - \varOmega \delta \phi (x) / B )} \sim \sqrt {\delta \phi }$ (
$x_{b} \simeq 0$), and can have both positive and negative values of 
$v_x$. These ions contribute to a term in the ion density proportional to 
$\sqrt {\delta \phi }$ (Geraldini et al. Reference Geraldini, Parra and Militello2018), heuristically owing to the size of the additional integration region in 
$v_x$. As no term in the electron density is proportional to 
$\sqrt {\delta \phi }$, type I ion orbits must be absent, requiring
\begin{equation} f_{\textrm{DSE}} (v_x= 0) = 0. \end{equation}
Recall from § 4 that all ions with 
$\bar {x} > \phi '_{\textrm {DSE}} / (\varOmega B )$, corresponding to a sufficiently large value of 
$\mu = \mu _{\textrm {op}} ( \bar {x})$, are in type I orbits. For there to be a complete absence of type I orbits, 
$\phi '(x)$ must be divergent at the Debye sheath entrance on the magnetic presheath scale, 
$\phi '_{\textrm {DSE}} \rightarrow \infty$.Footnote 2 As shown in the next subsection, this divergence also causes the asymptotic distribution function 
$f_{\textrm {DSE}}$ to decay exponentially for 
$v_x \rightarrow 0$ provided 
$F (\mu , U)$ decays exponentially for 
$U \rightarrow \infty$.
Excluding the presence of type I orbits, the ion density near the Debye sheath entrance is obtained by following ion characteristics backwards from the Debye sheath entrance. To lowest order in 
$\alpha$, the orbit parameters 
$\bar {x}$ and 
$U_{\perp }$ are constant; in addition, the total energy 
$U$ is exactly constant. Consider (4.10), (4.11) and (4.12) for the ion velocity in the magnetic presheath. The quantities 
$v_z$, 
$v_y + \varOmega x$ and 
$-\sqrt {v_x^{2} + 2\varOmega \delta \phi (x) / B - 2\varOmega ^{2} \bar {x} x + \varOmega ^{2} x^{2} }$ are constant and, from (4.19), (4.21) and (4.23), are equal to the components of the velocity at the Debye sheath entrance. Thus, the ion density at a distance 
$x$ from the wall, near the Debye sheath entrance, is
\begin{equation} n_{i} (x) \simeq \int f_{\textrm{DSE}}\left( -\sqrt{v_x^{2} + \frac{2\varOmega \delta \phi (x) }{ B } - 2\varOmega v_y x }, v_y + \varOmega x, v_z \right) \,\textrm{d}^{3} v. \end{equation}
Here, we have neglected the term 
$\varOmega ^{2} x^{2} \ll 2\varOmega v_y x$. The quasi-neutrality equation (1.3) to lowest order in 
$e\delta \phi (x) /T_{e} \ll v_x^{2} / v_{B}^{2}$ andFootnote 3 
$x \ll v_x^{2} / (\varOmega v_y ) \sim v_x^{2} / (\varOmega c_{s} )$ gives an equation for 
$\phi _{\textrm {DSE}}$,
\begin{equation} n_{\textrm{DSE}} \equiv \int f_{\textrm{DSE}} ( \boldsymbol{v} )\,\textrm{d}^{3} v = n_{\textrm{MPE}} \exp \left( \frac{e \phi_{\textrm{DSE}} }{T_{e}} \right). \end{equation}
In (6.4) we have denoted the lowest-order ion density at the Debye sheath entrance as 
$n_{\textrm {DSE}}$.
Considering the exponential decay of 
$f_{\textrm {DSE}}$ for 
$|v_x | \rightarrow 0$, the first argument of 
$f_{\textrm {DSE}}$ in (6.3) can be expanded in 
$e\delta \phi (x ) / T_{e} \ll v_x^{2} / v_{B}^{2}$ and 
$x \ll v_x^{2} / (\varOmega c_{s} )$ to give
\begin{equation} n_{i} (x) \simeq \int f_{\textrm{DSE}}\left( v_x + \frac{\varOmega \delta \phi }{ B v_x } - \frac{ \varOmega v_y x }{ v_x } - \frac{\varOmega^{2} \delta \phi^{2} }{ 2B^{2} v_x^{3} } , v_y + \varOmega x, v_z \right) \,\textrm{d}^{3} v. \end{equation}The result of Taylor expanding the integrand in (6.5) and subsequently integrating by parts is
\begin{align} n_{i} (x) &\simeq \int f_{\textrm{DSE}} ( \boldsymbol{v} ) \,\textrm{d}^{3} v + \frac{\varOmega \delta \phi}{ B } \int \frac{ f_{\textrm{DSE}} ( \boldsymbol{v} ) }{ v_x^{2} }\,\textrm{d}^{3}v - \varOmega x \int \frac{ v_y f_{\textrm{DSE}} ( \boldsymbol{v} ) }{ v_x^{2} } \,\textrm{d}^{3}v \nonumber\\ &\quad + \frac{3}{2} \left( \frac{ \varOmega \delta \phi }{ B } \right)^{2} \int \frac{f_{\textrm{DSE}} ( \boldsymbol{v} ) }{ v_x^{4} }\,\textrm{d}^{3} v. \end{align}
An alternative derivation of the same result is obtained by integrating the top-hat function in 
$v_x$ first and then expanding the resulting expression (Geraldini et al. Reference Geraldini, Parra and Militello2018). Note that the Taylor expansion of the second argument of 
$f_{\textrm {DSE}}$ in (6.3), 
$v_y + \varOmega x$, about 
$v_y$ did not give a variation in 
$x$. Collecting terms that are higher order than (6.4) in the quasi-neutrality equation gives an equation relating electrostatic potential variation and position,
\begin{align} &\frac{e\delta \phi}{T_{e}} \left( \int f_{\textrm{DSE}} ( \boldsymbol{v} ) \,\textrm{d}^{3} v - v_{B}^{2} \int \frac{f_{\textrm{DSE}} ( \boldsymbol{v} ) }{ v_x^{2} }\,\textrm{d}^{3}v \right) \nonumber\\ &\quad + \frac{1}{2} \left( \frac{e\delta \phi}{T_{e}} \right)^{2} \left( \int f_{\textrm{DSE}} ( \boldsymbol{v} )\, \textrm{d}^{3} v - 3 v_{B}^{4} \int \frac{f_{\textrm{DSE}} ( \boldsymbol{v} ) }{ v_x^{4} }\,\textrm{d}^{3} v \right) + \varOmega x \int \frac{f_{\textrm{DSE}} ( \boldsymbol{v} ) v_y }{ v_x^{2} }\,\textrm{d}^{3}v \simeq 0. \end{align}
As was concluded in the previous paragraph, because the electric field must diverge for 
$x \rightarrow 0$, the appropriate balance of terms in (6.7) is 
$\delta \phi ^{2} \propto x$. Therefore, the term linear in 
$\delta \phi$ must be set to zero, and we obtain the marginal form of the kinetic Bohm condition (Geraldini et al. Reference Geraldini, Parra and Militello2018),
\begin{equation} I_{\textrm{Bohm}} \equiv v_{B}^{2} \int \frac{f_{\textrm{DSE}} ( \boldsymbol{v} )}{v_x^{2}}\,\textrm{d}^{3} v = n_{\textrm{DSE}} . \end{equation}
In (6.8) we have defined the Bohm integral, 
$I_{\textrm {Bohm}}$, and we have used the definition of 
$n_{\textrm {DSE}}$ in (6.4).
The condition (6.8) applies to the lowest-order distribution function in the region 
$\lambda _{D} \ll x \ll \rho _{s}$. It does not apply to the exact distribution function measured near a target in an experiment (Baalrud & Hegna Reference Baalrud and Hegna2012; Riemann Reference Riemann2012). There are small corrections to the asymptotic distribution function 
$f_{\textrm {DSE}} ( \boldsymbol {v} )$ in the region 
$\lambda _{D} \ll x \ll \rho _{s}$. With a finite but large electric field, 
$\phi '(x)$, the distribution function in this region does not exactly satisfy 
$f(x, \boldsymbol {v})= 0$ for 
$v_x = 0$. One reason for this is the presence of a small number of very high-energy ions whose bottom turning point is only a few Debye lengths from the target, 
$x_{b} \sim \lambda _{D}$. A very small number of ion collisions or reflections from the target, both neglected, would also cause 
$f(x,\boldsymbol {v}) \neq 0$ for 
$v_x \geqslant 0$. If the exact distribution function, 
$f (x, \boldsymbol {v} )$, were used instead of the asymptotic one, 
$f_{\textrm {DSE}} (\boldsymbol {v})$, in the kinetic Bohm condition (6.8), then the left-hand side would diverge, 
$\int (\,f (x, \boldsymbol {v} ) / v_x^{2} )\,\textrm {d}^{3} v \rightarrow \infty$, and the condition could not even be approximately satisfied. Nonetheless, 
$f_{\textrm {DSE}} (\boldsymbol {v})$ is, within the validity of the underlying orderings, an approximation of the true distribution function in the region 
$\lambda _{D} \ll x \ll \rho _{s}$.
\begin{equation} \left( \frac{e\delta \phi}{T_{e}} \right)^{2} \left( \int f_{\textrm{DSE}} ( \boldsymbol{v} ) \,\textrm{d}^{3} v - 3 v_{B}^{4} \int \frac{f_{\textrm{DSE}} ( \boldsymbol{v} ) }{ v_x^{4} }\,\textrm{d}^{3} v \right) + 2 \varOmega x \int \frac{f_{\textrm{DSE}} ( \boldsymbol{v} ) v_y }{ v_x^{2} } \,\textrm{d}^{3}v \simeq 0. \end{equation}The electrostatic potential variation in the magnetic presheath, near the Debye sheath entrance, is thus given by
\begin{equation} \frac{e\left( \phi ( x ) - \phi ( 0 ) \right) }{T_{e}} \simeq \frac{ \sqrt{2 \bar{x}_{\textrm{av}} x} }{\rho_{B}} , \end{equation}
with 
$\bar {x}_{\textrm {av}}$, denoting a kinetic average of 
$\bar {x} = v_y / \varOmega$, given by
\begin{equation} \frac{ \bar{x}_{\textrm{av}} }{\rho_{B}} = \frac{ v_{B} \int \left( v_y f_{\text{DSE}} \left( \boldsymbol{v} \right) / v_x^{2} \right)\,\textrm{d}^{3} v }{ \int f_{\text{DSE}} \left( \boldsymbol{v} \right) \left( 3 v_{{B}}^{4} / v_x^{4} - 1 \right) \,\textrm{d}^{3} v} \sim \frac{\sqrt{1+\tau}}{(1+ \alpha \tau ) } . \end{equation}
Here, 
$\rho _{B} = v_{B} / \varOmega$ is referred to as the Bohm gyro-radius. As 
$f_{\textrm {DSE}}$ is exponentially small near 
$v_x=0$, the integral in the denominator of (6.11) is convergent. The ordering in (6.11) can be obtained as follows. Consider the smallest value of 
$|v_x |$ in the range (4.23) at the Debye sheath entrance (
$\phi (x) = \phi _{\textrm {DSE}}$),
\begin{equation} V_{x,\textrm{slow}} (\bar{x}) = \sqrt{2\left( \chi_{M} (\bar{x}) - \frac{1}{2} \varOmega^{2} \bar{x}^{2} - \frac{\varOmega \phi_{\textrm{DSE}}}{B} \right) }. \end{equation}
Ions with 
$v_x \simeq - V_{x, {\rm slow}}$ are referred to as ‘slow’ ions. From (6.4) and the ordering 
$|v_x| \sim v_{B} \sqrt {1+\alpha \tau }$ for typical values of 
$v_x$, the marginalized distribution function is ordered 
$\int f_{\textrm {DSE}}\,\textrm {d}v_y\,\textrm {d}v_z \sim n_{\textrm {DSE}} / ( v_{B}\sqrt {1 + \alpha \tau } )$. The kinetic Bohm condition (6.8) determines the size of slow ions, 
$\int (\,f_{\textrm {DSE}} / v_x^{2} ) \,\textrm {d}^{3} v \sim ( \int f_{\textrm {DSE}}\,\textrm {d}v_y\,\textrm {d}v_z ) / V_{x,\textrm {slow}} \sim n_{\textrm {DSE}} / (V_{x,\textrm {slow}} v_{B} \sqrt {1+\alpha \tau } )$. This gives the ordering 
$V_{x,\textrm {slow}} \sim v_{B} / \sqrt {1+\alpha \tau }$. Note that 
$V_{x, \textrm {slow}} \ll v_{B} \sqrt {1+\alpha \tau }$ only if 
$\alpha \tau \gg 1$, so that for 
$\alpha \tau \lesssim 1$ the normal velocity of slow ions is similar in size to the normal velocity of a typical ion. The size of 
$\bar {x}_{\textrm {av}}$ is obtained by considering the contribution of slow ions to the integrals in (6.11) and using also 
$v_y \sim c_{s}$, giving 
$\bar {x}_{\textrm {av}} / \rho _{B} \sim c_{s} V_{x,\textrm {slow}}^{2} / v_{B}^{3} \sim \sqrt {1+\tau } / (1+\alpha \tau )$.
The region of validity of (6.10) is obtained by investigating the validity of the expansion (6.7). For the expansion to be valid, the orderings 
$e\delta \phi (x ) / T_{e} \ll v_x^{2} / v_{B}^{2}$ and 
$x \ll v_x^{2} / (\varOmega c_{s} )$ must be satisfied. Using 
$x \ll V_{x, \textrm {slow}}^{2} / (\varOmega c_{s} )$, the ordering 
$x \ll \rho _{s} / [(1+\tau ) ( 1 + \alpha \tau ) ]$ for the region of validity of the expansion is obtained. The same ordering results from 
$e\delta \phi (x ) / T_{e} \ll V_{x,{\rm slow}}^{2} / v_{B}^{2}$ using (6.10) and (6.11).
6.2. Ion trajectories and ion distribution function for 
$\tau \gg 1$
In order to obtain 
$f_{\textrm {DSE}} (\boldsymbol {v})$ from (5.8), the electrostatic potential in the magnetic presheath is necessary to calculate: the function 
$\chi _{M} (\bar {x})$ from (4.17), the function 
$\mu _{\textrm {op}}(\bar {x})$ from (4.20), the quantity 
$\bar {x}_{c}$ from (4.18) and the quantity 
$\phi _{\textrm {DSE}}$. These quantities are calculated here using a model obtained by considering ion trajectories for 
$\tau \gg 1$ in the electrostatic potential of (6.10).
For 
$\tau \gg 1$, the thermal velocity of an ion is much larger than the Bohm velocity, 
$v_{t,i}^{2} \sim \tau v_{B}^{2} \gg v_{B}^{2}$. To calculate the adiabatic invariant, we can therefore neglect the small electrostatic potential variation throughout the orbit, 
$\varOmega \phi (x) / B \sim v_{B}^{2} \ll \varOmega \mu _{\textrm {op}}(\bar {x}) \sim v_{t,i}^{2}$, and using (4.20) obtain 
$\mu _{\textrm {op}} (\bar {x}) \simeq U_{\perp } / \varOmega \simeq \chi _{M} ( \bar {x} ) / \varOmega$. This does not specify the functional form of 
$\mu _{\textrm {op}} (\bar {x})$ and 
$\chi _{M} ( \bar {x} )$, but in relating them reduces the number of unknown functions from two to one. The approximate equivalence of 
$U_{\perp }$ and 
$\varOmega \mu$ and the conservation of 
$U$ and 
$\mu$ imply that 
$v_z = \sqrt {2(U - U_{\perp } )}$, has remained approximately unchanged from its value at the magnetic presheath entrance, 
$\sqrt {2(U - \varOmega \mu )}$. The quantity 
$\bar {x}_{c}$, defined in (4.18), corresponds to the orbit position of a gyro-orbit with adiabatic invariant equal to zero (since 
$x_{b} = x_{t} = x_{c}$), and thus 
$\bar {x}_{c}$ is obtained through 
$\mu _{\textrm {op}} (\bar {x}_{c}) =0$.
When an ion in a large gyro-orbit gets sufficiently close to the target, its gyro-motion is distorted as shown in figure 1(b). The net force away from the wall on an ion at a given instant is given by the effective potential gradient, 
$\chi '(x, \bar {x})$. The distortion of ion gyro-orbits is caused by a competition between the magnetic force pulling away from the wall and the electric force pushing towards the wall. As type I orbits are absent, 
$x_{M}$ is a stationary point where the electric force on the ion exactly balances the magnetic force. Its location can be obtained from (4.14) with 
$x_{\textrm {st}} = x_{M} < x_{c}$,
\begin{equation} \bar{x} = x_{M} + \frac{ \phi'(x_{M})}{\varOmega B}. \end{equation} In what follows, the electrostatic potential in (6.10) is used to approximate the electrostatic potential at distances from the wall corresponding to typical values of 
$x_{M}$. From (6.10) we obtain 
$\phi '(x_{M} ) = (T_{e} / e \rho _{B} ) \sqrt {\bar {x}_{\textrm {av}} / 2x_{M} }$. Using the ordering 
$x_{M} \ll \bar {x} \sim \phi '(x_{M})/\varOmega B \sim \rho _{s}$ in (6.13), we obtain
\begin{equation} x_{{M}} = \frac{\bar{x}_{\textrm{av}} \rho_{{B}}^{2}}{2 \bar{x}^{2}} . \end{equation}
By inserting (6.14) into 
$\chi _{M} ( \bar {x}) = \varOmega ^{2} (x_{M} -\bar {x} )^{2} / 2 + \varOmega \phi (x_{M}) / B$, neglecting the term 
$\varOmega ^{2} x_{M}^{2} / 2$ and remembering that 
$\mu _{\textrm {op}} (\bar {x}) \simeq \chi _{M} ( \bar {x} ) / \varOmega$, we obtain
\begin{equation} \varOmega \mu_{\textrm{op}} (\bar{x} ) \simeq \chi_{M}(\bar{x}) \simeq \frac{1}{2} \varOmega^{2} \bar{x}^{2} + \frac{v_{B}^{2} \bar{x}_{\textrm{av}}}{2\bar{x}} + \frac{\varOmega \phi_{\textrm{DSE}} }{ B } . \end{equation}
Imposing 
$\mu _{\textrm {op}} (\bar {x}_{c} ) = 0$ in (6.15) gives
\begin{equation} \bar{x}_{c} = \rho_{B} \sqrt{-\frac{2e \phi_{\textrm{DSE}}}{T_{e}} - \frac{v_{c}^{2}}{v_{B}^{2}} } , \end{equation}
where the quantity 
$v_{c} = v_{B} \sqrt { \bar {x}_{\textrm {av}} / \bar {x}_{c}}$, called the critical velocity, has been defined. With this definition, 
$\bar {x}_{\textrm {av}}$ is given by
\begin{equation} \bar{x}_{\textrm{av}} = \frac{ v_{c}^{2} }{v_{B} ^{2}} \bar{x}_{c}. \end{equation}From (6.15) and (4.23), large gyro-orbits at the Debye sheath entrance have a range of normal velocities given by
\begin{equation} \frac{v_{c}^{2} \bar{x}_{c}}{2\bar{x}} \leqslant \frac{v_x^{2}}{2} < \frac{v_{c}^{2} \bar{x}_{c}}{2\bar{x}} + 2 {\rm \pi}\alpha \mu_{\textrm{op}}'(\bar{x}) \sqrt{2\left( U - \varOmega \mu_{\textrm{op}}(\bar{x}) \right)} . \end{equation}Inserting the velocity spread (6.18) in the distribution function (5.8) the velocity distribution of ions in large gyro-orbits is
\begin{align} f_{\textrm{DSE}} (v_x, v_y, v_z) &\simeq F \left( \mu_{\textrm{op}}(\bar{x}) , U \right) \varTheta \left( \bar{x} - \bar{x}_{{c}} \right) \varTheta \left({-}v_x \right) \nonumber\\ &\quad\times \hat{\varPi} \left( \frac{1}{2}v_x^{2} - \frac{v_{c}^{2} \bar{x}_{c}}{2\bar{x}} , 0 , 2 {\rm \pi}\alpha \mu_{\textrm{op}}'(\bar{x}) \sqrt{2\left( U - \varOmega \mu_{\textrm{op}}(\bar{x}) \right)} \right) . \end{align} Despite being a useful analytical model for the ion distribution function, the large gyro-orbit model presented here is strictly not asymptotically self-consistent. For 
$\tau \gg 1$, 
$\bar {x} \sim \rho _{i}$ and 
$\bar {x}_{\textrm {av}} \sim \rho _{i} / (1+\alpha \tau )$, where 
$\rho _{i} = v_{t,i}/ \varOmega$ is the thermal ion gyro-radius. Using (6.14), it follows that 
$x_{M} \sim \rho _{i} /[ \tau ( 1+\alpha \tau ) ]$. Recall from the final paragraph of § 6.1 that the expansion used to derive (6.10) is valid, for 
$\tau \gg 1$, in the region 
$x_{M} \ll \rho _{i} /[ \tau ( 1+\alpha \tau ) ]$. Therefore, (6.14) is not valid for the majority of ions. There is, however, a minority of ions for which 
$\mu \gg v_{t,i}^{2} / \varOmega$ and 
$\bar {x} \gg \rho _{i}$, which have 
$x_{M} \ll \rho _{i} /[ \tau ( 1+\alpha \tau ) ]$. For these ions (6.14) is accurate. This can be used to derive the exponential decay of 
$f_{\textrm {DSE}}$ at 
$|v_x| \rightarrow 0$ as follows. The distribution function 
$F (\mu , U )$ is assumed to exponentially decay for 
$U \rightarrow \infty$ and consequently, because 
$U \geqslant \varOmega \mu$, for 
$\varOmega \mu \rightarrow \infty$, such that 
$F \sim \exp ( - 2\varOmega \mu / v_{t,i}^{2} )$. If follows from (6.15) and 
$\bar {x} \rightarrow \infty$ that 
$F \sim \exp ( - \varOmega ^{2} \bar {x}^{2} / v_{t,i}^{2} )$. The slowest value of 
$|v_x|$ in the top-hat function in (5.8) is given by the function 
$V_{x,\textrm {slow}}$ in (6.12), which in the model is
\begin{equation} V_{x,\textrm{slow}} (\bar{x}) = v_{c} \sqrt{\frac{ \bar{x}_{c} }{\bar{x}}}. \end{equation} For the top-hat function in (5.8) to be non-zero we require 
$V_{x,{\rm slow}} (\bar {x}) \leqslant |v_x|$ and so 
$\bar {x} \geqslant \bar {x}_{c} v_{c}^{2} / v_x^{2}$. Therefore, the largest value of 
$f_{\textrm {DSE}}$ for an ion with 
$v_x \rightarrow 0$ satisfies 
$f_{\textrm {DSE}} \sim \exp ( - \varOmega ^{2} \bar {x}_{c}^{2} v_{c}^{4} /(v_{t,i}^{2} v_x^{4} ) )$, which is exponentially small.
The critical velocity is the value of 
$| v_x |$ for an ion at the Debye sheath entrance with 
$\mu = 0$, which came from an infinitesimally small gyro-orbit, 
$v_{c} = V_{x, {\rm slow}} (\bar {x}_{c})$. These ions should have 
$\mu _{\textrm {op}}' ( \bar {x}_{c} ) = 0$ and thus 
$\varDelta _{M} (\bar {x}_{c}, U ) = 0$ for all values of 
$U$, which would give 
$v_x = - v_{c}$ as the only allowed value according to the velocity distribution (5.8). However, ions with 
$\mu = 0$ in the model have a finite range of velocities owing to the fact that 
$\mu _{\textrm {op}}' (\bar {x}_{c}) = 0$ is not imposed in order not to overconstrain the model. This could be concerning, because if 
$\mu _{\textrm {op}}'(\bar {x}_{c} ) < 0$ (and so 
$\varDelta _{M} (\bar {x}_{c}, U ) < 0$), the range of values of 
$v_x^{2}$ in (4.23) would allow for non-real values of 
$v_x$. Fortunately, 
$\mu _{\textrm {op}}'(\bar {x}_{c} ) = \varOmega \bar {x}_{c} - v_{c}^{2} / (2\varOmega \bar {x}_{c} )$ is always positive if 
$\alpha$ is sufficiently small that 
$e|\phi _{\textrm {DSE}}| / T_{e} \sim \ln \alpha \gg 1$, as (6.16) leads to 
$2 (\varOmega \bar {x}_{c})^{2} \sim |\ln \alpha | v_{B}^{2} \gg v_{c}^{2} \sim v_{B}^{2}$. In practice, 
$\mu _{\textrm {op}}'(\bar {x}_{c} ) > 0$ for all values of 
$\alpha \leqslant 5^{\circ }$ considered in this paper. With 
$\mu _{\textrm {op}}'(\bar {x}_{c} ) > 0$, ions with 
$\mu = 0$ (
$\bar {x} = \bar {x}_{c}$) have a non-zero range of values of 
$v_x$ according to (4.22) and (4.23). In the model, 
$|v_x | = v_{c}$ is therefore the smallest value of 
$|v_x|$ for an ion with 
$\mu = 0$. Although 
$\mu _{\textrm {op}}' (\bar {x}_{c}) \neq 0$ may look like a serious shortcoming of the model, for 
$\tau \gg 1$ the large discrepancy in the function 
$\mu _{\textrm {op}}' (\bar {x})$ is expected only for a small number of particles near 
$\bar {x} = \bar {x}_{c}$. In other words, the model does not correctly capture the small gyro-orbits, but there are assumed to be only a small number of them anyway.Footnote 4
6.3. Model closure: calculating $\phi _{\textrm {DSE}}$ and $v_{c}$
The only unknowns that specify the model distribution function (6.19) are the two constants 
$\phi _{\textrm {DSE}}$ and 
$v_{c}$. The value of 
$\phi _{\textrm {DSE}}$ is determined from quasi-neutrality at the Debye sheath entrance (6.4). The value of 
$v_{c}$ is determined by imposing the kinetic Bohm condition (6.8).
For numerical evaluation, it is best to re-express all velocity moments as
\begin{align} \int f_{\textrm{DSE}}(\boldsymbol{v}) v_x^{a} \,\textrm{d}^{3} v &= \int_{\bar{x}_{c}}^{\infty} \varOmega \,\textrm{d}\bar{x} \int_{\varOmega \mu_{\textrm{op}}(\bar{x})}^{\infty} F \left( \mu_{\textrm{op}}(\bar{x}), \varOmega \mu_{\textrm{op}}(\bar{x}) + \frac{1}{2} v_z^{2} \right) \nonumber\\ &\quad \times \frac{v_{c}^{a+1}}{a+1} \left( \left( \frac{ \bar{x}_{c} }{ \bar{x} } + \frac{ 4 {\rm \pi}\alpha \mu_{\textrm{op}}'(\bar{x}) v_z }{v_{c}^{2}} \right)^{(a+1)/2} - \left( \frac{ \bar{x}_{c} }{ \bar{x} } \right)^{(a+1)/2} \right)\,\textrm{d}v_z , \end{align}
obtained from (6.19) using the change of variables 
$v_y = \varOmega \bar {x}$ and 
$v_z = \sqrt {2( U - \varOmega \mu _{\textrm {op}}(\bar {x}) )}$, and substituting (6.17). In particular, to solve (6.4) and (6.8) for 
$\phi _{\textrm {DSE}}$ and 
$v_{c}$, we require the density,
\begin{align} n_{\textrm{DSE}} &= \int_{\bar{x}_{c}}^{\infty} \varOmega \,\textrm{d}\bar{x} \int_{\varOmega \mu_{\textrm{op}}(\bar{x})}^{\infty} F \left( \mu_{\textrm{op}}(\bar{x}), \varOmega \mu_{\textrm{op}}(\bar{x}) + \frac{1}{2} v_z^{2} \right) \nonumber\\ &\quad \times v_c \left( \left( \frac{ \bar{x}_{c} }{ \bar{x} } + \frac{ 4 {\rm \pi}\alpha \mu_{\textrm{op}}'(\bar{x}) v_z }{v_{c}^{2}} \right)^{1/2} - \left( \frac{ \bar{x}_{c} }{ \bar{x} } \right)^{1/2} \right)\,\textrm{d}v_z, \end{align}and the Bohm integral,
\begin{align} I_{\textrm{Bohm}} &= v_{B}^{2} \int_{\bar{x}_{c}}^{\infty} \varOmega \,\textrm{d}\bar{x} \int_{\varOmega \mu_{\textrm{op}}(\bar{x})}^{\infty} F \left( \mu_{\textrm{op}}(\bar{x}), \varOmega \mu_{\textrm{op}}(\bar{x}) + \frac{1}{2} v_z^{2} \right) \nonumber\\ &\quad \times \frac{1}{v_{c}} \left( \left( \frac{ \bar{x}_{c} }{ \bar{x} } \right)^{{-}1/2} -\left( \frac{ \bar{x}_{c} }{ \bar{x} } + \frac{ 4 {\rm \pi}\alpha \mu_{\textrm{op}}'(\bar{x}) v_z }{v_{c}^{2}} \right)^{{-}1/2} \right)\,\textrm{d}v_z . \end{align}
Note that the value of 
$I_{\textrm {Bohm}}$ decreases by increasing 
$v_{c}$, and vice versa.
Iterative expressions are used to determine 
$\phi _{\text {DSE}}$ from (6.4) and 
$v_{c}$ from (6.8). The first guesses, or zeroth iterates, are defined by 
$\phi _{\text {DSE},0} = (T_{e} / e )\ln \alpha$ and 
$v_{{c},0} = v_{B}$, and iteration values are denoted by 
$\phi _{\text {DSE},\nu }$ and 
$v_{{c},\nu }$. At each iteration, 
$n_{\text {DSE},\nu }$ and 
$I_{\text {Bohm},\nu }$ are evaluated from (6.22) and (6.23). The iterates 
$\phi _{\text {DSE},\nu +1}$ and 
$v_{{c},\nu +1}$ are obtained using
\begin{gather} \phi_{\text{DSE},\nu+1} = \frac{T_{e} }{e} \ln \left( \frac{ n_{\text{DSE},\nu} }{ n_{\textrm{MPE}} } \right), \end{gather}
\begin{gather} \left. \begin{aligned} v_{{c},\nu+1} & = \frac{v_{{c},\nu}}{ I_{\text{Bohm},\nu} } \left( I_{\text{Bohm},\nu} - n_{\text{DSE},\nu} \right), & \text{if } I_{\text{Bohm},\nu} > n_{\text{DSE},\nu} ,\\ & = \epsilon_{v_{c}}, & \text{otherwise}. \end{aligned} \right\} \end{gather}
Equation (6.24) originates from the rearranged form of (6.4), 
$\phi _{\text {DSE}} = (T_{e} / e ) \ln (n_{\text {DSE}} / n_{\textrm {MPE}} )$. Equation (6.25) is based on a Newton method with the approximations 
$\textrm {d} n_{\textrm {DSE}} / \textrm {d}v_{c} \approx 0$ and 
$\textrm {d}I_{\textrm {Bohm}} / \textrm {d}v_{c} \approx - I_{\textrm {Bohm}} / v_{c}$. The iteration is truncated when
\begin{gather} \frac{n_{\text{DSE},N} - n_{\text{DSE},N-1} }{ n_{\text{DSE},N} } < \epsilon_{n}, \end{gather}
\begin{gather} \left| \frac{ I_{\text{Bohm},N} }{ n_{\text{DSE},N} } - 1 \right | < \epsilon_{I} . \end{gather}
In the earliest iterations, it may happen that 
$v_{{c},\nu +1} \leqslant 0$, which is prevented by setting 
$v_{{c},\nu +1}$ to be a small number above zero (smaller than the solution 
$v_{{c}}$), denoted by 
$\epsilon _{v_{c}}$. The 
$N$th iteration values of 
$\phi _{\text {DSE}}$ and 
$v_{c}$, satisfying conditions (6.26) and (6.27), are considered to be acceptable numerical solutions of (6.4) and (6.8). The value of 
$\bar {x}_{\textrm {av}}$ is obtained from 
$v_{c}$ using (6.17). To obtain the results presented in the next section, 
$\epsilon _{n} = \epsilon _{I} = \epsilon _{v_{c}} = 10^{-10}$ was used.
Having solved (6.4) and (6.8) for 
$\phi _{\textrm {DSE}}$ and 
$v_{c} = v_{B} \sqrt { \bar {x}_{\textrm {av}} / \bar {x}_{c}}$, (4.19), (4.21), (6.15), (6.16) and (6.19) completely specify the large gyro-orbit model distribution function at the Debye sheath entrance, 
$f_{\textrm {DSE}}(\boldsymbol {v})$. The model distribution function at the wall is obtained by replacing (6.19) with
\begin{align} f_{W} (\boldsymbol{v}) &\simeq F \left( \mu_{\textrm{op}}(\bar{x}) , U \right) \varTheta \left( \bar{x} - \bar{x}_{{c}} \right) \varTheta \left({-}v_x \right) \nonumber\\ &\quad \times \hat{\varPi} \left( \frac{1}{2}v_x^{2} - \frac{v_{c}^{2} \bar{x}_{c}}{2\bar{x}}- \frac{\varOmega}{B} \left( \phi_{\textrm{DSE}} - \phi_{W} \right) , \ 0 ,\ 2 {\rm \pi}\alpha \varOmega \bar{x} \sqrt{2\left( U - \varOmega \mu_{\textrm{op}}(\bar{x}) \right)} \right) , \end{align}
where (5.7) determines the wall potential 
$\phi _{W}$.
To conclude this section, the application of the model to 
$\tau \lesssim 1$ is discussed. We have seen that the model is derived assuming 
$\tau \gg 1$, although it is not asymptotically self-consistent even in this limit. The Bohm condition closure (6.8) used in the model to obtain 
$v_{c}$ (and 
$\bar {x}_{c}$) is nonetheless valid for all 
$\tau$. Therefore, for 
$\tau \ll 1$ the model correctly recovers a distribution function that is centred around 
$v_x \simeq - v_{B}$, as expected from the fluid cold-ion result (Chodura Reference Chodura1982). This extends the applicability of the model to smaller values of 
$\tau$, though with less-accurate results. A measure of the accuracy of the model can be obtained by calculating the value of 
$\bar {x}_{\textrm {av}}$ from (6.11) and comparing it with the model value in (6.17). For 
$\alpha \tau \gg 1$, the two values are found to approach each other. For 
$\tau \ll 1$ the two values are found to differ approximately (with an 
$O(\alpha )$ error) by a factor of two: indeed, (6.11) results in 
$\bar {x}_{\textrm {av}} \simeq \bar {x}_{c} / 2$ upon using a cold-ion distribution function centred at 
$v_x = - v_{B}$ and 
$v_y = \varOmega \bar {x}_{c}$, whereas the model value from (6.17) is 
$\bar {x}_{\textrm {av}} \simeq \bar {x}_{c}$.
7. Numerical results
In this section, a comparison is presented of ion velocity distributions obtained from:
(i) equations (5.8), (5.10) and the full numerical solution
$\phi (x)$ of the quasi-neutrality equation in the magnetic presheath entrance; and(ii) equations (6.19), (6.28) and the closure equations of the large gyro-orbit model.
To obtain the solutions (i), the numerical scheme in Geraldini et al. (Reference Geraldini, Parra and Militello2018) is used. In § 7.1, the boundary conditions for the distribution function at the magnetic presheath entrance, as a function of 
$\tau$, are given. Then, in § 7.2, results for the distribution of the component 
$v_x$ of the ion velocity at the Debye sheath entrance, obtained using (i) and (ii), are presented. Finally, results for the energy–angle distributions of ions at the wall are presented in § 7.3 for some values of 
$\alpha$ and 
$\tau$. The possibility to extend the model for 
$\alpha \sim \sqrt {Zm_{e} / m_{i}}$ is briefly discussed in § 7.4.
7.1. Boundary conditions at the magnetic presheath entrance
The ion velocity distribution at the magnetic presheath entrance, 
$\rho _{s} \ll x \ll d_{c}$, is taken to be
\begin{equation} f_{\textrm{MPE}} \left( \boldsymbol{v} \right) = \begin{cases} \mathcal{N} n_{\textrm{MPE}} \dfrac{4 v_z^{2}}{{\rm \pi}^{3/2} v_{{t,i}}^{5}} \exp \left( - \dfrac{ \left| \boldsymbol{v} - u v_{{t,i}} \hat{\boldsymbol{e}}_z \right|^{2} }{v_{{t,i}}^{2}} \right) \varTheta \left( v_z \right), & \text{for } \tau \leqslant 1 ,\\ \mathcal{N} n_{\textrm{MPE}} \dfrac{ 4 v_z^{2} }{ {\rm \pi}^{3/2} v_{{t,i}}^{3} \left( v_{{t,i}}^{2} +r v_z^{2} \right)} \exp \left( - \dfrac{\left| \boldsymbol{v} \right|^{2} }{v_{{t,i}}^{2}} \right)\varTheta \left( v_z \right), & \text{for } \tau > 1 , \end{cases} \end{equation}
for any prescribed value of 
$\tau$, where 
$\varTheta$ is the Heaviside step function defined in (3.6) and 
$\hat {\boldsymbol {e}}_z$ is a unit vector in the 
$z$ direction. The family of velocity distributions (7.1) is the same as used in Geraldini et al. (Reference Geraldini, Parra and Militello2019) to study the dependence of the magnetic presheath solution on ion temperature, and is chosen to satisfy the marginal kinetic Chodura condition (Geraldini et al. Reference Geraldini, Parra and Militello2018)
\begin{equation} v_{{B}}^{2} \int \frac{ f_{\textrm{MPE}} \left( \boldsymbol{v} \right)}{v_z^{2}} \,\textrm{d}^{3}v = n_{\textrm{MPE}} . \end{equation}
The value of the normalization constant 
$\mathcal {N}$ is obtained from (5.3), giving
\begin{equation} \mathcal{N} = \begin{cases} \left[ \left( 1 + 2u^{2} \right) \left( 1 + \text{erf}(u) \right) + \dfrac{2u}{\sqrt{\rm \pi}} \exp({-}u^{2}) \right]^{{-}1}, & \text{for } \tau \leqslant 1 , \\ r^{3/2} \left[ 2\sqrt{r} - 2\sqrt{\rm \pi} \exp\left(\dfrac{1}{r}\right) \left( 1 - \text{erf} \left( \dfrac{1}{\sqrt{r}} \right) \right) \right]^{{-}1}, & \text{for } \tau > 1 . \end{cases} \end{equation}
The values of 
$u$ and 
$r$ are obtained by imposing (7.2), leading to
\begin{gather} 1 + \text{erf}(u) = \tau \left[ \left( 1 + 2u^{2} \right) \left( 1 + \text{erf} (u) \right) + \frac{2u}{\sqrt{\rm \pi}} \exp({-}u) \right], \end{gather}
\begin{gather} r \sqrt{\rm \pi} \exp\left(\frac{1}{r}\right) \left( 1 - \text{erf} \left( \frac{1}{\sqrt{r}} \right) \right) = \tau \left[ 2\sqrt{r} - 2\sqrt{\rm \pi} \exp\left(\frac{1}{r}\right) \left( 1 - \text{erf} \left( \frac{1}{\sqrt{r}} \right) \right) \right] . \end{gather}7.2. Narrowing of the wall-normal velocity distributions
The marginalized distribution function
\begin{equation} f_{x, \textrm{DSE}} (v_x) = \iint f_{\textrm{DSE}} (\boldsymbol{v}) \,\textrm{d}v_y \,\textrm{d}v_z, \end{equation}
is the distribution of wall-normal velocities 
$v_x$ of ions at the Debye sheath entrance. The numerical results obtained for 
$f_{x,\textrm {DSE}} (v_x)$ with the model and the theory for 
$\tau = 1$ and 
$\tau = 5$, for a number of angles 
$\alpha$, are shown in figure 3. The first thing to note is that the model distribution function (dashed lines) captures the essential features of the distribution function obtained from the full solution of the magnetic presheath electrostatic potential 
$\phi (x)$ (solid lines). Moreover, the agreement is better for the largest value of 
$\tau =T_{i} / ( ZT_{e} )$, 
$\tau = 5$, as expected.
(a,c) Wall-normal velocity distributions at the Debye sheath entrance from the numerical solution of 
$\phi (x)$ in the magnetic presheath (solid lines) and from the large gyro-orbit model (dashed lines), for 
$\tau = 1$ (a) and 
$\tau = 5$ (c) for angles 
$\alpha = 1^{\circ }, 3^{\circ }, 5^{\circ }$. (b,d). The variance 
$\langle \tilde {v}_x^{2} \rangle$ of the distributions from the numerical solution of 
$\phi (x)$ (circles) and from the model (crosses) for values of 
$\alpha$ between 
$0.5^{\circ }$ and 
$5^{\circ }$ for 
$\tau = 1$ (b) and 
$\tau = 5$ (d). The dotted lines are drawn to guide the eye, showing the linear scaling 
$\langle \tilde {v}_x^{2} \rangle / v_{t,i}^{2} \sim \alpha$ for 
$\alpha \gtrsim 1^{\circ }$. (Note: here 
$\alpha$ is measured in degrees.).
The width of the function 
$f_{x,\textrm {DSE}} (v_x)$ narrows as 
$\alpha$ decreases, a feature that was observed in Geraldini et al. (Reference Geraldini, Parra and Militello2018). The width of this function can be quantified using the variance 
$\langle \tilde {v}_x^{2} \rangle$, defined using the second moment of 
$f_{x,{\rm DSE}} ( v_x )$,
\begin{equation} \langle \tilde{v}_x^{2} \rangle = \sqrt{\frac{\int (v_x - u_{x,\textrm{DSE}})^{2} f_{x,\textrm{DSE}}(v_x) \,\textrm{d}v_x }{\int f_{x,\textrm{DSE}}(v_x) \,\textrm{d}v_x } } . \end{equation}Here
\begin{equation} u_{x,\textrm{DSE}} = \frac{ \int v_x f_{x,\textrm{DSE}}(v_x) \,\textrm{d}v_x }{ {\int f_{x,{\rm DSE}}(v_x) \,\textrm{d}v_x } } \end{equation}
is the average wall-normal velocity at the Debye sheath entrance. As can be seen in figure 3, the variance of the distribution function scales linearly with 
$\alpha$.
The scaling of the variance can be explained as follows. The ion velocity can be decomposed into two pieces: a piece coming from the electric field acceleration which depends only on 
$\bar {x}$ (or 
$\mu$), 
$V_{x,\textrm {slow}}(\bar {x} ) = \sqrt { 2( \chi _{M}(\bar {x}) - \varOmega ^{2} \bar {x}^{2} / 2 - \varOmega \phi _{\textrm {DSE}} / B ) }$, and an additional gyro-phase dependent piece which gives the velocity range in (4.23). In figure 4 the behaviour of 
$V_{x,\textrm {slow}} (\bar {x} )$ as a function of 
$\mu _{\textrm {op}} (\bar {x})$ is shown for some values of 
$\tau$ and 
$\alpha$. The slow decay of 
$V_{x,{\rm slow}}$ with 
$\mu$ is approximately captured by the model for 
$\tau = 5$, and for 
$(\tau , \alpha ) = ( 1, 5^{\circ } )$. For 
$(\tau , \alpha ) = (1, 1^{\circ } )$, the dependence of 
$V_{x,{\rm slow}}$ on 
$\mu$ is stronger than predicted by the model, but is nonetheless fairly weak. As 
$V_{x,\textrm {slow}}$ is only a weakly decreasing function of 
$\mu$, the distribution function sharply drops to zero around 
$|v_x | \approx V_{x,\textrm {slow}} ( \rho _{s} )$, a feature common to all velocity distributions in figure 3. The dominant contribution to the variance 
$\langle \tilde {v}_x^{2} \rangle$ therefore comes from the range of allowed values of 
$|v_x |$ in (4.23), instead of the dependence of 
$V_{x, \textrm {slow}}$ on 
$\mu$. For 
$\tau \gtrsim 1$, we order 
$\varOmega \mu _{\textrm {op}} \sim \varOmega ^{2} \bar {x}^{2} / 2 \sim v_{t,i}^{2}$ and 
$\sqrt {2( U - \varOmega \mu _{\textrm {op}}(\bar {x}) )} \sim v_{t,i}$, and obtain 
$2 {\rm \pi}\alpha \mu _{\textrm {op}}'(\bar {x}) \sqrt {2( U - \varOmega \mu _{\textrm {op}}(\bar {x}) )} \sim 2{\rm \pi} \alpha v_{t,i}^{2}$. Hence, the variance is 
$\langle \tilde {v}_x^{2} \rangle \sim \alpha v_{t,i}^{2} \sim \alpha \tau v_{B}^{2}$, as seen in the numerical results. The dependence of 
$V_{x,{\rm slow}}$ on 
$\mu$ does not cause a significant contribution to 
$\langle \tilde {v}_x^{2} \rangle$ unless 
$\alpha \tau$ is extremely small, seen in the numerical results of figure 3 as a saturation of the decrease of the variance for 
$\alpha \lesssim 1^{\circ }$.
The velocity of slow ions, 
$V_{x, \textrm {slow}} (\bar {x})$, is shown as a function of the adiabatic invariant 
$\mu _{\textrm {op}} (\bar {x})$ with (
$\tau , \alpha$) labelled. Solid lines are obtained from the numerical solution of 
$\phi (x)$; dashed lines are obtained from the large gyro-orbit model.
When deriving the scaling of (6.11), the typical value of 
$| v_x |$ of slow ions was found to be 
$V_{x,{\rm slow}} \sim v_{B} / \sqrt {1+\alpha \tau }$. From figure 4 it appears that the ordering 
$\alpha \tau \gtrsim 1$ is satisfied, as 
$V_{x,{\rm slow}}(\bar {x})$ is smaller than 
$v_{B}$ in most cases shown here. It may appear concerning that 
$V_{x,\textrm {slow}} / v_{B}$ is quite small also for 
$(\tau , \alpha ) = ( 1, 5^{\circ } )$, as this suggests that 
$\alpha \tau$ is large for 
$\tau = 1$ and for a value of 
$\alpha$ (
$=5^{\circ } \approx 0.09 \ \textrm {radians}$) which is considered small. This observation prompts a closer analysis of the validity of the asymptotic theory of the ion orbits, which assumes 
$\alpha \ll 1$. One of the consequences of this ordering is that the function 
$\varDelta _{M}(\bar {x}, U )$ is small. For 
$\tau \gtrsim 1$, the smallness of 
$\varDelta _{M}(\bar {x}, U )$ is measured relative to the kinetic energy of the ion,Footnote 5 estimated from the tangential components of the ion velocity, 
$(v_y^{2} + v_z^{2} ) / 2 \sim ( \varOmega ^{2} \bar {x}^{2} + v_{t,i}^{2} ) / 2$. The ratio 
$2\bar { \varDelta }_{M} / ( \varOmega ^{2} \bar {x}^{2} + v_{t,i}^{2} )$ is shown in figure 5 and highlights that, although for 
$\alpha = 5^{\circ }$ the validity of the asymptotic theory is not robust, the contribution of 
$\varDelta _{M}$ to the ion energy is smaller than the total kinetic energy for most ions, albeit by a factor of approximately two only. Note that 
$\alpha = 5^{\circ }$ corresponds to 
$2{\rm \pi} \alpha \approx 0.6 \ \textrm {radians}$, and so the factor of 
$2{\rm \pi}$ in (4.22) explains why the expansion in 
$\alpha$ starts to becomes inaccurate at 
$\alpha \approx 5^{\circ }$.
The quantity 
$2 \bar {\varDelta }_{M} / (\varOmega ^{2} \bar {x}^{2} + v_{t,i}^{2} )$, with 
$\bar {\varDelta }_{M} (\bar {x}) = 2{\rm \pi} \alpha \mu _{\textrm {op}}' (\bar {x} ) v_{t,i}$, is shown as a function of the adiabatic invariant 
$\mu _{\textrm {op}} (\bar {x})$ for labelled values of (
$\tau , \alpha$). Solid lines are obtained from the numerical solution of 
$\phi (x)$; dashed lines are obtained from the large gyro-orbit model. For 
$2 \bar {\varDelta }_{M} / (\varOmega ^{2} \bar {x}^{2} + v_{t,i}^{2} ) \ll 1$ the asymptotic theory in 
$\alpha \ll 1$ is valid.
Although this subsection presented ion distribution functions at the Debye sheath entrance, 
$f_{x,\textrm {DSE}} (\boldsymbol {v})$, the validity of the scaling 
$\langle \tilde {v}_x^{2} \rangle \sim \alpha v_{t,i}^{2}$ is expected to apply also to the ion velocity distribution at the wall, 
$f_{x,W} (\boldsymbol {v})$. In the next subsection, ion velocity distributions at the wall are considered for parameters where 
$\phi _{\textrm {DSE}} > \phi _{W}$, such that the assumption of Boltzmann electrons (recall (3.14)) remains at least approximately correct.
7.3. Energy–angle distributions at the target
As sputtering predictions depend on the distribution of kinetic energy and angle of impact of ions reaching the target, it is useful to calculate the energy–angle distribution of ions at the wall. To obtain our results, we considered a Deuterium plasma such that 
$\sqrt { m_{e} / m_{i}} = 0.0165$ and 
$Z=1$.
The kinetic energy of an ion at the wall is 
$E = U - \varOmega \phi _{W} / B$ and the angle of impact of an ion with the wall surface is 
$\sin \theta = |v_x| / \sqrt {2E }$. Thus, the components 
$v_z$ and 
$v_x$ of the ion velocity can be expressed as functions of 
$\bar {x}$, 
$E$ and 
$\theta$ via
\begin{gather} v_x ={-} \sqrt{2E} \sin \theta, \end{gather}
\begin{gather} v_z = \sqrt{2\left( E - \chi_{M} (\bar{x} ) + \frac{\varOmega \phi_{W} }{ B} \right) }. \end{gather}
The energy–angle distribution 
$\zeta _{W} (E, \theta )$ is calculated from 
$f_{W} ( \boldsymbol {v} )$ using the equation
\begin{equation} \zeta_{W} (E,\theta ) = \int_{ \bar{x}_{c}}^{\chi_{M}^{{-}1} ( E + \varOmega \phi_{W} / B )} \frac{\sqrt{2E}\cos \theta }{\sqrt{2\left( E - \chi_{M} (\bar{x} ) + \varOmega \phi_{W} / B \right) } } f_{W} ( \boldsymbol{v} ) \varOmega \,\textrm{d}\bar{x} , \end{equation}where the Jacobian
\begin{equation} \frac{\partial (v_x, v_z )}{ \partial ( E, \theta ) } = \frac{\sqrt{2E}\cos \theta }{\sqrt{2\left( E - \chi_{M} (\bar{x} ) + \varOmega \phi_{W} / B \right) } } \end{equation}
was used to change variables from 
$v_x$ and 
$v_z$ to 
$E$ and 
$\theta$. The inverse function of 
$\chi _{M } (\bar {x})$, denoted by 
$\chi _{M}^{-1}$, is used to obtain the maximum value of 
$\bar {x}$ for a given value of 
$E$, which is, from (7.10), the solution of 
$\chi _{M} (\bar {x}) = E + \varOmega \phi _{W} / B$.
The energy–angle distributions calculated from the numerical solution of the electrostatic potential in the magnetic presheath and from the large gyro-orbit model are shown for 
$\alpha = 3^{\circ }$ and 
$5^{\circ }$, for 
$\tau = 0.5$, in figure 6, and for 
$\tau =2$, in figure 7. The qualitative features of the distribution function obtained from the full electrostatic potential solution are, even for 
$\tau = 0.5$, adequately captured by the model, including the average angle of impact of ions with the wall. The model performs better at the largest of the two values of 
$\tau$ (
$\tau = 2$, figure 7), as expected.
Energy–angle distributions at the target, 
$\zeta _{W} ( E, \theta )$, obtained from the full electrostatic potential solution (full theory) and from the large gyro-orbit model for 
$\tau = 0.5$ and 
$\alpha = 3^{\circ }$ and 
$5^{\circ }$, are shown normalized to their peak value.
Energy–angle distributions at the target, 
$\zeta _{W} ( E, \theta )$, obtained from the full electrostatic potential solution (full theory) and from the large gyro-orbit model for 
$\tau = 2$ and 
$\alpha = 3^{\circ }$ and 
$5^{\circ }$, are shown normalized to their peak value.
7.4. Accounting for $\sqrt {Zm_{e} / m_{i}} \sim \alpha$
For some of the angles we have considered, the assumption of adiabatic electrons, 
$\sqrt {Zm_{e} / m_{i} } = 0.0165 \approx 1^{\circ } \ll \alpha$, is not well-satisfied. Once 
$\phi _{\textrm {DSE}} - \phi _{W} \leqslant 0$ our assumption that the Debye sheath repels most electrons back into the magnetic presheath is clearly incorrect. In fact, the Boltzmann distribution for the electron density becomes inaccurate when 
$\phi _{\textrm {DSE}} - \phi _{W}$ becomes sufficiently small that the ordering (3.14) is no longer satisfied. The critical value of 
$\alpha$ for which 
$\phi _{\textrm {DSE}} = \phi _{W}$ in the model increases slightly with 
$\tau$: for 
$\tau = 2$ it is 
$\alpha \approx 3^{\circ }$, whereas for 
$\tau = 10$ it is 
$\alpha \approx 5^{\circ }$. In order to solve for the self-consistent electrostatic potential across the magnetic presheath, a more accurate expression for the electron density must be used. In the context of the large gyro-orbit model, this is expected to change (6.4) and (6.8).
8. Conclusions
The velocity distribution of ions reaching a planar target when the angle between the magnetic field and the target is small, 
$\alpha \ll 1$, has been calculated using a model consisting of (4.19), (4.21), (5.7), (6.4), (6.8), (6.15), (6.16) and (6.28) (replaced with (6.19) at the Debye sheath entrance instead of the target). The model, like the asymptotic theory it is based on, has been argued to be valid for 
$\alpha \leqslant 5^{\circ }$. The advantage of the model is that the full solution of the quasi-neutrality equation in the magnetic presheath is bypassed, and replaced with constraints derived from quasi-neutrality near the Debye sheath entrance only. The treatment is more accurate for large ion gyro-orbits, 
$\tau = T_{i} / ZT_{e} \gg 1$. Yet, it can also be used for 
$\tau \sim 1$ and reproduces the main qualitative features of distribution functions obtained by solving the self-consistent electrostatic potential across the magnetic presheath (for 
$\alpha \ll 1$), as shown in figures 3, 6 and 7. As the sputtering yield of an ion striking a target depends on the ion's energy and angle of incidence with the target, calculations of energy–angle distributions (7.11) using the model, shown in figures 6 and 7, may be valuable for sputtering predictions.
The narrowing of the wall-normal velocity distribution with the angle 
$\alpha$, shown in figure 3 at the Debye sheath entrance, is explained from the model as follows. Ions reaching the Debye sheath have a minimum normal velocity, 
$V_{x,\textrm {slow}}$, which is related to the size of the gyro-orbit, and so to the adiabatic invariant 
$\mu$. Ions with smaller gyro-orbits have a smaller gyration velocity, and so a smaller magnetic force acts on them to maintain the gyro-motion. Consequently, a weaker electric force is needed to overcome the magnetic force and accelerate these ions towards the target. Ions in smaller gyro-orbits (smaller 
$\mu$) are thus accelerated towards the wall for a larger distance, as shown schematically in figure 1(b,c). However, the dependence of 
$V_{x,\textrm {slow}}$ on the adiabatic invariant 
$\mu$ is weak, as seen in figure 4. As the distribution function exponentially decays with 
$\mu$, the distribution function sharply drops to zero for 
$|v_x |$ below the typical values of 
$V_{x,\textrm {slow}}$, as seen in figure 3. The width of the wall-normal velocity distribution is therefore dominated by the gyro-phase dependence of 
$v_x$ at the target. This dependence is represented, schematically, by pairs of ion trajectories with the same value of 
$\mu$ and 
$U$ in figure 1(b,c). It results in the scaling 
$\langle \tilde {v}_x^{2} \rangle \sim \alpha v_{t,i}^{2}$ for the variance of 
$v_x$.
The orderings (2.1) and (2.2) are required in the asymptotic theory and in the large gyro-orbit model, and are typically well-satisfied in fusion devices except for 
$\alpha \gg \sqrt {Zm_{e} / m_{i}}$ and 
$\rho _{e} \ll \lambda _{D}$. Therefore, a kinetic model (instead of an adiabatic model) for the electrons should be used in the quasi-neutrality equation for 
$\phi (x)$ in the magnetic presheath. This would change the electron contribution to the closure (6.4) (quasi-neutrality) and (6.8) (kinetic Bohm condition) of the large gyro-orbit model.
Acknowledgements
The author is grateful to Felix Parra for stimulating discussions and feedback. This work was supported by the US Department of Energy through grant DE-FG02-93ER-54197.
Editor Paolo Ricci thanks the referees for their advice in evaluating this article.
Declaration of interests
The author reports no conflict of interest.
Appendix A. Change of $U_{\perp } - \chi _{{M}} (\bar {x})$ during the last ion gyro-orbit
In this appendix, the change in the quantity 
$U_{\perp } - \chi _{{M}} (\bar {x})$ during the last gyro-orbit of an ion is calculated. This quantity is denoted by 
$\varDelta _{M} (\bar {x}, U )$, and is responsible for the spread of values of 
$v_x$ in the ion distribution function at the Debye sheath entrance (5.8) and at the wall (5.10).
Recalling from the discussion after (4.8) that 
$\dot {\bar {x}} = - \alpha v_z$, we obtain 
$\dot {\chi }_{M} (\bar {x}) = \dot {x}_{M} \chi '(x_{M}, \bar {x} ) + \dot {\bar {x}} \partial \chi (x_{M}, \bar {x} ) / \partial \bar {x} = \alpha v_z \varOmega ^{2} (x_{M} - \bar {x} )$. Here we have used that 
$\dot {x}_{M} \chi '(x_{M}, \bar {x} ) = 0$ owing to 
$\dot {x}_{M} = 0$ for type I orbits (
$x_{M} = 0$) and 
$\chi ' (x_{M} , \bar {x}) = 0$ for type II orbits. Also recalling 
$\dot {U}_{\perp } = - \alpha \varOmega v_z v_y = \alpha \varOmega ^{2} (x - \bar {x})$, the rate of change of the quantity 
$U_{\perp } - \chi _{{M}} (\bar {x})$ is
\begin{equation} \frac{\textrm{d}}{\textrm{d}t} \left( U_{{\perp}} - \chi_{{M}} (\bar{x}) \right) = \alpha \varOmega^{2} V_{{\parallel}} (\chi_{M} (\bar{x} ), U ) \left( x - x_{{M}} \right) . \end{equation}
This is always positive for closed orbits which satisfy 
$x \geqslant x_{b} \geqslant x_{M}$. Consider an ion, at a position 
$x$, that has just reached values of 
$\bar {x}$ and 
$U_{\perp }$ such that 
$U_{\perp } = \chi _{M}(\bar {x})$. The time taken for the ion to reach 
$x=0$ is approximated by integrating the equation 
${\textrm {d}x}/\textrm {d}t = v_x \simeq \sigma _x \sqrt {2( \chi _{M} - \chi (x, \bar {x} ) )}$, where 
$\sigma _x = \pm 1$ is the sign of 
$v_x$, to obtain
\begin{equation} t = \left(\sigma_x +1 \right) \int_{x}^{x_{t}} \frac{\textrm{d}s}{\sqrt{2\left( \chi_{M} - \chi (s, \bar{x} ) \right)} } + \int_0^{x} \frac{\textrm{d}s}{\sqrt{2\left( \chi_{M} - \chi (s, \bar{x} ) \right)} } . \end{equation}
Here, and in the rest of this appendix, we denote the position 
$x$ by the symbol 
$s$ when under an integral if the symbol 
$x$ is already used for one of the limits of the integration. The problem with the approximation in (A 2) is that the second integral is logarithmically divergent for type II orbits owing to the form of the integrand for 
$s \rightarrow x_{M}$,
\begin{equation} \lim_{s\rightarrow x_{M}} \frac{1}{\sqrt{2\left( \chi_{M} - \chi (s, \bar{x} ) \right)}} = \frac{1}{\sqrt{\chi'' (x_{M} )} \left| s- x_{M} \right| }. \end{equation}
However, the time 
$t$ taken by an ion to reach the Debye sheath entrance from a point in its last gyro-orbit would only be infinite if 
$v_x = \sigma _x \sqrt {2( \chi _{M} - \chi (x, \bar {x} ) )}$ was exactly true. In practice, the quantity 
$U_{\perp } - \chi _{M}(\bar {x})$ is not exactly zero. To calculate this quantity, the time evolution of 
$U_{\perp } - \chi _{M}(\bar {x})$ is estimated in the same way the time 
$t$ was estimated (incorrectly): we replace the time derivative in (A 1) with a spatial derivative using the substitution 
$\textrm {d}/\textrm {d}t = v_x \textrm {d}/{\textrm {d}x}$, and the approximation 
$v_x \simeq \sqrt {2( \chi _{M} - \chi (x, \bar {x} ) )}$ to obtain
\begin{equation} \frac{\textrm{d}}{\textrm{d}x} \left( U_{{\perp}} - \chi_{{M}} \right) = {\pm}\alpha \varOmega^{2} V_{{\parallel}} \left( \chi_{{M}} (\bar{x}), U \right) \frac{ x - x_{{M}} }{ \sqrt{2\left( \chi_{M} - \chi (x, \bar{x} ) \right)} }. \end{equation}This equation is then integrated in the same way as before to obtain
\begin{align} U_{{\perp}} - \chi_{{M}}(\bar{x}) &= \alpha \varOmega^{2} V_{{\parallel}} \left( \chi_{{M}} (\bar{x}), U \right) \left[ \left(\sigma_x +1 \right) \int_{x}^{x_{t}} \frac{s - x_{{M}}}{\sqrt{2\left( \chi_{{M}} - \chi (s, \bar{x} ) \right) }}\,\textrm{d}s \right. \nonumber\\ &\quad + \left. \int_0^{x} \frac{s - x_{{M}}}{\sqrt{2\left( \chi_{{M}} - \chi (s, \bar{x} ) \right) }} \,\textrm{d}s \right] . \end{align}
The second integral in (A 5) is not divergent near 
$x=x_{M}$ because the integrand tends to
\begin{equation} \lim_{s\rightarrow x_{M}} \frac{s-x_{M}}{\sqrt{2\left( \chi_{M} - \chi (s, \bar{x} ) \right)}} = \frac{s-x_{M}}{\sqrt{\chi'' (x_{M} )} \left| s- x_{M} \right| } = \varTheta \left( s-x_{M} \right) \frac{1}{\sqrt{\chi'' (x_{M} )} } , \end{equation}
which is always finite (moreover, the contribution from the region near 
$s=x_{M}$ in the integral (A 5) vanishes because the integrand changes sign there). Considering (A 5), 
$U_{\perp }$ is only ever exactly equal to 
$\chi _{M} ( \bar {x})$ at an instant, and at all other times it is different. Therefore, the time estimated in (A 2) is incorrect, and the divergence in (A 3) does not occur. In practice, ions cross the effective potential maximum in a time 
$t \sim 2{\rm \pi} |\ln \alpha | / \varOmega$ (Geraldini et al. Reference Geraldini, Parra and Militello2018).
Upper and lower bounds for the values of 
$U_{\perp } - \chi _{M}(\bar {x})$ of ions reaching 
$x=0$ can be obtained using the fact that these ions must have past trajectories with a bottom bounce point 
$x_{b}$. We consider the following two limiting cases: (i) an ion crossing the maximum 
$x=x_{M}$ towards the sheath with 
$U_{\perp } = \chi _{M}(\bar {x}) + \epsilon$; (ii) an ion bouncing back (for the last time) from 
$x=x_{M}$ with 
$U_{\perp } = \chi _{M}(\bar {x}) - \epsilon$, where 
$\epsilon$ is an energy difference so small it can be neglected. The minimum value of 
$U_{\perp } - \chi _{M}$ of an ion entering the Debye sheath is calculated from case (i),
\begin{equation} U_{{\perp}} - \chi_{{M}}(\bar{x}) ={-} \varDelta_+ (x, \bar{x}, U ), \end{equation}where
\begin{equation} \varDelta_+ (x, \bar{x}, U ) = \alpha \varOmega^{2} V_{{\parallel}} \left( \chi_{{M}}, U \right) \int_0^{x_{M}} \frac{ x_{{M}} - s}{\sqrt{2\left( \chi_{{M}} - \chi (s, \bar{x} ) \right) }}\ \textrm{d}s \end{equation}
is a positive quantity. Here, we have added to 
$U_{\perp } - \chi _{{M}}(\bar {x}) = 0$ the amount obtained by integrating (A 1) from 
$x_{M}$ to the Debye sheath entrance (
$x\simeq 0$ here). The maximum value of 
$U_{\perp } - \chi _{M} (\bar {x})$ is calculated from case (ii),
\begin{equation} U_{{\perp}} - \chi_{{M}}(\bar{x}) = \varDelta_{M} (\bar{x}, U ) - \varDelta_+ (x, \bar{x}, U ) , \end{equation}where
\begin{equation} \varDelta_{M} (\bar{x}, U ) = 2 \alpha \varOmega^{2} V_{{\parallel}} \left( \chi_{{M}}, U \right) \int_{x_{M}}^{x_{t}} \frac{x - x_{{M}}}{\sqrt{2\left( \chi_{{M}} - \chi (x, \bar{x} ) \right) }} \,{\textrm{d}x}. \end{equation}
Here, we have added to 
$U_{\perp } - \chi _{{M}}(\bar {x}) = 0$ the amount obtained by integrating (A 1) from 
$x_{M}$ to 
$x_{t}$, then back again all the way to the Debye sheath entrance (
$x\simeq 0$). The quantity 
$\varDelta _{+}$ was shown to be negligible when calculating 
$v_x$ from (4.12), as it is always small relative to either 
$\chi _{M} (\bar {x}) - \chi ( x, \bar {x})$ or 
$\varDelta _{M}(\bar {x}, U )$ (Geraldini et al. Reference Geraldini, Parra and Militello2018). Thus, we can consider 
$0 \leqslant U_{\perp } - \chi _{M}(\bar {x}) < \varDelta _{M} (\bar {x}, U )$ for ions reaching the Debye sheath entrance.
Equation (4.22) for 
$\varDelta _{M} (\bar {x}, U )$ follows from (A 10) and from the equality
\begin{equation} \mu_{\textrm{op}}'(\bar{x}) = \varOmega^{2} \int_{x_{M}}^{x_{t}} \frac{x - x_{{M}}}{\sqrt{2\left( \chi_{{M}} - \chi (x, \bar{x} ) \right) }}\,{\textrm{d}x}, \end{equation}which can be verified from (4.20).
Appendix B. Ion conservation
The ion distribution function at the Debye sheath entrance (5.8) is proved here to be consistent with ion conservation in the magnetic presheath. Equation (5.4) gives the current flowing normal to the wall at the magnetic presheath entrance. In steady state, the current flowing normal to the wall at the Debye sheath entrance should be the same. At the Debye sheath entrance, the ion density is small in 
$\alpha$ and the ion current flowing normal to the wall is owing to the component 
$v_x$ of the velocity of all ions,
\begin{align} \frac{j_{{i},x}}{Ze} &={-} 2{\rm \pi}\int_{ \bar{x}_{c}}^{\infty} \varOmega\, \textrm{d} \bar{x} \int_{\varOmega \mu}^{\infty} \frac{ F \left( \mu_{\textrm{op}} (\bar{x}) , U \right)\,\textrm{d}U}{ V_{{\parallel} } ( \chi_{M} (\bar{x} ), U ) }\nonumber\\ &\quad \times \int_{-\infty}^{\infty} \hat{\varPi} \left( \frac{1}{2}v_x^{2} - \chi_{M} ( \bar{x} ) + \frac{1}{2} \varOmega^{2} \bar{x}^{2} + \frac{\varOmega \phi_{\textrm{DSE}}}{B} , 0 , \varDelta_{M}(\bar{x}, U ) \right) v_x\,\textrm{d}v_x. \end{align}
The last integral in 
$v_x$ is taken by replacing 
$v_x \,\textrm {d}v_x = \,\textrm {d}(v_x^{2}/2)$, and the result is 
$\varDelta _{M} (\bar {x}) = 2\alpha {\rm \pi}\mu _{\textrm {op}}'(\bar {x}) V_{\parallel } ( \chi _{M} (\bar {x} ), U )$,
\begin{equation} \frac{j_{{i},x}}{Ze} ={-} 2 \alpha {\rm \pi}\int_{0}^{\infty} \varOmega \,\textrm{d}\bar{x} \mu_{\textrm{op}}'(\bar{x}) \int_{\varOmega \hat{\mu}}^{\infty} F \left(\mu_{\textrm{op}}(\bar{x}) , U \right) \,\textrm{d}U . \end{equation}
Using 
$\mu _{\textrm {op}}'(\bar {x}) = \textrm {d}\mu / \textrm {d}\bar {x}$ and changing integration variable to 
$\mu = \mu _{\textrm {op}}(\bar {x})$ leads to (5.4). The same argument applies to the ion distribution function at the wall (5.10) and to the large gyro-orbit model distribution functions (6.19) and (6.28).
