Empirically extending 1D Child-Langmuir theory to a finite temperature beam

Numerical solutions to the 1D steady-state Vlasov-Poisson system are used to develop a straightforward empirical formula for the electric current density transmitted through a vacuum diode (voltage gap) as a function of gap distance, gap voltage, the injected current density, and the average velocity and temperature of injected particles, as well as their charge and mass. This formula generalizes the 1D cold beam Child-Langmuir law (which predicts the maximum transmitted current for mono-energetic particles in a planar diode as a function of gap voltage and distance) to the case where particles are injected with a finite velocity spread. Though this case is of practical importance, no analytical solution is known. Found by a best-fit to results from particle-in-cell (PIC) simulations, the empirical formula characterizes the current transmitted across the diode for an injected velocity distribution of a drifting Maxwellian. It is not meant to yield a precise answer, but approximately characterizes the effect of space charge on transmitted current density over a large input space. The formula allows quick quantitative estimation of the effect of space charge in diode-like devices, such as gate-anode gaps in nanoscale vacuum channel transistors.


I. INTRODUCTION
The space charge limited (SCL) current between two electrodes is a problem of great interest with applications to nanoscale vacuum channel transistors, 1 semiconductor diodes, 2,3 ion beams, 4 x-ray detectors in fusion devices, 5 magnetron sputtering, 6 and ion thrusters. 7The simplest version of the SCL current problem is the 1D planar geometry where charged particles are injected with zero velocity at one electrode and are subsequently accelerated by the applied gap voltage to the other electrode.There exists a maximum steady-state current (the SCL current) beyond which the self-fields of the particles in the gap prevent injected particles from entering the gap, precluding a steady-state solution.The SCL current for the zero injection velocity case is given by the well known Child-Langmuir (CL) law: where J CL is the SCL current density, ϵ 0 is the permittivity of free space, q is the charge of the injected particles, m is the mass of the injected particles, Φ D is the gap potential, and D is the gap distance.(We note that Eq. 1 only gives the CL law for qΦ D < 0, that is, an accelerating gap potential, but we use the nonstandard notation including absolute values for later generalization when qΦ D ≥ 0 for the case of finite injection velocity.)For an injected current density magnitude (|J i |) less than |J CL |, a) jesse.snelling@colorado.edu the magnitude of the electric field at the injection electrode is reduced by space charge, but particles are still accelerated into the gap.For J i = J CL , the electric field at the injection electrode goes to zero and no particles enter the gap.This problem was first examined by Child 8 and Langmuir. 9Since then, it has been studied extensively to account for associated instabilities, 10,11 higher dimensional geometries, 12,13 relativistic regimes, 14 quantum regimes, 15 and time-dependent boundary conditions. 16here exists a simple analytical extension 17 for the SCL current density for the case of particles injected with the fixed velocity v i > 0 (a cold beam): where J CL is given in Eq. 1 and β ≡ mv 2 i /2|qΦ D | 1 2 .We note that J cb is monotonic in β.Additionally, if v i → 0 (β → 0) then J cb → J CL ; thus β quantifies the effect of injection velocity on the SCL current.When mv 2 i /2 < qΦ D , particles are energetically forbidden from crossing the gap and we see that Eq. 3 yields a complex value.
A natural next step in extending the 1D diode SCL current problem is to include a thermal velocity spread, v th , in the description of the injected velocity distribution, which allows the case where only a portion of the particles are reflected by space charge.However, a general analytical solution to this problem has not yet been presented.In this paper, we utilize simulation to pro-vide an approximate formula for the transmitted current density as a function of particle and gap characteristics, taking into account the effect a thermal particle distribution has on space charge.
We will begin with a precise definition of the finite temperature (warm) beam 1D diode problem in Section II.In Section III, we reduce the large parameter space to just three dimensions, which makes a numerical scan over the input space feasible.In Section IV, we describe the individual simulations.The procedure for numerically scanning over the parameter space and fitting the simulation results is presented in Section V.The model parameters along with the model's performance are presented in Section VI.Finally, we provide a brief summary of the paper in Section VII.The appendix gives a recipe for straightforward implementation of the empirical model in general and in the simplified special case of a "nearly" cold beam.

II. PROBLEM DEFINITION
The aim of this paper is to find steady-state solutions of the one-species 1D Vlasov-Poisson system for phase space density f (x, v, t) and electric potential Φ(x, t) with the distribution of injected current given by j i (v) in a gap of length D with applied voltage Φ D .The equations describing such a system are: where the final equation gives the current density boundary condition (BC).We impose particle absorption boundaries at the electrodes at x = 0 or x = D. Nominally, this problem has inputs q, m, D, Φ D , and j i (v) with the output J t , the steady-state current density that is transmitted to x = D.
While the Vlasov-Poisson system above describes the problem for arbitrary j i (v), for the sake of practicality, we restrict our study to a single distribution function.A Maxwellian distribution is common in many applications and is a standard choice to capture the first and second moments of a distribution when modeling.In this paper, we consider only the constant injected current density distribution function j i (v) of the form of a truncated drifting Maxwellian with drift velocity v D > 0, thermal velocity v th > 0, and injection current density J i : where v norm is given by: Because the Maxwellian is truncated, the mean injection velocity v i differs from v D .In principle, v D can be zero or negative, though we do not explore these cases.Regardless, we find it more convenient to work with v i instead of v D , so we employ the mapping between the two to formulate the problem in terms of the mean injection velocity, which is given by: With this current density distribution function, there are a total of seven scalar input parameters, q, m, D, Φ D , J i , v i , and v th .However, in the following section, the seven inputs are nondimensionalized to just three dimensionless parameters.This problem is not known to be analytically solvable, so we utilize time-domain particle-in-cell (PIC) simulation to find steady-state solutions.

III. REDUCTION TO THREE DIMENSIONLESS PARAMETERS
In general, the input space for this 1D diode SCL current transmission problem comprises seven parameters: q, m, D, Φ D , J i , v th , and v i .These seven dimensional parameters involve four basic units (charge, mass, length, and time), so by dimensional analysis, solutions to Eqs. 4 can be characterized by just 7 − 4 = 3 independent dimensionless parameters. 18I.e. by choosing normalizations for the four basic units, Eqs. 4 can be nondimensionalized so that q, m, D, Φ D , J i , v th , and v i do not appear, except in combinations of the three dimensionless parameters.Doing so vastly simplifies the empirical description, but maintains full generality since any desired result can be obtained from a function of just the three dimensionless parameters.While any reasonable normalization yields a similar simplification to a 3D parameter space, choosing normalizations that yield values near unity can help to provide physical insight and simplify empirical fitting.We normalize charge, mass, and length to q, m, and D, respectively; we normalize velocities to v i , or equivalently, time to D/v i .
Unfortunately, the straightforward normalization for current density, qv i /D 3 , is not very meaningful.Instead, we normalize J i to J norm , which is carefully chosen to reduce to the cold beam SCL current in Eq. 2 (i.e.lim v th →0 J norm = J cb ) while being well defined for all v th > 0.
To construct J norm , we first define the average injection velocity of particles that will traverse the gap in the absence of space charge (v it ) and the final velocity of a typical particle crossing the gap in the absence of space charge (v f ): where the reflection velocity v r = 0 for qΦ D < 0 and Comparing Eq. 10 to Eq. 2 shows significant similarities.
In particular, if v it were simply replaced with v i , then these two equations would be identical.Consequently, since lim v th →0 v it = v i , J norm = J cb in the cold beam limit.Thus J norm has been defined so Eq. 10 reduces to Eq. 2 as v it → v i , which occurs as v th → 0. With this built in, the model presented later in Sec.V inherits properties in the v th → 0 limit corresponding to the cold beam solution.
Now we construct the three dimensionless parameters: Every possible scenario for the warm steady-state 1D SCL current problem with the injected current density distribution given in Eq. 5 is represented by some combination of these three quantities.
We nondimensionalize the solution output J t with J max (the current able to cross the gap in the absence of space charge), resulting in a measure of the fractional current transmitted across the gap, F: |J max | is less than |J i | for decelerating voltages (qΦ D > 0).F → 1 corresponds to the weakly SCL regime (i.e.emission limited regime where J t ≈ J max ).On the other hand, F → 0 corresponds to the strongly SCL regime where most of the current is reflected rather than transmitted.For the cold beam case (v * = 0), F is just a step function (at j * = 1): For a warm beam (v * > 0), however, the transition of F between 1 and 0 is gradual, reflecting the smoothness of the velocity distribution.Ultimately we have reduced the problem to determining the dimensionless function F(j * , ϕ * , v * ) of three dimensionless parameters.
If F(j * , ϕ * , v * ) is known, then trivial scaling will yield the physical output J t (q, m, D, Φ D , J i , v i , v th ).Since this work only addresses distributions in the form of Eq. 5 with a positive drift velocity, there is a straightforward mapping between v i and v D to give the equivalent output IV. SOLUTION VIA PIC SIMULATION FOR PARTICULAR (j * , ϕ * , v * ) We performed standard 1D electrostatic PIC simulations with the VORPAL 19 code distributed in VSim-12. 20First, physical parameters were determined from (j * , ϕ * , v * ).Without loss of generality, we simulated electrons with q = −e and m = m e , a gap distance D = 150 × 10 −6 m, and drift velocity v D = 4.19 × 10 6 m/s.Then, v th , Φ D , and J i were uniquely determined from chosen j * , ϕ * , and v * .Due to nondimensionalization, the particular choice of physical parameters that give the dimensionless parameters is essentially an arbitrary choice of units.Only cases with different dimensionless parameters are meaningfully distinct in this context.Boundary conditions were specified by Φ(0, t) = 0 and Φ(D, t) = Φ D .
The simulation grid cell size ∆x was chosen based on an estimation of the minimum Debye length.First, a conservative nominal number density was calculated as n nom = 100J i /qv D .
Then, ∆x = min (ϵ 0 mv 2 th )/(q 2 )(n nom ), D/150 .(The upper bound of D/150 was determined by convergence studies of multiple cases.)After the simulation was run, the actual minimum Debye length λ D = (ϵ 0 mv 2 th )/(q 2 )n max was calculated from the maximum steady state number density n max to ensure the minimum Debye length was re-solved.If it wasn't, an additional simulation was run to resolve it.For all but nine simulations used in this study, λ D /∆x > 2.5.
Due to computational expense, these other nine simulations were unable to be resolved to the same degree and have 0.44 < λ D /∆x < 1.88.Fortunately, these simulation points are not clustered together in the (j * , ϕ * , v * ) parameter space, meaning the majority of their nearest neighbors are well resolved.These few unresolved simulations yield results for F that differ from the linearly interpolated value from their nearest neighbors by less than 0.1.Therefore, these points would not significantly affect the fit, and we chose to leave them in the analysis.
All simulations were run for 50 drift velocity crossing times, i.e.T = 50D/v D .(We observed that an apparent steady state was reached by time 4T /5 in all considered cases.)At each timestep, a number, n particles/step , of equal weight macroparticles were emitted with velocities randomly sampled from the truncated Maxwellian current density distribution j i (v) in Eq. 5 further truncated to four standard deviations about v D , |v − v D | < 4v th .Note that this neglects the low velocity (v < v D − 4v th ) and high velocity (v > v D + 4v th ) tails of the velocity distribution.The velocity distribution is numerically sampled at each time step and particles are injected at the plane of the boundary at a uniform, random time within ∆t.We ensured that sufficient macroparticles were emitted per cell per step such that neglecting each tail is equivalent to omitting at most a single macroparticle representing 2% of the transmitted current.There is then a maximum macroparticle velocity within the simulation, , where v hi = v D + 4v th is the largest injected velocity.The timestep was chosen to prevent particles with velocity v max from crossing a cell, i.e. ∆t = ∆x/v max .The transmitted current density was averaged over time interval [4T /5, T ] to yield a measurement for J t in steady state with associated standard deviation.
For cold beams above the SCL, it is known that the transmitted current density can oscillate in the time domain.However, simulations have shown that introducing a small thermal velocity can heavily damp these oscillations. 21If present in our simulations, any oscillations cannot have amplitudes larger than the standard deviation of J i , which was on average 0.06J i and does not exceed 0.11J i for the measurements used in this pa-per.Differentiating oscillations from particle noise and subsequently characterizing them requires further study.The ultimate result of the simulation is J t , or equivalently, at F sim (j * , ϕ * , v * ) = J t /J max .
We observed that, for small v * , the time to reach steady state can increase dramatically when j * ≈ 1.Although we have not explored this in quantitative detail, we speculate that this happens because, as v * becomes small, F approaches a step function at j * = 1, where it is highly sensitive to small changes in j * .Thus F is sensitive to small fractions of the injected distribution, including the particles on the verge of being reflected, which slow down to nearly zero velocity.Reaching a self-consistent equilibrium may take many crossing times of these slowed particles.To avoid unconverged measurements in this near-SCL regime, simulations were not run within the window j * ∈ [0.9, 1.1].This choice is supported by the observation that all simulations outside this window yield the same F sim measurement to within 0.1 with half the simulation time (T = 25D/v D ), while test simulations within this window change F sim by more than 0.1.This j * window corresponding to uncertainty in F model could be shrunk by increasing the simulation run time, though further study would be needed to characterize the time evolution.
After all simulations were run, we fit the resultant data F sim (j * , ϕ * , v * ) to the following model, which is essentially a sigmoid function in j * , decreasing monotonically from 1 to 0 with a varying center j c (ϕ * , v * ) and two widths j w1 (ϕ * , v * ) and j w2 (ϕ * , v * ), each dominating in a different j * limit (j w1 for large j * and j w2 for small j * ): where j c (ϕ * , v * ), j w1 (ϕ * , v * ), and j w2 (ϕ * , v * ) are given by Eq. 15 is defined such that  I.
The error was minimized to 0.17.
We note that this model automatically incorporates limiting behavior known from the cold beam SCL in Eq. 3. Particularly, in the v th → 0 limit, or v * → 0, we have (j c , j w1 , j w2 ) → (1, 0, 0).This yields a step function at j * = 1, or, equivalently, J i = J norm .Due to the definition of J norm in Eq. 10, J norm goes to J cb in the cold beam limit, so Eq. 15 is consistent with Eq. 14, as desired.Additionally, F(j * → 0) = 1 exactly.
Though Eq. 15 is monotonically decreasing in J i (and j * ), we note that the expected increasing monotonic behavior in Φ D is not strictly maintained.Holding all other parameters fixed, varying Φ D within the studied parameter space yields, at worst, a single non-monotonic dip in F model of magnitude 0.086.So long as accuracy greater than ∼ 0.10 is not required, this effect can be ignored.

VI. RESULTS
Before describing the primary result, F sim , over the thousands of simulations performed, we first describe the final steady state of a few selected simulations shown in Figure 1.Each panel shows the phase-space density in (x, v)-space of electrons using the blue-to-yellow color scale and the left vertical axis; overlaid in red is the potential Φ(x), using the right axis.These cases range from minimally space charge limited (Figure 1A-B), to moderately space charge limited (Figure 1C-D), to heavily space charge limited (Figure 1E-F).
In Figure 1A, electrons are emitted from x = 0 with velocities v > 0 centered around v i .The applied potential slows emitted electrons, and the curvature of Φ(x) indicates the effect of space charge.A small population of reflected particles are indicated by the faint green/yellow region near x = 0 with v < 0. Since F sim is not equal to 1, we know that some of these particles are reflected due to space charge.
In Figure 1B, the applied potential is zero; any space charge, no matter how small, causes Φ(x) to have a minimum below Φ D = 0, and that will reflect some electrons.However, because v * is smaller than in 1A, there are not many electrons emitted with a small enough v to be reflected by the small (negative) Φ(x), and we see that a significant space charge potential forms before even 1% of the current is reflected.Emitted electrons are visibly slowed by the potential until they are either reflected back to x = 0, or pass the potential minimum near x ≈ 0.5D and are accelerated to x = D.
Figures 1C-F exhibit significantly more space charge effects and show the characteristic ">"-shaped reflection pattern of slow electrons being decelerated, stopped, and then accelerated back towards the emitting plane.These four cases all have accelerating applied potentials.We see that as we consider larger j * , space charge increases, and the minimum Φ(x) decreases.This reflects more particles and lowers F sim .
The model parameters that best fit F sim (j * , ϕ * , v * ) are given in Table I.With these parameters, Eqs. 15 and 16 yield an approximate solution F model (j * , ϕ * , v * ) for a range of 1D diode problems with ϕ * ∈ [−0.9, 4] and v * ∈ [0.022, 0.39].This corresponds to 1D SCL current problems to 0.022 < v th /v i < 0.39 and −qΦ D < 2mv 2 i .(In the parameter space considered, v i differs from v D by at most 1%.)Multiplying F model by J max (see Sec. III) yields J t .The appendix summarizes explicitly how to find J t (q, m, D, Φ D , J i , v D , v th ) and also provides a muchsimplified calculation for beam-like cases with v th /v D < 0.05.Python implementations of these recipes are provided in supplementary material.
Figures 2:A-D show sample sweeps over j * for fixed (ϕ * , v * ) pairs with typical maximum errors, ranging from 0.07 to 0.09. Figure 2:E shows an example where F sim is fit particularly well, with a maximum error of 0.05.Finally, we show in Figure 2:F a sweep with large error, corresponding to a maximum error of 0.16.(The ϕ * and v * choices for these panels correspond to those found in Figure 1 discussed below.) We note a general trend between the relationship between v * and the derivative of F model (red curves in Figure 2) with respect to j * ; a greater v * results in a shallower slope.Physically, this is because a greater v * corresponds to a larger velocity spread, which means that an even greater change to the potential is required to cause the same change in F.
As stated in Sec.
V, the pointwise error, |F sim (j * , ϕ * , v * ) − F model (j * , ϕ * , v * )|, was reduced to a maximum of 0.17 over the studied parameter space.Figure 3 shows the maximum error, max j * |F sim − F model |, of each sweep over j * as a function of ϕ * and v * .This figure also shows the covered ϕ * and v * in black dots.

VII. SUMMARY
We have developed an empirical model (Eq.15) that accurately estimates the fraction, F, of steady-state current passing through a 1D planar diode system as a function of particle and gap characteristics, with an error of less than |Fsim − Fmodel| = 0.17 for the transition from emission-limited regime (F ≈ 1) to SCL regime (F ≈ 0).This formula empirically generalizes the 1D cold beam Child-Langmuir law to the case where particles enter the diode with a finite velocity spread and applies broadly to systems within a large parameter space if the injected particle velocity distribution resembles a truncated Maxwellian (Eq.5).Future research will extend the model's applicability to wider parameter ranges for ϕ * and v * , expand the input space to include in-0.0 0.2 0.4 0.6 0.8 1.0 x/D   jected Maxwellian distributions with negative drift velocities, and investigate oscillatory time-dependent behavior above the SCL.

SUPPLEMENTARY MATERIAL
A Python implementation of the warm model described in this paper and detailed in the appendix is included as supplementary material.Additionally, the inputs and transmitted current density outputs of all 2764 simulations used in this study are also provided.

FIG. 1 :FFFFFFFIG. 2 :
FIG.1:Each subfigure (A-F) shows two quantities for a simulation with the provided dimensionless parameters: the particle phase space density (velocity and distance) from one time step at the end of the simulation (background color map) and the electric potential Φ(x) (red lines).All quantities are normalized, with distance (x) normalized by the gap distance (D), velocity normalized to the mean injection velocity (v i ), and the potential (Φ) normalized in the same manner that Φ D is normalized to get ϕ * .

TABLE I :
List of fitting parameters used in Eq. 16.
FIG.3:The maximum error in F model over the simulated j * for fixed (v * , ϕ * ).The black dots show the simulated values of (v * , ϕ * ).The worst error is max j * |F sim − F model | = 0.17.