{'fraction_non_alphanumeric': 0.03316293092148381, 'fraction_numerical': 0.009973903406856612, 'mean_word_length': 4.269370840135627, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sheath collapse at critical shallow angle due to kinetic effects

The Debye sheath is known to vanish completely in magnetised plasmas for a sufficiently small electron gyroradius and small angle between the magnetic field and the wall. This angle depends on the current onto the wall. When the Debye sheath vanishes, there is still a potential drop between the wall and the plasma across the magnetic presheath. The magnetic field angle corresponding to the predicted sheath collapse is shown to be much smaller than previous estimates, scaling with the electron-ion mass ratio and not with the square root of the mass ratio. This is shown to be a consequence of the kinetic electron and finite ion orbit width effects, which are not captured by fluid models. The wall potential with respect to the bulk plasma at which the Debye sheath vanishes is calculated. Above this wall potential, it is possible that the Debye sheath will invert.


Introduction and orderings
At the boundary between plasmas and absorbing surfaces, non-neutral layers, called sheaths, are known to form [1][2][3]. The nature of these sheaths is controlled by the difference between ion and electron mobility and the physics of their interaction with the wall. Plasma wall interaction is of great importance in many devices such as spacecrafts [4], Hall thrusters [5], plasma probes [6,7] and magnetic filters [8]. It is of particular importance for magnetic confinement fusion, where fluxes of electrons and ions from the plasma are directed by magnetic fields onto the divertor or limiter targets, which can lead to overheating and damage [9]. In order to reduce damage to solid targets, a common method is to incline the magnetic field at a small angle, α, relative to the wall, thus spreading heat flux over a larger area [10].
In this paper we consider the effect of magnetisation on the plasma-wall interaction. We focus on the magnetic presheath, a quasineutral layer with a size on the order of the ion gyroradius that appears on top of the typically much thinner Debye sheath in magnetised plasmas [11]. Within the magnetic presheath, we assume that the magnetic field is at a small angle α to the plane parallel to the wall of the device, B = − sin αê x + cos αê z .
Here we are using a set of orthonormal coordinates, {ê x ,ê y ,ê z }, such thatê x points in the direction perpendicular to the wall andê z points in the direction of the component of B parallel to the wall. The electric field is electrostatic and points solely in the x direction, such that where φ is the electrostatic potential. There is assumed to be no dependence on the y or z coordinates and no time dependence. The first treatment of the magnetic presheath was due to Chodura [11,12]. This initial study used fluid equations to describe the electrons and ions, complemented by numerical simulations with a "particle-in-cell" kinetic model. Chodura found a solution to the potential throughout the presheath and a condition on the entrance ion distribution function, known as the Chodura condition, for the existence of this solution. The Chodura condition specifies that for cold ions, the fluid speed of ions entering the presheath along the magnetic field lines must be at least equal to the Bohm speed where Z is the ion charge number, T e is the electron temperature and m i the ion mass (we call this speed the Bohm speed to distinguish it from the sound speed c s = (ZT e + γT i )/m i ; where the numerical factor γ depends on the particular ion physics). Between the magnetic presheath and the wall lies the Debye sheath.
The Debye sheath has a typical thickness of several Debye lengths and, unlike the magnetic presheath, is not expected to be quasineutral. Analogous to (and predating) the Chodura condition, the Debye sheath has a 'Bohm condition' [13], which states that ions entering the Debye sheath must be moving at least at the Bohm speed towards the wall (in contrast to along the field line). Subsequent fluid treatments of the magnetic presheath have built upon the work of Chodura. Fluid treatments are, however, only well motivated if the ion temperature T i is much less than the electron temperature T e [14][15][16]. For ion temperatures comparable to the electron temperature, the strong distortion of the ion orbits means that an accurate description of the the presheath requires a kinetic treatment. The difficulty of kinetic theory and the complication of an absorbing boundary has meant that most approaches have been numerical, with further particle-in-cell simulations and Eulerian-Vlasov codes [17][18][19][20][21], although some analytical progress has been made [22][23][24]. This paper will build on the recent work of Geraldini et al [14][15][16], whose kinetic model of the ions takes into account not only the distortion of the ion orbits, but also the contribution of ions that have recently broken free from their orbits due to the strong electric field near the Debye sheath. The electron density assumed by Geraldini et al, and a swathe of the literature, is Boltzmann distributed. This model is, in general, only valid in the limit of a very large potential drop in the Debye sheath, when most electrons are reflected. A common model to give the electrons a net flux is to truncate the entrance electron distribution in the velocity component along the magnetic field line [25][26][27]. The truncation velocity is determined by the wall potential because electrons considered energetic enough to reach the wall are assumed to be absorbed instead of reflected. In this paper, we will provide a drift kinetic derivation for such a model. We will treat the wall potential as a parameter, denoted φ cut , referred to as the 'cutoff potential' or simply the 'cutoff'. It is assumed that the potential and its first two derivatives are monotonic converging to zero as x → ∞ throughout this paper. Therefore, the imposed cutoff potential also sets a lower bound for the potential throughout the magnetic presheath. However, this lower bound is not usually reached within the presheath. Indeed, when the potential in the presheath does not reach the cutoff potential, this implies the remainder of the potential drop occurs in the Debye sheath (see figure 1).
The width of the presheath will be set by the distance at which ion orbits begin to contact the wall, which is of the order of the ion gyroradius, ρ i . The potential drop must be sufficient to repel a considerable fraction of the electrons. Hence, we arrive at the following ordering for the potential and its first derivative, and where φ denotes the derivative of φ with respect to x. One could imagine having a dial on the wall of the device that controlled the cutoff potential. Since the typical velocity of electrons is ordered as cutoffs corresponding to energies much greater than the electron thermal energy would result in only very fast electrons reaching the wall, of which there will be exponentially few. As such, it is cutoffs at ordering (4) that we expect to significantly alter the electron distribution function. The typical velocity of ions is ordered as Throughout this article we will use an 'i' subscript to refer to quantities associated with the ions and an 'e' subscript for the electrons. When we wish to refer to both electrons and ions we will use the subscript 's' for species. Due to the Bohm and Chodura conditions, the fluid velocity of the ions in the x direction goes from αv B to v B : the potential drop accelerates the ions towards the wall. By conservation of particle flux, the density must therefore decrease. This allowed Chodura to arrive at his prediction for the ion density at the Debye sheath entrance of [11] n i (0) ≈ αn i,∞ , where n i,∞ is the ion density at x → ∞. Since the system is quasineutral and the electron density and the potential are approximately Boltzmann related (n e ≈ n e,∞ exp(eφ/T e )), the corresponding potential drop in the presheath is eφ(0) T e ≈ ln α.
Typical estimates of the total potential drop (φ cut in our notation) over the presheath and Debye sheath come from considering the case of an electrically floating wall accepting equal and opposite electron and ion currents. The typical scaling of such cutoff potentials was found to be [9,11,28] eφ cut T e = 1 2 ln m e m i g T i T e .
Here g is a function of the ion to electron temperature ratio T i /T e which is independent of α. The important feature of this scaling for the total potential drop is its dependence on mass ratio and its lack of α dependence. Stangeby [28] noted that, for sufficiently shallow magnetic field angle, the potential drop in the presheath predicted by Chodura could equal or exceed the total potential drop in the combined presheath and Debye sheath. This would imply that there is a magnetic field angle at which the Debye sheath would collapse. By direct comparison of equations (9) and (10), Stangeby predicted sheath collapse at an angle that scales as the square root of mass ratio, and had value 3.35 • for a deuterium plasma.
The assumptions that led to this ordering for sheath collapse are, however, incorrect as they fail to account for the alteration to Chodura's fluid model due to the loss of electrons onto the wall and the resulting breakdown of the assumption of Boltzmann electrons. We will prove that, in fact, there is no Debye sheath collapse for the fluid model. Stangeby's central idea that an increasing potential drop across the presheath can cause the Debye sheath to collapse is nonetheless correct, but in fact this collapse owes its origins to finite ion orbit widths and is therefore an intrinsically kinetic effect. Furthermore, we will show that it occurs at angles which scale linearly with mass ratio, much smaller than those predicted by Stangeby.
The discrepancy can be understood by first noting that the approximation of Boltzmann electrons breaks down for potentials φ near the cutoff potential φ cut . At such potentials, the electron that just grazed the wall and was reflected has a small parallel velocity and a small change in the potential therefore results in a large change in this cutoff velocity of returning electrons. This makes the derivative of electron density with respect to potential diverge as one approaches the wall potential since This divergence occurs in the presheath precisely when the Debye sheath has collapsed and we will show that this implies the electric field does not require a singularity on the presheath scale to satisfy quasineutrality. As a result, equation (8) for the drop in density in the presheath is incorrect exactly in the case when the Debye sheath vanishes. Without the singularity, the electric force cannot overcome the magnetic force for ions with a large enough perpendicular velocity. Therefore, a fraction of ions reach the wall travelling almost tangentially to it, much like the circular orbit shown in figure 2. This gives rise to a density at the wall (or equivalently the Debye sheath entrance) that scales like α 1/2 , instead of like α [15,23]. To see this, we appeal to the geometrical picture of a circular orbit near the wall shown in figure 2. Because the ions are on circular orbits, the ion reaching the wall with the largest component of velocity perpendicular to the wall is the ion that just grazed it on its previous cycle. If such an ion is moving at typical thermal velocities along the magnetic field line, in the time it takes for the ion to make a single gyration, it will have drifted closer to the wall by a distance of order αρ i . The ion velocity then makes an angle θ to the wall, shown in figure 2, which we argue scales as √ α.
This means the velocity onto the wall scales as √ αv t,i . Since the ion flux at infinity scales like αv B n i,∞ , by particle flux conservation, the density at the wall must scale like √ αn i,∞ . Therefore, the ordering (8) for the ion density at the Debye sheath entrance used to recover the ordering (11) for the angle at which the Debye sheath vanishes is incorrect precisely when the Debye sheath vanishes, if finite ion orbit width is included. We will prove this scaling algebraically in section 3.1, but this simplified picture serves to gain an estimate and also shows why this could not be predicted by a fluid model.
Equipped with this estimate for the ion density, we see that the new estimate for the potential drop in the presheath is Direct comparison between equations (14) and (10) now implies that the Debye sheath will vanish at an angle α that scales with mass ratio rather than square root mass ratio. This motivates a different, smaller, ordering for α, given by The rest of the paper is structured as follows. In section 2 we will discuss the single particle motion of electrons and ions within the presheath, showing both to be quasi-periodic. We will build on this with the electron model in section 2.2, detailing the simple model of the electron orbits within the presheath which we will use to find the electron density in section 2.3. We give a more rigorous derivation of this electron model from drift kinetics in Appendix A.1 which will prove the validity of our model for the small values of α required to see the collapse of the Debye sheath. For the ion density, we review the ion model developed in [14] [15] in section 3.1. We will use this model to derive the equation for the ion density in section 3.2 which we will show can give rise to the √ αn i,∞ we predict, provided the ion orbits are sufficiently undistorted. In section 4 we will expand the quasineutrality equation close to and far from the wall analytically. Far from the wall, this will give rise to the modified kinetic Chodura condition, while close to the wall we will prove that, for the majority Wall Figure 2: A sketch of a circular ion orbit drifting towards the wall. The ions approaching the wall with the largest velocity perpendicular to the wall are those that just grazed the wall on the previous gyration (dashed line). In a time ∼ 1/Ω i the ion orbit drifts a distance αρ i towards the wall, causing these ions to make an angle of order √ α to the wall. The angles in this diagram have been exaggerated for clarity.
of the solutions, the presheath will satisfy a marginal kinetic Bohm condition. The truncation of the electron distribution will lead to a subsonic Bohm condition, which agrees qualitatively with the predictions of [29]. In order to gain insight into these specific solutions, we will appeal to the limits of hot and cold ion temperatures in section 5. We will see that there exist cutoff potentials that do not admit solutions, and that this lack of solutions is tied to finite orbit effects. These effects manifest themselves most obviously in the hot ion limit. We will show they are not present in the limit of cold ions. In section 6 we will detail how our solutions are reached numerically and where this procedure differs from [15]. In section 7 we will present our numerical solutions to the presheath for a range of temperature ratios, angles α, and cutoff potentials. We will show that there are indeed ranges of cutoffs for which no solution exists. Furthermore, we will show that this lack of solutions implies that there are angles of the magnetic field for which the current from the plasma cannot be zero unless the electric field of the Debye sheath reverses. These angles are of order m e /m i .

Trajectories within the presheath and the electron model
In this section we will discuss the model for the electron orbits within the presheath. This will begin with a general discussion for both ions and electrons of the constants of the motion at α = 0, and a justification of treating the electron and ion motion as periodic when α = 0. In subsection 2.2 we will use this insight to discuss the relevant variables for the electron model. Finally in subsection 2.3 we will use this model to calculate the electron density as a function of potential. We will see that this may deviate strongly from the Boltzmann response near the Debye sheath entrance.

The α = 0 problem and quasi-periodicity
The equations of motion for ions and electrons in the presheath arė andv z = sgn(Z s )Ω s v y sin α.
Here Ω s = |Z s |eB/m s is the gyrofrequency of a particular species. The factors of sgn(Z s ) have appeared so we may refer to Ω s as a positive quantity for both electrons and ions. The trajectories defined by equations (16)(17)(18), as in [14,15], can be solved for α = 0, giving the constants of the motions U ⊥ , U andx defined bȳ and Herex is a constant of the motion that will later be useful for characterising orbits, and U ⊥ and U are the perpendicular and total energy, respectively. For α = 0,x and U ⊥ are no longer constants of the motion. In this case we take equations (19) and (20) to be the definitions of these quantities, with the understanding that they are no longer constants of the motion. Indeed, we may now use our equations of motion to determine the rates of change ofx and U ⊥ , andU Hence, whilex and U ⊥ are no longer constants of the motion for α = 0, they are still useful as they vary slowly in α.
We may use our definition (19) to rewrite U ⊥ as Here we have defined the effective potential χ(x,x) (which will be most relevant for the ion model) as When α = 0, U ⊥ andx are constants and thus the motion is truly periodic in x for certain values of U ⊥ andx. Equation (19) implies that the range of x sampled in one period will be ∼ ρ s = v t,s /Ω s and thus the period of the motion will be ∼ 1/Ω s by equation (24). When α = 0, the change in U ⊥ over a time 1/Ω s will be a factor of α smaller than the value of U ⊥ by equation (23). Hence the motion is no longer truly periodic, but we will refer to it as 'quasi-periodic' since the particles are moving in a potential well that is changing slowly compared to the time taken to cross the well. While we have thus far used U ⊥ to discuss the quasi-periodic motion of the electrons and ions, it is possible for this quasiperiodic motion to have another approximate conserved quantity. Such a variable is known as an adiabatic invariant. We will give an expression for the adiabatic invariant for the electrons, µ, in the next section and adiabatic invariant for the ions, µ gk , in section 3.2.

The electron model
As we have motivated in section 2.1, the motion of both the electron and ion orbits is split into two timescales: one rapid oscillatory timescale and another, much longer, timescale over which the parameters of the orbits can change. With the electrons the case is further simplified by the small orbits of size ∼ ρ e . These orbits, unlike the ions, sample only a small change in the electric field over one period and as such are only weakly deformed by the potential. This allows us to give an anayltical expression for another variable, the adiabatic invariant, which is approximately constant over the electron's motion. This is µ given by where b = B/|B| is the normalised magnetic field and v E is the E × B drift velocity given by Together with U , the total energy of the electron, we will use this to specify a given electrons velocity. Of course, we still need one variable that will capture the rapid variation of the electron motion. This is given by the gyrophase, ϕ, defined by Hereê n = cos αê x + sin αê z denotes the unit vector perpendicular to both the magnetic field and the y direction. Assuming that the electron distribution function at each point x is independent of the rapidly varying gyrophase, we can cast the electron motion in terms of only the constants of the motion U and µ, f e (x, v) ≈ f e (U µ). Therefore, in these variables the electron distribution will be constant along the trajectory of the electrons and to determine the electron density at a given position one only has to ask what values of U and µ could still be present in this position. This will be done in the next section, where we will derive the electron density as a function of position.
In fact, the adiabatic invariant for the electrons is only approximately constant, as is the independence of the electron distribution function on gyrophase. Therefore, while the simple discussion above suffices to fully describe the electron model, it is unclear how accurate it is. Since we wish to consider the small angle ordering (15) as well as (11), it is necessary to ensure that the electron model is still valid where the electron and ion density are small in α. This is done in Appendix A.1, Appendix A.2 and Appendix A.3 using a drift kinetic derivation for the evolution of the electron distribution function. It will turn out that, despite various complications, the errors are of sufficiently small magnitude to be ignored in both ordering (11) and ordering (15).

The electron density
The electron density at any position x is given by over suitable limits. The aim now is to change variables to those given in section 2.2. The Jacobian for the coordinate transformation to {U, µ, ϕ} is given by So we have where the sum is over electrons moving into (σ = 1) and out of (σ = −1) the presheath. Here, we have exploited the decomposition of f e in equation (A.19) into the part independent of gyrophase, f e ϕ , and a correction which we find to be small in Appendix A.1. When written in the variables U and µ, f e ϕ is independent ofx (the distinction betweenx and x will turn out to be unimportant for the electrons) and as such we may use the distribution function evaluated at infinity. All that remains is to determine for which values of U and µ a given electron will be present at position x. The effect of the potential drop within the presheath is to decelerate electrons in the direction parallel to the magnetic field. In Appendix A.2, we write the total energy U as Since the total energy and the adiabatic invariant are constant within the presheath, the parallel velocity must decrease to compensate for the potential drop (the E × B drift will not compensate for the potential drop as it is smaller than the electron thermal speed by a factor of m e /m i ). Thus, an electron that enters the presheath at a finite velocity may at some point be reflected. As such, the minimum value of the total energy for a particle at position x with adiabatic invariant µ and gyrophase ϕ corresponds to a particle that has zero parallel velocity. Thus, Note that the term v y v E in equation (33) makes U min a function the gyrophase ϕ. However, this correction is small (of size m e /m i ), and when integrated over ϕ will be of size m e /m i , so will turn out to be unimportant. The maximum possible parallel velocity for an electron entering the presheath is infinite since particles can enter the presheath at any speed. For electrons exiting the presheath, however, a different consideration must be made. We assume that an electron that reaches the wall of the device is absorbed onto the wall and does not return to the presheath. Then, if an electron is leaving the presheath (i.e. has negative parallel velocity), it must have been reflected before reaching the wall. There must be a maximum energy of these particles leaving the presheath that depends on the adiabatic invariant, which we denote by U max (µ). To determine the value of U max (µ), we must consider which particles are lost to the wall within the Debye sheath. Energy is conserved for particles entering the Debye sheath. We may parametrise U max by the function Φ(µ, ∆φ DS ) (where we have denoted the potential drop in the Debye sheath, φ(0) − φ cut , by ∆φ DS ), declaring where the choice to isolate the Ω e µ is arbitrary, but driven by the form of our expressions for U , which in Appendix A.2 we write as with errors small by mass ratio (note that φ is evaluated atx and not x). In general the physics of absorption onto the wall of a plasma is complicated by the fact that the length scale of the Debye sheath, which is a few Debye lengths λ D , may be of the order of the electron gyroradius, which would lead to a complex form of Φ(µ, ∆φ DS ), depending implicitly on the form of the electron and ion distribution functions at the Debye sheath entrance as well as on µ and ∆φ DS . Inserting the relevant limits on our variables, we arrive at the following expression for the electron density at position x (36) Here, the integral over U has been explicitly split into the incoming and outgoing particles, contained in the first and second terms, respectively. The factors of v y v E in equation (36) are due to the fact that we aim to evaluate the potential at a given position x rather than atx. However, these factors are small by m e /m i and the integral can be Taylor expanded in them. After expanding, the lowest order correction due to v y v E vanishes when integrating over ϕ leaving a higher order correction of size m e /m i , which we neglect. We therefore arrive at an expression for the electron density that depends only on the input distribution function, the local potential, and the cutoff potential Important features of this density can be seen easily by making the change of variables from U tov defined byv This changes equation (37) to  Here we see an example entrance electron distribution (the truncated Maxwellian) that has a discontinuity. Note that as the φ cut is taken to infinity the distribution tends to a Maxwellian and the density profile to a Boltzmann distribution, as equation (37) would imply.
In the limit ρ e /λ D 1, we expect that the electron orbit width can be neglected. Neglecting the electron orbit width is equivalent to assuming that the electrons are absorbed when they reach a value ofx corresponding to the wall, such that φ(x) = φ cut . Appealing to equations (34) and (35), this corresponds to the assumption of Φ(µ, ∆φ DS ) = 0. This justifies our convention on Φ(µ, ∆φ DS ) in equation (34) as it isolates the potential required to reach the wall in this limit. When the electron gyroradius is comparable to the Debye length, the form of Φ(µ, ∆φ DS ) is unknown. As such, for the remainder of this paper we will only consider the limit where ρ e /λ D 1, and Φ = 0. In this particular limit, it will prove very important to consider the derivative of the electron density with respect to the potential at the origin. This is given by It follows that the derivative of the electron density with respect to potential at x = 0 will diverge as φ(0) → φ cut . This can be seen in figure 3, where we plot n e as a function of φ with a truncated Maxwellian input electron distribution function. In section 4 we will prove this infinite derivative in n e to be an important distinction from the Maxwell-Boltzmann response.
3. The ion model and density

The ion model
In this section we will review the ion model presented by [14,15]. Like the electrons, the aim is to change to convenient variables to describe the ion orbits. However, unlike the electrons, these orbits will be strongly distorted by the electric field and will be classified into two groups: 'closed orbits' where the ion will not contact the wall within the next period of the motion and 'open orbits' where, at some point within the next period of the motion, the ion on this orbit will contact the wall. Both must be included to give a correct expression for the density of ions near the wall. (25) defines the effective potential in which an ion with a givenx and U ⊥ oscillates. We refer to an orbit as being closed if in the future it will have a bounce point, which is a point at which the x directed velocity of the particle changes sign. Since we have defined closed orbits in relation to bounce points, and bounce points are characterised by Figure 4: Curves ξ(x), defined in equation (43), for two possible potentials in the presheath: one with finite gradient at x = 0, the other with infinite gradient. For a given value ofx, the effective potential χ(x,x) will have a stationary point wherever ξ(x) =x. The effective potential will have negative gradient whereverx > ξ(x), and positive gradient whereverx < ξ(x). Thus, the effective potential gains a well forx >x c . Forx m,I >x >x c (shaded dark grey), the effective potential will then be a well with two extrema, creating a type 2 orbit. Forx >x m,I (shaded light grey), the potential has only one extremum; a minimum, creating a type 1 orbit. Note that, in the right hand plot, the gradient of the potential is infinite at x = 0. Hence, there are always two solutions to equation (42) meaning allx give type 2 orbits. One possible ξ(x) curve not shown is one where ξ(x) is an increasing function of x, implying x c = 0. Such a ξ(x) would give rise to only type 1 orbits.

Closed orbits. Equation
where the extrema of the effective potentials are. The locations of the extrema are given by the equation Thus,x Since we consider potentials with positive and monotonic decreasing, convex up, derivative (see figure 1), we expect equation (42) to have zero, one or two solutions.
This can be shown by defining the function ξ(x), Whenever ξ(x) =x, the effective potential will have a stationary point. Furthermore, the effective potential will have positive gradient when ξ(x) >x and negative gradient when ξ(x) <x. From figure 4, it is evident that, at lowx, there will be no solutions to equation (42) and the gradient of the effective potential will be positive for all x.
The firstx for which equation (42) has a solution isx c , defined bȳ Forx >x c it is possible that equation (42) will have one or two solutions. Having two solutions implies the effective potential will have two stationary points (one maximum and one minimum), creating a potential well in which an ion could oscillate, provided its energy did not exceed that of the effective potential maximum. This is not the only type of potential well in which ions could oscillate, as all that is required for a potential well is a local minimum. Appealing once more to figure 4, we see that equation (42) will only have one solution whenx exceedsx m,I , defined bȳ The parameter regionx >x m,I corresponds to the effective potential decreasing initially from x = 0 before reaching a minima and increasing. Thus, this also defines a potential well in which the ion can oscillate. We now introduce the nomenclature used in [14,15,23] to classify orbits. An orbit is referred to as being of 'type 1' if it has a single local minimum, and it is referred to as being of 'type 2' if it has a maximum and a minimum. Effective potentials corresponding to these orbits are shown in figure 5 and we see from the above discussion thatx >x m,I gives rise to type 1 orbits whilex m,I >x >x c gives rise to type 2 orbits. Note that the derivative of the potential at x = 0 need not be finite and in such a case all orbits will be of type 2.
An ion associated with any of these effective potentials will be in a closed orbit provided the perpendicular energy, U ⊥ , is not sufficient to allow the ion to reach Orbits with perpendicular energy U ⊥ 1 will be closed orbits in both cases while orbits at U ⊥ 2 will by open orbits. Note that the typical size of U ⊥ 2 has been exaggerated for this sketch as open orbits will only have perpendicular energies marginally higher than χ M (x).
the wall. In other words, the perpendicular energy must be less than the effective potential maximum, χ M (x), defined by Here we have also defined x M (x), the x coordinate of the effective potential maximum for a givenx. Note that for type 1 orbits, χ M (x) = χ(0,x) and x M (x) = 0. In section 2.1, we showed thatx slowly decreases in time. As such, for a given ion, χ M (x) will also decrease slowly, eventually falling below U ⊥ (which is also decreasing, but more slowly). At this point, the ion will be free to travel towards x = 0 and is no longer considered to be on a closed orbit. To account for this effect, it is necessary not only to consider closed orbits, but also 'open orbits' 3.1.2. Open orbits. While closed orbits dominate throughout most of the presheath, in the area close to x = 0, there are much fewer closed orbits. Clearly we expect more orbits to become unbound as the ions drift slowly towards the wall. This motivates the fact that the presheath cannot be understood using only closed orbits and instead we need to include open orbits. We define an open orbit to be an ion trajectory which at a past time had bounce points and has no bounce points at future times. These orbits are of type 1 or type 2 depending on the type of closed orbit they last 'escaped' from. In reference [15], it was shown that the perpendicular energies of these open orbits was strongly concentrated around the effective potential maximum. This implies that the v x distribution of the ions with a particular value ofx near the wall is strongly concentrated around V 0 , defined by Physically, this makes sense as the orbits in consideration are those that have only recently become unbound and therefore their perpendicular energy will be only a small amount above χ M (x) because U ⊥ does not increase appreciably over a single orbit, by equation (23). Moreover, Geraldini et al [15] found that the velocity spread of ions in the x direction in open orbits is where the function ∆ M (x, U ) is defined by and x t and x b are the top and bottom bounce points: the locations where the velocity in the x direction reverses in the quasiperiodic motion. Therefore Although ∆ M is, for most cases, not analytically tractable, its scaling will be sufficient for our purposes. We first examine the integral in equation (49). The denominator of the integral vanishes at both end points of the integral but it is otherwise finite. Dealing with the lower bounce point first, we see that the integrand is always finite since χ M (x) − χ(s, x) is at worst quadratic in s − x b as s → x b (this occurs for type 2 orbits, see figure 5) and so the denominator χ M (x) − χ(s, x) is at worst linear in s − x b , which is cancelled by the numerator. At the upper bounce point, χ M (x) − χ(s, x) vanishes linearly as s → x t for both type 1 and 2 orbits. Therefore, as the integrand contains an integrable singularity (which is present in regular orbits), the integral in equation (49) will scale as v t,i /Ω 2 i regardless of the type of orbit, and we arrive at the overall scaling of ∆ M = O αv 2 t,i . As a result of the size of ∆ M , when the ion is near the effective potential maximum, such that (48)). However, away from the effective potential maximum, we

The ion density
As with the electrons, to compute the ion density, we need the integral over the regions that will not fictitiously include ions that ought to have encountered the wall in their prior trajectory. We have found these regions in the variablesx, U ⊥ and U in section 3.1 and so just need to change coordinates. The Jacobian of the (19), (20) and (21) and is The ion distribution function is most convenient when written in the conserved variables of the ions which are the energy U and the adiabatic invariant for the ions, µ gk , given in [14] as Here x b and x t are the top and bottom bounce points respectively, as defined in section 3.1. In line with [15], we may now denote the closed orbit ion distribution function as a function of the conserved variables {µ gk , U } to be F cl (µ gk , U ). Hence the density integral (50) is transformed to The sum is over the particles with positive and negative z directed velocities. A similar sum exists over the positive and negative x directed velocities but, since the distribution function is gyrophase independent at lowest order, this simply becomes a factor of 2.
All that remains is to find the appropriate limits on the coordinates to ensure that only true closed orbits are being integrated over. Consider an ion at position x. For this ion to be in a closed orbit, two conditions must be met: first, there must be a point in the ions effective potential at s < x such that χ(s,x) > χ(x,x). This implies 1 2 which rearranges to givex The minimum possible value ofx for which this condition could be satisfied isx m (x), defined byx giving the conditionx >x m (x). Having chosen a possiblex, the second condition for a closed orbit is that the particle must not have sufficient perpendicular energy to pass the maximum of the effective potential. This amounts to an upper bound on U ⊥ of χ M (x). The only condition on the parallel energy is that the velocity in the z direction must be real valued. The resulting ion density integral is given by Here we have also removed the sum over parallel directions due to that fact that no ions are reflected back from the sheath into the presheath. As well as closed orbits, we must also consider open orbits. Just as with closed orbits, we make the same change of variables from {v x , v y , v z } to {x, U ⊥ , U }. However, we must impose different limits on the allowed values of these coordinates at a given x.
To do so, we note that an ion on an open orbit will hit the wall on a timescale of order 1/Ω i . Thus we may exploit the smallness of equations (22) and (23) to argue that the value ofx and U ⊥ that this ion had a time of order 1/Ω i before will not change appreciably. An orbit becomes open when U ⊥ grows minutely larger than the value of the effective potential maximum χ M (x) that had stopped the ion from reaching the wall. This small increase will however leave the upper bounce point  For these simulations, α = 0.05 radians, and we use the boundary conditions specified in section 6 with T e = T i could not come from thisx: the perpendicular energy required to reach x from thatx would imply the orbit had been open for much longer than one period of the motion and thus no ions could exist on such an orbit. If x < x c , however, all open orbits are free to move towards x = 0 (by definition) and hence an ion with x < x c could have come from any valid closed orbit:x >x c . We denote the minimumx for an open orbit at position x to bex m,o given bȳ Once we have established the value ofx that an ion in an open orbit previously had, we can determine the range of possible v x . As argued above, it is only a minute change in U ⊥ that causes an orbit to become open so v x will be tightly concentrated around a single value. This is made precise by (48). Thus, While the presence of ∆v x ensures that the open orbit density is everywhere small, it is still the dominant contribution to the ion density near x = 0, where the closed orbit ion density tends to zero. We can see from figure 6 that this is indeed true and that in all cases it would be a mistake to neglect the contribution of the open orbits to the ion density.
As discussed in the introduction, the scaling of n i,op (x) with α near x = 0 is of special importance due to its relation to the total potential drop in the presheath. By definition, type 1 orbits have their effective potential maximum at x = 0, while the locations of the maxima of type 2 effective potentials occupy a continuous range of x. Recalling the discussion at the end of section 3.1, we see that, away from x = 0, the vast majority of ions are not near their effective potential maxima, and so ∆v x = O (αv t,i ) which therefore implies n i,op (x) = O(αn i,∞ ). As we approach x = 0, however, if there is a significant number of type 1 orbits present, these orbits will be near their effective potential maxima, and as such will have ∆v . This is the algebraic counterpart to the argument based on circular orbits in figure 2 as orbits that are not significantly deformed by the potential are naturally of type 1. By conservation of particle flux, we know that the product of the ion density and fluid velocity is constant. Therefore, when only type 2 orbits are present, the fluid velocity towards the wall is much larger than when type 1 and type 2 orbits are present. In section 4.2 we will see that this is an early sign of how the Bohm condition will emerge from the theory, and how the presence of type 1 orbits can have a startling impact on the presheath.

Quasineutrality
The solutions to the collisionless magnetic presheath are expected to be quasineutral and thus satisfy the equation The problem of solving the magnetic presheath amounts to finding a potential throughout the presheath such that equation (60) is satisfied. In this section we will study the tractable asymptotic limits of this. These limits are at large and small x where the Chodura and Bohm conditions are expected to manifest, respectively. Compared to section 5 of [15] only the electron physics has been changed and therefore the expansions of the ion density will be identical. We will find in section 4.1 that the new electron physics creates only a small modification to the Chodura condition in the large x expansion. However, in section 4.2, we will take more care in deriving the modification to the Bohm condition, which, despite looking only superficially different, we will see can be qualitatively different depending on φ(0) and φ cut . Indeed, we will show in section 4.2.2 how the breakdown of the Bohm condition for φ(0) = φ cut (a collapsed Debye sheath) arises naturally as the selfconsistent culmination of the Bohm condition's weakening as φ(0) approaches φ cut .

Expansion at x → ∞
At large x, ion and electron orbits are far from the wall. This means that the value of x itself becomes less relevant and the more relevant quantity is the relation between the potential and the density of the plasma. In [15] this relation was found in the asymptotic limit of a small slowly varying field, defined by the two equations and Equipped with the small parameters φ and φ one can then expand the ion and electron density. The calculation can therefore be ported across with the only difference being the more exact electron density expression as the extension to non-Boltzmann electrons. The central result of the expansion is then the "kinetic chodura condition" which is a condition for there to be non-oscillatory solutions of the potential at x → ∞. This alters the result of [15] to include the effect of kinetic electrons at the presheath entrance. Despite this being an inequality, it is expected that in most cases the plasma beyond the magnetic presheath self-consistently marginally satisfies the kinetic Chodura condition just as unmagnetised plasma beyond a Debye sheath usually marginally satisfies a kinetic Bohm condition [1]. We will later impose this condition numerically in section 6.

Expansion around x → 0
While at infinity we were safe to assume that the potential had a Taylor expansion in x, care must be taken as we approach the Debye sheath entrance. Indeed, it was shown in [15] that an infinite gradient at the origin is always present in a magnetic presheath with Boltzmann electrons.
With this in mind, in this section we will first expand quasineutrality in the parameter δφ defined by Thus we will not make any assumptions about the nature of δφ save for the fact that it is finite (equivalent to the statement that the potential is finite at the origin). We conduct this expansion for two cases: first, in the usual case of φ(0) = φ cut , and secondly, in the special case of φ(0) = φ cut . (39), we see that the electron density will have a Taylor expansion in δφ, giving (66) n e (φ) − n e (φ(0)) = dn e dφ δφ + 1 2 where the derivatives of density with respect to potential are evaluated at the origin. The ion density at the origin has contributions from both open and closed orbits.
In [15], it was shown that the leading order term of n i,cl is due to closed type 1 orbits and gives a contribution of where δχ is defined as Thus, all that remains to consider is the contribution from open orbits. To this end, we rewrite equation (59) as where we have defined ∆v x 0 as If type 2 orbits are present, then x c > 0. Hence, sufficiently close to x = 0, we will have x < x c and the second integral in equation (69) will be zero by equation (58). If no type 2 orbits are present, then the second integral scales likex m (x), defined in equation (56), which in turn scales like x or δφ and is small compared to √ δφ. Hence, only the first term in equation (69) can give a contribution comparable to that of n i,cl . The size of the first term in equation (69) is set by ∆v x − ∆v x 0 which is given by Depending onx, ∆v x −∆v x 0 will be evaluated for a type 1 orbit or a type 2 orbit. For a type 1 orbit, χ M (x) = χ(0,x) so the last term in equation (71) is zero and the first and third terms differ by a term of order δχ, while the second term, of order δχ 1/2 , would be the largest. If instead we are on a type 2 orbit, since χ M (x) ≥ χ(0,x), we may choose x sufficiently small that the above terms cancel leaving a piece of size at most δχ (first term cancels with third, second cancels with fourth). As such, to lowest order, we need only consider the part of the integral for which we are on a type 1 orbit, corresponding to a lower limit onx ofx m,I . After this, the lowest order term in the quasineutrality equation is Given that the integrand is always positive, this equation is only satisfied ifx m,I → ∞, which by equation (45) implies the gradient of the potential is divergent at the origin.
This divergence offers a considerable simplification to the quasineutrality equation: with an infinite gradient at the origin, all effective potentials are of type 2 and so all ion trajectories will be open a finite distance before x = 0. By examining equation (56), we see thatx m (x) (the minimum allowedx for a closed orbit to be present) diverges at small x (as we now expect). From this, we see thatx must be increased greatly to allow for closed orbits closer to the origin. However, asx is increased, the value of U ⊥ necessary to reach the origin also increases. Thus, at sufficiently small x, only exponentially few trajectories will be present that are still in quasi-closed orbits. Since the majority of ions are therefore far from their effective potential maximum, we may expand equation (71) in the limit of Collecting what are now the lowest order terms in the quasineutrality equation gives Here the dn e /dφ is evaluated at the entrance to the Debye sheath. This could be solved in the form δφ = qx, but this solution would break our assumption that the potential has infinite gradient at the origin. Instead, the terms of order δφ must cancel. This cancellation gives the marginal kinetic Bohm condition, modified compared to that derived in [15] dn e This recovers [15]'s kinetic Bohm condition in the limit of φ cut → −∞. Furthermore, we may note that this condition can be recast in the following form by using the following identity Here, the top hat functionΠ is defined bŷ Equation (77) is the marginal (equality) form of equation (79) of [1]. The inequality form of equation (77) was originally derived (with adiabatic electrons) in reference [30] as a generalisation of the Bohm condition to kinetic ions. As such, what we call the 'kinetic Bohm condition' is often referred to as the 'generalised Bohm condition'. We note that the kinetic Chodura and Bohm conditions (equations (64) and (77)) are both similar in form, as the kinetic Chodura condition demands that ions enter the presheath on average above the Bohm speed while the kinetic Bohm condition says that the presheath must accelerate these ions in the x direction into the Debye sheath. In the context of the magnetic presheath however, the kinetic Chodura condition arises as a constraint on the boundary conditions (in this case the entrance ion distribution function) while the marginal kinetic Bohm condition is an analytical property of the solution to the magnetic presheath. We further note that, while the derivative of density with respect to potential plays a small role in the kinetic Chodura condition, it can have a large impact on the marginal kinetic Bohm condition since, as discussed after equation (40), the derivative of the electron density with respect to potential at x = 0 will diverge as φ(0) → φ cut . This means that as φ(0) → φ cut the marginal kinetic Bohm condition becomes weaker, in the sense that the divergence of dn e /dφ makes the ion distribution move into lower regions of v x so the ions are being accelerated less strongly in agreement with figure 7b.
Despite giving important information about the distribution function at x = 0, the marginal kinetic Bohm condition does not tell us about the form of the potential. To find the true dependence of δφ on x, we must expand to a further order, keeping only terms that scale like δφ 2 or x. Thus, and this is solved by with p defined as Together with the modified kinetic Bohm condition given in equation (76), we may exhibit some important features of the factor p. First, in Appendix B we find a sufficient condition for p to be finite. This condition is where the derivatives are evaluated at φ(0). We further prove that the entrance electron distribution function we choose in section 6 satisfies this condition, and hence that we do not need to further expand the quasineutrality condition to compute the asymptotic behaviour of the potential at x = 0. We can also show that p tends to zero as φ(0) → φ cut . To see this, we differentiate equation (40) with respect to φ giving Provided that the distribution function f e ϕ (Ω e µ − Ω e φ cut /B, µ) is non-zero for some range of µ, it is guaranteed that the first derivative of n e with respect to φ (given in equation (40)) will diverge to +∞ as φ(0) → φ cut while the second derivative will diverge to −∞. This implies that the denominator of p diverges like 1/(φ(0) − φ cut ) 3/2 . While the numerator in the bracket in p will also diverge, it will diverge with the same divergence as the term in the marginal kinetic Bohm condition (76), since the additional presence ofx does not change the form of the divergence. Thus the numerator diverges only as 1/(φ(0) − φ cut ) 1/2 . It follows that p will tend to zero as φ(0) tends to φ cut . Therefore, the divergence of the electric field vanishes as the magnitude of the potential drop across the Debye sheath decreases. This is in agreement with the numerical results shown in figure 7a which shows the weakening of the electric field as the magnitude of the potential drop in the Debye sheath is decreased. While the vanishing divergence of the electic field is consistent with the absence of the Debye sheath implied by φ(0) = φ cut , it indicates the electron density must be expanded with more care in this regime. This is done in the following section.
4.2.2. φ(0) = φ cut . From equation (39), we see that when φ(0) = φ cut , it will no longer be possible to Taylor expand the electron density in δφ and instead we must expand in √ δφ, Since the ion physics has remained unchanged, the only alteration to equation (72) is that now there is another term of size √ δφ on the left hand side, giving the square root order quasineutrality equation Having assumed the potential to be finite at the origin, this equation is solved by a finite gradient in potential at x = 0. Therefore, we are motivated to (for sufficiently small x) rewrite δφ as φ (0)x. This makes equation (86) an implicit equation forx m,I , related to the gradient at the origin through equation (45), Thus, for φ cut = φ(0), the electric field at the origin is finite and there exist type 1 closed orbits. The transition from a divergent electric field to a finite one can be seen clearly in figure 7a.

The limits of hot and cold ion temperature
In section 4.2 we considered the distinction between the cases φ(0) = φ cut and φ(0) > φ cut . Given that the cutoff potential is a parameter of the problem, it is natural to ask whether solutions exist for all φ cut and whether they have φ(0) = φ cut or φ(0) > φ cut . We can answer these questions by appealing to the hot ion theory developed in [23] and [16]. The hot ion theory models the limit T e T i in which the potential cannot substantially distort the ion orbits. In this section we will first show that, for hot ions, there exists a 'critical cutoff': a cutoff potential above which quasineutral solutions that adhere to our monotonicity constraints on the elecrostatic potential do not exist. At this critical cutoff, we find φ(0) = φ cut , while for all solutions below this cutoff, we will have φ(0) > φ cut .
The existence of critical cutoffs is not unique to hot ions and it is possible to show that they apply to all finite ion temperatures. For zero ion temperature, we prove that no critical cutoffs exist due to the lack of finite ion orbit widths. This has an important implication: fluid models of the magnetic presheath fail at very shallow magnetic field angles

The hot ion limit
As the ion temperature is increased, the distortion of ion orbits due to the potential becomes weaker. In the limit of a very large ion temperature, the effect of this is that the ion motion at lowest order becomes independent of the potential. This was investigated in [16] with the input ion distribution function where Θ is the Heaviside function, For an electric field negligible for the ions, it is purely the finite orbit sizes that result in the depletion of ions in the presheath. It was shown in [16] that the density at x = 0 is given by where the gamma function Γ(x) is defined by Γ(x) = ∞ 0 t x−1 e −t dt. Note that equation (90) is predictably only a function of α. This is a manifestation of the fact that the ions are unaltered by the potential in the presheath.
Since the ion density is fixed regardless of potential, and the presheath is quasineutral, the potential will take the required value to make the electron density equal to the ion density. The electron density at x = 0 is If the input electron distribution function is a decreasing function of energy, then the lower bound of equation (91) is reached when φ(0) = φ cut (see figure 3). We denote the value of n e (0) when φ(0) = φ cut by n e,min , − Ω e B φ cut , µ dv .
If n e,min (φ cut ) < n i,op (0), then there is a quasineutral solution with φ(0) > φ cut such that n e (0) = n i,op (0). If n e,min (φ cut ) = n i,op (0), this implies that the quasineutral solution has φ(0) = φ cut . For such critical solutions, the Debye sheath has collapsed as the wall potential is reached at x = 0 on the magnetic presheath scale. Noting that equation (92) increases to n e,∞ as φ cut approaches zero, there is a 'critical cutoff' φ crit for which n e,min (φ crit ) = n i,op (0). For φ cut > φ crit , the electron density is constrained, sufficiently close to the wall, to be greater than the ion density and thus no quasineutral solutions can exist. Heuristically, the reason for the existence of the critical cutoffs is as follows. The electron density within the presheath is controlled entirely by the potential. A greater potential drop in the presheath results in more electrons being repelled and the density of electrons at x = 0 lowering. Raising the cutoff potential φ cut limits the potential drop in the presheath and thus limits the minimum electron density as fewer electrons can be repelled. By raising the cutoff, the minimum electron density can be made arbitrarily close to the electron density far from the wall. The ion density, however, is not only reduced by the presence of an electric field but also by the finite width of the ion orbits. This finite width means that the density will reduce as we let x → 0 as there is a smaller chance an ion could be found close to the wall without having previously touched the absorbing wall. Thus, the cutoff can be raised to a point where the minimum electron density is larger than the fixed ion density at x = 0 and, as a result, there can be no quasineutral solution within our assumption of a monotonic increasing  potential profile. The critical cutoffs can be recovered numerically for finite T e /T i (see section 6 for the numerical method). Figures 8a and 8b show the critical cutoffs as a function of α for two electron to ion temperature ratios. Even in the case of equal electron and ion temperatures, the solutions lie close to the theoretical hot ion limit of equation (90). This is a manifestation of the reduced distortion of the ion orbits near the wall due to the weaker electric field present in the critical solutions as these solutions are not longer forced to satisfy the marginal kinetic Bohm condition. In section 5.2 we show that this effect is indeed only present when modelling finite orbit widths, as modifying Chodura's fluid model does not give rise to any critical cutoffs.

The cold ion limit
In the limit of T i /T e → 0, the ions can be considered mono-energetic upon entering the magnetic presheath. As such, the ions will follow a single trajectory. This means that all ions move at the same velocity u, which is a function of the position x. This allows us to rewrite equations (16), (17) and (18) for the ions as equations for the fluid velocities u x , u y and u z as a function of position and u x du z dx = Ω i u y sin α.
As well as this, we make use of the fluid continuity equation Here, we have invoked quasineutrality in order to relate the flow equation to our much simpler equations involving the electron density. We can then integrate equation (96) to give where we have further specified that the electron density is a function of only φ locally, as in equation (39). Using equation (93), we find an expression for the y directed velocity in terms of the potential and its first derivative, This then allows us to integrate equation (95), We note that equation (98) is a function of the first derivative of φ. Together with the energy conservation equation (21), this allows us to write a first order ODE inφ, A salient property of this ODE is that the potential gradient need not be finite. Indeed, motivated by the discussion of the gradient in the general case, we see that the infinite gradient should mark the transition from the quasineutral presheath into the Debye sheath. This occurs when which can be recast, using equation (97), as a Bohm-like condition, consistent with equation (77) on the ions entering the Debye sheath, In the case of Maxwellian electrons, equation (101) recovers the result of Chodura [11] for the potential drop across the magnetic presheath. Note that equation (101) has a solution for any α and φ cut . For smaller α the density at the origin will reduce itself so as to solve the equation. Thus, lower α results in a more negative potential at the Debye sheath entrance as in [11,15]. In our electron model, there is a minimum density at the origin since the potential at the origin cannot be reduced below the cutoff potential. Yet, equation (101) still has solutions for all α because the derivative of the density with respect to potential can be made arbitrarily large by making φ(0) arbitrarily close to φ cut , as was shown explicitly in equation (40). Thus, unlike with the finite temperature ions, there is no critical cutoff. Solutions exist for every φ cut and α, implying that the fluid ion model will be unable to predict the collapse of the Debye sheath. This is in contradiction to the results of [28] where purely Boltzmann electrons were assumed. Physically, while the large circular orbits of the hot ions meant that the reduction in density was inevitable, the cold ions are influenced entirely by the electric field. Hence, if the cutoff potential is raised towards zero, the cold ions are accelerated less and the density remains high. Thus, there is no critical cutoff for cold ions. This explains why the critical cutoff is much higher for the colder ions in figure 9, since the cutoff must be raised much further before the finite orbit widths cause the electron density to exceed the ion density.

Numerical method
In the following section, we will outline how the numerical method used here differs from that used in reference [15]. We will first detail the entrance distribution functions for electrons and ions. After this we will discuss the alterations to the iteration scheme necessary to handle solutions with cutoffs and to find the critical cutoffs.

Electron boundary conditions
In place of the Maxwellian electrons that were used by [14][15][16], our model accounts for the loss of electrons to the wall. Without a theory of the collisional presheath to provide an appropriate boundary condition on the electrons, the choice of entrance electron distribution functions is arbitrary, just as the original choice of a Maxwellian input was. We choose a distribution function that will tend to the Maxwellian in the limit of φ cut → −∞: a truncated Maxwellian where erf(x) is the error function defined by erf(x) = (2/ √ π) x 0 exp (−u 2 ) du. To write equation (103) we have furthermore made the simplest choice of the function Φ(µ, ∆φ DS ) as described in section 2.3, and set Φ(µ, ∆φ DS ) = 0. In our conserved variables {U, µ}, the distribution function is Using equation (39) we arrive at the density as a function of φ(x) (see figure 3), To allow for different possible electron distribution functions at the entrance of the presheath, instead of using equation (105), the code takes as input the electron distribution function integrated over µ on an evenly spaced velocity grid to determine the electron density using the procedure detailed in Appendix C. For the simulation results shown in this paper, we use a spacing in parallel velocity of 0.03v t,e although resolutions as low as as 0.1v t,e would not alter the results greatly. Since we have now included a cutoff potential, it is possible to calculate the electron current density along the magnetic field line. This is given by This expression can also clearly be determined numerically for a general entrance electron distribution function.

Ion boundary conditions
As in [15,16], without a scheme for determining the proper input distributions for the ions, there is an infinite number of possible distribution functions. We propose a family of possible distribution functions parametrised byτ , This parameter deviates from T i /ZT e by only a small amount, but we use this parameter as it is more convenient for the numerical handling of the kinetic Chodura condition. Equipped with this, we choose identical input distributions functions to those chosen in [16], In the conserved variables for the ions, this gives The parameter N ensures that the distribution functions are normalised such that The resulting value of N is We motivated in section 4.1 that the entrance distribution function must at least satisfy the modified kinetic Chodura condition (64), and we further choose to marginally satisfy it. As such, we select the parameters u and r so that the distribution functions satisfy the modified Chodura equation (64). The resulting equations determining u and r are thus forτ > 1. The fluid velocity in the z direction of the distribution functions is then a function of u or r, respectively, and is given by and where we have introduced the exponential integral

The iteration scheme
The iteration procedure from [15] remains largely unchanged with the alteration of the electron model. We continue to use the numerical method detailed in [15] to calculate the ion density as a function of position for a given entrance ion distribution function and potential in the presheath. We discretise the presheath into a grid (x n ) in x and will calculate the potential and electron and ion densities at each x n . Let φ m (x n ) denote the potential at position x n on the m th iteration. Likewise n e,m (x n ) and n i,m (x n ) will denote the electron and ion densities at position x n on the m th iteration. We will see that there are steps in the iteration in which using the desired cutoff (denoted φ cut,desired ) would produce an error. As such, the cutoff potential is also iterated and φ cut,m will denote the cutoff at the m th iteration. At each step in the iteration, we verify how close the system is to fulfilling quasineutrality condition (60), and we deem quasineutrality to be satisfied if  Here N + 1 is the number of points in x n , and was taken to be 0.007 for our numerical solution. The electron density is calculated from (39) using the integration scheme detailed in Appendix C. For a given cutoff, this gives the electron density as a function of potential which we will denoteñ e,φcut (φ), The electron density (105) is a monotonic increasing function of potential for the distribution functions we consider, as well as others. As a result, we can define an inverse function denoted byφ φcut (n), This is achieved numerically by tabulating the electron density as a function of potential at a given cutoff and then defining a piecewise linear inverse. This tabulation used a potential grid evenly spaced between the cutoff potential and zero. Similarly, we further defineφ cut (n min ) to be the inverse function of n e,min (φ cut ), defined in equation (92) (i.e. this returns the cutoff potential that would make n min the minimum possible electron density in the presheath). Unlike the tabulation of the electron density as a function of potential for a particular cutoff, which must be recalculated each time the cutoff is changed, this function only needs to be calculated once. As such, for our simulations, the typical spacing for theφ φcut grid was ∼ 0.004T e /e while the grid spacing forφ cut was 0.001T e /e. We may now present the iteration scheme. After specifying an initial potential profile φ 1 (x n ), we begin the iteration: (ii) set φ cut,m = min φ m (0), φ cut,desired ,φ cut (n i,m (0)) ; (iii) set n e,m (x n ) =ñ e,φcut,m (φ m (x n )); (iv) stop the iteration if quasineutrality satisfies (117); (v) set φ m+1 (x n ) =φ φcut,m (wn e,m (x n ) + (1 − w)n i,m (x n )); (vi) update the entrance ion distribution function to ensure that the modified kinetic Chodura condition (64) is satisfied; (vii) go to step (i).
Here w is a weight factor (which took values between 0.8 and 0) that quantifies how much the potential will change at each iteration step. After the first potential guess, each subsequent iteration is calculated by inverting the electron density relation on a weighted average of the electron and ion density. The procedure (i) used to calculate the ion density for a given potential φ m (x n ) remains unchanged from [15]. The remainder of the procedure is devoted to finding the electron density and then inverting the weighted average of the densities to achieve the potential for the next iteration. This procedure is analogous to the one used in [15], but is now made more complicated by the presence of a cutoff potential. This is because the cutoff potential places a lower limit on both the electron density and the potential itself. If φ m (x n ) < φ cut,desired , then it would not be possible to find n e,m (x n ) usingñ e,φ cut,desired since φ m (x n ) would not be in the domain ofñ e,φ cut,desired . Similarly, if the ion density is lower than the electron density, then the weighted density wn e,m (x n ) + (1 − w)n i,m (x n ) can be outside the domain ofφ φ cut,desired orφ φ(0) , which would make step (v) impossible. Thus, in step (ii), we choose the cutoff to be the minimum of φ m (0), φ cut,desired andφ cut (n i,m (0)). For the majority of solutions, the iteration will terminate with φ cut,m = φ cut,desired , indicating that a solution has been found at the required cutoff. However, it is possible for the iteration to terminate with φ cut,m = φ m (0) =φ cut (n i,m (0)). These correspond to the critical cutoff solutions discussed in section 4.2.2. It is also evident from the iteration scheme that setting the desired cutoff any higher will result in the same critical cutoff solution.

Figure 7a
shows the potential throughout the presheath for a range of cutoff potentials. By plotting this against √ x, we see the linear dependence near the origin that characterises a diverging electric field at the entrance to the Debye sheath. This is what we predicted. As the cutoff potential is raised, the gradient of the potential with respect to √ x decreases. This corresponds to the factor p in equation (82) decreasing, which we predicted in section 4.2. Finally, when φ cut = φ(0), we see that the potential becomes quadratic in √ x. Thus the divergence of the electric field vanishes along with the Debye sheath for this solution.
The shrinking of the region with diverging electric field is closely related to the weakening of the modified kinetic Bohm condition discussed in section 4.2.1. Figure 7b shows this directly by displaying the marginalised ion distribution function in v x , f 0 (v x ) = f i (x = 0, v x , v y , v z )dv y dv z , at the entrance to the Debye sheath (x = 0) for the same solutions displayed in figure 7a. For solutions with φ cut = φ(0), the distribution falls rapidly to zero as v x → 0 which is in agreement with equation (77). As the cutoff potential is raised towards the critical cutoff, we find the density of particles with low v x increases. This is a manifestation of the electric field being weaker and ions being accelerated less strongly into the Debye sheath. Finally, when we reach the critical solution (φ cut = φ(0)), there is a non-zero density of particles with v x = 0, indicating that there are particles at a turning point as they reach x = 0, as was expected from section 4.2.2. In the nomenclature of section 3.1, this is equivalent to the statement that type 2 and type 1 orbits are present for solutions with φ cut = φ(0), whereas only type 2 orbits are present for solutions with φ(0) > φ cut (as type 1 orbits have turning points at x = 0). Figure 6 shows the ion densities throughout the presheath for two solutions, one at the critical cutoff and one at φ cut = −∞. The dashed and dotted lines show the closed and open orbit densities, respectively. In both cases, the ion densities drop significantly as we approach x = 0. However, there are qualitative differences between the two solutions. In solutions with φ(0) > φ cut , there is a region near x = 0 where the closed orbit density becomes very close to zero and the open orbit density makes up almost the entire portion of the ion density. This is an identical result to the one found in [15] and confirms that the open orbit density is necessary to correctly understand a kinetic treatment of the ions in the presheath. The steep gradients in the potential are causing all the closed orbits to be open before reaching the wall. Again, this can be phrased in the language of type 1 and type 2 orbits discussed in section 3.1: the steep gradient of the electric field makes all orbits type 2. Sufficiently close to x = 0, the only closed orbits that can be found are those with high energy on very large orbits (so that the restoring force of the magnetic field can compete with the large electric field) of which there are exponentially few. In contrast, when φ(0) = φ cut , the closed orbit density is linear near the origin. Closed orbits are present arbitrarily close to the wall. A subtle but important feature of figure 6 which we have not discussed yet is the increase in open orbit density at the origin when going from the non-critical to critical regimes. Of particular importance is the change of scaling of this density with α. In [15,23], it was shown that the density would scale proportionally with α in the shallow-angle limit for non-critical solutions. It was further demonstrated that, if there were a significant number of type 1 open orbits, the density would scale like √ α. This feature was not realisable outside of the hot ion theory (studied in [23]) as the diverging electric field did not permit type 1 orbits to exist. At the critical cutoff, however, the electric field is finite and no such limitations on type 1 orbits are present. Figure 8b shows the plot on logarithmic scales of density at the origin against magnetic field angle for two separate electron to ion temperature ratios. The solution with T e /T i = 0.1 lies close to the prediction made by the hot ion theory given in equation (90), although the agreement worsens as α is reduced since the approximation of section 5.1 that the ions are independent of the potential worsens as the presheath potential becomes more negative. The fact that the critical cutoff becomes more negative as α is decreased can be seen in figure 8a which shows the critical cutoffs as a function of α for the same two temperature ratios.
The existence of critical cutoffs, or equivalently maximum ion densities, presents a further interesting feature relating to the wall current. Raising the cutoff raises The value of α below which the ion current onto the wall will exceed the electron current as a function of the electron to ion mass ratio. The dashed line represents the lower limit in α that the code can solve numerically. Below this value of α the data has been extrapolated using a power law fit.
the electron current onto the wall, as fewer electrons are repelled when the potential is less negative. However, since the cutoff cannot be raised above the critical value, for each α there is a maximum electron current onto the wall. Using equation (106), we find that the maximum electron current density is [21] J e,max = −e v t,e √ π n e,min (φ crit )α.
The ion current onto the wall is simply Zeαn i,∞ u z,∞ , with u z,∞ given in equations (114-115), leading to the scaling of the ion current density where A is some factor of order unity and c s = (ZT e + T i )/m i the characteristic ion speed. Noting that the maximum ion density at the origin falls with α, we see by quasineutrality that the maximum electron current is guaranteed to fall faster in α.
There is therefore an α below which the maximum electron current cannot equal the ion current. At this α, if the system is driven to have zero total current, φ(0) is forced to equal φ cut (to achieve maximum electron current) and the Debye sheath will collapse. Numerically, this value of α could be found by scanning through the critical cutoff solutions at different α and finding the solutions for which the total current is zero. Equivalently, and for ease of calculation, we find these α values by fitting the ion density at x = 0 at the critical cutoff (figure 8b) to a power law in α and solving J i = J e,max for α. This critical angle is shown in figure 10 as a function of the electron to ion mass ratio. The linear behaviour is exactly what would follow from a √ α dependence of the density, predicted in [15], when type 1 orbits are present. Note that for a Deuterium plasma with equal ion and electron temperature, the critical angle for which the Debye sheath collapses in ambipolar conditions is 0.3-0.4 degrees. This is an order of magnitude smaller than previous estimates of about 3 degrees [28]. We show in figure 11 that the Debye sheath does not vanish at this angle for realistic values of the temperature ratio within our model.

Conclusion
We have extended the model of Geraldini et al [15] to account for the effects of absorption of electrons with sufficient energy to reach the wall. Drift kinetics was used to find the kinetic equation satisfied by the electron distribution function. The density of electrons within our model was derived in equation (39), which we simplified by neglecting electron orbit size. The kinetic Bohm and Chodura conditions were derived in section 4, where we showed that the truncation of the electron distribution leads to subsonic sheaths, which agrees qualitatively with the results of [26]. It was shown in section 5.1 that the kinetic Bohm condition fully breaks down at a certain value of the wall potential, referred to as the critical cutoff. At this critical cutoff, the Debye sheath vanishes entirely, bringing the quasineutral presheath into direct contact with the wall. It was further shown in section 5.2 that a Chodura-like fluid theory cannot capture the existence of a critical cutoff due to the lack of finite ion orbit widths (equation (101)). The numerical method of Geraldini et al [15] was altered to include the revised electron model, and the existence of critical cutoffs for finite ion temperatures was shown. Finally, it was argued that, should the wall be at a floating potential, there is a critical magnetic field angle at which the cutoff potential is equal to the critical cutoff. sheath is present. Below this angle of incidence, or equivalently for a wall potential above the critical cutoff, the behaviour of the magnetic presheath and Debye sheath is unknown. The solution to the plasma wall boundary above the critical cutoff, therefore, must break one of the assumptions on the nature of our ion or electron model. The reversal of the sheath potential, φ cut > φ(0), which is excluded by our assumptions, is a possible scenario. In section 2.2, we argued that by using variables which are approximately constant for the electron motion, the electron distribution function became independent of x. As a result, finding the electron density became a one dimensional integral where one only had to take into account the deceleration of the electron along the magnetic field line. However, several approximations have been used in this argument. Firstly, the adiabatic invariant given in (26) is not a constant of the motion, but varies slowly. Secondly, the distribution function at a given position need not be independent of gyrophase as changes in the electric field can distort orbits. While these are only small corrections, the electron and ion densities associated with the Debye sheath entrance are also very small and could be of the same size as these corrections. In this appendix therefore, we will provide a drift kinetic derivation of the electron model, which will allow us to track the size of the errors in our approximation.

Appendix A.1. Drift kinetic electrons
We begin with the Vlasov equation for the electrons: The aim of this section is to rewrite equation (A.1) to achieve an expression for how the distribution function changes as we move into the presheath. Equation (22) implies that, despite their fast oscillations with timescale 1/Ω e , the electrons approach the wall on a timescale ρ i /v t,e α. Furthermore, the electrons oscillate on length scales of order ρ e , much smaller than the characteristic width of the presheath. We will therefore be able to expand the electron motion in small ρ e /ρ i ∼ m e /m i as well as small angle α. Given such conditions, the classic approach to this problem is that of drift kinetics [31]. Hence we will change variables to those quantities that vary slowly over a single particle oscillation. Of course, the gyrophase given in equation (28) will capture the fast motion and can be alternately written as Hereê n = cos αê x + sin αê z denotes the unit vector perpendicular to both the magnetic field and the y direction. The quantity v n = v ·ê n is the component of velocity in the direction ofê n while v E denotes the E × B drift velocity, We have chosen to use velocities suited to the basis {ê n ,ê y ,b} as it is in the plane perpendicular to the magnetic field that the electrons oscillate. Using this basis removes certain terms small in α that would complicate the algebra.
Clearly the other variables we will choose are the energy U and the adiabatic invariant given in equation (26) which we can rewrite as This variable was chosen over U ⊥ as it is designed to vary on average zero over a single period and is also very convenient for expressions involving ϕ.
Using equations (19), (21) where we have removed the t, y and z related terms as the distribution is independent of these (as explained in the introduction), and we have removed the termU (∂f e /∂U ) sinceU = 0. Here, the dot operator is defined bẏ which we will compute for our variablesx, µ and ϕ. First noting that equations (A.2)) and (A.4) give and we find the dot operator of our variables equal to Our choice of variables has made the characteristic frequency of the operator L acting on G much smaller than the characteristic frequency of particle gyration. Note that presently the choice of ordering (11) or (15) does not change the overall scaling of L (G), which is dominated by the largest terms, of order m e /m i . Now, for functions G = G(x, U, µ, ϕ), we further define the gyroaverage of G, G ϕ , by This is an integral over the gyrophase at constantx, µ and U . The gyroaverage allows us to decompose G into its gyrophase dependent and independent pieces, Decomposing f e into its gyrophase dependent and independent pieces, we get This equation tells us that the gyroaveraged distribution function, written as a function of the total energy U and adiabatic invariant µ, is independent ofx. However, the left hand side of equation (A.26) scales with α. Therefore, for small α, we must carefully justify that the neglected terms of equation (A.25) are smaller than the term we retain. The second neglected term on the right hand side of equation (A.25) includes a factor of α, hence, is trivially smaller by a factor of mass ratio than the left hand side for all α. A naive glance at the other neglected term, L(f e,1 ) ϕ , would indicate that it has size (m e /m i ) 3/2 . Since this does not scale with α, this would imply that, in our small α ordering in equation (15), the errors in our equation for f e ϕ would have size m e /m i . This is however, not the case, as we show in Appendix A.3 that the offending term which would appear to have size (m e /m i ) 3/2 in fact cancels, leaving terms that are smaller in mass ratio in both orderings (11) and (15).

Appendix A.2. Gyroaverages
Here we calculate the gyroaverages of the quantitiesẋ,μ andφ. To do this we must rewrite equations (A.9), (A.10) and, (A.11) in terms of the variables that we hold constant for the gyroaverage, and ϕ. Since all our variables contain v , it is first helpful to show that its ϕ dependence is small by a factor of m e /m i . We start by rewriting our equation for U Here we have Taylor expanded φ aroundx since we require our terms to be functions only of variables that we hold constant (this is a procedure that will be employed frequently in the gyroaveraging process). Now analysing the last three terms of equation (A.27), we have Here the error in equation (A.33) comes from the term μ ϕ ∂ f ϕ /∂µ.

Appendix B. Proof that p is finite
We will first show that p defined by equation (82) is finite for the boundary conditions defined in section 6. To do this we need only show that the denominator in equation (82) will always be non-zero. We rewrite the first term of the denominator as follows x dv x . (B.1) Hence it will be sufficient to prove that To do this we use the Bohm condition (equation (77)) and Schwarz's inequality to write We see that we manifestly satisfy our sufficient inequality (83). Thus, p is finite for all solutions.
The terms involving square roots can be solved by hyperbolic substitution. This gives Here, the first and second derivatives of the function are calculated using finite difference methods to the required level of accuracy, We will now address the size of the errors of our approximation to the integral. A common feature of strip based methods of integration is the scaling of the error with the strip size. Due to our choice of evenly spaced strips in x, however, the strips in z are not even. Therefore, we will split our discussion into one of regular strips in z with typical width ∼ h and irregular strips of width ∼ √ h. In both cases, the size of the neglected term is of order h 3 from two contributions: the truncation of the Taylor expansion of g(x), and the finite difference expansion of the derivatives g(x). There are an order unity number of irregular width strips so the total error from the irregular strips scales as h 7/2 . There are of order 1/h regular strips of width h and therefore the total error of these strips scales like h 3 .