A model for a Langmuir sheath in a stagnating dense plasma with secondary ion formation

This simplified model provides solutions for the current-voltage characteristics of a sheath in a dense flowing plasma when surface chemistry contributes secondary ions. The problem is motivated by the recent discovery that strong transient signals in industrial ion current sensors are caused by chemical reactions with carbon in the steel being cut or welded by oxyfuel processes. The one-dimensional model considers a quasi-uniform dense plasma flowing towards and stagnating on an absorbing surface, above which there is a source of secondary ions. Because the secondary ions are formed directly in the plasma sheath, they have strong impacts on the current-voltage characteristic. With ionic Reynolds number, R, and integral length scale, α, secondary ion formation rate, Ω, and length scale, β, saturation currents are simply R + βΩ until β ≪ 1, at which point, new electrons cannot escape the sheath, and secondary ions have no effect. Floating potential, ϕ ∞, scales like exp(ϕ∞)∝R−3/4 , and secondary ions have little impact unless β 2Ω > 1. Even then, floating potential is only weakly affected by secondary ion formation. The integral length scale, α, is not found to strongly affect the results.


( )
, and secondary ions have little impact unless β 2 Ω > 1.Even then, floating potential is only weakly affected by secondary ion formation.The integral length scale, α, is not found to strongly affect the results.

Introduction
This simplified model provides solutions for the current-voltage characteristics of a sheath in a dense flowing plasma when surface chemistry contributes secondary ions.The problem is motivated by the recent discovery that strong transient signals in industrial ion current sensors [1] are caused by chemical reactions with carbon in the steel being cut or welded by oxyfuel processes [2,3].In these systems, computer controlled torches issue an oxygen-fuel-gas flame to heat a steel work piece like in Figure 1.Ion current sensors interrogate the currentvoltage characteristic of the flame between the torch and work piece to deduce process parameters, so poorly understood variations in the observed electrical properties are deeply problematic.So far, it is apparent that carbonaceous molecules evaporated from the steel result in chemions in the region immediately above the solid surface, but how they interact with the sheath there is not yet clear.Figure 1 is a representation of the oxyfuel flame over a steel work surface.Conical flame fronts anchor on the torch tip, and the luminous outer cone stagnates on the work surface.Carbonaceous molecules are supposed to diffuse through the steel to its surface, where they sublime.Despite more than a half-century of research on the topic, carbon sublimation remains nuanced in materials far simpler than steel alloys [4,5], but it is known that carbonaceous molecules can be found in the vapor phase at temperatures typical in the oxyfuel flame (3000 K to 3300 K).This provides some justification for the idea that ion formation is in the gas phase, but there is an alternative view that the reaction begins at the surface with carbon still in the solid phase.In either case, for there to be an asymmetrical transport of positive and negative charge, at least one of the charged products must enter the gas phase, which will affect the space charge in the sheath.
In a previous work [6], the authors give an analytical solution for ion transport in one dimension throughout the domain.Though the original model includes some simple treatment of the sheath at the metal surface, it ignores secondary ion generation and the non-uniform velocity distribution above the work surface.A preliminary numerical investigation also predicts enhanced current signals with surface reactions [7], but even with computational tools, it is difficult to thoroughly study the effect without a theoretical basis to reduce the parameter space into dimensionless groups.The present study is entirely devoted to the sheath at the metal surface with the goals of (1) establishing dimensionless groups, and (2) establishing corresponding predictions for the phenomena observed in the sheath.
There are three natural length scales in these problems [8], pp.406, with a fourth unique to this particular problem: (1) the mean free path between molecular collisions, l ; (2) the Debye length, λ; (3) the integral length scale, L, representing the macro-scale physical features of the system; and (4) the length over which secondary ions are formed above the surface.In low-pressure plasmas common to gas discharge tubes, l l  , so individual charged particle trajectories may be modeled without regard for collisions.In atmospheric flames, l L l   , so in addition to diffusion and electrical mobility, it is essential to model advection with the bulk flow as well.
There has been tremendous progress establishing chemical mechanisms for the formation of these ions in flames starting in the mid twentieth century [9,10], so the pathways by which carbon leads to chemions are now well understood.Contemporary works in combustion literature focus on chemical kinetic mechanisms, modeling of secondary and tertiary charged species like - O 2 , and ionic wind (high voltage) [11, 12].Probably the most appropriate parallel experimental work was a 2019 study for ion current sensing in furnaces by Chien et al. [13].
In the sheath, transport of hydronium (H 3 O + ), the free electron (e − ), and the any other relevant charged species must be resolved at the Debye length, with longer length scales dominating in the free stream.As we will examine in detail, the Debye length is on the order of 10 μm, and the equations are highly nonlinear and numerically stiff near metal surfaces.The complexity of this problem is uniquely compounded by the secondary source of ions very close to the lower boundary [1], which can vary independently of the source of ions in the flame.Our model is specifically informed by the works of Su and Lam [14], Clements and Smy [15], and similar works that followed.
Though little is known about the exact chemical mechanisms for this particular solid-phase process, we will adopt a two-species model that considers only a representative heavy positive ion and the free electron.These models are popular in flames because the H 3 O + and e − pair usually dominates other charge carriers by orders of magnitude.Even if the positive charge carriers are more diverse, it will be safe to assert that they are heavy compared with the free electron.The identity of the negative charge carrier is more nuanced since it is the severe imbalance in properties between positive and negative charge carriers is critical to so many of the phenomena we study.As this model is scrutinized in future works, it may become clear that a three-species model is essential, but for now, we limit ourselves to two species.

Model
We imagine a flow in which a quasi-uniform plasma flows downwards towards a flat surface.At the surface, the velocity is zero, and the flow begins to decelerate some distance above.Meanwhile, carbonaceous molecules sublime and diffuse into the plasma from the surface.For the purposes of the present model, we neglect horizontal spatial information, and we are not concerned with the details of the chemical pathways by which secondary ions are generated.
If horizontal ion concentration gradients are small in the region of interest, the problem can be reduced to a single vertical dimension despite the non-uniform horizontal velocity components.The focus of the model is on the fine details of the sheath very close to the work surface, so we will regard the flame front, in which primary ions originate, as far away in the "free stream." For a species k, the upward-traveling flux density, Γ k , along an axis, x, in number per area per time is These terms are transport due to the vertical bulk velocity component, U, diffusion due to mobility, μ k , and drift due to electrical forces on the ion's charge, q k .Each species is assigned its own non-equilibrium temperature, T k .Finally, k b is the Boltzmann constant.The conservation of mass dictates that the gradient of flux density of a species is equal to the volumetric formation rate, w¢¢¢ f  , When there is only a single kind of cation and anion, they are formed in equal portions, so there is no need to distinguish between individual species' formation rates.The Poisson equation governing the electric potential, V, is when the number densities are due to the positive ion, n i , and the free electron, n e .The permittivity is ò 0 , and fundamental charge is e.

Nondimensionalization
Near the metal surface, there is a length scale, λ, over which the essential sheath physics occur.We choose to adopt a dimensionless length, z, that maps z ä [0, 1) to x ä [0, ∞ ) while resolving the essential length scales.The infinite free stream (where gradients vanish and transport is uninteresting) is compressed to the domain near z = 1.Meanwhile, near z = 0, the dominant length scale is well resolved, and the dimensionless coefficient, γ, determines the z scales at which the transition occurs.When γ > 1, the center of the z-domain reaches into the free stream where x > λ, and when γ < 1, the center of the domain is shifted deep into the sheath where x < λ.
Most of the remaining physical scales of the problem are established by the free stream fluid properties.Here, we introduce dimensionless velocity, u, formation rate, w, number densities, η and ν, field strength, ψ, and voltage, f.
The free-stream velocity and number densities, U ∞ and n ∞ , are constrained by distant aspects of the flow field, and are treated as external parameters.The peak secondary ion formation rate, w¢¢¢ f ,0  , and voltage scale, V 0 , are also treated as external parameters, but these will be shown to be determined by other aspects of the flow.Finally, note that ψ is actually proportional to the derivative of voltage, which is the negative of field strength.
When these dimensionless parameters are substituted into (3) and (2) written for ions and free electrons, we obtain three governing equations, The ionic Reynolds number, R, electron temperature ratio, τ, electron mobility ratio, μ, and dimensionless secondary ion formation rate, Ω, are defined as ( ) To obtain this result, the voltage scale is the electron thermal energy,

( ) 
It is important to note that (11), (12), and (13) include the assumption that ion temperatures and mobility are uniform.While not strictly true in the presence of a cold surface, when nonuniformities in the temperature field are relatively slight, they merely have the effect of exaggerating or diminishing local transport length scales without substantively changing the vital physics of the solution.

Velocity and ion formation functions
At the work surface, the vertical bulk velocity is zero, and very far from the work surface, it will be equal to the free stream velocity, U ∞ .The exact function and scale over which the velocity declines are determined by the details of the flow field.So, to build a model as generally applicable as possible, we approximate the transition with an exponential with characteristic length scale α Debye lengths thick.So, As we will see in later sections, the model is not very sensitive to this function, so an approximation with a length scale that can be selected for a given flow field is justified.
Once nondimensionalized, we obtain Note that u has been formulated with the assertion that U ∞ is in the negative direction.
The secondary reaction rate is due to the reaction of carbonaceous molecules that have diffused from the surface.If this were presumed to be a simple reaction, the length scales over which they are found will be a balance of their diffusion against their reaction rate.If w¢¢¢ f 0  and n 0 were the reaction rate and carbonaceous molecule number density at the surface, secondary ion formation would decline proportionally with carbonaceous concentration, n c , so This represents a formation region that decays away over a scale β Debye lengths thick.When nondimensionalized, When γ/β or γ/α is unity, these curves are simply linear on z, but as β and α shrink, these curves compress to form a rapid gradient very close to the metal surface.These are, therefore, dimensionless measures of the secondary ion formation and velocity length scales.

Reasonable parameter values
Oxyfuel flame temperatures are on the order of 3000K or higher, and ion densities in the outer cone of the flame have been measured to be on the order 10 17 m −3 [16].Under these conditions, the Debye length is on the order λ ≈ 10 μm.
The electrical mobility of the positive ion has been approximated to be μ i e ≈ 10 −3 m 2 /Vs [17], and based solely on their dissimilar masses the electron mobility is at least 200 times that.Velocities in the outer cone of the oxyfuel flame have been approximated to be on the order of 50m/s [7].Under these conditions, the practical maximum for the ionic Reynolds number is about 2.
The major flow features will vary on the integral length scale, but the decay to zero velocity will occur at viscous length scales.By adopting μ i k b T i in place of the kinematic viscosity in the formulation of the Reynolds number in (14), we have already implicitly asserted that viscous length scales are the same as diffusive length scales, or that the Schmidt number is approximately 1. Therefore, α should be not far from 1.
Less is known about the carbonaceous molecule and its properties, so it is more difficult to make substantive statements about reasonable limits on β.However, it can be immediately observed that When β 2 Ω is near unity, the carbonaceous molecules responsible for forming secondary ions are diffusing into the sheath at roughly the same rate as primary ions from the free-stream.So, when β 2 Ω = 1, secondary ion formation is probably irrelevant.On the other hand, when β ? 1, secondary ion formation is essentially uniform throughout the domain.In this case, the ion concentration above the work surface would diverge since there is no mechanism to reach a steady-state concentration.We disregard this case not because it is impossible, but because (1) its behaviors are already addressed in our prior paper [6], and (2) only ion formation inside the sheath will cause physical effects that are of interest in the present study.

Net electric current
In the study of a sheath, the natural limiting transport scale is diffusion of the positive ion, so the dimensionless flux density is best written Electrical current density in Amperes per unit area, is Since gradients of F i and F e are equal to the dimensionless ion formation rate, the difference of their gradients, ¢ J , must be equal to zero.Therefore, current density is the same everywhere in the domain, though positive and negative flux are free to change in tandem.
Because total ion current density is constant everywhere in the domain, when we are ready to calculate it from a solution, we may choose a place in the domain where the calculation is simplest.In the free stream (z → 1), there is no ion formation and velocity is uniform, so the flux of ions is The net electrical current is their difference, e i e i 2

Sheath potential
It is a deeply intuitive result that any net current conducted through an infinite dense plasma requires a nonzero electric field to overcome the plasma's bulk resistivity.Unfortunately, this also means that there is no such parameter as the "plasma potential."In classical probe theory, the probes are spheres, cylinders, or plates of finite length, so the current density is allowed to vanish in the free stream, leading to a uniform voltage field.That is not true in a 1D problem, in which net electric currents conduct through an infinite medium.In the real system, the plasma between electrodes is a finite length with a finite electrical resistance, but this model has deliberately neglected these details.
To isolate voltage drops across the sheath from ohmic voltage drops, we define the sheath potential as For all solutions where ψ approaches ψ ∞ faster than (1 − z) approaches zero, this integral converges.

Boundary conditions
Very close to the metal surface, slow diffusion starves the plasma of ions there, so they are rapidly depleated by collision with the surface.Therefore, their density at z = 0 is commonly approximated to be zero [14].
Meanwhile, far from the surface (z = 1), the ion density is equal to the free stream density.Therefore, There are two types of electrical boundary conditions that might be applied: constant voltage or constant current.A current condition is easily enforced by specifying Following the discussion above, constant voltage conditions can be established by specifying a value for the integral in (34).Fortunately, in a numerical approach to the problem, this is no more challenging than specifying any other boundary condition.

Numerical solution
The model was evaluated using a finite difference approach with a quadratic interpolation scheme.While the model was coded to permit a non-uniform grid spacing, it was found that using γ = 10 with a uniform 200-node grid adequately resolved the important features.Grid sensitivity studies showed that repeating solutions with a 1000-node grids did not appreciably change the results.Once discretized, the governing equations can be expressed in a quadratic form.When X  is a serialized vector of all node values for η, ν, and ψ, the problem is expressible with constant, linear, and quadratic terms, Here, we are using tensor notation so that a i b i is a shorthand for the sum, ∑ i a i b i .The tensors are error, E  , a vector of constants, C  , a linear tensor coefficient, L, and a quadratic coefficient tensor, Q.It was found that a simple Newton scheme consistently converged in about ten or fewer iterations.Here, the change in the solution vector, D


, is calculated by linearizing about the current guess, Note that L is symmetrical, and because of the transpose operation on Q, the resulting matrix inversion problem will also be symmetrical.
The only remaining nuance is implementing a discretized form of the boundary condition in (36).When using trapezoidal integration with a uniform step size, δ, the contribution of each node to the first integral is γδ/ (1 − z).We presume that the integral converges, so that the singular value at z = 1 can be neglected.Meanwhile, the negative of δ/(1 − z) is integrated over the entire grid, which becomes the coefficient for ψ at z = 1.At low sheath voltages, the rapid diffusion of electrons into the surface drives strong negative currents, requiring a uniform electric field in the free stream (on the right) to overcome the plasma resistivity.As the sheath voltage rises, saturation occurs quickly, marked by the utter depletion of electrons (dashed line) in the vicinity of the metal surface (z = 0).The small positive currents that still flow are marked by a weak electric field in the free stream.Starved of a source of electrons, additional applied voltage merely strips away electrons from a wider sheath, intensifying the electric field there.

Sheath Formation
In Figure 3, it is clear that the addition of secondary ion formation enhances the ion density in the sheath.It is less obvious, but the electron density in the sheath is also enhanced.This observation is important for understanding results presented later in this section.Note that solutions for f ∞ at 0 and 1 are omitted in Figure 3 to keep the plot well scaled.These curves are pushed more strongly negative by the accumulation of charge in the sheath.

Current-Voltage Characteristic
Figure 4 shows the current-voltage characteristic for the sheath by plotting current density, J, against the applied sheath voltage, f ∞ , for various dimensionless ion formation rates.Saturation can clearly be seen beginning around f ∞ > 8 (above roughly 2V in a flame).Intuitively, saturation currents increase with the addition of secondary ions, but it also appears that the floating potential (voltage when J = 0) decreases.

Saturation Current
When the applied voltage is high enough to evacuate all electrons from the sheath, additional positive ions formed there merely add to the ions convected from the free-stream.From (1), the ion flux density at the work surface is

( ) 
Recall that Γ will be negative when transport is down into the surface, so the negative before the integral has the effect of enhancing ion flux.If the flux due to ion mobility is small compared with convection, the dimensionless form of this is simply ( ) The first term, R, is the dimensionless rate of ion delivery in the free stream.The multiple of the formation intensity, Ω, and thickness, β, is the total secondary ion formation rate.
Figure 5 shows J at f ∞ = 40 plotted against βΩ for a wide range of values for Ω and β.When the formation thickness is much wider than the sheath thickness (β ?1), J ≈ − F i,0 , confirming that electrons do not contribute significantly to the current at the surface.However, when secondary ions are formed in a layer very close to the surface, a dimensionless sheath potential of 40 is still not sufficient to prevent electrons from diffusing into the surface, inhibiting the electrical current the same formation rate would have otherwise produced.
Simulations in which α were varied show no impact on the saturation current.Once ions are delivered by the free stream, it does not matter how far from the surface the fluid decelerates-they will eventually diffuse to the surface, where they will contribute to the saturation current.

Floating Sheath Potential
When no net electrical current flows, the potential that develops between the surface and the free-stream is the floating potential.Figure 6 shows the floating potential versus the Reynolds number for various values of α and with no secondary ion formation.
Even without an explicit analytical solution for f ∞ , it is possible to find a form that approximates the behavior in Figure 6.Because the net electric current is zero, there are three special constraints that may be added  Though it does not deserve separate treatment here, it is important to note that this is written in a linear domain instead of the logarithmic domain used throughout the rest of this paper.Using an integration factor, f f -¥ exp( ) , this can be manipulated to yield This approximation is no longer valid as soon as the electric field vanishes, and convection of electrons resumes its importance in the free stream.At the metal surface, ν = f = 0.At the sheath's edge, z = z s , ν = 1 and f = f ∞ , so Near the metal surface, f exp( )in the integral begins at 1 and shrinks as f grows to the free-stream value.The value of the integral, then, is mostly determined by the sheath's thickness.Figure 7 shows the spatial solutions for the case when α = 1, which shows increasing Reynolds numbers "blowing" the sheath closer to the metal surface.It can be determined, therefore, there is an inverse relationship between R/μ and the value of this integral.
Using this line of reasoning, the integral might reasonably be approximated by A(R/μ) −m , which leads to an empirical form where m may be between 0 and 1.The line in Figure 6 is drawn using A=0.16 and m=0.25.The approximation is worst when α is small and R is large-the case where neglecting convection of electrons in the sheath is most suspect.
As with saturation current, the velocity length scale, α, has little impact, while the velocity, R, is paramount.It may appear alarming that the floating potential tends to infinity at low velocities, but that is a peculiarity of the semi-infinite sheath in which ion formation is neglected.Without a source of fresh ions, the surface would preferentially absorb electrons until the sheath grew to be infinitely large, causing an infinite floating potential.This mechanism for determining the size of the sheath is reflected in our choice for an empirical approximation to floating potential.
The impact of ion formation on the floating potential can be seen in Figure 8.Here, f ∞ is shown with Ω ranging from .025 to 1000 while Reynolds number and α are both held at 1. Despite the wide range of formation rates, the decrease in floating potential is very little-only about -2 (d-less) or -0.5V.
When the horizontal axis is scaled to be β 2 Ω, these curves collapse from a vast scatter into a relatively tidy cluster.As predicted, when β 2 Ω = 1, the order of secondary ions in the sheath will drop below primary ions, so the floating potential resumes its value at R = 1 in Figure 6.Note that the scatter between the curves is tightest when β is small.As the secondary ion formation reaches out into the free-stream, only the portion near the sheath actually impacts the floating potential.
It is potentially counter-intuitive that secondary ions should reduce the floating potential, but just as we see in Figure 5, electrons that spawn in the sheath can have a dramatic impact on the IV characteristic.The positive ions are primarily transported by convection, but electron mobility is roughly 200 times that of positive ions, so the electric field in the sheath must collapse dramatically to bring their flux into balance if the number of ions there is increased.

Figure 1 .
Figure 1.A diagram of the oxyfuel torch flame geometry.
the length scale, λ, is the familiar Debye length

Figures 2 and 3
show superimposed solutions for how the sheaths grow with applied voltages.In each figure, conditions are held constant while incrementally increasing the sheath voltage, f ∞ .Since these figures use γ = 10, values at z=0.1 correspond to approximately one Debye length.

Figure 4 .
Figure 4. Current-voltage characteristic for the sheath with various secondary generation rates (Ω).
to the problem: (1) the flux of electrons must be equal to the flux of positive ions, (2) the free stream electric field is zero, (3) and the plasma potential converges without the special dispensation of (34).The flux of positive ions is simply F i = − R everywhere.The flux of electrons is dominated by diffusion and electrical mobility, so

Figure 5 .
Figure 5. Current at f ∞ = 40 versus secondary formation rate and length scale.The asymptotic thick formation rate solution is shown as a solid line.