Equilibrium probability distribution of a conductive sphere's floating charge in a collisionless, drifting Maxwellian plasma

A dust grain in a plasma has a fluctuating electric charge, and past work concludes that spherical grains in a stationary, collisionless plasma have an essentially Gaussian charge probability distribution. This paper extends that work to flowing plasmas and arbitrarily large spheres, deriving analytic charge probability distributions up to normalizing constants. We find that these distributions also have good Gaussian approximations, with analytic expressions for their mean and variance.


I. INTRODUCTION
A dust grain in a plasma acquires an electric charge by collecting electrons and ions that land on it. Because electrons and ions arrive at the grain at random, the grain's charge fluctuates, and because this fluctuating charge affects the dust's physical behaviour, the probability distribution of the charge is of practical interest.
Deriving this probability distribution for an arbitrary grain away from equilibrium in a plasma of arbitrary collisionality is very difficult, so we consider here a simpler case: a spherical, conductive grain at equilibrium in a collisionless plasma. This problem has been solved before [1,6,7,15] but these solutions have two major limitations. Firstly, they assume a stationary plasma. Secondly, most of these solutions apply the OML (orbital motion limited) grain charging theory, but OML is limited by its requirement that the grain be small relative to the Debye length [12].
In this paper we go beyond previous work by relaxing these assumptions. Instead of OML we use the more general charging model SOML (shifted OML), which allows for plasma flow by assuming a shifted Maxwellian ion velocity distribution far from the sphere [3,13,14]. We also circumvent the small sphere requirement that orthodox SOML inherits from OML [13, pp. 94-95], by applying SOML differently for arbitarily large spheres. Ultimately we derive two sets of results which, between them, cover spheres of all sizes (except those so small that their average charge is of order e). One set of results applies to "mid-sized" spheres, the other to "large" spheres. A "mid-sized" sphere is one with a λ D and a "large" sphere is one with a ≫ λ D , where a is the sphere's radius and λ D the Debye length. We use these definitions of "mid-sized" and "large" throughout the paper.
We start by building two stochastic models of sphere charging, one based on ordinary SOML and the other on the alternative large sphere model. We then solve both stochastic models for their equilibrium probability distributions, and derive Gaussian approximations to both. * Electronic address: dmt107@imperial.ac.uk † Electronic address: m.coppins@imperial.ac.uk

II. ELECTRON AND ION CURRENTS
Before building our stochastic charging models we need the electron and ion currents onto a sphere. In this paper we assume the electrons conform to the Boltzmann relation, and hence that a sphere of radius a and surface electric potential φ has an electron current I e = −4πa 2 n 0 e k B T e 2πm e exp eφ k B T e (1) where n 0 is the electron density far from the sphere and T e the electron temperature. The Boltzmann relation works well regardless of plasma flow because flow is invariably much slower than the electron thermal speed.
To assume the Boltzmann relation we require φ < 0. This is sometimes false for spheres so tiny that random fluctuations can render their charge positive, so we exclude those spheres from our calculations. This leaves us with mid-sized (a λ D ) and large (a ≫ λ D ) spheres.
SOML gives good estimates for ion currents onto the former. According to SOML [10,14], the ion current is when φ ≤ 0, where Z is the ions' charge state, T i is the ion temperature, m i is the ion mass, v is the flow velocity normalized by 2k B T i /m i , and are auxiliary functions of the normalized flow velocity. Both functions tend to 1 in the no flow limit (v → 0). SOML also gives I i for φ > 0 but the resulting stochastic model is intractably complicated -see appendix A. For large spheres the ion current is less than that obtained from eq. 2 because a larger a/λ D brings into existence absorption radii that undermine the SOML model [14]. To avoid this problem we follow Willis et al. in assuming that these absorption radii are "at or within the sheath", and therefore apply SOML at the sheath edge instead of the sphere's surface [14]. I i is then given by equation 2, but with the sheath edge potential φ s substituting for φ . To find φ s , we suppose that ions cross the sheath edge perpendicularly at the Bohm speed u B (irrespective of flow) and make the thin sheath (a ≫ λ D ) assumption that all ions crossing the sheath hit the sphere. This implies a second expression for I i , where the ion density at the sheath edge follows from assuming quasineutrality at the sheath edge. Equating this with eq. 2 and setting φ = φ s gives Estimating the Bohm speed as where γ is the heat capacity ratio [14], and substituting it into eq. 6 gives where Θ ≡ T i /T e is the ion-to-electron temperature ratio. Solving this equation for φ s gives where c ≡ Θs 1 (v) and W (x) is the principal branch of the Lambert W special function [2]. The expression for φ s is unwieldy but has the nice property of being independent of the sphere's charge. We assume that the plasma contains only one species of singly charged positive ion, so we can simply take Z = +1 and c = Θs 1 (v)/s 2 (v) from this point on.

III. BUILDING THE STOCHASTIC MODELS
Over a sufficiently short period of time δ , it is extremely unlikely that the sphere has time to collect multiple particles. Hence during a small enough δ effectively only three events may happen: the sphere absorbs nothing; the sphere absorbs an electron; or the sphere absorbs an ion. The chances of these events happening during a given δ depend (to a good approximation) on only the sphere's net charge Ne at the start of that period. As such we can model the sphere's charge fluctuations as a Markovian "one-step process" [11, p. 134], where the sphere's state is its net charge (in elementary charges) N, and N changes only in sporadic increments of ±1.
A one-step process is characterized by its rate coefficients r N , the probability per unit time of a shift from state N to state N − 1, and g N , the probability per unit time of a shift from state N to state N + 1 [11, p. 134]. In our model these rates correspond to the electron collection rateṄ e and ion collection rateṄ i .
The electron collection rate iṡ where is the collection rate for a neutral sphere of the same size.
The ion collection rate depends on which potential we use in equation 2. For a λ D we use the surface potential φ , but for a ≫ λ D we use the sheath edge potential φ s , as explained above. To accommodate both options we writė where µ ≡ m i /m e is a normalized ion mass, and φ ′ is either φ (if a λ D ) or φ s (if a ≫ λ D ). This will lead to two different models, one for mid-sized grains and one for large grains.
For both models the rate coefficients r N and g N are (14) and To complete the stochastic models, one must define φ in terms of N. Assuming the sphere is conducting, and neglecting polarization from nearby plasma particles, φ is is a dimensionless characteristic parameter analogous to the Coulomb coupling parameter Γ c for a simple plasma; in fact α = (d/a)Γ c , where d is the average inter-particle distance. The parameter α is key to the results that follow, with many of our approximations relying on its being small (α ≪ 1 or α ≪ Θ 1). Substituting eq. 16 into equation 14, For a λ D we likewise substitute eq. 16 for φ ′ in eq. 15: which we rewrite as to streamline the next section's algebra. For a ≫ λ D , φ = φ s and so eq. 15 is independent of N, being where we drop the N subscript to emphasize the independence from N for large spheres. We now have one complete stochastic model for mid-sized grains (comprising eqs. 18 & 20) and one for large grains (comprising eqs. 18 & 21). The next step is solving them for the charge probability distribution f N . Normally one would solve each model's master equation [11, passim], but these models' master equations aren't exactly solvable. For an exact solution we use a more direct approach.

IV. SOLVING THE STOCHASTIC MODELS
At equilibrium detailed balance holds [11, p. 142]. That is, in a huge ensemble of sphere-in-plasma systems at equilibrium, just as many should be going from charge state N to N − 1 as are going from charge state N − 1 to N (on average). As such which is a recurrence relation that has f N as its solution. We solve it for the a λ D case first. Substituting eqs. 18 and 20 into eq. 22, for N ≤ 0. To solve this equation, note that for N < 0. By inspection the first product is and the second is equivalent to where (x) + n ≡ is the rising factorial (or "Pochhammer symbol"), defined as Putting together eqs. 24, 25 & 26, which is, unfortunately, not analytically solvable. However, equation 28 permits numerical calculation of N's probability distribution in practice; one can simply compute ∑ f N / f 0 for those N where f N / f 0 is non-negligible, and set equation 28's f 0 to that sum's reciprocal. As one may rewrite (x) + n as a ratio of gamma functions (except when x or x + n is a negative integer), eq. 28 constitutes an analytic definition of f N in terms of elementary functions and the gamma function, lacking only f 0 's value: We now turn to the a ≫ λ D case. Substituting eqs. 18 and 21 into eq. 22, where g * ≡ g/(χδ ) is a more convenient form of g. The implied product of exponentials is readily solvable by inspection: The normalization condition does not appear to be analytically solvable for this distribution, either.

V. THE MODAL CHARGE AND THE STOCHASTIC MODEL'S VALIDITY
Deriving a closed form expression for f N 's mode is conceptually straightforward. Given the mode at M, We derive M for mid-sized grains first. Setting equation 23's left hand side to 1, substituting M for N, and solving, where W (x) is again the Lambert W special function's principal branch. For small α (i.e. large T e a), M's dependence on the Lambert W term is weak and M's dependence on α goes approximately as O(1/α). Equation 33 leads to an obvious precondition for the stochastic model's validity. As the model assumes the sphere never has a positive charge, a positive value of M indicates that the model has broken down and become self-inconsistent. Therefore M ≤ 0 is a necessary condition for model validity.
When is M ≤ 0? Rearranging eq. 33, M ≤ 0 when By exploiting W (z)'s definition, monotonicity and positivity for positive arguments, one can rewrite the inequality to remove the right hand side's dependence on α: This sets an upper bound on α, above which the inequality is unsatisfied and the model fails. This is not surprising, as the model assumes the sphere is not very tiny, which implies a small α ∝ 1/(T e a). Substituting in c (eq. 10), Evidently, in the cold ion limit (Θ → 0), the inequality reduces to α 0 and is never satisfied, indicating model failure. This is also unsurprising, as a vanishing Θ requires either an infinite T e (and hence an infinite electron current, from eq. 1) or a zero T i (and hence an infinite ion current whenever φ < 0). Inversely, for vanishing T e (Θ → +∞), the inequality's right hand side tends to −∞, in which case the inequality is again never fulfilled and the model fails.
The inequality also shows that flow affects the model's validity. Because Θ > 0 and s 1 (v) increases with v (figure 1), equation 36's right hand side becomes negative for large v, and the model eventually fails. Fortunately this only occurs at truly huge flow velocities, as shown by solving the trivial sub-inequality Because ions are rarely hotter than electrons, and the lightest ions are protons, µ/ √ Θ m p /m e = 42.85, which s 1 (v) is always less than for any reasonable flow speed (v < 48.34).
More sedate flow velocities have the effect of increasing eq. 36's right hand side, as shown by calculating that and substituting into equation 36: Inevitably µ √ Θ > 2Θ, so flow's effect (to second order) is to increase the right hand side, raising the chance of satisfying the inequality and loosening the validity constraint. The second order effect can also compensate for a low Θ; as Θ shrinks, µ √ Θ gets greater relative to the negative O(Θ) terms, enhancing flow's beneficial effect on the model's validity.
Deriving M for large grains is trivial. Setting equation 31's left hand side to 1 immediately leads to The direct 1/α dependence dovetails with eq. 33's approximate 1/α dependence for small α.
The same M ≤ 0 validity condition applies here. From eq. 41, M ≤ 0 if and only if g * ≤ 1, i.e. if We can get some insight from this knotty expression by considering limiting cases. For example, we can expand about Θ = 0 to obtain a validity inequality for a cold ion plasma: Neglecting higher order terms, this inequality is always true, since the LHS is always less than √ 2π, and √ 2π < µ. With a large grain and cold ions the stochastic model is always selfconsistent, regardless of flow. Another limiting case is that of vanishing flow velocity. With v = 0, s 1 (v) = s 2 (v) = 1, and eq. 42 becomes By exploiting W (z)'s definition, monotonicity and positivity for positive arguments once again, For ease of solution, we replace this with a more stringent validity condition. Specifically, when the inequality is clearly an even tighter bound on Θ than eq. 45. This tighter condition quickly reduces to Θ µ 2 . This bound, together with eq. 46, implies the condition when γ ≥ 1/(2π), which is always true because γ is between 1 and 3 [12]. Even when µ 2 is as small as realistically possible (m p /m e = 1836) and γ as large as realistically possible (3), Θ must be ludicrously high (at least 97) to violate even this conservative validity condition.
Thus the large grain model is always valid when at least one of Θ or v is small. To find a regime where the model breaks down, we now consider the large v limit. For large v, and eq. 42 becomes Taking the inverse Lambert W function of both sides and cancelling common terms, Taking logarithms and rearranging, The right hand side is smallest when µ is smallest and γ and Θ are largest. Realistically, µ ≥ 42.85, γ ≤ 3 and Θ 1, so the RHS is at least 1.6. As such, setting the RHS to zero gives the tighter inequality Like the mid-sized grain model, the large grain model breaks down only in the face of exceptional flow (v ∼ 50).

VI. THE MASTER EQUATION AND GAUSSIAN APPROXIMATIONS
Although the stochastic models' master equations have no exact, analytic solution, we can follow Matsoukas, Russell & Smith [5][6][7] in finding approximate solutions by treating the models as if continuous. A one-step process has the master equation [11, p. 134] ∂ This one-step master equation is approximated well by the following Fokker-Planck equation when r N and g N are smooth, slowly varying functions of N [11, pp. 197-198 & 207-208]: For our models, the F-P approximation's conditions are satisfied when the equilibrium N is large. Except for the tiniest grains, the equilibrium N is approximately M, which is of order 1/α for mid-sized grains satisfying α ≪ 1 and of order (ln g * )/α for large grains. Thus the F-P approximation is a good one for mid-sized grains when α ≪ 1, and for large grains when α ≪ − ln g * . At equilibrium, eq. 57's left hand side is nil. This banishes f N (t)'s time dependence, so we write the equilibrium probability distribution as f N as before. Integrating both sides with respect to N, where s is a constant of integration corresponding to the relative probability current between charge states [9, p. 72]. At equilibrium this current is a constant, and for this system must be zero because N is bounded [9, p. 98]. Applying the boundary condition s = 0 and then the product rule, The integral is insoluble for both models, but approximate solutions are possible by linearizing y(N) about an N where most of f N 's probability density is concentrated. We could use the mode M but the algebra is tidier if we use the value Substituting into equation 62, where y ′ 0 ≡ y ′ (N 0 ) for brevity. This is a Gaussian probability distribution with mean N 0 and variance −1/y ′ 0 , so the final approximate probability distribution is Because y(N 0 ) = 0, the mean N 0 is implicitly defined by Inserting r N and g N for mid-sized grains (eqs. 18 & 20), This has the solution which is of similar form to the expression for M (eq. 33), and again of order 1/α for small α. For a Gaussian distribution the mean equals the mode, so it makes sense that N 0 ≈ M algebraically.
Substituting the large grain r N and g N (eqs. 18 & 21) into eq. 67 gives which has the solution similar to the large grain formula for M (eq. 41). Even simpler expressions for N 0 arise when the derivatives in eq. 67 are small compared to g N 0 and r N 0 , which occurs when α ≪ Θ 1. In that regime eq. 67 reduces to giving the solutions for mid-sized grains (becoming close to eq. 33 for very small α) and for large grains, which matches M (eq. 41) and makes the asymptotic O(1/α) dependence very explicit. Eq. 72 amounts to equating I i and I e , so eq. 73 implies the same normalized electric potential as the SOML equation [14, eq. 15], which comes from explicitly taking I i = I e . The more exact N 0 given by eq. 69 implies a more negative electric potential. The same simplification allows a concise approximation for y ′ 0 and so the distribution's variance. Neglecting derivatives, and so, applying the quotient rule and simplifying, Solving for y ′ 0 for mid-sized grains is tedious but feasible. Substituting in r N and g N , and applying eq. 73 eventually gives The variance is then for α ≪ Θ, so σ 2 ∝ 1/α ∝ T e a in this regime, consistent with Matsoukas and Russell's finding that σ 2 ∝ a [6, p. 4288] (α being proportional to 1/a). Also consistent is the implication that the normalized standard deviation σ /N 0 ∝ √ α ∝ 1/ √ a. While the dependence of σ 2 on α goes as O(T e a), the dependence on v is more complicated. nonlinearly increases with v. Gentle flows (v ≪ 1) have a negligible effect on σ 2 , but as the flow speed exceeds the ion thermal speed (v ∼ 1) σ 2 rises appreciably with v, plateauing at a higher value for very rapid flow. Eq. 79 systematically underestimates σ 2 , but not appreciably. For still faster flows the model breaks down because the sphere's chance of acquiring a positive charge becomes non-negligible. For large grains, Applying eq. 74, the variance is which aligns with the asymptotic 1/α dependence for midsized grains. Notice that σ 2 depends only on α for large grains. This remains the case even if we use the more exact value of N 0 given in equation 71: Using the more exact N 0 has the sole effect of making σ 2 slightly greater than 1/α; it does not reveal any dependence on µ, Θ, or v. We therefore reach the interesting conclusion that for a given large grain, the mean charge is sensitive to the values of the plasma parameters µ, Θ, v, and γ, while the charge's variance is not. The variance depends only on T e .

VII. SKEWNESS OF THE CHARGE DISTRIBUTION
The Gaussian approximation to f N roughly matches f N 's mean and variance, but ignores f N 's higher order moments. This may be problematic if f N has appreciable skew. A priori it is likely that f N is negatively skewed, since its support excludes the positive integers, cutting off f N 's upper tail. However, the degree of negative skew may be negligible for realistic parameter values. We explore this possibility for mid-sized grains first. To assess f N 's skewness, we compute it for a hydrogenic plasma with 1 eV electrons, as a function of the normalized flow velocity v and the sphere's radius a ( figure 3). Numerical experiments reveal that the skewness is less with heavier ions (figure 4), so the hydrogenic plasma results we discuss here are a worst-case scenario. Figure 3 shows decreasing skewness with increasing radius and flow speed. With equal ion and electron temperatures (Θ = 1) skewness is consistently small, except in the limit of vanishing sphere radius, but in that limit the model becomes invalid anyway as α ≫ 1. With cooler ions, skewness is greater and less affected by flow. The large grain model's f N is even less skewed for realistic parameter values. For large grains the skewness depends on the four parameters α, Θ, v and γ, but as it increases only marginally with γ over the range 1 ≤ γ ≤ 3 we may ignore γ.
Like the mid-sized grain model, the skewness decreases with increasing µ. Unlike the mid-sized grain model, the large grain model's skewness decreases with Θ (except when α ∼ 1, v is huge, and Θ is already very small). This fits our finding above that for cold ions the mid-sized grain model breaks down, while the large grain model improves its selfconsistency.
The foregoing means that f N 's skewness is highest for large grains when γ, Θ, α and v are high. Figure 6 presents numerical calculations of f N 's skewness where Θ and γ take on their highest realistic values (1 and 3 respectively). Even in this most pessimistic case, apprecible skewness is only a risk when v is enormous or the grain's size tends towards the nanometre scale, and is merely a symptom of the large grain model failing as its validity conditions are progressively violated. When those conditions are instead satisfied, the skewness of the predicted charge distribution is negligible, and the large grain model's Gaussian approximation appears to be even more robust than the mid-sized grain model's.

VIII. CONCLUSION
We have derived equilibrium probability distributions of a spherical grain's charge in a flowing, collisionless plasma, using stochastic models based on the SOML charging theory. It transpires that these distributions are expressible in closed form in terms of exponential and gamma functions. The modal grain charge is proportional to T e a for large grains, and remains approximately proportional to T e a for mid-sized grains with small α.
When the grain is large enough (and the ions are no hotter than the electrons, which is usually true) Gaussian distributions approximate the exact distributions well, affirming Matsoukas et al.'s demonstration that particle charge "fluctuations are Gaussian, regardless of the detailed form of the charging currents" [7]. For mid-sized grains, the Gaussian distribution's variance increases with the normalized flow velocity v, with the dependence on v strongest for transonic flows. For large grains there is no v dependence.
One possible use of the Gaussian approximation is estimating the likelihood of various deviations from the mean charge. For example, in a low-pressure argon discharge with n e = 10 16 m −1 and T e = 100T i = 4 eV, a 10 µm sphere (a ≪ λ D = 149 µm) has a 53% chance of being within 0.1% of its equilibrium charge, while the same sphere in a tokamak edge deuterium plasma (n e = 10 20 m −1 , T e = T i = 100 eV, and a > λ D = 7 µm) has a >99% chance of being within 0.1% of its equilibrium charge.
There are many ways in which this paper's results could be extended. For example, we assume our spheres are conducting yet remain unpolarized in the face of approaching electrons and ions. Physically this is a contradiction in terms, but it simplifies our calculations and has a substantial effect only on tiny spheres (a 10 nm) in plasmas with cool electrons (k B T e 1 eV) [7]. It might nonetheless be useful to incorporate electrical polarization into our models. Our models might also be generalizable to collisional plasmas, grains out of equilibrium, non-spherical grains, plasmas with non-Maxwellian electron velocity distributions, grains in a sheath, and grains that emit electrons (whether by photoelectric, thermionic, or field emission). As our models stand, however, they should remain applicable to a broad range of close-to-spherical grains in typical flowing plasmas. For a negatively charged sphere in a singly ionized plasma, the SOML ion current is which is equation 2 above with Z = 1. (There seems to be a typo, incidentally, in Willis et al.'s presentation [14] of this result. Their equation 13 lacks an exponent of 1/2 for the parenthetical fraction.) However, a tiny sphere has a non-negligible chance of being positively charged. Under SOML a positively charged sphere has the ion current [10,14] where the new velocity-dependent auxiliary functions are with v 0 ≡ eφ /(k B T i ) as a normalized speed cutoff. These auxiliary functions are complicated functions of φ via v 0 , and this blocks an exact analytic solution for f N where N > 0. We can and do ignore this for larger spheres as they're almost always negatively charged, but for spheres small enough to attain a positive charge it is an insuperable obstacle. Not all is lost: one can derive a complete solution for f N with plain OML, i.e. when there is no flow. In this special case, one can solve a recurrence relation for f N when N ≤ 0, solve another recurrence relation for f N when N ≥ 0, and graft the two solutions together with the detailed balance condition r 1 f 1 = g 0 f 0 . We have not done this here as the mechanics of that calculation are little different to those of section IV, and Draine & Sutin have already given an analogous solution (albeit assuming T i = T e ) for small grains and low-temperature plasmas [1, p. 808].

Appendix B: Unimodality of the charge probability distribution
Here is a demonstration that f N is unimodal in the sense of Medgyessy [8], i.e. that the sequence has exactly one change of sign after discarding zero terms. In intuitive terms, this asserts that f N has only one peak, though that peak may spread across multiple adjacent abscissae with the same maximal ordinate. As this paper's stochastic model assumes a non-positive charge a priori, f N = 0 ∀ N > 0. Hence the term f 0 − f 1 = f 0 in sequence B1, and the terms before it are zero and dispensable. Therefore we need only show that has one change of sign after discarding zero terms. We now prove this explicitly for mid-sized grains; the same basic logic applies for large grains. Consider f N−1 / f N , given in equation 23. Because α > 0, exp(αN) strictly increases in N. Similarly, because s 1 (v), s 2 (v), and α are always positive, eq. 23's parenthetical divisor strictly decreases in N. As the parenthetical divisor must always be non-negative, it follows that eq. 23 is strictly increasing in N for N ≤ 0. Because f N−1 / f N strictly increases in N for N ≤ 0, f N−1 / f N − 1 and hence f N−1 − f N can change sign at most once for N ≤ 0. From the second term onwards, then, sequence B2 has at most one sign change.
Suppose there were such a sign change. This requires that f N−1 / f N − 1 changes sign as N becomes more negative, which means f N−1 / f N must go from being more than 1 to being less than 1, because f N−1 / f N decreases as N becomes more negative. This implies that f −1 / f 0 > 1, implying f −1 − f 0 > 0, which in turn implies no sign change in sequence B2's first two terms (because f 0 > 0). As such, if sequence B2 has a sign change after the second term, it is the only sign change.
Suppose there were instead no sign change after the second term. Then either f N−1 / f N > 1 for N ≤ 0, or f N−1 / f N < 1 for N ≤ 0. The former is impossible, because it asserts that f N becomes ever larger as N → −∞, which would render f N unnormalizable. The latter implies f −1 < f 0 , and so a sign change between sequence B2's first two terms. Thus, were there no sign change after the sequence's second term, there would have to be a sign change between the first two terms.
The last two paragraphs mean that sequence B2 has exactly one sign change, completing the proof that f N is unimodal.