A new simple algorithm for space charge limited emission
Peter H. Stoltz, John W. Luginsland, Andrew M. Chap, David N. Smithe, John R. Cary
NNew Emission Algorithm
AIP/123-QED
A new simple algorithm for space charge limited emission
P. H. Stoltz, a) J. W. Luginsland, A. Chap, D. N. Smithe, and J. R. Cary Tech-X Corporation Confluent Sciences, LLC. Tech-X Corporation (Dated: 7 July 2020)
Many high power electronic devices operate in a regime where the current they draw islimited by the self-fields of the particles. This space-charge-limited current poses partic-ular challenges for numerical modeling where common techniques like over-emission orGauss’ Law are computationally inefficient or produce nonphysical effects. In this paperwe show an algorithm using the value of the electric field in front of the surface insteadof attempting to zero the field at the surface, making the algorithm particularly well suitedto both electromagnetic and parallel implementations of the PIC algorithm. We show howthe algorithm is self-consistent within the framework of finite difference (for both elec-trostatics and electromagnetics). We show several 1D and 2D benchmarks against boththeory and previous computational results. Finally we show application in 3D to highpower microwave generation in a 13 GHz magnetically insulated line oscillator. a) Electronic mail: [email protected] a r X i v : . [ phy s i c s . p l a s m - ph ] J u l ew Emission Algorithm I. INTRODUCTION
High power microwave (HPM) devices frequently operate at the highest possible voltage inorder to produce the maximum rf power. The maximum power typically is coincident with themaximum current, and frequently the maximum current is limited by the self fields (or spacecharge) of the emitted electrons. This situation serves to maximize the input power to the highpower electromagnetic source. This space charge limiting (SCL) current is well understood boththeoretically and experimentally, however computer models can struggle to reproduce the exactphysical limiting current. This stems primarily from the formal singularity that exists for the spacecharge limiting current with electrons emitted with zero emission velocity . One computationalapproach is to over-emit electrons and let the self-consistent electric field push the excess chargeback to the cathode. This approach has a number of disadvantages, including being computation-ally inefficient and resulting in a nonphysical virtual cathode that must be resolved by the mesh atthe point where the electrons are pushed back. Another popular approach is to use Gauss’ Law toonly emit enough charge to exactly zero out the electric field at the emission surface. While thisis an improvement on over-emission, the Gauss’ Law approach is still quite inefficient computa-tionally. For instance, this approach requires either a particle sort and potentially an additionalelectrostatic solve. Especially in electromagnetic applications, the extra sort and solve are burden-some. Furthermore, in a parallel particle-in-cell implementation, both the sort and electrostaticsolve can be less efficient. The goal of this work is to find a new approach to SCL emission that isboth physically accurate and computationally efficient.In this paper we show an algorithm using the value of the electric field in front of the surfaceinstead of attempting to zero the field at the surface. Using the well-established particle-in-cell(PIC) code VSim , we show how the algorithm is self-consistent within the framework of finitedifference (for both electrostatics and electromagnetics). We show several 1D and 2D benchmarksagainst both theory and previous computational results for the electrostatic case. We also showa benchmark 1D electromagnetic case, showing that the voltage rise time can impact transientbehavior in a manner consistent with the literature on electromagnetic effects on space-chargelimited flow. Finally we show application in 3D to high power microwave generation in a 13 GHzmagnetically insulated line oscillator. 2ew Emission Algorithm II. DERIVATION OF THE ALGORITHM
The motivation for this algorithm is to use the finite difference electric field just above a surfaceto inform the choice of emitted current. Previous techniques typically relied instead on trying tointerpolate the field to the surface, inject charge and then use Gauss’ Law to zero out the surfacefield. However, we assume that the finite difference electric field at the midpoint of the cell edges(see E x in Fig. 1) is the field resulting from space charge limited emission, and we use that fact toinform our choice of particle velocities (or equivalently, current densities).In particular, we make the assumption that near the emission surface, the system is emittingsteady-state space charge limited current. For a one-dimensional, Cartesian geometry, the steady-state limiting current in terms of the electric field and cell size is given by the Child-Langmuirlaw: J SCL = ε (cid:114) qm E / ∆ x / . (1)In the above, ∆ x is the cell size and E is the normal component of the electric field interpolatedto the position of a new particle being added as part of the emission. In VSim, this interpolationinvolves first interpolating the edge fields to the nodes (labeled E i , j , etc, in Fig. 1), and theninterpolation from the nodal fields to the particle position. By using the same interpolation schemefor evaluation of the space charge limited emission as we do for the particle push, we are buildingin a self-consistency to the scheme. We believe this consistency is one reason for the performanceof this algorithm despite the potential challenge of applying the Child-Langmuir solution to justthe first cell rather than the more traditional procedure of applying the law to the entire gap .We can rewrite Eq. 1 in several useful ways, specifically we can write the density and velocityas a function of position : ρ ( x ) = (cid:18) m ε J SCL q (cid:19) / x − / . (2) v ( x ) = (cid:18) qJ SCL m ε (cid:19) / x / . (3)This scaling of density and velocity with position will be important when we examine the con-vergence of our algorithm. While Eq. 1 paradoxically appears to depend on the cell size, theelectric field in the 1D Child-Langmuir law scales as x / , so J SCL remains constant as the cellsize decreases.For the cases shown here, we use a variable weight approach, with a fixed number of emittedcomputational particles per step, placed randomly within the cell, and the weight of the particles3ew Emission Algorithmvaried to match the current density given by Eqs. 1. We give a list of the steps we follow here:1. Pick a number of particles to emit per step (usually based on computational resources andaccuracy requirements)2. For each particle, pick a location within one cell of the emission surface3. Evaluate the electric field at the particle location, using the same interpolation as for theforce on the particle4. Take the component of electric field normal to the emission surface and use that normalelectric field and the cell size in Eq. 1 to get a J SCL
5. Multiply that J SCL by the time step and the emission area, and divide it by the number ofparticles from Step 1, to obtain the weight for this particle.For a fixed-weight approach, one would add particles with a probability scaled by the densityprofile, keeping track of the current added, until the space charge limit is reached. For our imple-mentation, we added all particles with zero initial velocity, however one could use Eq. 3 to obtainpotentially even more self-consistency. The advantages of this approach are that it requires noextrapolation of the electric field to the surface, requires no additional evaluation of Gauss’ Lawand no associated sorting of particles, and does not rely on numerical virtual cathode formation.Thus, the method is both easy to implement in parallel PIC codes and highly effective as we showbelow.
III. VERIFICATION AND VALIDATION OF 1D ELECTROSTATICS
The initial test case for our new emission algorithm is a 1D planar diode. We chose unitvoltage and gap distance for this test, and we chose computational parameters that resulted inroughly 100,000 simulation particles in the steady state. We show in Fig. 2 the current densityversus time at both the anode and cathode. The current density is normalized to the theoretical1D limit, and the time is normalized to the time for an electron to cross the gap in the absenceof space charge. For the simulations in this figure, we used 2000 cells across the diode gap.The new algorithm exhibits a transient period where current does return to the cathode, but ona time scale short compared to the vacuum gap transit time, the cathode current drops to zero4ew Emission Algorithmand remains zero, demonstrating that this new algorithm results in no numerical virtual cathodeformation. Furthermore, by three transit times, the anode current density is steady at the expectedspace charge limit. We show the scaling of the error with resolution in Fig. 3. For this plot, weheld the time step and number of particles in the simulation fixed, and we varied only the numberof cells. We found that error scaled like the number of cells to the -2/3 power. This is consistentwith the fact that the theoretical Child-Langmuir electron density scales as x − / , as shown in Eq.2, indicating that the current is only as accurate as the simulated density near the cathode. As onefurther test of the 1D electrostatic behavior, we tested how this algorithm behaves in the presenceof particles moving in the opposite direction, for example due to elastic reflection of secondaryelectrons. We wanted to test how rapidly the algorithm can adapt to this non-standard case. Wechose as a simple test to reflect elastically all electrons at the anode. In this case, the charge ateach position is roughly doubled (the particles moving toward the anode and the particles movingtoward the cathode), so one expects the steady-state cathode-to-anode current to be roughly halfthe Child-Langmuir value. In Fig. 4, we show the longitudinal phase space for this case afterroughly ten crossing times, showing the presence of the return current. In Fig. 5, we show thecurrent reaching the anode in dark blue and the current returning to the cathode in light blue. Therough estimate of 0.5 J SCL is slightly under the simulated anode current in the steady state. Thisis due to the finite transit time across the gap and back, during which a virtual cathode forms(seen in Fig. 4). This results in return current to the physical cathode. The amount of this returncurrent is the amount by which the anode current exceeds 0.5 J SCL . This test demonstrates thatthe algorithm handles not just pure vacuum cases, but also cases where charge or current (evenoppositely directed) is already present in a computational cell.
IV. VERIFICATION AND VALIDATION OF 2D ELECTROSTATICS
A second benchmark test with widely accepted results is a 2D planar diode with a finite widthemission region . This case illustrates one key difference from the 1D diode, namely the ap-pearance of enhanced current density at the edges of the emission region. The exact amount ofenhancement depends on the width of the emission region, the size of the diode gap, and thestrength of any focusing magnetic field, but is on the order of three times the 1D Child-Langmuirvalue for emission width and diode gap close to the same size with a strong focusing field. There-fore, we set up in a 2D planar diode with dimensions (following previous work ), d=1cm, w=1cm,5ew Emission Algorithml=8cm, as shown in Fig. 6. Further following previous work, we simulated V=30kV across thegap, d, and applied a static magnetic field, By = 0.5T to keep the flow essentially purely in they-direction (this field is roughly 10x larger than the field required to keep the Larmor radius under1cm, even at the largest electron velocities). For these parameters, the transit time is on the order of0.2ns. We ran our simulation for several transit times (1.6 ns) to allow for an equilibrium to form.We chose computational parameters such that there were roughly 100,000 simulation particles atsteady state. We selected electrons within 2mm vertically of the emission region and measuredthe current density, both as a function of time and space. In Fig. 7 we show the current densityat the center of the emission region compared to the 1D Child-Langmuir versus time. Becausespace charge limited flow is a steady state condition, one expects effects like the initial spike overtime scales short compared to the transit time (this effect has a physical analog in the short pulsecommunity, where researchers can generate short electron pulses with current density greater thanthe Child-Langmuir value if the pulse duration is shorter than the transit time ). After a few transittimes, the current density is steady to within 10% of the 1D Child-Langmuir value. In Fig. 8, wehave plotted the current density versus transverse position across the emission region. The factorof three enhancement in the edge is a well-established characteristic of finite width 2D emission and further establishes that this new algorithm is behaving as expected. V. VERIFICATION AND VALIDATION OF 1D ELECTROMAGNETICS
The above tests demonstrated how the algorithm behaves for systems with time-independentfields, and in those tests we used an electrostatic field solver. However, a main goal of this work isthat this algorithm also performs well for electromagnetic systems. Consequently, we repeated ourtest of a 1D diode using an electromagnetic field solver, and allowing for a transient voltage in thegap. We changed the values of voltage and gap size to better suit an electromagnetic system, sowe chose Vfinal = 50kV and d = 0.01m. We chose computational parameters of 200 cells acrossthe gap and roughly 100,000 particles in the simulation in the steady state. To demonstrate theperformance of the algorithm in a time-dependent setting, we performed two tests: one with thevoltage varied quickly compared to the transit time, and one with the voltage varied slowly. Weshow in Fig. 9 the results for the voltage (green), current density emitted from the cathode (lightblue), and current density absorbed at the anode (dark blue). In the first case (top), we varied thevoltage linearly over one transit time. In the second case (bottom), we varied the voltage more6ew Emission Algorithmgradually over two transit times. In both cases, the simulation reaches a steady state after 2.5-3transit times, and in both cases, the simulation had zero return current to the cathode, indicating nonumerical virtual cathode formation. These metrics indicate that the algorithm performs well forelectromagnetic simulations. The algorithm also captures time-dependent transients such as theinductive increase in electric field for quick voltage rise, seen in the top plot as the spike of emittedcurrent to 1.5 times the steady-state Child-Langmuir value at around 0.5 transit times. The moregradual waveform used in the second case reduces the inductive effect. This behavior is consistentwith the physics of electromagnetic space-charge limited flow . VI. APPLICATION TO THE PHYSICS OF A MILO
With the new algorithm well benchmarked, we apply it to a problem of interest in the HPMcommunity, namely the physics of a magnetically insulated line oscillator (MILO), a promisingHPM source for high frequency. We chose to simulate a MILO similar to one from the literature to allow us to compare our results with known values. In simulation, these researchers saw 45kA at 475 kV in the oscillator region, for an input power of just over 20 GW. We show a crosssection of the full 3D geometry in Fig. 10. For simplicity, we chose to omit the extractor andmodel only the oscillator region. The slow wave structure has a choke region with two discs ofslightly larger area, followed by six discs of smaller area. An additional reason for choosing thisparticular geometry is the high frequency, which results in smaller bunching, and provides an evenmore strict test of the new emission algorithm.The MILO is a 3D cylindrical structure that does not conform to the rectangular grids used in astandard finite difference scheme. VSim uses an embedded boundary scheme to more accuratelyrepresent complex shapes within a rectangular grid. As an example, we show in Fig. 11 how ametal surface that does not conform to the grid lines might be represented. We have implementedan emission algorithm that follows these same embedded boundaries. We use the same five-stepalgorithm shown in Sec. II, with the same approach of interpolating from the nodal fields forself-consistency (even if the nodes involved are within the metal). In general, with the embeddedboundaries, the nodal fields that are inside the metal (for example, the nodes labeled E i , j and E i + , j in Fig. 11) will have values potentially different from their values in the absence of the embeddedboundary . This is because adjacent cells may not contribute to those nodes as they would inthe case of no embedded boundaries (for example, one can imagine the cell below the one shown7ew Emission Algorithmwould also be within the metal and therefore would not have any current to contribute to the nodes E i , j and E i + , j ). These modifications to the nodal fields are the way the emission algorithm ismodified for embedded boundary emission.To isolate the performance of the emission algorithm, we proceed by removing the extractorfrom the MILO in the literature. Without the extractor, in order to match the power of the realdevice, we chose a longitudinal cathode emission region to produce roughly the same amount oftotal current as seen in the previous simulations. The SCL current per unit length (when correctedfor the coaxial geometry and for relativistic effects ) for this geometry and voltage is roughly500kA/m. Therefore, we choose a longitudinal emission region of 0.1m (seen in Fig. 12). Ifthe new algorithm is emitting at the SCL limit, we expect this to result in approximately 50kA.Figure 13 (top) shows that the current is approximately 55 kA, within 10% of the expected value(because the geometry is not perfectly coaxial, we expect some difference). The bottom plot inFig. 13 shows the applied voltage of 475kV and the rf cavity voltage with amplitude of nearly thesame value. This strong coupling of input voltage to cavity voltage is an indication of fully spacecharge limited flow. Finally in Fig. 14, we show the FFT of the rf cavity voltage, showing a strongpeak at 13.5 GHz, within a few percent of the previously seen value of 13 GHz (the lack of theoutput coupler may explain this small difference). Also in the FFT, one sees the second harmonicclearly present, another indication of the strong coupling between the electrons and the rf, possibleonly with space charge limited flow. VII. CONCLUSIONS
In this paper we showed a new algorithm for simulating space charge limited emission. Thealgorithm uses the value of the electric field in front of the surface instead of attempting to zerothe field at the surface. The algorithm uses the same interpolation scheme when evaluating thecurrent as for evaluating the force on the particle, leading to a natural self-consistency. We showedapplication of the algorithm to both electrostatic and electromagnetic test cases in 1D and 2Dbenchmarks against both theory and previous computational results. We also showed applicationin 3D to high power microwave generation in a 13 GHz magnetically insulated line oscillator.Finally, while we have implemented and tested this algorithm particularly for a PIC code, thesame technique could work as a boundary condition for a fluid code.8ew Emission Algorithm
ACKNOWLEDGMENTS
This work was supported by the Office of Naval Research under Contract N68335-18-C-0060.The authors thank Christine Roark for help with several details of the simulations. The authorsalso thank Keith Cartwright, Kristian Beckwith, Thomas Gardiner, and Allen Robinson for helpfuldiscussions. 9ew Emission Algorithm
REFERENCES P. Zhang, A. Valfells, L. Ang, J. Luginsland, and Y. Y. Lau, Ap. Phys. Reviews , 011304(2017). C. Nieter and J. R. Cary, “Vorpal: A versatile plasma simulation code,” J. Comput. Phys. ,448–473 (2004). S. Mahalingam, C. Nieter, J. Loverich, D. Smithe, and P. Stoltz, “Space charge limited currentscalculations in coaxial cylindrical diodes using particle-in-cell simulations,” The Open PlasmaPhysics Journal , 63–69 (2009). J. W. Luginsland, Y. Y. Lau, R. J. Umstattd, and J. J. Watrous, “Beyond the child–langmuir law:A review of recent results on multidimensional space-charge-limited flow,” Physics of Plasmas , 2371–2376 (2002), https://doi.org/10.1063/1.1459453. R. Umstattd, C. Carr, C. Frenzen, J. Luginsland, and Y. Lau, Am. J. Phys , 160 (2005). J. J. Watrous, J. W. Lugisland, and G. E. Sasser, “An improved space-charge-limited emis-sion algorithm for use in particle-in-cell codes,” Physics of Plasmas , 289–296 (2001),https://doi.org/10.1063/1.1329152. J. J. Watrous, J. W. Luginsland, and M. H. Frese, Phys. Plasmas , 4202 (2001). W. S. Koh, L. Ang, and T. Kwan, “Multidimensional short-pulse space-charge-limited flow,”Physics of Plasmas (2006), 10.1063/1.2208086. J. Luginsland, S. McGee, and Y. Lau, IEEE Trans. Plasma Sci , 901 (1998). T. Jiang, J. Zhang, J. He, Z. Li, and J. Ling, “Experimental research on ku-band mag-netically insulated transmission line oscillator,” Physics of Plasmas , 102112 (2015),https://doi.org/10.1063/1.4933351. C. S. Meierbachtol, A. D. Greenwood, J. P. Verboncoeur, and B. Shanker, “Conformal elec-tromagnetic particle in cell: A review,” IEEE Transactions on Plasma Science , 3778–3793(2015). Y. Zhang, G. Liu, Z. Yang, Q. Xing, H. Shao, R. Xiao, H. Zhong, and Y. Lin, “Simple solutionsfor relativistic generalizations of the child–langmuir law and the langmuir–blodgett law,” Physicsof Plasmas , 044511 (2009), https://doi.org/10.1063/1.3124135.10ew Emission Algorithm FIG. 1. A typical finite difference scheme, showing the electric field ( E x ) one half of a cell offset from thesurface and the interpolated nodal fields ( E i , j and others). An electrostatic scheme would include chargedensity on the cell corners. An electromagnetic scheme would include current density co-located with theedge electric field.FIG. 2. The anode (dark blue) and cathode (light blue) return current versus time in a 1D planar diodesimulated with 2000 cells. The new algorithm results in steady current at the expected theoretical spacecharge value after about two transit times after current arrives at the anode. FIG. 3. The convergence of the current to the theoretical value with number of cells. The error convergeslike the number of cells to the -2/3 power, consistent with the fact the theoretical electron density scales as2/3 power.FIG. 4. The electron phase space for the case where the anode is a perfect elastic electron reflector. Thiscase shows that the new algorithm works not just for simple cases, but also for cases where charge or currentalready resides in a cell.
FIG. 5. The current density at the anode (dark blue) and re-absorbed at the cathode (light blue) for thecase of an anode that emits elastic secondary electrons. The simple estimate of steady-state current at theanode of 0.5 J SCL is a slight under-estimate due to the finite transit time that allows a virtual cathode to form,resulting in some return current to the cathode.FIG. 6. The geometry for the finite width 2D test case.
FIG. 7. The current density at the center of the emitter versus time. The algorithm we have implementedtakes only a few transit times to find a steady state that has less than 10% variation from the expected value.FIG. 8. The current density as a function of space after several transit times for a 2D planar diode witha finite emission region of width w. The factor of three enhancement matches previous results, furtherdemonstrating that this new algorithm is working properly.
FIG. 9. The voltage (green), emitted current density (light blue) and absorbed (dark blue) current densityfor an electromagnetic model of a 1D diode. We show two voltage waveforms, a rapid linear rise (top) anda more gradual rise (bottom). The new algorithm is able to capture both time-dependent (inductive spikes)and steady-state behavior.
FIG. 10. The cross sectional geometry of a 13 GHz MILO, showing the anode, cathode, and slow wavestructure. We monitor the voltage between the anode and cathode and in one of the cavities of the slowwave structure.FIG. 11. A finite difference electromagnetic scheme with an embedded boundary. Nodal fields within themetal ( E i , j and E i + , j in this example) will have in general smaller values than nodes in the vacuum, andthis is how the space charge limited current then is modified in the cut cells. FIG. 12. The simulated electron distribution at 25ns (showing the electrons only in the upper y plane forsimplicity). The strong bunching (seen in the vertical striations) is one sign of space charge limited emissionworking well.FIG. 13. The emitted current and current at the anode for the VSim model of the 13 GHz MILO (top),and the voltage input and measured in between one set of teeth in the slow wave structure (bottom). Theamplitude of the cavity voltage being equal to the input voltage is an indication that the space charge limitedemission is working properly.