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

TL;DR

This paper investigates how ion temperature influences the structure and behavior of shallow-angle magnetic presheaths with adiabatic electrons, relevant for fusion device boundary conditions, using kinetic and fluid models.

Contribution

It introduces a kinetic model for magnetic presheaths that captures ion temperature effects and compares it with fluid models, highlighting differences at various temperature regimes.

Findings

Kinetic and fluid models agree at low ion temperature, τ << 1.

Ion gyro-orbits occupy large regions at low τ, affecting potential.

At high ion temperature, τ >> 1, ions reach the wall at very shallow angles.

Abstract

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

Equations

$$\tau =\frac{T_{\text{i}}}{Z{T}_{\text{e}}}\text{,}$$

$$\mathbf{B}=B\cos \alpha \hat{\mathbf{z}}-B\sin \alpha \hat{\mathbf{x}}\text{.}$$

$$\mathbf{E}=-{\varphi }^{\prime }(x)\hat{\mathbf{x}}\text{,}$$

$$e\varphi \sim T_{\text{e}}\text{.}$$

$$v_{\text{B}}=\sqrt{\frac{Z{T}_{\text{e}}}{m_{\text{i}}}}\text{.}$$

$$v_{\text{t,i}}=\sqrt{\frac{2T_{\text{i}}}{m_{\text{i}}}}\text{.}$$

$$|\mathbf{v}|\sim c_s=\sqrt{v_{\text{B}}^2+\frac{1}{2}{v}_{\text{t,i}}^2}=\sqrt{\frac{Z{T}_{\text{e}}+T_{\text{i}}}{m_{\text{i}}}}\text{.}$$

$${\rho }_{\text{s}}=\frac{c_s}{\Omega }\text{,}$$

$${\rho }_{\text{B}}=\frac{v_{\text{B}}}{\Omega }\text{.}$$

$$v_{\text{t,e}}=\sqrt{\frac{2T_{\text{e}}}{m_{\text{e}}}}\text{,}$$

$$\sqrt{\frac{m_{\text{e}}}{m_{\text{i}}}}\sqrt{\tau +1}\ll \alpha \text{.}$$

$${\dot{v}}_x=-\frac{\Omega {\varphi }^{\prime }}{B}+\Omega v_y\cos \alpha \text{,}$$

$${\dot{v}}_y=-\Omega v_x\cos \alpha -\Omega v_z\sin \alpha \text{,}$$

$${\dot{v}}_z=\Omega v_y\sin \alpha \text{.}$$

$${\dot{v}}_x\simeq -\frac{\Omega {\varphi }^{\prime }}{B}+\Omega v_y\text{,}$$

$${\dot{v}}_y\simeq -\Omega v_x-\Omega v_z\alpha \text{,}$$

$${\dot{v}}_z\simeq \Omega v_y\alpha \text{.}$$

$$\bar{x}=x+\frac{1}{\Omega }{v}_y\text{,}$$

$$U_{\perp }=\frac{1}{2}{v}_x^2+\frac{1}{2}{v}_y^2+\frac{\Omega \varphi (x)}{B}\text{,}$$

$$U=\frac{1}{2}{v}_x^2+\frac{1}{2}{v}_y^2+\frac{1}{2}{v}_z^2+\frac{\Omega \varphi (x)}{B}\text{.}$$

$$v_x$$

$$v_y=\Omega \left(\bar{x}-x\right)\text{,}$$

$$v_z=V_{\parallel }\left(U_{\perp },U\right)=\sqrt{2\left(U-U_{\perp }\right)}\text{,}$$

$$\chi (x,\bar{x})=\frac{1}{2}{\Omega }^2{\left(x-\bar{x}\right)}^2+\frac{\Omega \varphi (x)}{B}$$

$${\chi }^{\prime }(x_{\text{m}},\bar{x})={\Omega }^2\left(x_{\text{m}}-\bar{x}\right)+\frac{\Omega {\varphi }^{\prime }(x_{\text{m}})}{B}=0$$

$${\chi }^{\prime \prime }(x_{\text{m}})={\Omega }^2+\frac{\Omega {\varphi }^{\prime \prime }(x_{\text{m}})}{B}>0\text{.}$$

$${\chi }_{\text{m}}(\bar{x})=\chi (x_{\text{m}},\bar{x})=\frac{1}{2}{\left(\frac{{\varphi }^{\prime }(x_{\text{m}})}{B}\right)}^2+\frac{\Omega \varphi (x_{\text{m}})}{B}\text{.}$$

$$\mu ={\mu }_{\text{gk}}(\bar{x},U_{\perp })=\frac{1}{\pi }{\int }_{x_{\text{b}}}^{x_{\text{t}}}{V}_x\left(x,\bar{x},U_{\perp }\right)dx\sim \frac{v_{\text{t,i}}^2}{\Omega }\text{,}$$

$${\rho }_x=x-x_{\text{m}}\sim x_{\text{t}}-x_{\text{b}}\text{,}$$

$$w_x={\dot{\rho }}_x\sim \sqrt{2\left(U_{\perp }-{\chi }_{\text{m}}\right)}\text{.}$$

Full text

1 Introduction

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

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

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

[TABLE]

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

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

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

2 Orderings

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

[TABLE]

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

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

[TABLE]

where $\phi$ is the electrostatic potential and a prime denotes differentiation with respect to $x$. The electrostatic potential is assumed to monotonically converge to some value at $x\rightarrow\infty$, and this value is set to be $\phi=0$. Moreover, it has been shown that $\phi(x)-\phi(0)\propto\sqrt{x}$ at $x\rightarrow 0$ (see equations (141) and (142) in Geraldini et al. (2018)), so that the magnetic presheath electric field diverges at the Debye sheath entrance111This is not a real divergence of the electric field, but is rather a large electric field satisfying $T_{\text{e}}/e\rho_{\text{s}}\ll\phi^{\prime}(0)\ll T_{\text{e}}/e\lambda_{\text{D}}$. See, for example, Riemann (1991) for detailed explanations on the use of asymptotic methods for Debye sheaths and for certain types of presheath.. The co-ordinate system and the geometry are depicted in figure 1.

Since the electric field is present to repel electrons from the wall, the characteristic size of the electrostatic potential $\phi$ is given by

[TABLE]

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

[TABLE]

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

[TABLE]

From equations (5) and (6) it follows that, in general, the ion’s speed has a characteristic size equal to the ion sound speed222Our definition of the ion sound speed is not the most general one, as in fluid treatments this quantity is often defined with an adiabatic constant multiplying the ion temperature. Since the adiabatic constant is normally of order unity in size, the discrepancy in these definitions does not matter. $c_{\text{s}}$,

[TABLE]

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

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

[TABLE]

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

[TABLE]

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

The assumption of an electron-repelling wall is not valid for any value of $\alpha$ and $\tau$. We proceed to obtain the condition on $\alpha$ and $\tau$ for this assumption to be valid. We expect electrons to travel at characteristic velocities equal to their thermal speed,

[TABLE]

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

[TABLE]

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

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

3 Kinetic ion model

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

When written as a function of $\mu$ and $U$, the distribution function is constant across the magnetic presheath, to lowest order in $\alpha$. This is used, in Section 3.2, to write an expression for the ion density. In Section 3.3, we write the quasineutrality equation and summarize the main equations of the shallow-angle kinetic model.

3.1 Ion trajectories in terms of slowly changing orbit parameters

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

[TABLE]

[TABLE]

[TABLE]

Expanding equations (12)-(14) in $\alpha\ll 1$ and neglecting second order terms, we obtain

[TABLE]

[TABLE]

[TABLE]

We introduce three orbit parameters: the orbit position

[TABLE]

the perpendicular energy

[TABLE]

and the total energy

[TABLE]

The orbit parameters vary over a timescale which is longer by a factor of $1/\alpha$ than the timescale $1/\Omega$ over which $x$, $v_{x}$ and $v_{y}$ vary, $\bar{x}/\dot{\bar{x}}\sim U_{\perp}/\dot{U}_{\perp}\sim 1/\alpha\Omega\gg|\mathbf{v}|/|\dot{\mathbf{v}}|\sim 1/\Omega$. The total energy $U$ is exactly constant, $\dot{U}=0$. The instantaneous particle velocities can be expressed in terms of the instantaneous position $x$ and the orbit parameters:

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

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

For times comparable to the typical ion gyroperiod, $2\pi/\Omega$, the orbit parameters are constant to lowest order in $\alpha$ and equations (21)-(24) can be used to infer the approximate particle trajectory. From equation (21), the ion motion is periodic to lowest order in $\alpha$ if, for some $\bar{x}$ and $U_{\perp}$, turning points $x_{\text{b}}$ (bottom) and $x_{\text{t}}$ (top) exist such that: (i) $V_{x}\left(x_{\text{b}},\bar{x},U_{\perp}\right)=V_{x}\left(x_{\text{t}},\bar{x},U_{\perp}\right)=0$ and (ii) $\chi(x,\bar{x})\leqslant U_{\perp}$ in the interval $x_{\text{b}}\leqslant x\leqslant x_{\text{t}}$. Then, the ion will move back and forth between $x_{\text{b}}$ and $x_{\text{t}}$ with period $\sim 2\pi/\Omega$. In order to satisfy (ii), the turning points must lie on either side of an effective potential minimum $x_{\text{m}}$ which, by definition, satisfies

[TABLE]

and

[TABLE]

The value of $\chi$ evaluated at the effective potential minimum is, using equations (24) and (25),

[TABLE]

The ion motion in the $x$ direction (normal to the wall) is exactly periodic for $\alpha=0$, with constant orbit parameters (in the $y$ direction, the motion is a sum of an exactly periodic motion and a constant $\mathbf{E}\times\mathbf{B}$ drift, as explained in Geraldini et al. (2017)). The small angle $\alpha$ perturbs the periodic motion by a small amount, since the orbit parameters become slowly changing in time. Thus, the ion velocity can be approximately decomposed into a periodic piece, with period $\sim 2\pi/\Omega$, and a piece that is approximately constant over the timescale of the periodic motion. Under such circumstances, there is a quantity related to the underlying periodic motion, called an adiabatic invariant, which is a constant of the overall quasi-periodic motion to lowest order in the perturbation parameter333In fact, adiabatic invariants can usually be corrected at every order in such a way that they are conserved to all orders in the perturbation parameter.. The adiabatic invariant in this system is given by

[TABLE]

and is constant to lowest order in $\alpha$. The ordering $\mu\sim v_{\text{t,i}}^{2}/\Omega$ on the far right is obtained in the following way. We define the quantities

[TABLE]

and

[TABLE]

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

3.2 Ion density

Treating the ion motion as periodic to lowest order in some expansion parameter is akin to conventional gyrokinetics (Rutherford & Frieman, 1968; Taylor & Hastie, 1968; Catto, 1978; Antonsen & Lane, 1980; Frieman & Chen, 1982). At every point, the ion’s trajectory can be approximated to lowest order by a periodic orbit whose period is faster than any other timescale of interest. As in gyrokinetic theory, the ion distribution function can be shown to be independent of the fast timescale to lowest order in $\alpha$. Moreover, since $\mu$ and $U$ are both constants of the perturbed motion (at least to lowest order in $\alpha$), the distribution function written in terms of the variables $\mu$ and $U$, $F(\mu,U)$, can be shown to be constant across the magnetic presheath (Cohen & Ryutov, 1998; Geraldini et al., 2017). Therefore, the function $F(\mu,U)$ is completely determined by ions entering the magnetic presheath at $x\rightarrow\infty$. In order to write $F(\mu,U)$ from the distribution function at $x\rightarrow\infty$ expressed in terms of $\mathbf{v}$, denoted $f_{\infty}(\mathbf{v})$, we use the equations

[TABLE]

and

[TABLE]

These equations are obtained by setting $\phi=0$ in equations (20) and (28), and are thus valid at $x\rightarrow\infty$. Note that the self-consistent form of $f_{\infty}(\mathbf{v})$ should be independent of the gyrophase angle, which at $x\rightarrow\infty$ is $\tan^{-1}(v_{x}/v_{y})$.

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

[TABLE]

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

[TABLE]

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

[TABLE]

In equations (33) and (3.2), we have introduced several quantities which are derived and explained in detail in Geraldini et al. (2017) and Geraldini et al. (2018), and we have assumed that $\phi(x)$, $\phi^{\prime}(x)$ and $\phi^{\prime\prime}(x)$ are all monotonic functions of $x$. The minimum allowed orbit position $\bar{x}_{\text{m}}$ for an ion at position $x$ to be in an orbit that is periodic to lowest order in $\alpha$ is

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

where

[TABLE]

Equation (41) is derived in Appendix B from the expression for $\Delta_{\text{M}}$ given in Geraldini et al. (2018).

3.3 Quasineutrality and summary of equations

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

[TABLE]

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

[TABLE]

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

[TABLE]

Equation (45) is used to determine the self-consistent electrostatic potential $\phi(x)$ across the magnetic presheath. A condition that must be satisfied in order for equation (45) to have a solution is (Geraldini et al., 2018)

[TABLE]

which we refer to as the kinetic Chodura condition.

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

[TABLE]

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

[TABLE]

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

[TABLE]

where $\bar{x}=v_{y}/\Omega$ at $x=0$, and $\hat{\Pi}$ is the top-hat function defined by

[TABLE]

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

[TABLE]

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

[TABLE]

Equation (3.3) is obtained by integrating (49) over $\bar{x}$ and $U$, without integrating over $v_{x}$. Equation (3.3) is obtained by integrating (49) over $v_{x}$ and re-expressing the distribution as a function of $v_{y}$ and $v_{z}$, using $v_{y}=\Omega\bar{x}$ (valid at $x=0$) and $v_{z}=\sqrt{2\left(U-\chi_{\text{M}}(\bar{x})\right)}$. In Geraldini et al. (2018) it was shown that the equation

[TABLE]

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

4 The limits of small and large ion temperature

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

4.1 Cold ions ($\tau\ll 1$)

In this subsection, we argue that our kinetic model is equivalent to Chodura’s fluid model, which is valid for $\tau=0$, in an appropriate limit for $\tau\ll 1$.

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

[TABLE]

The choice (54) corresponds to flow parallel to the magnetic field satisfying the Chodura (or Bohm-Chodura) condition (Chodura, 1982) with the equality sign. Using equation (47) and (54), the ion fluid velocity at every position can be written in terms of the electrostatic potential at that position,

[TABLE]

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

[TABLE]

Equation (4.1) was originally derived in Chodura (1982) (and later in Riemann (1994)), in terms of $u_{x}$ instead of $\phi$, and is valid for all values of $\alpha$ (provided that $\alpha\gg\sqrt{m_{\text{e}}/m_{\text{i}}}$ as discussed in section 2).

For $\alpha\ll 1$, the relationship between electrostatic potential and fluid velocity, equation (55), simplifies to

[TABLE]

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

[TABLE]

In Appendix C.2, we expand equation (4.1) for $\alpha\ll 1$, thus obtaining an equation for the electrostatic potential for $\tau=0$ and $\alpha\ll 1$,

[TABLE]

In order to obtain a correspondence between our kinetic model and Chodura’s fluid model, we define a new expansion parameter, $\epsilon\equiv 1/|\ln\alpha|$, and take the ordering

[TABLE]

The ordering (60) restricts $\alpha$ (and $\tau$) to be exponentially small, $\alpha=\exp(-1/\epsilon)$. However, with the ordering (60), the kinetic model is asymptotically equivalent to the fluid model with an exponentially small error in $\epsilon$ (or, equivalently, a small error in $\alpha$ and $\tau$). Hence, since the error is so small, in practice $\alpha$ need not be excessively small (we require $\alpha<0.1$).

Before analyzing the kinetic model for $\tau\ll 1$ using the ordering (60), it is instructive to solve equation (59) explicitly for $\epsilon\ll 1$. Using the boundary condition (58), we order $e\phi/T_{\text{e}}\sim 1/\epsilon$ in the magnetic presheath. Then, we expand equation (59) in $\epsilon\ll 1$ to obtain

[TABLE]

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

[TABLE]

From equation (62), the length scale of the magnetic presheath is $\sim\rho_{\text{B}}/\sqrt{\epsilon}$.

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

• •

a region far from the wall,

[TABLE]

where all ions are in small approximately periodic orbits (closed orbits);

• •

a region close to the wall,

[TABLE]

where all ions are in open orbits;

• •

an intermediate region,

[TABLE]

where ions moving towards the wall transition from small closed orbits to larger, distorted closed orbits, and finally to open orbits.

In subsections 4.1.1-4.1.3, we study the three regions in the order listed above. Instead of taking the limit $\tau\ll 1$ of equations (33) and (3.2) directly, which we leave to Appendix E, in subsections 4.1.1 and 4.1.2 we derive the flow velocity of ions in closed and open orbits. For ions in closed orbits, the flow velocity is much smaller than the particle velocity, as most of the particle velocity is periodic in time and gives no contribution to the flow (because the periodic motion is averaged over); hence, the flow velocity is equal to the drift velocity of the ion gyro-orbits. However, for ions in open orbits sufficiently close to the wall the motion has no periodic piece, and so the flow velocity is equal to the individual particle velocity. From the flow velocity $u_{x}$ and equation (57), we derive equations for the electrostatic potential $\phi$ in the regions (63) and (64). To lowest order in $\alpha$ and $\tau$, the solution for the electrostatic potential in the part of these two regions that overlaps with the intermediate region (65) is a parabola. Therefore, we assume that the lowest order solution for $\phi(x)$ in the whole intermediate region (65) is a parabola, and use this to write an approximate kinetic quasineutrality equation for the region (65). Finally, in subsection 4.1.4 we write an approximate differential equation for the electrostatic potential, whose solution is (59), and show that it is equivalent to the equations describing the electrostatic potential in the three regions.

4.1.1 Far from the wall

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

[TABLE]

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

[TABLE]

Consider an ion moving in an effective potential given by (66). The turning points $x_{\text{b}}$ and $x_{\text{t}}$ are solutions of the equation $U_{\perp}=\chi(x,\bar{x})$, and so

[TABLE]

Recalling the definitions and orderings in (29)-(30), equation (68) corresponds to $w_{x}^{2}\sim\chi^{\prime\prime}(x)\rho_{x}^{2}$. Thus, we obtain the ordering $\rho_{x}\sim w_{x}/\sqrt{\chi^{\prime\prime}(x)}$ relating the typical spatial extent in the $x$ direction (normal to the wall) of the ion orbit, $\rho_{x}$, to the typical orbital velocity component in the same direction, $w_{x}$. Note that $w_{x}/\rho_{x}\sim\sqrt{\chi^{\prime\prime}(x_{\text{m}})}$ is the characteristic gyrofrequency of the approximately periodic motion of the ion. This is consistent with the elliptical gyro-orbits studied in the Appendix of Geraldini et al. (2017). Moreover, from equation (28) we have the relationship $\mu\sim w_{x}\rho_{x}\sim v_{\text{t,i}}^{2}/\Omega\sim\tau v_{\text{B}}^{2}/\Omega$, from which we obtain the estimates

[TABLE]

and

[TABLE]

The electrostatic potential $\phi$ given in equation (62) has a discontinuous second derivative: for $x/\rho_{\text{B}}\geqslant\sqrt{2/\epsilon}$, we have $\phi^{\prime\prime}(x)\simeq 0$ and $\chi^{\prime\prime}(x)\simeq\Omega^{2}$, while for $x/\rho_{\text{B}}<\sqrt{2/\epsilon}$ we have $\phi^{\prime\prime}(x)\simeq-\Omega^{2}$ and $\chi^{\prime\prime}(x)\simeq 0$. Hence, to lowest order in $\epsilon$, the second derivative of the electrostatic potential is not determined (note that equation (4.1) does not specify $\phi^{\prime\prime}(x)$). However, the abrupt jump in the value of $\phi^{\prime\prime}(x)$ and $\chi^{\prime\prime}(x)$ occurring at $x/\rho_{\text{B}}=\sqrt{2/\epsilon}$ is a reflection of a decrease of $\phi^{\prime\prime}(x)$ and $\chi^{\prime\prime}(x)$ in going from $x\rightarrow\infty$ to $x/\rho_{\text{B}}<\sqrt{2/\epsilon}$. From equation (70), the size of ion orbits is $\rho_{x}\sim\sqrt{\tau}\rho_{\text{B}}\sim\rho_{\text{i}}$ when $\chi^{\prime\prime}(x_{\text{m}})\simeq\Omega^{2}$. Conversely, when $\chi^{\prime\prime}(x_{\text{m}})\ll\Omega^{2}$, the spatial extent of the ion orbits is larger, $\rho_{x}\gg\rho_{\text{i}}$. The growth of the ion orbit as it approaches the wall in the magnetic presheath is shown in figure 3. Note that as $\rho_{x}$ becomes larger, the typical orbital velocity $w_{x}\simeq v_{x}$ becomes smaller (see equation (69)). When $\rho_{x}$ grows so large that $\rho_{x}\sim l$, equations (66) and (68)-(70) cease to be valid as the effective potential can no longer be Taylor expanded near its minimum. This happens when the ion reaches the shaded region in figure 2.

We proceed to solve for the ion motion by assuming that $\rho_{x}$ is small,

[TABLE]

From equation (25) and (71), we obtain the value of $v_{y}=\Omega(\bar{x}-x)$,

[TABLE]

Indeed, since the orbital velocity is small, by equation (69), and the angle between the magnetic field and the wall is shallow, the motion of the ion is approximately parallel to the wall. Thus, the magnetic force away from the wall, $ZeBv_{y}$, is approximately equal to the electric force towards the wall, $Ze\phi^{\prime}(x)\simeq Ze\phi^{\prime}(x_{\text{m}})$. Using equations (19) and (72), the perpendicular energy of an ion at a position $x\simeq x_{\text{m}}$ is given by

[TABLE]

The first error in (73) is a combination of the orbital component of $\frac{1}{2}v_{y}^{2}$, $\Omega^{2}\rho_{x}^{2}/2$, and the quadratic term, $\Omega\phi^{\prime\prime}(x_{\text{m}})\rho_{x}^{2}/2B$, of the Taylor expansion of $\Omega\phi(x)/B$ near $x=x_{\text{m}}$. Note that the term $O\left(\phi^{\prime}(x_{\text{m}})\Omega\rho_{x}/B\right)$ arising from taking the square of equation (72) has cancelled with the linear term of the Taylor expansion of $\Omega\phi(x)/B$. The second error in (73) comes from neglecting $\frac{1}{2}v_{x}^{2}$. From equations (69) and (70), the two errors have the same size. From equation (54), the total energy of an ion at $x\rightarrow\infty$ is given by $U=v_{\text{B}}^{2}/2$. Hence, the $z$-component of the ion velocity is, using equations (23) and (73) with $U=v_{\text{B}}^{2}/2$,

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

Taking the derivative of (72), we obtain

[TABLE]

Inserting equations (74), (76), and (79) into (78), the terms proportional to $w_{x}$ on the left and right hand side cancel and we obtain an implicit equation for $v_{\text{d}}$,

[TABLE]

The right hand side of equation (80) consists of the small component of parallel streaming in the $x$ direction, approximately given by $-\alpha v_{z}$, a polarization drift, approximately given by $-v_{\text{d}}\phi^{\prime\prime}/\Omega B$, and the error term coming from the error in $v_{z}$. By manipulating (80), we obtain

[TABLE]

An alternative procedure to derive equation (81) is to obtain the time derivative of $x_{\text{m}}$ by using the chain rule, $v_{\text{d}}=\dot{x}_{\text{m}}=\dot{\bar{x}}dx_{\text{m}}/d\bar{x}$, as shown in Appendix D. Equation (81) is divergent for $\phi^{\prime\prime}(x_{\text{m}})=-\Omega B$, but approximating the ion motion as a closed orbit becomes invalid close to the divergence, as it requires $v_{\text{d}}$ to be small by equation (77).

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

[TABLE]

The $O(w_{x}^{2}/v_{\text{B}}^{2})$ error in equation (81) is neglected in (82) as it is smaller than the $O\left(\rho_{x}^{2}/l^{2}\right)$ error.

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

[TABLE]

where the error has been rewritten using the ordering (70). This result can also be derived by taking the limit $\tau\ll 1$ in equation (33), which is a more direct though perhaps less intuitive approach (see Appendix E). We can substitute either of (82) or (83) into (57) or (45), respectively, to obtain a differential equation for the electrostatic potential,

[TABLE]

Recalling that $1+\Omega\phi^{\prime\prime}(x)/B=\chi^{\prime\prime}(x)/\Omega^{2}$, the ordering that results from equation (84) is $\chi^{\prime\prime}(x)/\Omega^{2}\sim\exp\left(e\phi/T_{\text{e}}\right)$, which also leads to the ordering $1/l\sim\chi^{\prime\prime\prime}/\chi^{\prime\prime}\sim e\phi^{\prime}/T_{\text{e}}$. Moreover, using the fact that balancing the terms in the denominator of equation (84) gives $\rho_{\text{B}}^{2}\left(e\phi^{\prime}/T_{\text{e}}\right)^{2}\sim e\phi/T_{\text{e}}$, we obtain the ordering

[TABLE]

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

[TABLE]

Multiplying equation (86) by $e\phi^{\prime}/T_{\text{e}}$, integrating once and using the boundary condition $\phi=\phi^{\prime}=0$ at $x\rightarrow\infty$, gives

[TABLE]

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

[TABLE]

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

For $-e\phi/T_{\text{e}}\gg 1$, all the terms on the right hand side of equation (88) become small except for the $O(\tau)$ term which approaches a constant, and the solution approaches the parabola

[TABLE]

Here, $C$ is a constant determined by boundary conditions, and we denoted the constant $O(\tau)$ error coming from the right hand side of (88) as $\kappa\tau$, where $\kappa$ is an unknown constant of order unity. The electrostatic potential at the wall is large, $-e\phi(0)/T_{\text{e}}=|\ln\alpha|=1/\epsilon\gg 1$, and so we expect equation (89) to become valid closer to the wall. If we assume that (89) is valid at $x=0$ to lowest order in $\epsilon$ and impose $-e\phi(0)/T_{\text{e}}=|\ln\alpha|=1/\epsilon$, we obtain

[TABLE]

To lowest order in $\epsilon$, equation (89), with $C$ given by (90), is equivalent to equation (62), which was obtained from the fluid model. However, note that the non-constant piece of the error in equation (88) becomes comparable to the first term on the right hand side when $\exp\left(e\phi/T_{\text{e}}\right)\sim\tau\left(e\phi/T_{\text{e}}\right)\exp\left(e\phi/2T_{\text{e}}\right)$. Hence, equation (88) fails to correctly determine the potential when $\exp\left(e\phi/T_{\text{e}}\right)\sim\tau^{2}|\ln\tau|^{2}\sim\tau^{2}/\epsilon^{2}$. From equation (70), with the ordering $\chi^{\prime\prime}/\Omega^{2}\sim\exp(e\phi/T_{\text{e}})$, and equation (85), this value of $\exp\left(e\phi/T_{\text{e}}\right)$ corresponds to $\rho_{x}\sim l\sim\sqrt{\epsilon}\rho_{\text{B}}$, which is the point at which the approximation in (66) ceases to be valid. The validity of (88) is thus restricted to

[TABLE]

Note that, from (89), (90) and (91), the validity region is given by (63) to lowest order in $\epsilon$. Ion gyro-orbits grow in size as they approach the wall, as shown in figure 3, making the treatment of this section invalid for $x/\rho_{\text{B}}\leqslant\sqrt{2/\epsilon}-\sqrt{4|\ln\tau|}$, where $\rho_{x}$ is no longer small.

4.1.2 Near the wall

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

[TABLE]

where we used $d\left(\chi_{\text{M}}-\chi(x,\bar{x})\right)/d\bar{x}=\Omega^{2}\left(x-x_{\text{M}}\right)\sim\Omega^{2}\rho_{x}$ and estimated $\left(d\mu/d\bar{x}\right)_{\text{open}}\sim\Omega^{2}\rho_{x}^{2}/w_{x}\sim\Omega\rho_{x}^{3}/\tau\rho_{\text{B}}^{2}$ from equation (158). Since typical ion orbits have values of $\mu$ differing by $O(\tau v_{\text{B}}^{2}/\Omega)$, the values of $\chi_{\text{M}}(\bar{x})-\chi(x,\bar{x})$ of such orbits change by $O(\tau^{2}v_{\text{B}}^{2}\rho_{\text{B}}^{2}/\rho_{x}^{2})$. Recall that an open orbit has $U_{\perp}-\chi(x,\bar{x})=\chi_{\text{M}}(\bar{x})-\chi(x,\bar{x})+O(\Delta_{\text{M}})$: from equation (41) and the previous estimate for $\left(d\mu/d\bar{x}\right)_{\text{open}}$, we obtain the scaling $\Delta_{\text{M}}\sim\alpha\Omega v_{\text{B}}\rho_{x}^{3}/\tau\rho_{\text{B}}^{2}$. If the condition

[TABLE]

is satisfied, with $A$ a positive constant, the $O(\Delta_{\text{M}})$ error term is small in $\alpha$ compared to the $O(\tau^{2}v_{\text{B}}^{2}\rho_{\text{B}}^{2}/\rho_{x}^{2})$ error. Then, using equation (21) with $U_{\perp}-\chi(x,\bar{x})\simeq\chi_{\text{c}}-\chi(x,\bar{x}_{\text{c}})+O(\tau^{2}v_{\text{B}}^{2}\rho_{\text{B}}^{2}/\rho_{x}^{2})$, the velocity of an ion in an open orbit near the wall is

[TABLE]

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

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

[TABLE]

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

[TABLE]

In Appendix E.2, we derive equation (96) by expanding $n_{\text{i,op}}(x)$ in equation (3.2) to lowest order in $\tau\ll 1$. In order for equation (96) to be valid, we require the error term to be small, implying $n_{\text{i,op}}\ll\left(\alpha n_{\infty}/\tau\right)\left(\rho_{x}/\rho_{\text{B}}\right)$. Moreover, since $\rho_{x}$ here quantifies the characteristic size of closed orbits while the ion is transitioning from a closed to an open orbit, the ion density changes from $n_{\text{i,cl}}\gg n_{\infty}\tau^{2}/\epsilon^{2}$ (recall equation (91) and the fact that $n_{\text{i,cl}}\sim n_{\infty}\exp\left(e\phi/T_{\text{e}}\right)$) to $n_{\text{i,op}}\ll\left(\alpha n_{\infty}/\tau\right)\left(\rho_{x}/\rho_{\text{B}}\right)$ over a length scale of $\rho_{x}$. If the validity condition (93) is satisfied, this drop in density corresponds to a decrease of $\ln\left(\tau^{3}/\alpha\right)\sim 1/\epsilon$ in the normalized electrostatic potential $e\phi/T_{\text{e}}$. In order to be consistent with the lowest order electric field obtained from (89) and (90), $e\phi^{\prime}/T_{\text{e}}=-(x-C)/\rho_{\text{B}}^{2}\sim 1/\sqrt{\epsilon}\rho_{\text{B}}$, transitioning ions must have $\rho_{x}\sim\rho_{\text{B}}/\sqrt{\epsilon}$.

Inserting equation (96) into (45) with $n_{\text{i,cl}}(x)=0$, or inserting equation (95) into (57), we obtain

[TABLE]

The constants

[TABLE]

and $\bar{x}_{\text{c}}$ are to be determined; they are related by the boundary condition $e\phi(0)/T_{\rm{e}}=\ln\alpha$, giving

[TABLE]

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

[TABLE]

For equation (100) to be valid, we require $\alpha^{2}\exp\left(-2e\phi/T_{\text{e}}\right)\sim\alpha^{2}\exp\left((x-\bar{x}_{\text{c}})^{2}/\rho_{\text{B}}^{2}\right)\sim\exp\left(-x\bar{x}_{\text{c}}/\rho_{\text{B}}^{2}+x^{2}/2\rho_{\text{B}}^{2}\right)\ll 1$, where we have used $\bar{x}_{\text{c}}\simeq\rho_{\text{B}}\sqrt{2/\epsilon}$; hence, the electrostatic potential becomes well-approximated by (100) for $x\gg\rho_{\text{B}}\sqrt{\epsilon}$. The derivation of equation (97) fails when $n_{\text{i,op}}\sim\alpha n_{\infty}/\sqrt{\epsilon}\tau$ (recall equation (96) with $\rho_{x}\sim\rho_{\text{B}}/\sqrt{\epsilon}$), and so the validity of (97) is restricted to

[TABLE]

From equations (100) and (101), we obtain the estimate (64) for the region where (97) is valid. Outside of the validity region (64), the velocity of a typical ion is of the order of the gyration velocity, $\sqrt{2\left(\chi_{\text{M}}(\bar{x})-\chi(x,\bar{x})\right)}\sim w_{x}$, and so the assumption that all ions are in open orbits is invalid.

4.1.3 Intermediate region

With the ordering (93), there is a finite region where equations (88) and (97) are not valid: from equations (63) and (64), this region is

[TABLE]

However, the solution of equation (88) tends to (89) for $\sqrt{2/\epsilon}-x/\rho_{\text{B}}\gg 1$ and (97) tends to (100) for $x/\rho_{\text{B}}\gg\sqrt{\epsilon}$. Hence, we proceed by assuming that in the intermediate region (65), which includes the region (102), the electrostatic potential is simultaneously given by the parabolas in equations (89) and (100) to lowest order in $\alpha$ and $\tau$. This provides the value of $K$, $K=\kappa\tau\ll 1$, and the equality $C=\bar{x}_{\text{c}}$. Using equation (99) with $K\simeq 0$, the value of $C$ and $\bar{x}_{\text{c}}$ is

[TABLE]

The neglected term $\kappa\tau$ causes a small constant correction to the value of $C$, as we had claimed in the discussion following equation (87).

From equation (24), the effective potential curves associated with the parabolic electrostatic potential of equation (89) are a set of straight lines,

[TABLE]

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

[TABLE]

must be retained. In equation (105), $\phi(x)$ is the solution to the quasineutrality equation (45) for a given value of $\tau$ and $\alpha$. In the intermediate region we take $\phi(x)\simeq\phi_{\text{p}}(x)$ and calculate $\phi_{\text{np}}(x)$ as a higher order asymptotic correction from the following equation:

[TABLE]

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

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

[TABLE]

where we re-expressed the lowest order inequality $x/\rho_{\text{B}}<\sqrt{2/\epsilon}$ to the more precise form $1\ll\sqrt{2/\epsilon}-x/\rho_{\text{B}}$ in order to emphasize the necessity of the ordering $|\ln\tau|\sim 1/\epsilon\gg 1$. We proceed to calculate $\phi_{\text{np}}(x)$ in this region. Inserting $\phi=\phi_{\text{p}}+\phi_{\text{np}}$ in (84) and rearranging the error term, we obtain

[TABLE]

where we have used $\phi\simeq\phi_{\text{p}}$ in the denominator to get $\sqrt{v_{\text{B}}^{2}-\left(\phi^{\prime}(x)/B\right)^{2}-2\Omega\phi(x)/B}\simeq 2v_{\text{B}}$. From the definition of $l$ in (67) and using $\chi^{\prime\prime}=\Omega\phi_{\text{np}}^{\prime\prime}/B$ with equation (108), we obtain $l\sim\rho_{\text{B}}^{2}/(C-x)\sim\rho_{\text{B}}\sqrt{\epsilon}$, consistent with our previous estimate for $l$ in this region (before equation (91)). Integrating (108) twice and imposing $\phi^{\prime}_{\text{np}}(x)=\phi_{\text{np}}(x)=0$ at $(C-x)/\rho_{\text{B}}\rightarrow\infty$ (where the electrostatic potential becomes more parabolic) gives

[TABLE]

where we have used that the double integral of the term $O\left((\tau\rho_{\text{B}}/l^{2})\sqrt{e\phi_{\text{np}}^{\prime\prime}/T_{\text{e}}}\right)$ is $O\left((\tau\rho_{\text{B}}/l)\sqrt{e\phi_{\text{np}}/T_{\text{e}}}\right)$. We proceed to consider the part of the intermediate region that is closest to the wall. Equation (97) is valid in the region (64) near the wall and equation (106) is valid in the intermediate region (65). Hence, the region in which these two equations are both valid is the overlap of (64) and (65),

[TABLE]

From equation (97) and using $e\phi_{\text{np}}/T_{\text{e}}\ll 1$ (the assumption behind equation (106)), we extract

[TABLE]

In the region (102), equation (106) cannot be simplified further. Ion orbits are large, $\rho_{\text{B}}\sqrt{\epsilon}\lesssim\rho_{x}\lesssim\rho_{\text{B}}/\sqrt{\epsilon}$, and the non-parabolic piece of the electrostatic potential is small, $-\tau^{2}\epsilon\lesssim e\phi_{\text{np}}(x)/T_{\text{e}}\lesssim\tau^{2}/\epsilon$ (consistent with the errors in (4.1.3) and (111), noting that $\phi_{\text{np}}$ is positive in (4.1.3) and negative in (111)). From the previous paragraph we deduce that $l\sim\sqrt{\epsilon}\rho_{\text{B}}$, and so the double derivative of the effective potential is small, $-\tau^{2}\lesssim\chi^{\prime\prime}(x)/\Omega^{2}=\rho_{\text{B}}^{2}e\phi_{\text{np}}^{\prime\prime}/T_{\text{e}}\lesssim\tau^{2}/\epsilon^{2}$. The characteristic size of the periodic piece of the ion velocity is obtained from the relation $\mu\sim\rho_{x}w_{x}\sim\tau v_{\text{B}}^{2}/\Omega$ (which holds provided the adiabatic invariant is still conserved), giving $\tau\sqrt{\epsilon}\lesssim w_{x}/v_{\text{B}}\lesssim\tau/\sqrt{\epsilon}$. The characteristic size of the drift of closed ion orbits is obtained from equation (81), giving $\alpha\epsilon^{2}/\tau^{2}\lesssim v_{\text{d}}/v_{\text{B}}\lesssim\alpha/\tau^{2}$. Ignoring the dependences on $\epsilon$ and focusing only on how quantities scale with $\alpha$ and $\tau$ (which, from equation (60), are both exponentially small in $\epsilon$), we obtain $\rho_{x}\sim l\sim\rho_{\text{B}}$, $w_{x}\sim\tau v_{\text{B}}$ and $v_{\text{d}}\sim\alpha v_{\text{B}}/\tau^{2}$. Hence, the condition $v_{\text{d}}/l\ll w_{x}/\rho_{x}$ (equation (77)), which is necessary for the ion motion to be approximately periodic, implies that $v_{\text{d}}/w_{x}\sim\alpha/\tau^{3}\ll 1$. In equation (93) $\alpha/\tau^{3}$ is required to be small in $\alpha$ to ensure that the motion remains periodic to lowest order in $\alpha$. The scaling with $\tau^{3}$ in (93) implies that the kinetic model is not valid for relatively large values of $\tau\ll 1$. This unfortunate scaling arises because of the growth of the ion orbits: if small ion orbits reached $x=0$, we would expect $l\sim\rho_{x}\sim\rho_{\text{i}}$, $\chi^{\prime\prime}(x)\sim\Omega^{2}$, $w_{x}\sim\sqrt{\tau}v_{\text{B}}$, $v_{\text{d}}\sim\alpha v_{\text{B}}$, and so the condition (77) would give the weaker requirement $\alpha\ll\sqrt{\tau}$ for the model to be valid near $x=0$. Hence, the orbit growth and the associated large polarization drift have a strong negative effect on the condition for validity of the model, multiplying the power by which $\tau$ is raised by a factor of six. It is for this reason that, as we will see in section 5, we do not obtain numerical solutions of equation (45) for values of $\tau$ lower than $\tau=0.2$.

4.1.4 Uniformly valid solution

We proceed to obtain an expression for $\phi(x)$, equation (59), that is uniformly valid across the whole magnetic presheath to lowest order in $\alpha$ and $\tau$. In order to do this, we first make a change of variables: guided by the form of (97), we introduce the function

[TABLE]

The term $-\alpha^{2}/2$ is small but is included in the definition (112) in order to have the desirable exact property that $\phi=0$ when $\psi=0$: then, far from the wall, where $-e\phi/T_{\text{e}}\ll 1/\epsilon$, the relation $\psi=\left(e\phi/T_{\text{e}}\right)\left(1+O(\alpha^{2})\right)$ is satisfied.

We proceed to show that the equation

[TABLE]

is equivalent to the equations describing the electrostatic potential in the three regions of the magnetic presheath. All the errors on the right hand side of (113) are exactly equal to zero at $\psi=0$. Moreover, the $O(\tau)$ error tends to a constant for $-\psi\gg 1$, and this constant does not have an effect on the functional form of $\psi(x)$. Hence, we retain the smaller errors $O(\tau\psi\exp\left(\psi/2\right),\tau^{2})$. First, we compare equation (113) with equation (88), valid in the region (63). Since, from (91), $\alpha^{2}\exp(-2e\phi/T_{\text{e}})\ll\alpha^{2}/\tau^{4}$ in this region, it follows that $e\phi/T_{\text{e}}=\psi+O(\alpha^{2}/\tau^{4})$. Hence, equation (113) directly follows from (88), after noting that the $O(\alpha^{2}/\tau^{4})$ error term is smaller than the $O(\tau^{2})$ error term (the smallest in (113)) because of the validity condition (93). Next, we compare the solution to equation (113) with equation (97) (recall that $K=\kappa\tau$, where $\kappa$ is an unknown constant of order unity, and $\bar{x}_{\text{c}}=C$), valid in the region (64) close to the wall. From (112), in this region $-e\phi/T_{\text{e}}=-\psi+O(1)\sim 1/\epsilon$ and so, from (101), $\exp(\psi/2)\ll\sqrt{\alpha/\tau}\ll\tau$. Therefore, all terms on the right hand side of (113), except for the constant $O(\tau)$ term, become smaller than the $O(\tau^{2})$ term. Re-expressing the constant $O(\tau)$ error term as $\kappa\tau$ and integrating gives

[TABLE]

The size of the error terms in (114) is obtained as follows. From $\psi\sim 1/\epsilon$ and $l\sim\psi^{\prime\prime}/\psi^{\prime\prime\prime}\sim\rho_{\text{B}}^{2}/(x-C)\sim\sqrt{\epsilon}\rho_{\text{B}}$, the error terms in the expression for $\psi^{\prime}$ are smaller by a factor of $\sqrt{\epsilon}$ compared with the error terms in the expression for $\psi^{\prime 2}$, equation (113). Upon integrating the expression for $\psi^{\prime}$ to obtain $\psi$, we multiply these error terms by a factor of $l\sim\sqrt{\epsilon}\rho_{\text{B}}$. Using (112), and remembering that in the region (64) the size of the largest error term in (114) is $O(\epsilon\tau^{2})$, observe that equation (114) is equivalent to (97).

It only remains to be shown that (113) is valid in the part of the intermediate region (102) where neither (88) (the equation determining $\phi$ far from the wall) nor (97) (the equation determining $\phi$ close to the wall) are valid. From the equation $\phi=\phi_{\text{p}}+\phi_{\text{np}}$, equation (89) for $\phi_{\text{p}}$, and the fact that $e|\phi_{\text{np}}|/T_{\text{e}}\lesssim\tau^{2}/\epsilon$, the electrostatic potential in this region is given by $e\phi/T_{\text{e}}=-3/2+\kappa\tau-\left(x-C\right)^{2}/2\rho_{\text{B}}^{2}+O(\tau^{2}/\epsilon)$. From equation (91) for the size of the electrostatic potential in the region far from the wall, and equation (101) for the size of the electrostatic potential in the region close to the wall, we obtain the ordering $\alpha/\sqrt{\epsilon}\tau\lesssim\exp(e\phi/T_{\text{e}})\lesssim\tau^{2}/\epsilon^{2}$ for the size of the electrostatic potential in this region. Hence, the first error term in equation (114) becomes $O(\tau^{2}/\epsilon)$. Moreover, from $\epsilon^{4}\alpha^{2}/\tau^{4}\lesssim\alpha^{2}\exp(-2e\phi/T_{\text{e}})\lesssim\epsilon\tau^{2}$ we obtain $e\phi/T_{\text{e}}=\psi+O(\epsilon\tau^{2})$. Hence, from equation (114) and the associated $O(\tau^{2}/\epsilon)$ error in this region, we obtain $e\phi/T_{\text{e}}=-3/2+\kappa\tau-\left(x-C\right)^{2}/2\rho_{\text{B}}^{2}+O\left(\tau^{2}/\epsilon\right)$. Equation (113) is thus a good approximation also in the region (102).

Using the definition (112), and equation (58), the boundary condition at the wall is $\psi(0)=\ln\alpha+1/2-\alpha^{2}/2$. This can be used to integrate equation (113) (neglecting the error terms) and obtain the approximate electrostatic potential solution, as in equation (59).

4.2 Hot ions ($\tau\gg 1/\alpha$)

In the limit of very hot ions, $\tau\gg 1/\alpha$, we assume that the ion distribution function is a half-Maxwellian at the magnetic presheath entrance,

[TABLE]

where we introduced the Heaviside step function

[TABLE]

Since $\phi(\infty)=0$, $U=|\mathbf{v}|^{2}/2$ and we re-express (115) to

[TABLE]

Equation (115) is one of many choices that could be made. The reason we choose this distribution function is that it was also used in Cohen & Ryutov (1998) in the equivalent limit of small electron temperature. We consider the limit $\tau\ll\alpha^{2}m_{\text{i}}/m_{\text{e}}$ in order to be consistent with condition (11) for an electron repelling sheath.

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

[TABLE]

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

[TABLE]

The minimum value of $\bar{x}$ necessary for an ion at position $x$ to be in a closed orbit or an open orbit is, using equations (36) and (37) with $\phi(x)=0$,

[TABLE]

Moreover, the adiabatic invariant is $\mu=U_{\perp}/\Omega$.

Inserting the distribution function (117) into equation (33), the closed orbit density is

[TABLE]

Changing variables to $\tilde{v}_{y}=\left(\bar{x}-x\right)/\rho_{\text{i}}$, $\tilde{U}_{\perp}=m_{\text{i}}\left(U_{\perp}-\frac{1}{2}\Omega^{2}(x-\bar{x})^{2}\right)/T_{\text{i}}$ and $\tilde{U}=m_{\text{i}}\left(U-U_{\perp}\right)/T_{\text{i}}$ gives

[TABLE]

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

[TABLE]

where we introduced the error function,

[TABLE]

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

[TABLE]

Re-expressing the integral over $\tilde{v}_{y}$ in terms of the standard Gamma function $\Gamma$,

[TABLE]

we simplify (125) to

[TABLE]

The density of open orbits is given by

[TABLE]

Note that, in equation (4.2), we have used

[TABLE]

for the adiabatic invariant of ions with $U_{\perp}=\chi_{\text{M}}(\bar{x})=\chi(0,\bar{x})=\Omega^{2}\bar{x}^{2}/2$. Using equation (41), we obtain

[TABLE]

Then, using the dimensionless integration variables $\tilde{v}_{z}=\sqrt{m_{\text{i}}\left(U-\Omega^{2}\bar{x}^{2}/2\right)/T_{\text{i}}}$ and $\tilde{\bar{x}}=\bar{x}/\rho_{\text{i}}$, equation (4.2) reduces to

[TABLE]

Equation (4.2) does not simplify further for general values of $x$, but can be simplified for $x\ll\alpha\rho_{\text{i}}$. Evaluating the ion density at $x=0$ using equation (4.2), we obtain

[TABLE]

We then proceed to evaluate $n_{\text{i,op}}(x)-n_{\text{i,op}}(0)$ for $x\ll\alpha\rho_{\text{i}}$. In this ordering, $\left(x/\rho_{\text{i}}\right)\left(2\tilde{\bar{x}}-x/\rho_{\text{i}}\right)\ll 4\alpha\pi\tilde{\bar{x}}\tilde{v}_{z}$ and so

[TABLE]

Note that there is a small integration region, $\tilde{v}_{z}\lesssim x/\alpha\rho_{\text{i}}$, where $\left(x/\rho_{\text{i}}\right)\left(2\tilde{\bar{x}}-x/\rho_{\text{i}}\right)\gtrsim 4\alpha\pi\tilde{\bar{x}}\tilde{v}_{z}$, but the contribution to the integral from this region is higher order in $x/\alpha\rho_{\text{i}}\ll 1$. In equation (133) we take $\sqrt{\left(x/\rho_{\text{i}}\right)\left(2\tilde{\bar{x}}+x/\rho_{\text{i}}\right)}\simeq\sqrt{2\tilde{\bar{x}}x/\rho_{\text{i}}}$, since this is accurate everywhere except where $\tilde{\bar{x}}\sim x/\rho_{\text{i}}\ll\alpha$, and evaluate the integrals over $\tilde{v}_{z}$ and $\tilde{\bar{x}}$ to obtain

[TABLE]

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

[TABLE]

for $x\ll\alpha\rho_{\text{i}}$.

The ion density profile for $\tau\rightarrow\infty$ is, according to (35), the sum of equations (123) and (4.2). The potential profile is obtained by imposing quasineutrality and inverting the Boltzmann relation for the electron density, to find

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

5 Numerical results

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

5.1 Boundary conditions

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

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

[TABLE]

The assumption behind equation (140) is that, just as the magnetic presheath solution self-consistentely satisfies the kinetic Bohm condition with the equality sign, the collisional presheath will self-consistentely satisfy the kinetic Chodura condition with the equality sign. In order to be consistent with the models in section 4 in the limits $\tau\rightarrow 0$ and $\tau\rightarrow\infty$, we also choose a set of distribution functions that:

• •

for $\tau\rightarrow 0$ is a Maxwellian that peaks at $v_{z}=v_{\text{B}}$;

• •

for $\tau\rightarrow\infty$ is a half Maxwellian that peaks at $v_{z}=0$.

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

[TABLE]

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

[TABLE]

Note that, from equations (31), (32) and (141), we can write the distribution function in the form $F(\mu,U)$,

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

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

[TABLE]

and

[TABLE]

In equation (148), we have introduced the exponential integral,

[TABLE]

Using equations (145)-(148), in figure 4 we plot the value of $u_{z\infty}$ as a function of $\tau$. Equations (144)-(148) are derived in Appendix F.

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

5.2 Numerical solutions

The numerical scheme presented in Geraldini et al. (2018) is used to obtain numerical solutions to the quasineutrality equation (45) for values of $\alpha$ and $\tau$ in the range $0.01\leqslant\alpha\leqslant 0.2$ (roughly corresponding to $0.57^{\circ}\leqslant\alpha\leqslant 11^{\circ}$) and $0.2\leqslant\tau\leqslant 10$. We define a quantity

[TABLE]

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

[TABLE]

where $E$ is a small number. An iteration scheme, outlined in Geraldini et al. (2018), is performed to find the numerical electrostatic potential solution $\phi_{\mu}=\phi(x_{\mu})$ for a given value of $\alpha$ and $\tau$. The solution numerically satisfies the quasineutrality equation with an error $E=0.7\%$ for all values of $\tau$ except for $\tau=0.2$, where $E=1.2\%$.

The electrostatic potential drop across the magnetic presheath is shown on the left in figure 5 as a function of $\alpha$ and $\tau$. The numerical results approaching $\tau=0.2$ and $\tau=10$ are consistent with the results obtained using equation (58) (valid for small $\tau$, $3/|\ln\alpha|<1/|\ln\tau|\ll 1$) and using equation (137) (valid for $\alpha\tau\gg 1$), shown with dashed lines. The shaded region is where we expect the assumption of an electron-repelling wall not to be suitable for Deuterium ions, $\alpha\lesssim\sqrt{1+\tau}\sqrt{m_{\text{e}}/m_{\text{i}}}\sim 0.02\sqrt{1+\tau}$. Considering the unshaded region in figure 5, the potential drop with finite ion temperature is up to $10-15\%$ smaller than the cold ion ($\tau=0$) potential drop. For a fixed angle, $\alpha=0.05\text{ rad}\approx 3^{\circ}$, the electrostatic potential profiles for different values of $\tau$ are shown on the right in Figure 5. The blue dashed curve labelled “0” in Figure 5 is obtained from equation (59), while the red dashed curve marked “$\infty$” is obtained from equation (136). The numerical profiles are consistent with the limits $\tau=0$ and $\tau=\infty$.

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

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

• •

for $\tau\gg 1$, the velocity components tangential to the wall, $v_{y}$ and $v_{z}$, remain unaffected while the velocity component normal to the wall, $v_{x}$, becomes of the order of whichever is largest between $\sqrt{\alpha}v_{\text{t,i}}$ and $v_{\text{B}}\sim v_{\text{t,i}}/\sqrt{\tau}$;

• •

for $\tau\lesssim 1$, all velocity components are affected by the magnetic presheath electric field and become of order $v_{\text{B}}$ (ignoring factors of $|\ln\alpha|$).

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

6 Conclusion

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

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

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

This work was supported by the US Department of Energy through grant number DE-FG02-93ER-54197. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Appendix A Glossary of notation

Here, we provide a glossary of some of the notation used in this paper. For each symbol, we give a brief description and a reference to the equation where the symbol first appears.

[TABLE]

Appendix B Derivation of equation (41)

In Geraldini et al. (2017) the quantity $\Delta_{\text{M}}$ appearing in the open orbit density (3.2) was expressed as

[TABLE]

We proceed to show that equations (152) and (41) for $\Delta_{\text{M}}$ are equivalent.

Open orbits have $U_{\perp}=\chi_{\text{M}}(\bar{x})$ to lowest order. Hence, their orbit position $\bar{x}$ determines the perpendicular energy $U_{\perp}$. Every ion in an open orbit must have come from a closed orbit which had an adiabatic invariant equal to $\mu=\mu_{\text{gk}}(\bar{x},\chi_{\text{M}}(\bar{x}))$, where $\mu_{\text{gk}}$ is defined in equation (28). Taking the total derivative of $\mu$ with respect to $\bar{x}$ leads to

[TABLE]

Using equation (21), we obtain the partial derivatives $\partial V_{x}/\partial U_{\perp}=1/V_{x}$, $\partial V_{x}/\partial\bar{x}=\Omega^{2}\left(x-\bar{x}\right)/V_{x}$. Then, differentiating equation (28) under the integral sign (which is possible because the limits of integration are points where the integrand vanishes), we get

[TABLE]

and

[TABLE]

To obtain $d\chi_{\text{M}}/d\bar{x}$, we first write

[TABLE]

As was argued in Geraldini et al. (2018), one of the two terms in $d\chi_{\text{M}}/d\bar{x}$ is $\chi^{\prime}(x_{\text{M}},\bar{x})dx_{\text{M}}/d\bar{x}=0$, because $\chi^{\prime}(x_{\text{M}},\bar{x})=0$ if the maximum is a stationary point of $\chi$, and $dx_{\text{M}}/d\bar{x}=0$ if the maximum is the non-stationary point $x_{\text{M}}=0$. Hence, only one term is left when differentiating equation (156),

[TABLE]

Inserting (154), (155) and (157) into (153), we obtain

[TABLE]

Then, equation (41) follows from (158) and (152).

Appendix C Chodura’s fluid model

In this appendix, we first recap Chodura’s fluid model, valid for any angle $\alpha$, and derive the differential equation (4.1). We then proceed to expand the fluid model to lowest order in $\alpha$ using the ordering $\alpha\ll 1$. We thus derive equation (59), which coincides with the solution of the kinetic model in the ordering (60) to lowest order in $\alpha$ and $\tau$.

C.1 General oblique angles: derivation of equation (4.1)

In this appendix subsection, we consider general oblique angles, $\alpha\sim 1$ (in radians). For $\tau=T_{\text{i}}/T_{\text{e}}=0$, all ions have the same velocity, the ion fluid velocity $\mathbf{u}=(u_{x},u_{y},u_{z})$, and thus the ion equations of motion (12)-(14) reduce to

[TABLE]

[TABLE]

[TABLE]

Here, ′ indicates differentiation with respect to $x$. The fluid equations (159)-(161) follow from the particle equations of motion (12)-(14) by setting $\mathbf{v}=\mathbf{u}$ and using $u_{x}=\dot{x}$ to write $\dot{\mathbf{u}}=u_{x}\mathbf{u}^{\prime}$ (thus changing the time derivative of every velocity component to a spatial derivative).

Adding equations (159)-(161) multiplied by $u_{x}$, $u_{y}$ and $u_{z}$ respectively, dividing by $u_{x}$ and integrating leads to

[TABLE]

where we used $\phi(\infty)=0$ and the boundary condition (54). We proceed to obtain a differential equation for $\phi(x)$ from equation (162), following the derivation in Riemann (1994)555In Riemann (1994) (originally in Chodura (1982)) the corresponding differential equation for $u_{x}(x)$ was derived.. Differentiating (55) gives

[TABLE]

Inserting (163) in (159) and re-arranging gives

[TABLE]

Equations (55) and (164) are substituted in equation (161) to obtain

[TABLE]

Using the boundary conditions in (54), equation (165) integrates to

[TABLE]

Substituting equations (55), (164) and (166) into the energy equation (162) results in equation (4.1), which is solved by imposing a boundary condition at $x=0$, the Debye sheath entrance.

We proceed to discuss this boundary condition. First, we note that equation (4.1) has a singularity at $|u_{x}|/v_{\text{B}}=\sin\alpha\exp\left(e\phi/T_{\text{e}}\right)=1$ and that our boundary condition at $x\rightarrow\infty$ imposed $|u_{x}|/v_{\text{B}}=\sin\alpha\exp\left(e\phi/T_{\text{e}}\right)=\sin\alpha<1$. Since a crossing of the singularity in equation (4.1) would not be physical, it follows that the quantity $|u_{x}|/v_{\text{B}}=\sin\alpha\exp\left(e\phi/T_{\text{e}}\right)$ should stay below unity or reach unity at $x=0$, $|u_{x}(0)|/v_{\text{B}}\leq 1$. However, the Bohm condition for a stationary Debye sheath requires that $|u_{x}(0)|/v_{\text{B}}\geq 1$. Therefore, the only way to match the magnetic presheath with the Debye sheath is by using the boundary condition $u_{x}(0)/v_{\text{B}}=\sin\alpha\exp\left(e\phi(0)/T_{\text{e}}\right)=1$. The electrostatic potential profile in the magnetic presheath can then be obtained by numerically integrating equation (4.1) using $e\phi(0)/T_{\text{e}}=\ln\left(\sin\alpha\right)$ as a boundary condition.

C.2 Shallow angles: derivation of equation (59)

We proceed to expand equation (4.1) for $\alpha\ll 1$, with the aim of obtaining $e\phi(x)/T_{\text{e}}$ correct excluding terms that are small in $\alpha$. The electrostatic potential $\phi$ in equation (4.1) changes from $\phi(\infty)=0$ to $e\phi(0)/T_{\text{e}}\simeq\ln\left(\alpha\right)$ at $x=0$. Neglecting terms that are small in $\alpha$ over the entire range of values of $\phi$, equation (4.1) becomes

[TABLE]

By substituting the definition of $\psi$ in equation (112), equation (C.2) becomes

[TABLE]

Notice that equation (168) satisfies $\psi^{\prime}=0$ for $\psi=0$ and therefore also satisfies $\phi^{\prime}=0$ for $\phi=0$, a condition which is satisfied by the exact equation (4.1)666As Riemann (1994) showed, the derivative of the right hand side of equation (4.1) evaluated at $\phi=0$ is equal to zero, and the second derivative is equal to zero when the Chodura condition is marginally satisfied (which is the case we consider). Equation (168) has both of these properties, while equation (C.2) has neither of them. but is not exactly satisfied by equation (C.2). There are two terms of equal size that give rise to the error in (168). One is the error in (C.2), and the other arises from the non-equivalence of $\psi$ and $e\phi/T_{\text{e}}$, giving

[TABLE]

Hence, equation (59) gives the uniformly valid magnetic presheath electrostatic potential in Chodura’s fluid model to lowest order in $\alpha$.

In figure 8, we plot the electrostatic potential, $\phi(x)$, that results from solving (4.1) (exact) and (59) (approximate) for four different values of $\alpha$: the approximate solution is different from the exact solution for $\alpha=0.4$, is very close to the exact solution for $\alpha=0.2$ and almost overlaps with the exact solution for $\alpha=0.1$ and $\alpha=0.05$.

Appendix D Alternative derivation of drift velocity of closed ion orbits into the wall

The drift velocity $v_{\text{d}}=\dot{x}_{\text{m}}$ can be obtained using the relation

[TABLE]

From equation (25) we have

[TABLE]

which can be differentiated to obtain

[TABLE]

Therefore, the drift velocity is

[TABLE]

Inserting $\dot{x}=v_{x}$ and equation (16) into $\dot{\bar{x}}=\dot{x}+\dot{v}_{y}/\Omega$, we obtain the relation $\dot{\bar{x}}=-\alpha v_{z}$. As a final step, we insert $\dot{\bar{x}}=-\alpha v_{z}$ into equation (173), and we use equation (74) for $v_{z}$ to recover equation (81).

Appendix E Alternative derivation of closed and open orbit ion density for small $\tau$

For $\tau\rightarrow 0$, the ion distribution function $F$ tends to

[TABLE]

We proceed to use this distribution function to derive equations (83) and (96).

E.1 Closed orbit density

Using equation (66), the adiabatic invariant of an ion in a closed orbit is given by

[TABLE]

with $x_{\text{b}}=x_{\text{m}}-\sqrt{2\left(U_{\perp}-\chi_{\text{m}}(\bar{x})\right)/\chi^{\prime\prime}(x_{\text{m}})}$ and $x_{\text{t}}=x_{\text{m}}+\sqrt{2\left(U_{\perp}-\chi_{\text{m}}(\bar{x})\right)/\chi^{\prime\prime}(x_{\text{m}})}$. In equation (175), the $O\left(\rho_{x}^{2}/l^{2}\right)$ error comes from the fourth order term of the Taylor expansion of $\chi$ around $x_{\text{m}}$, since the third order term integrates to zero. Equation (175) thus reduces to

[TABLE]

Inserting the distribution function of equation (174) into the closed orbit integral (33) and changing from $U_{\perp}$ to $\mu$ using equation (176) gives

[TABLE]

In equation (E.1), the upper limit of integration in $\mu$ was extended to $\infty$ because $\delta(\mu)$ is zero for $\mu\neq 0$ (in practice, $F(\mu,U)$ is small for orbits with $\mu\gg\tau v_{\text{B}}^{2}/\Omega$).

To calculate the integral in equation (E.1), we change variable from $\bar{x}$ to $x_{\text{m}}$ and change the order of integration so that the integral over $x_{\text{m}}$ is carried out first. By using the relation (172) for $d\bar{x}/dx_{\text{m}}$, and taking $\chi^{\prime\prime}(x_{\text{m}})=\chi^{\prime\prime}(x)\left(1-\rho_{x}/l+O\left(\rho_{x}^{2}/l^{2}\right)\right)$, equation (E.1) becomes

[TABLE]

Note that, when Taylor expanding the integrand, the terms linear in $\rho_{x}=x-x_{\text{m}}$ coming from the correction to $\chi^{\prime\prime}(x_{\text{m}})\simeq\chi^{\prime\prime}(x)$ integrate to zero. Hence, the size of the relative error has remained $O(\rho_{x}^{2}/l^{2})$. The rightmost integral evaluates to $2\pi$, and thus equation (E.1) becomes

[TABLE]

The straightforward integrals over Dirac delta functions give the density of closed orbits in (83).

E.2 Open orbit density

Expanding the integrand in equation (3.2) gives

[TABLE]

By changing variable from $\bar{x}$ to $\mu$, substituting (180) and inserting $\chi_{\text{M}}(\bar{x})-\chi(x,\bar{x})=\chi_{\text{c}}-\chi(x,\bar{x}_{\text{c}})+O(\tau^{2}\epsilon v_{\text{B}}^{2})$ (recall the discussion preceding equation (94)), where $\chi_{\text{c}}=\chi(x_{\text{c}},\bar{x}_{\text{c}})$, the integral (3.2) simplifies to

[TABLE]

Inserting the relation (41) into (E.2) gives

[TABLE]

Using (174) for the distribution function, the density of open orbits becomes (96).

Appendix F Integrals of distribution functions (141)

We proceed to derive equations (144)-(148) for the values of $\mathcal{N}$, $u$, $r$ and $u_{z\infty}$ associated with the distribution functions in (141). Integrating (141) over $v_{y}$ and $v_{z}$, we obtain the functions

[TABLE]

All the integrals in this appendix are carried out using the dimensionless variables $\tilde{w}_{z}=v_{z}/v_{\text{t,i}}-u$ and $\tilde{v}_{z}=v_{z}/v_{\text{t,i}}$.

Using (183), the normalization condition (144) is

[TABLE]

Evaluating equation (184) for $\tau\leqslant 1$, and changing integration variable to $\tilde{w}_{z}$ gives

[TABLE]

Thus,

[TABLE]

The integral in equation (186) evaluates to

[TABLE]

Hence, equation (144) for $\tau\leqslant 1$ follows.

Evaluating equation (184) for $\tau>1$, one finds, after changing the integration variable to $\tilde{v}_{z}$,

[TABLE]

The last integral in equation (188) is calculated in the following way. First, one can obtain the integral of the function $\exp(-\tilde{v}_{z}^{2})/(1+r\tilde{v}_{z}^{2})$ (which will be useful when imposing the kinetic Chodura condition (140) in the next paragraph). Re-expressing $1/(1+r\tilde{v}_{z}^{2})=\int_{0}^{\infty}\exp\left(-\eta\left(1+r\tilde{v}_{z}^{2}\right)\right)d\eta$, one has

[TABLE]

Changing the integration variable to $\xi=\sqrt{\eta+1/r}$ gives

[TABLE]

Then, using the relation

[TABLE]

the integral

[TABLE]

is obtained. Inserting this integral into (188), we obtain the expression for $\mathcal{N}$ in (144).

Equation (140) is used to obtain the values of the positive constants $u$ and $r$. For $\tau\leqslant 1$, one inserts the distribution function (183) into (140) and changes variable to $\tilde{w}_{z}=v_{z}/v_{\text{t,i}}-u$ to obtain

[TABLE]

Rearranging equation (193) and inserting the value of $\mathcal{N}$ gives equation (145). For $\tau>1$, one changes variable to $\tilde{v}_{z}=v_{z}/v_{\text{t,i}}$ in the integral (140) to obtain

[TABLE]

Inserting the value of $\mathcal{N}$ and the integral in equation (F) gives equation (146).

The ion fluid velocity is evaluated using

[TABLE]

For $\tau\leqslant 1$ one has

[TABLE]

The integrals in (196) evaluate to

[TABLE]

giving (147). For $\tau>1$, one has

[TABLE]

The integral in equation (198) is calculated, as before, by expressing $1/(1+r\tilde{v}_{z}^{2})$ as a definite integral,

[TABLE]

Then, integrating by parts and changing the integration variable to $\xi=\eta+1/r$ gives

[TABLE]

Using the definition of the exponential integral in equation (149), we obtain

[TABLE]

leading to equation (148).