Particle-in-cell simulations of an RF emission mechanism associated with hypervelocity impact plasmas

Hypervelocity projectiles contain enough energy to form a plasma upon impact. A projectile striking a hard target produces a shockwave provided that the impact speed is greater than the sound speed of the target (typically 5–10 km/s). The shock expands into the solid target, vaporizing and ionizing the material behind it. As the crater forms, a combination of dust particles, gas, and plasma is forced into the surrounding vacuum environment. In the early phases of expansion, the plasma is near solid density, and thus strongly coupled and highly collisional. The density drops as the plasma expands into vacuum and becomes collisionless, freezing the ionization state. Upon expansion, the temperature is ∼2 eV and the plasma flow speed is on the order of the original projectile speed.¹ Figure 1 illustrates the process. Under certain circumstances, these hypervelocity impacts are observed to produce radio frequency (RF) emission. There have been many attempts to understand this RF emission, both experimentally and theoretically. There is not yet a consensus on its source.

Radio frequency (RF) emission associated with hypervelocity impacts was first observed by Bianchi et al.² in 1984. Bianchi et al. used 1 g aluminum projectiles and hard, rocky targets in a chamber at atmospheric pressure and at 0.1–6 Torr. The detected RF emission frequency was ∼300 kHz, which suggests that the source is unlikely to be coupling from electrostatic to electromagnetic modes in the plasma since the plasma frequency is on the order of gigahertz. Bianchi et. al. also argue that the source is unlikely to be a dipole-like mechanism in the plasma since the scale of the plasma is too small to produce the observed frequency. Crawford and Schultz^3,4 measured charge separation, currents, and field strengths in a series of experiments using a light gas gun at the NASA Ames Vertical gun range. The projectile speed was 5–7 km/s and mass was 0.16–0.37 g. Foschini⁵ argued that electromagnetic radiation from impacts is capable of disrupting satellite electronics⁶ observed microwave emission from rail gun impacts at projectile speeds of 2–6.7 km/s and attributed the emission to micro-cracking.⁷ Starks⁸ did not find conclusive evidence of emissions from impacts at 11 km/s.

Close et al.⁹ measured RF from high-speed impacts under biased conditions. In contrast to previous experiments seeking RF emission, these projectiles were faster (>30 km/s) and orders of magnitude less massive (1–100 fg), while the impact occurred in a vacuum chamber with sufficiently low pressure (≲10⁻⁶ Torr) to allow a collisionless expansion of plasma. Broadband emission was measured by patch antennas at 315 MHz and 916 MHz and the impact-produced plasma was detected with retarding potential analyzers. The measured emission was found to depend strongly on impact velocity but not on impact mass. There is a velocity threshold of 15–20 km/s under which no RF emission is observed, while for higher speeds (up to 66 km/s) the detection rate reaches 60% for biased targets. In impact continuum dynamics simulations, Fletcher et al.¹ found that the same threshold occurs for the transition from weakly ionized plasma to fully ionized plasma, again independent of projectile mass. This suggests that the RF emission mechanism in these faster impacts in a sufficient vacuum is associated with the plasma, which is possibly different than the emission mechanisms (e.g., micro-cracking) in the slower impacts.

In addition to these ground-based impact experiments, there have been many reported electrical signatures associated with hypervelocity impacts on spacecraft, such as in-situ radio frequency (RF) emissions. The solar wind, which is the continuous stream of plasma emanating from the sun, accelerates dust particles (∼1 nm) up to 450 km/s. Measurements from the Cassini spacecraft orbiting Saturn measured RF signatures temporally associated with impacts from these particles.¹⁰ Larger particles (∼1 μm) traveling roughly 10 km/s originating from Saturn's rings have also produced RF measurable by Cassini.¹¹ Ten meter dipole antennas that are part of the Plasma Wave Science instrument^12–14 observed the RF through broadband electric field measurements. Kelley et al.¹⁵ argued that these signals correspond to an electrostatic discharge with decaying oscillations. The S/WAVES and SECCHI instruments aboard the two STEREO spacecraft observed electrical signals which were then correlated with impacts with optical measurements.^16,17

This paper presents simulations of an RF emission mechanism due to coherent electrostatic oscillations coupling to electromagnetic waves. A model of this mechanism was first developed by Close et al.¹⁸ The idea is that there is an extremely rapid transition from a collisional to a collisionless state, allowing the electrons to separate en masse from the ions due to their higher mobility. This sets up an ambipolar field which pulls the electrons back and the ions forward. The coherent motion of the electrons is the source of the emission. In this paper, we assume that at the point of collisional-collisionless transition, the electrons gain a bulk velocity greater than the ions by a factor of $mi/me$⁠, where m_s is the mass of species s. Under this assumption, we test the rest of the model and examine whether this is a reasonable mechanism for the observed RF emission. Section II covers the physical governing equations and a summary of the methodology used in the simulations, Section III contains the initial conditions results and discussion, and finally concluding remarks are in Section IV.

In order to simulate the rapid expansion of plasma into the vacuum surrounding the target material, the coupled Maxwell-Vlasov system of equations must be used. At this point, the plasma has reached a collisionless and frozen-in state, and kinetic effects are important. The assumptions underlying continuum models used for early stages of impact are no longer valid,¹ particularly the neglect of electromagnetic forces and the assumption of the electrons acting as a perfectly neutralizing fluid. Maxwell's equations are required over Poisson's equation since the free-space electromagnetic radiation is the phenomenon of interest. The Maxwell-Vlasov system of equations is

where s = e, i is the species, q_s is the charge, m_s is the mass, E is the electric field, B is the magnetic field, v is the velocity, f_s is the velocity distribution function of species s, ϵ₀ is the permittivity of free space, and c is the speed of light. The magnetic field in Equations (1)–(5) is the self-generated magnetic field from the plasma. The charge density and current density are, respectively,

To simulate the hypervelocity impact plasma, we use a high-order Discontinuous Galerkin (DG) particle-in-cell (PIC) method. For details on the implementation of PIC, DG, and DG-PIC, we refer the reader to Birdsall and Langdon;¹⁹ Hesthaven and Warburton;^20,21 and Jacobs and Hesthaven,²² respectively. We briefly summarize the DG-PIC method here.

DG refers to the algorithm used to solve Maxwell's equations on an unstructured mesh. It is a method for spatially discretizing the Maxwell curl equations (Equations (1) and (2)), leaving a system of ordinary differential equations in time which are integrated using a low storage, strong stability preserving Runge Kutta method.²³ Gauss' Law (Equation (3)) is a global constraint that is enforced by way of an auxiliary scalar field and an accompanying hyperbolic differential equation, a technique called hyperbolic divergence cleaning.²⁴ Equation (4) is enforced at t = 0 and is maintained within acceptable error throughout the simulation without correction. We use the transverse electric (TE) formulation of Maxwell's equations in two dimensions; this yields a current density $j=jxx̂+jyŷ$⁠, an electric field $E=Exx̂+Eyŷ$⁠, and a magnetic field $B=Bzẑ$⁠.

The Vlasov equation (Equation (5)) is effectively solved via sampling the distribution function of N_p computational macro-particles that move under the standard Lorentz force $mx¨=q(E+v×B)$⁠. It is necessary to couple the Lagrangian particles and the Eulerian mesh. We use a shape function S(r, R) to represent the ith particle p at location $xp,i$⁠. Each particle has a finite radial extent to a distance R, which is chosen to be on the order of the grid spacing Δh. Lower values of R/Δh yield faster but noisier solutions. The value at mesh node n of the charge density and current density are

The shape function is the Wendland C⁶ kernel,²⁵ which is given by (for dimensions D = 2, 3 only)

where N₂ = 78/7π and N₃ = 1365/64π. In testing of which shape function comes closest to reproducing a constant charge density from uniformly distributed particles, the Wendland C⁶ outperforms all tested kernels except the “polynomial1” shape function from Jacobs and Hesthaven.²² 

Noting that the summations in Equations (8) and (9) are over all particles, but that the shape function in Equation (10) is compactly supported, only particles within a radius R of each node contribute to the source terms at that node. To find the nodes with which each particle interacts, we build a kd-tree from all nodes at t = 0, without regard to which element they belong. To determine the forces on each particle, one must know the electromagnetic field values at the particle locations. We interpolate the fields using the multivariate Legendre polynomial basis functions L_i within each element.²² Figure 2 shows a computational particle and unstructured mesh in a DG-PIC simulation.^26,27

Particles are located via a barycentric coordinate transformation

where x_i and y_i are the coordinates of vertex i of an element and b_p,j is the jth barycentric coordinate of particle p. A particle is in the element if Equations (11)–(13) are all positive. If any b_p,j is negative, the minimum one dictates the next element to check via the element-to-element connectivity information associated with the mesh. For example, if b_p,1 is the minimum barycentric coordinate (and is negative), then the next element to check is the one connected to the edge between vertices 2 and 3.

The simulated plasma expands with a bulk fluid velocity of 30 km/s, corresponding to an impact speed of approximately 30 km/s, and an initial temperature of 2.5 eV. The impact angle is 0°, i.e., normal to the target. The initial density, scale of the plasma, and mass of the projectile are functions of the normalization. An impact plasma in the very early stages of expansion has a peak density of ∼10²⁷ m⁻³ and a length scale of ∼10 μm. An impact plasma with a plasma frequency near the RF emission frequency that was measured in experiments⁹ has a peak density of ∼10¹⁶ m⁻³ and a length scale of ∼40 mm. In normalized units, the initial peak plasma frequency in the simulations is 400. There is no background magnetic field as the gryofrequency due to the Earth's magnetic field is much less than the plasma frequency. The particles would only complete a small fraction of a cyclotron orbit over the course of the simulation.

We consider two cases

where v_d is the initial drift velocity, v_t is the initial thermal velocity, and the cold case represents the impact plasma with temperature ∼2.5 eV. In both cases, the electrons initially drift with a velocity greater than the ions by a factor of $mi/me$⁠, where m_s is the particle mass of species s.

Particles are distributed according to the radial distribution function shown in Figure 3 and uniformly in θ. This distribution function is meant to imitate results from hydrocode simulations of impact plasmas.¹ The impact point is at (x, y) = (0, 0) with the y = 0 plane as the impact surface. The expanding plume subtends an angle of π/8 and initially expands normal to the impact surface.

These simulations contain 1.15 million particles in a domain comprised of 32 258 fourth order triangles. Note that the particles are focused in about 15% of the elements, resulting in up to ∼400 particles per element. The radius of each particle is $3/2A$⁠, where A is the average area of the triangles (in this particular case, each triangle is nearly the same size). There are absorptive boundary conditions on all sides for fields and particles, although in practice the particles never approach the boundary. The reason the particles are concentrated in a small part of the domain is that we also wish to study the free space radiation produced by the plasma.

The charge density and current density for the warm case are found in Figure 4. The plasma is confined to the conical area. One can see the effect of the density gradient in both the charge and current densities. For the charge density, the Debye length, and thus the average wavelength of electrostatic fluctuations, becomes smaller near the peak density at r = 0.475. The Debye length is not entirely resolved by the mesh, but no appreciable numerical heating is observed due to the high-order nature of the method.²² The current density in $ŷ$ shows the striations that produce coherent RF emission. These are still present despite the thermal velocity exceeding the initial drift velocity by a factor of 10.

Figure 5 shows the electromagnetic fields for both the cold and warm case at the same instance in time as Figure 4. Note that the effect of the effect of the plasma is seen in the conical region in the center. The free space transverse electromagnetic waves are clearly seen in all three components for the cold case. These waves still exist in the warm case but are not perfectly clear anymore. The produced waves are similar to those produced by a dipole. There is a clear preferred direction (⁠$x̂$⁠) that can most easily be seen in B_z of the cold case. It appears that adding a significant temperature destroys this directionality. Striations similar to those in j_y in Figure 4 can also be seen in E_y.

The maximum amplitude of the wave at the simulation boundary is ∼0.1 in both cases, despite the difference in field strengths within the plasma itself. We can redimensionalize this value by multiplying by a factor $e/ϵ0L02$⁠, where L₀ is the length scale mentioned previously. We find electric field values of 18 V/m for a small scale (L₀ = 10 μm) dense impact plasma down to 1.1 × 10⁻⁶ V/m for a large scale (L₀ = 40 mm) lower density impact plasma. Spatially, we measure the wave to decay as ∝d^−0.73, where d is the distance from the plasma. We expect this decay rate to be between that of cylindrical (d^−1∕2) and spherical (d⁻¹) waves, which it is. In reality, the wave likely decays with a factor greater than 0.73 (but still less than 1) since there is a relieving effect in three dimensions. Scaling the 18 V/m value out to 30 cm, we find a predicted field of 9.8 mV/m, which again is an overestimate. In experiments, the patch antenna measured an electric field of 1.9 mV/m at the same distance.⁹ 

The waves that develop in and around the plasma can be seen most easily in the $ŷ$ component of the electric field in the cold case. The power spectra of E_y are shown in Figure 6 for both k_y and k_x. For k_y, one can see plasma oscillations for frequencies ω < 400 with the peak power near the peak plasma frequency (ω_p = 400). For small k, it follows the expected dispersion relation for electrostatic plasma oscillations, $ω2=ωp2$⁠, but for large k the frequency decreases due to interactions with the shape function (Equation (10)). The line with phase velocity $vϕ≈4$ is an artifact of the divergence cleaning algorithm; the error in the divergence is effectively propagated out of the domain at this phase velocity. For k_x, the free space electromagnetic wave with dispersion relation ω² = c²k² is visible as the diagonal line, as is the ordinary electromagnetic mode in the plasma with dispersion relation $ω2=ωp2+c2k2$⁠. The dispersions shown in Figure 6 are different due to the bulk expansion being in the $ŷ$ direction, which set up electrostatic waves with a wave vector primarily in the $ŷ$ direction and couple to electromagnetic modes with a wave vector primarily in the $x̂$ direction.

To compare with experiment, we measure the electric field at two hypothetical sensor locations: beside and above the plasma; this is shown in Figure 7. Figure 8 shows the power spectral density of the fields that reach any boundary of the simulation domain. Again, the peak power is at a frequency near the peak plasma frequency of ω = 400.

The simulations presented here show that under the assumption of a small but coherent bulk separation between ions and electrons being created by a difference in electron and ion expansion velocities, RF emission can be produced by a hypervelocity impact plasma. This result is consistent with the theory by Close et al.,¹⁸ but produces emission at a frequency higher than that detected in experiments. For radiation coherence, it is critical that the electrons expand faster than the ions en masse rather than naturally due to their significantly higher mobility.

This assumption of bulk separation is critical and must be examined more carefully. To this end, we are currently building electrostatic simulations of expanding plasma that includes a Coulomb collision algorithm. The hope is that these electrostatic simulations will show whether the sudden transition from a collisional to a collisionless state can produce the necessary bulk separation. Furthermore, these future simulations will allow us to study other possible mechanisms that could create the separation. It is also critical to examine the effect of a biased target. In experiments, the detection rate of RF emission rises significantly for biased targets.⁹ It is possible a diamagnetic bubble forms when a background magnetic field is included; investigating this, as well as the effect of impact angle, will be the goal of future simulations.

This work was supported by Department of Energy Grant No. DE-SC0010390, Air Force Office of Scientific Research Grant No. FA9550-14-1-0290 DEF, NSF-AGS Postdoctoral Research Fellowship Award No. 1433536, and the Engineering Risk Assessment team at NASA Ames Research Center.