A new simple algorithm for space charge limited emission

High power microwave (HPM) devices frequently operate at the highest possible voltage in order to produce the maximum rf power. The maximum power typically is coincident with the maximum current, and frequently, the maximum current is limited by the self-fields (or space charge) of the emitted electrons. This situation serves to maximize the input power to the high power electromagnetic source. This space charge limiting (SCL) current is well understood both theoretically and experimentally; however, computer models can struggle to reproduce the exact physical limiting current. This stems primarily from the formal singularity that exists for the space charge limiting current with electrons emitted with zero emission velocity.¹ One computational approach is to over-emit electrons and let the self-consistent electric field push the excess charge back to the cathode. This approach has a number of disadvantages, including being computationally inefficient and resulting in a nonphysical virtual cathode that must be resolved by the mesh at the point where the electrons are pushed back. Another popular approach is to use Gauss' Law to only emit enough charge to exactly zero out the electric field at the emission surface. While this is an improvement on over-emission, the Gauss' Law approach is still quite inefficient computationally. For instance, this approach requires either a particle sort and potentially an additional electrostatic solve. Especially in electromagnetic applications, the extra sort and solve are burdensome. Furthermore, in a parallel particle-in-cell implementation, both the sort and electrostatic solve can be less efficient. The goal of this work is to find a new approach to SCL emission that is both physically accurate and computationally efficient.

In this paper, we show an algorithm using the value of the electric field in front of the surface instead of attempting to zero the field at the surface. Using the well-established particle-in-cell (PIC) code VSim,^2,3 we show how the algorithm is self-consistent within the framework of finite difference (for both electrostatics and electromagnetics). The approach we present is based on successful work in other PIC codes.^4,5 The goal of our work is to put this approach on a more firm theoretical ground through analysis and benchmarking. We show several 1D and 2D benchmarks against both theory and previous computational results for the electrostatic case. We also show a benchmark 1D electromagnetic case, showing that the voltage rise time can impact transient behavior in a manner consistent with the literature on electromagnetic effects on space charge limited flow. Finally, we show the application in 3D to high power microwave generation in a 13 GHz magnetically insulated line oscillator.

The motivation for this algorithm is to use the finite difference electric field just above a surface to inform the choice of the emitted current. Previous techniques typically relied instead on trying to interpolate the field to the surface, inject charge, and then use Gauss' Law to zero out the surface field. 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 to determine the particle weights (or equivalently, current densities).

In particular, we make the assumption that near the emission surface, the system is emitting steady-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–Langmuir law:

In the above, $Δx$ is the cell size and E is the normal component of the electric field interpolated to the position of a new particle being added as part of the emission. In VSim, this interpolation involves first interpolating the edge fields to the nodes (labeled $Ei,j$⁠, etc., in Fig. 1) and then interpolation from the nodal fields to the particle position. By using the same interpolation scheme for the evaluation of the space charge limited emission as we do for the particle push, we are building in a self-consistency to the scheme. We believe this consistency is one reason for the performance of this algorithm despite the potential challenge of applying the Child–Langmuir solution to just the 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 velocity as a function of position:⁷ 

This scaling of density and velocity with position will be important when we examine the convergence of our algorithm. While Eq. (1) paradoxically appears to depend on the cell size, the electric field in the 1D Child–Langmuir law scales as⁷ $x1/3$⁠, so J_SCL remains constant as the cell size decreases.

For the cases shown here, we use a variable weight approach, with a fixed number of emitted computational particles per step, placed randomly within the cell, and the weight of the particles varied to match the current density given by Eq. (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 and accuracy requirements).

2. For each particle, pick a location within one cell of the emission surface.

3. Evaluate the electric field (including space charge) at the particle location, using the same interpolation as for the force on the particle.

4. Take the component of electric field normal to the emission surface and use that normal electric 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 of particles 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 density profile, keeping track of the current added, until the space charge limit is reached. For our implementation, we added all particles with zero initial velocity; however, one could use Eq. (3) to obtain potentially even more self-consistency. The advantages of this approach are that it requires no extrapolation of the electric field to the surface, requires no additional evaluation of Gauss' Law and no associated sorting of particles, and does not rely on the numerical virtual cathode formation. Thus, the method is both easy to implement in parallel PIC codes and highly effective as we show below.

The initial test case for our new emission algorithm is a 1D planar diode. We chose unit voltage and gap distance for this test, and we chose computational parameters that resulted in roughly 100 000 simulation particles in the steady state. We show in Fig. 2 the current density vs time emitted and returned to the cathode, as well as the current transmitted to the anode. The current density is normalized to the theoretical 1D limit, and the time is normalized to the time for an electron to cross the gap in the absence of space charge [see Eq. (2) of Ref. 5]. 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 on a timescale short compared to the vacuum gap transit time, the cathode current drops to zero and remains zero, demonstrating that this new algorithm results in no numerical virtual cathode formation. Furthermore, by three transit times, the anode current density is steady at the expected space charge limit. We show the scaling of the error with resolution in Fig. 3. For this plot, we held the time step and number of particles in the simulation fixed, and we varied only the number of cells. We found that the error scaled like the number of cells to the −2/3 power. This is consistent with the fact that the theoretical Child–Langmuir electron density scales as $x−2/3$⁠, as shown in Eq. (2), indicating that the current is only as accurate as the simulated density near the cathode. As one further test of the 1D electrostatic behavior, we tested how this algorithm behaves in the presence of particles moving in the opposite direction, for example, due to the elastic reflection of secondary electrons. We wanted to test how rapidly the algorithm can adapt to this non-standard case. We chose as a simple test to reflect elastically all electrons at the anode. In this case, the charge at each position is roughly doubled (the particles moving toward the anode and the particles moving toward the cathode), so one expects the steady-state cathode-to-anode current to be roughly half the Child–Langmuir value. In Fig. 4, we show the longitudinal phase space for this case after roughly 10 crossing times, showing the presence of the return current. In Fig. 5, we show the current reaching the anode in the dark blue and the current returning to the cathode in light blue. The rough estimate of 0.5 J_SCL is slightly under the simulated anode current in the steady state. This is due to the finite transit time across the gap and back, during which a virtual cathode forms (seen in Fig. 4). This results in the return current to the physical cathode. The amount of this return current is the amount by which the anode current exceeds 0.5 J_SCL. This test demonstrates that the algorithm handles not just pure vacuum cases, but also cases where charge or current (even oppositely directed) is already present in a computational cell.

A second benchmark test with widely accepted results is a 2D planar diode with a finite width emission region.^8,9 This case illustrates one key difference from the 1D diode, namely, the appearance of enhanced current density at the edges of the emission region. The exact amount of enhancement depends on the width of the emission region, the size of the diode gap, and the strength of any focusing magnetic field, but is on the order of three times the 1D Child–Langmuir value for emission width and diode gap close to the same size with a strong focusing field. Therefore, we set up a 2D planar diode with dimensions (following previous work⁸) d = 1 cm, w = 1 cm, l = 8 cm, as shown in Fig. 6. Further following previous work, we simulated V = 30 kV across the gap, d, and applied a static magnetic field, by = 0.5T to keep the flow essentially purely in the y-direction (this field is roughly 10× larger than the field required to keep the Larmor radius under 1 cm, even at the largest electron velocities). For these parameters, the transit time is on the order of 0.2 ns. 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 at steady state. We selected electrons within 2 mm vertically of the emission region and measured the current density, as a function of both time and space. In Fig. 7, we show the current density at the center of the emission region compared to the 1D Child–Langmuir vs time. Because the space charge limited flow is a steady state condition, one expects effects like the initial spike over time scales short compared to the transit time (this effect has a physical analog in the short pulse community, where researchers can generate short electron pulses with current density greater than the Child–Langmuir value if the pulse duration is shorter than the transit time¹⁰). After a few transit times, the current density is steady to within 10% of the 1D Child–Langmuir value. In Fig. 8, we have plotted the current density vs transverse position across the emission region. The factor of three enhancements in the edge is a well-established characteristic of finite width 2D emission⁸ and further establishes that this new algorithm is behaving as expected.

The above tests demonstrated how the algorithm behaves for systems with time-independent fields, and in those tests, we used an electrostatic field solver. However, the main goal of this work is that this algorithm also performs well for electromagnetic systems. Consequently, we repeated our test of a 1D diode using an electromagnetic field solver, and allowing for a transient voltage in the gap. We changed the values of voltage and gap size to better suit an electromagnetic system, so we chose V_final = 50 kV and d = 0.01 m. We chose computational parameters of 200 cells across the gap and roughly 100 000 particles in the simulation in the steady state. To demonstrate the performance of the algorithm in a time-dependent setting, we performed two tests: one with the voltage varied quickly compared to the transit time and another with the voltage varied slowly. We show in Fig. 9 the results for the voltage (green), current density emitted from the cathode (light blue), and current density absorbed at the anode (dark blue). In the first case (top), we varied the voltage linearly over one transit time. In the second case (bottom), we varied the voltage more gradually over two transit times. In both cases, the simulation reaches a steady state after 2.5–3 transit times, and in both cases, the simulation had a zero return current to the cathode, indicating no numerical virtual cathode formation. These metrics indicate that the algorithm performs well for electromagnetic simulations. The algorithm also captures time-dependent transients such as the inductive increase in the electric field for quick voltage rise, seen in the top plot as the spike of emitted current to 1.5 times the steady-state Child–Langmuir value at around 0.5 transit times. The more gradual waveform used in the second case reduces the inductive effect. This behavior is consistent with the physics of the electromagnetic space charge limited flow.¹¹ 

With the new algorithm well benchmarked, we apply it to a problem of interest in the HPM community, namely, the physics of a magnetically insulated line oscillator (MILO), a promising HPM source for high frequency. We chose to simulate a MILO similar to that from the literature¹² to allow us to compare our results with the known values. In simulation, these researchers saw 45 kA at 475 kV in the oscillator region, for an input power of just over 20 GW. We show a cross section of the full 3D geometry in Fig. 10. For simplicity, we chose to omit the extractor and model only the oscillator region. The slow wave structure has a choke region with two disks of slightly larger area, followed by six disks of smaller area. An additional reason for choosing this particular geometry is the high frequency, which results in smaller bunching, and provides an even more strict test of the new emission algorithm.

The MILO is a 3D cylindrical structure that does not conform to the rectangular grids used in a standard finite difference scheme. VSim uses an embedded boundary scheme to more accurately represent complex shapes within a rectangular grid. As an example, we show in Fig. 11 how a metal surface that does not conform to the grid lines might be represented. We have implemented an emission algorithm that follows these same embedded boundaries. We use the same five-step algorithm shown in Sec. II, with the same approach of interpolating from the nodal fields for self-consistency (even if the nodes involved are within the metal). VSim has machinery in place to extract the normal components to any arbitrary cut. In general, with the embedded boundaries, the nodal fields that are inside the metal (for example, the nodes labeled $Ei,j$ and $Ei+1,j$ in Fig. 11) will have values potentially different from their values in the absence of the embedded boundary.¹³ This is because adjacent cells may not contribute to those nodes as they would in the case of no embedded boundaries (for example, one can imagine the cell below the one shown would also be within the metal and therefore would not have any current to contribute to the nodes $Ei,j$ and $Ei+1,j$⁠). These modifications to the nodal fields are the way the emission algorithm is modified for embedded boundary emission.

To isolate the performance of the emission algorithm, we proceed by removing the extractor from the MILO in the literature. Without the extractor, in order to match the power of the real device, we chose a longitudinal cathode emission region to roughly produce the same amount of the total current as seen in the previous simulations. The SCL current per unit length (when corrected for the coaxial geometry and for relativistic effects¹⁴) for this geometry and voltage is roughly 500 kA/m. Therefore, we choose a longitudinal emission region of 0.1 m (seen in Fig. 12). If the new algorithm is emitting at the SCL limit, we expect this to result in approximately 50 kA. Figure 13 (left) 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 right plot in Fig. 13 shows the applied voltage of 475 kV and the rf cavity voltage with the amplitude of nearly the same value. This strong coupling of the input voltage to cavity voltage is an indication of the fully space charge limited flow. Finally in Fig. 14, we show the FFT of the rf cavity voltage, showing a strong peak at 13.5 GHz, within a few percent of the previously seen value of 13 GHz (the lack of the output coupler may explain this small difference). Also in the FFT, one sees the second harmonic clearly present, another indication of the strong coupling between the electrons and the rf, possible only with the space charge limited flow.

In this paper, we showed a new algorithm for simulating space charge limited emission. The algorithm uses the value of the electric field in front of the surface instead of attempting to zero the field at the surface. The algorithm uses the same interpolation scheme when evaluating the current as for evaluating the force on the particle, leading to a natural self-consistency. We showed the application of the algorithm to both electrostatic and electromagnetic test cases in 1D and 2D benchmarks against both theory and previous computational results. We also showed the application in 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, the same technique could work as a boundary condition for a fluid code.

This work was supported by the Office of Naval Research under Contract N68335–18-C-0060. The authors thank Christine Roark for her help with several details of the simulations. The authors also thank Keith Cartwright, Kristian Beckwith, Thomas Gardiner, and Allen Robinson for their helpful discussions.

The data that support the findings of this study are available from the corresponding author upon reasonable request.