Nonlinear saturation of the Weibel instability

TL;DR

This paper investigates the nonlinear saturation mechanisms of the Weibel instability in plasmas, revealing the significant role of electrostatic potentials in colder electron populations, which was previously underexplored.

Contribution

It demonstrates that electrostatic potentials significantly influence Weibel instability saturation in colder plasmas, expanding understanding beyond magnetic trapping and skin depth effects.

Findings

Electrostatic potential develops in colder electron populations during WI saturation.

Electrostatic fields can create potential wells that influence saturation dynamics.

Secondary instabilities can strengthen electrostatic fields affecting WI saturation.

Abstract

The growth and saturation of magnetic fields due to the Weibel instability (WI) have important implications for laboratory and astrophysical plasmas, and this has drawn significant interest recently. Since the WI can generate a large magnetic field from no initial field, the maximum magnitudes achieved can have significant consequences for a number of applications. Hence, an understanding of the detailed dynamics driving the nonlinear saturation of the WI is important. This work considers the nonlinear saturation of the WI when counter-streaming populations of initially unmagnetized electrons are perturbed by a magnetic field oriented perpendicular to the direction of streaming. Previous works have found magnetic trapping to be important and connected electron skin depth spatial scales to the nonlinear saturation of the WI. 2 Results presented in this work are consistent with these findings for a high-temperature case. However, using a high-order continuum kinetic simulation tool, this work demonstrates that, when the electron populations are colder, a significant electrostatic potential develops that works with the magnetic field to create potential wells. The electrostatic field develops due to transverse flows induced by the WI, and in some cases is strengthened by a secondary instability. This field plays a key role in saturation of the WI for colder populations. The role of the electrostatic potential in Weibel instability saturation has not been studied in detail previously.

Tables

Table 1. Table 1: Overview of the simulations

	$u_d/c$	$v_{th}/c$	$k_0{\lambda }_D$
High-temperature beams	$\pm$0.3	0.3	0.04
Low-temperature beams	$\pm$0.3	0.031	0.04

Equations

$$\frac{\partial f}{\partial t}+\text{v}\cdot \frac{\partial f}{\partial \text{x}}+\frac{q}{m}\left(\text{E}+\text{v}\times \text{B}\right)\frac{\partial f}{\partial \text{v}}=0,$$

$$\begin{gathered} \frac{\partial }{\partial t}(\frac{1}{2}\iint m{v}^{2}f(\mathbf{x},\mathbf{v},t)d\mathbf{v}d\mathbf{x}+ \\ +\frac{1}{2}\int {ϵ}_0{E}^{2}d\mathbf{x}+\frac{1}{2}\int \frac{B^2}{{\mu }_0}d\mathbf{x})=0. \end{gathered}$$

$$B_z(x)=B_{z,0}\sin (k_{0}x),$$

$$\begin{gathered} -\mathrm{i}\omega f_{1_s}+\mathrm{i}{v}_x{k}_x{f}_{1_s}+ \\ \frac{q_s}{m_s}\left[(E_x+v_y{B}_z){\partial }_{v_x}{f}_{0_s}+(E_y-v_x{B}_z){\partial }_{v_y}{f}_{0_s}\right]=0, \end{gathered}$$

$$-\mathrm{i}{k}_x{B}_z={\mu }_{0}q\left[{\int }_{\mathcal{V}}{v}_y{f}_{1_1}{\mathrm{d}}^2\mathbf{v}+{\int }_{\mathcal{V}}{v}_y{f}_{1_2}{\mathrm{d}}^2\mathbf{v}\right]-\frac{\mathrm{i}\omega }{c^2}{E}_y.$$

$$\frac{1}{2}=\frac{{\omega }_0^2}{c^2{k}^2}\left[\zeta Z(\zeta )\left(1+\frac{u_d^2}{v_{th}^2}\right)+\frac{u_d^2}{v_{th}^2}\right]+\frac{v_{th}^2}{c^2}{\zeta }^2,$$

$$Z(\zeta )=\frac{1}{\sqrt{\pi }}{\int }_{-\infty }^{\infty }\frac{e^{-x^2}}{x-\zeta }\mathrm{d}x.$$

$$\frac{{\omega }^4}{2k^2}-\left(\frac{1}{2}+\frac{1}{k^2}\right){\omega }^2-u_d^2=0,$$

$$\frac{\omega }{{\omega }_{pe}}\approx \frac{2\pi }{82}=0.076,$$

$$\frac{{\omega }_B}{{\omega }_{pe}}=\sqrt{k\frac{q}{m}{u}_y{B}_z}\approx 0.077,$$

$$q\mathbf{u}\times \mathbf{B}=-\nabla {\varphi }_{uB}.$$

$${\epsilon }_0\frac{\partial E_x}{\partial t}=-j_x=e\left[n^{+}{u}_x^{+}+n^{-}{u}_x^{-}\right].$$

$$E_x(t)$$

$$=-\frac{e}{{\epsilon }_0}\int \left[n^{+}\int \frac{e}{m}\left(E_x+u_y^{+}{B}_z\right)dt+n^{-}\int \frac{e}{m}\left(E_x+u_y^{-}{B}_z\right)dt\right]dt$$

$$\approx -\frac{e^2{u}_{y0}^{+}}{{\epsilon }_{0}m}\int \left[\left(n^{+}-n^{-}\right)\int B_{z}dt\right]dt-\frac{e^2}{{\epsilon }_{0}m}\int \left[\left(n^{+}+n^{-}\right)\int E_{x}dt\right]dt$$

$$\approx -\frac{e^2{u}_{y0}^{+}}{{\epsilon }_{0}m\gamma }\int \left[\left(n^{+}-n^{-}\right)C{e}^{\gamma t}\right]dt-\frac{e^2}{{\epsilon }_{0}m}\int \left[\left(n^{+}+n^{-}\right)\int E_{x}dt\right]dt$$

$$|E_x(t)|$$

$$1-\frac{1}{4k^2{\lambda }_D^2}\left[Z^{\prime }\left({\zeta }_1\right)+Z^{\prime }\left({\zeta }_2\right)\right]=0,$$

$${\zeta }_{1,2}=\frac{\frac{\omega }{{\omega }_0}}{\sqrt{2}k{\lambda }_D}\pm \frac{u_d}{\sqrt{2}{v}_{th}}$$

Full text

Plasma physics; Continuum kinetic simulation; Weibel instability; Nonlinear saturation

††preprint: AIP/123-QED

I Introduction

The Weibel instability (WI)Fried (1959); Weibel (1959) has been studied as a leading mechanism for the origin and growth of magnetic fields in a number of laboratoryFox et al. (2013); Okada and Ogawa (2007); Silva et al. (2002); Califano, Pegoraro, and Bulanov (1997) and astrophysical plasmaLazar et al. (2009); Ghizzo, Sarrat, and Del Sarto (2017) applications. Note that especially in the regime when drift velocities are larger than thermal velocities, this instability is also referred to as the current filamentation instability (CFI). WI can generate a large magnetic field from no initial field and can amplify a small existing field by many orders of magnitude. Hence, the WI has generated a significant amount of interest in the laboratory and astrophysics communities in recent years and a comprehensive study of the growth and nonlinear saturation of the WI is critical to estimate the saturated magnetic field magnitudes that may be achieved. Previous works have emphasized the role of magnetic trapping in the nonlinear saturation of the WIDavidson et al. (1972) and relate saturation to when the effective electron gyroradius in the generated magnetic field becomes of the order of the electron collision-less skin depth.Califano et al. (1998) In the work presented here, using fully kinetic simulations, it is shown that for relatively cold beams, in addition to the magnetic potential, an electrostatic potential develops and plays a critical role in saturating the WI. Note that in this context, the magnetic potential does not refer to the vector potential, A, but rather to the integral of the magnetic part of the Lorentz force, $\int\left(\mathbf{v}\times\mathbf{B}\right)_{x}dx$.

The WI is studied in this work using two counter-streaming populations of electrons perturbed by a magnetic field perpendicular to the beam longitudinal direction of both beams. The thermal velocity may be varied with respect to the drift velocity to understand how the saturation of the instability changes across these regimes. Results from two regimes are presented here. The first is when the counter-streaming velocity is smaller than the thermal velocity, the system is analogous to a single population with anisotropic temperature; the second is for the case when the counter-streaming velocity is larger than the thermal velocity of two distinct counter-streaming populations. The magnetic component of the Lorentz force (referred to as the “filamentation force” in this work) causes the two populations to repel each other resulting in exponential growth of the magnetic field. For hot populations, the saturation of the magnetic field occurs due to magnetic trapping, consistent with previous work Davidson et al. (1972). The work presented here shows that the saturation of the magnetic field for cold populations occurs due to the formation of potential wells that counter the filamentation force and halt the growth of the WI. The wells are caused by a combination of the magnetic and electrostatic potentials.

II Problem Description

To understand the nonlinear physics of the WI, continuum kinetic simulations of the WI in one configuration space dimension and two velocity dimensions (1X2V) are performed. The continuum kinetic model uses the discontinuous Galerkin (DG)Cockburn and Shu (2001) scheme with serendipity basis setArnold and Awanou (2011) to directly discretize the Vlasov-Maxwell equations. The Vlasov equation is solved for the electron species,

[TABLE]

where $f$ is the distribution function, $q/m$ is the charge-to-mass ratio, $\mathbf{E}$ and $\mathbf{B}$ are the electric and magnetic fields evolved using Maxwell’s equations. Ions are considered stationary in the time scales of the interest and are used only as a non-evolving neutralizing background. Extensive benchmarks are presented in a companion numerics paperJuno et al. (2017). The base method conserves energy exactly, meaning

[TABLE]

Tests performed to evaluate the energy conservation properties of this algorithmJuno et al. (2017) show that the relative energy change is on the order of $10^{-11}$. However, the limiter, which is included to ensure positivity of the distribution function, leads to small (of order one percent) energy conservation errors.

The WI simulations are initialized using two electron streams of uniform density and temperature along the $x$-axis and uniform drift along the $y$-axis. This initial uniform, but unstable, equilibrium is disturbed with a perturbation in $B_{z}$ given by,

[TABLE]

where the $k_{0}$ is the initial perturbation wave-number and $B_{z,0}=10^{-3}$ in the dimensionless units. Note that $B_{z,max}/B_{z,0}\approx 200$. The configuration space ($x$-axis) is periodic and ranges from 0 to $2\pi/k_{0}$, therefore, the initial perturbation is exactly one sine wave. In this work, results for $k_{0}\lambda_{D}=0.04$ are presented. Simulations using higher $k$ have been performed to verify that the results shown here are consistent for the range of unstable mode numbers. (For sufficiently high mode numbers the WI is stable.)

The results presented here summarize findings that rely on the subtle interplay between the magnetic and electric fields that leads to nonlinear saturation of the WI. Hence, the ability to obtain a smooth, noise-free phase-space distribution function is critical.

III Linear theory

In order to obtain the kinetic dispersion relation, the Vlasov equation (1) is linearized

[TABLE]

where $f_{0}$ is the equilibrium distribution function and $f_{1}$ is the distribution function perturbation.

$f_{0}$ is the Maxwellian distribution function. The perturbation is combined with the linearized Ampere’s law,

[TABLE]

Elimination of the fields from the previous equation then leads to the following kinetic dispersion relation,

[TABLE]

where $\omega_{0}$ is the plasma oscillation frequency, $c$ is the speed of light, $k$ is the instability wave-number, $u_{d}$ is the drift speed of each population (assuming symmetric drift velocities $u_{d}$ with respect to zero), and $v_{th}$ is thermal speed. $\zeta=\omega/(\sqrt{2}v_{th}k)$, where $\omega=\omega_{r}+\mathrm{i}\gamma$. $Z(\zeta)$ is the plasma dispersion function defined as

[TABLE]

In the cold fluid limit, using the asymptotic expansion of $Z(\zeta)$ for large $\zeta$, the cold fluid Weibel dispersion relation is obtained

[TABLE]

where $\omega$ is normalized to $\omega_{0}$, $k$ is normalized to $\omega_{0}/c$ and $u_{d}$ is normalized to $c$. Eq. 12 in Ref [Califano, Pegoraro, and Bulanov, 1997] obtains the same cold fluid dispersion relation described by Eq. 8 for the case of two counter-streaming, but otherwise identical electron beams. The cold fluid dispersion relation predicts a larger growth rate compared to the growth rate that is obtained from the kinetic dispersion relation in Eq. 6.

Previous workCagas et al. (2017) presents some initial results of the linear growth of the WI benchmarked to kinetic theory using the same continuum kinetic framework used here. The linear growth rates obtained from the kinetic simulations are in reasonably good agreement with theory for lower wavenumbers but begin to differ from theory for higher wavenumbers. The difference between theory and simulations for high wavenumbers is likely due to plasma heating which occurs as the instability grows and phase-space mixing occurs. This heating is not accounted for in the linear dispersion relation where the temperature is assumed constant. This is in agreement with the results discussed in this work in Fig. 8. Numerical values of the growth rates for the two cases used in this work together with the simulation results are presented in Fig. 1.

The growth rates from simulation are calculated by fitting the simulation data to an exponential function (with 2 free parameters). This is done by gradually increasing the region (in time) that is used to calculate the fit. The fit with the highest coefficient of determination, $R^{2}$, is then selected and used to determine the simulation growth rate.

IV Continuum kinetic simulation results

Simulation results are presented for two cases, a high-temperature conuter-streaming population of electrons and a low temperature version. Using the problem description outlined in Sec. II, the first case uses $v_{th}=u_{d}$, which is relevant to the classical Weibel instability configuration, while the second features distinct electron streams with $v_{th}<u_{d}$.

Overview of the simulations is in the Table 1. Note that the speed of light is set to be artificially lower to overcome time-step limitations. Relativistic effects are not considered in this work.

IV.1 High-temperature beams

Simulations are performed for the high-temperature case using the parameters described in Table 1 for the hot electron beam. A magnetic field grows significantly from the initial perturbation. Figure 2(a) presents the magnetic field energy and the electric field energy as a function of time. Note the exponential growth of the magnetic field energy followed by nonlinear saturation. There is negligible growth of the electric field energy in comparison to the magnetic field energy for this case. Also note the distinct periodic behavior of the magnetic field energy in the nonlinear phase of the instability. The period of the oscillations in the nonlinear part of the instability is

[TABLE]

which compares very well to the theoretical magnetic bounce period Davidson et al. (1972)

[TABLE]

when the values of $k=0.4$, $u_{y}=0.3$, and $B_{z}=0.5$ are used. These values are consistent with the normalized parameters presented in Table 1 for the hot electron case. The high-temperature case provides excellent agreement with previous results Davidson et al. (1972) showing that magnetic trapping is the primary mechanism for saturation of this instability. Figures 2b) and c) present the magnetic field and its Fourier transform as a function of the 1-dimensional space (in the y-axis) and time (in the x-axis). Note that the magnetic field spatio-temporal profile shows that only a single-mode grows throughout, and this is the same as the initialized mode. No other modes are growing in this case.

IV.2 Low-temperature beams

Simulations are performed for the low-temperature case using the parameters described in Table 1 for the cold electron beam. The evolution of the magnetic field energy and electric field energy are presented in Fig. 3(a). There are key differences in these results compared to the high-temperature case. Firstly, there is significant growth of the electric field energy, which is of the order of the magnetic field energy at saturation. Secondly, the magnetic field energy contains a similar low frequency periodic behavior in the nonlinear phase, but there is also a higher frequency perturbation superimposed on top of it (note that the continuum kinetic simulations used do not produce statistical noise like particle-in-cell codes do). Figures 3b) and c) present the spatio-temporal profile of the magnetic field and its Fourier transform similar to the high-temperature case. An additional key difference from the high-temperature case observed from the spectrogram includes the presence of higher-k modes (compared to the single initialized mode) that occur approximately around the time of nonlinear saturation.

The filamentation force is introduced in the $x$-direction owing to the drift $u_{y}$ (note that this work distinguishes bulk velocity $u$ and local velocity $v$; the bulk velocities for populations with positive $u_{y}$ and negative $u_{y}$ are denoted as $u_{y}^{+}$ and $u_{y}^{-}$) and the magnetic field perturbation $B_{z}$. This force, $qv_{y}^{\pm}B_{z}$, results in a transverse flow, $u_{x}$, in opposite directions for each of the electron populations. This filamentation leads to an exponential growth of the magnetic field and the corresponding magnetic field energy as seen in Figs. 2 and 3, which show magnetic and electric field energies. Both energies are converted from the free kinetic energy of the electrons. It is worth noting that due to the dimensions (1X2V) and problem setup, only the $B_{z}$ component of the magnetic field grows. For electric fields, both $E_{x}$ and $E_{y}$ grow, but the energy corresponding to the $E_{y}$ component is over an order of magnitude lower than the energy corresponding to $E_{x}$.

The evolution of $E_{x}$, $B_{z}$, density, and the $u_{x}$ velocity is presented in Fig. 4 as a function of space and time. More precisely, the bulk velocities and densities for each of the populations with positive $u_{y}$ and negative $u_{y}$ are computed separately. Figure 4 presents the results for the population with $u_{y}>0$, and spatial symmetry of the populations can be used to understand the profile of the population with $u_{y}<0$. During the time of linear growth, a $u_{x}$ velocity develops (Fig. 4d) and density flows from one part of the domain to the other (Fig. 4c). However, unlike the high-temperature case, there is a significant increase in the electric field (Fig. 4b) (mechanisms for this will be discussed later in more detail) at the time of instability saturation which rapidly stops the flow and consequently saturates the growth of the instability. Without the $x$-flow, the electric field decays, and the $u_{x}$ flow is reintroduced by the filamentation force. The second saturation occurs soon after. Figure 4d also provides information about the nonlinear periodic phase. When $u_{x}^{+}>0$ (red) in the region of $0.5<x\cdot k_{0}/2\pi<1.0$ and negative in the other half, the particles flow from the low-density to the high-density region increasing the filamentation. Note that the second population behaves in the opposite manner, i.e. $n^{-}$ has a maximum where $n^{+}$ has a minimum and vice versa. Because of this, the currents of the counter-streaming populations do not cancel out and $B_{z}$ is growing. On the other hand, when $u_{x}^{+}<0$ (blue) in the region of $0.5<x\cdot k_{0}/2\pi<1.0$ and positive in the other half, the gradients of density are decreasing, currents cancel out, and magnetic field decreases as well. Comparison of Fig. 3a and Fig. 4d shows that this direction of $u_{x}$ corresponds well to regions where the magnetic field is increasing and decreasing. Particles are moving from the lower density region to the higher density region until the first saturation ($t\cdot\omega_{pe}\approx 50$) and the flow is stopped. At $t\cdot\omega_{pe}>60$, the flow is reintroduced and magnetic field is rising again. After the second saturation, the flow is reversed and magnetic field is decreasing until $t\cdot\omega_{pe}\approx 75$, which is the local minimum of the magnetic field energy, and the process repeats. This behaviour is also visible in Fig. 7.

Figure 5 provides additional insight into the saturation through electromagnetic forces and the corresponding potentials. Note that, as was mentioned previously, the potential does not refer to the magnetic vector potential, A, but rather to the integral of the forces. In other words,

[TABLE]

Figure 5a presents only the magnetic part, $\mathbf{u}\times\mathbf{B}$, of the Lorentz force, whereas Fig. 5b shows the full Lorentz force including the non-negligible electric field contribution as well. Figures 5c and d present the potentials (integrals of the forces) corresponding to Figs. 5a and b, respectively. As expected, the filamentation force potential, $\phi_{uB}$, creates a potential well over one half of the domain ($x\cdot k_{0}/2\pi<0.25$ and $x\cdot k_{0}/2\pi>0.75$; note that the domain is periodic in $x$), which is consistent with the magnetically trapped particles bouncing between magnetic field extremes (see Fig. 5c). However, for this cold-temperature case, the electric field is significant in the Lorentz force. This is seen in Fig. 5d) which describes modifications to the overall potential well located at the boundary between the two populations (maxima of the electric field). The wells in these boundary regions are narrower in comparison to the potential due to the filamentation force alone, but with comparable depth (as noted from the magnitudes). If magnetic trapping was the sole mechanism in this instability, the potential described by Fig. 5c would represent the trapping potential. However, note the presence of a $\phi_{E}$ potential in Fig. 5d due to the $E_{x}$ that develops and grows. The net result of these two potentials shows the net regions of particle trapping as a function of time. Hence, the electric field trapping plays a significant role along with the magnetic trapping.

The electric field periodically rises and decays (see either Fig. 4b or 5d) which is consistent with the high frequency oscillations in Fig. 3a.

IV.3 Origin of the electric field

The transverse flows introduced by the magnetic field coupled with the nonuniform densities of the populations ($n^{+}(x)\neq n^{-}(x)$) are the main source of the electric field. Since $B_{y}=B_{x}=0$ and $\partial B_{z}/\partial y=0$, Ampere’s law in the $x$-direction reduces to

[TABLE]

The growth of the electric field can be estimated with

[TABLE]

where $\gamma$ is the magnetic field growth rate. Using the continuity equation and that $u_{x}\propto e^{\gamma t}$, one can estimate the maximum $n^{+}-n^{-}$ to be of order $e^{\gamma t}$. The contribution of the second term on the right-hand-side of Eq. 15 is bounded and cannot be an exponential. Using these approximations, an upper limit of the electric field is obtained

[TABLE]

and $C_{1}$ and $C_{2}$ are constants of integration. This shows that the electric field should be constrained by twice the exponential growth of the magnetic field.

Fig. 6 plots the growth of the electric field energy as a function of time for the low-temperature case. Note that right before saturation (before $t\omega_{pe}\sim 50$), there is an enhancement in the growth of the electric field energy.

The enhancement of electric field energy growth right before saturation could be due to a secondary instability. In this work a two-stream-like instability is explored. The green shaded region in Fig. 6 represents the regime where the transverse velocities ($u_{x}$), which vary with time, are in a regime that is unstable to the electrostatic two-stream instability. The transverse velocities change rapidly in the $x$-direction hence, there is not a single classical two-stream growth rate that is relevant to this regime. The growth rate of the secondary instability changes rapidly with time. The dispersion relation of the classical two-stream instability is given as

[TABLE]

where

[TABLE]

and $Z^{\prime}$ is the first derivative of the plasma dispersion function. The time period over which Equation (17) shows a growing two-stream instability is denoted in green in Figs. 6 and 7.

IV.4 Phase-space and temperature evolution

The trapping of particles in the potential wells near the magnetic field peaks is also seen directly in the phase-space plots of the distribution function in Fig. 8 for the low-temperature case. In order to present 2D descriptions of the 3D (1X2V) distribution function, $f(x,v_{x},v_{y})$ is integrated in $v_{x}$ to give $\hat{f}(x,v_{y})$ and in $v_{y}$ to give $\hat{f}(x,v_{x})$. The first row of Fig. 8 shows the initial conditions, the second-row plots these quantities at the time of kinetic saturation ($t\cdot\omega_{pe0}=50$), and the third row is at the end of the simulation ($t\cdot\omega_{pe0}=150$). Particle trapping and phase space mixing are clearly seen in panel (e). The bright spots are separatrices between the trapped/passing regions. The last column shows the 1D $v_{x}$ profiles of the distribution function integrated over all $v_{y}$ and averaged over $x$ between the magnetic peaks (i.e., from $1/4$ to $3/4$ of the domain). Late in time, the distribution function has significantly broadened due to phase space mixing, consistent with the earlier description plasma heating during instability evolution. The broadening of the distribution function is also seen in the high-temperature case (not shown here).

Changing the wavelength of the initial perturbation leads (for unstable wave numbers) to the same qualitative behavior, i.e growth and saturation of electric and magnetic fields and quasi-periodic nonlinear behavior due to transverse flow polarity reversal. However, higher wavenumbers display a much shorter period. This indicates that the saturation mechanisms observed here are likely universal for the WI particularly in regimes where the plasma is relatively cold compared to the drift velocity. Furthermore, the inclusion of temperature anisotropy in the counter-streaming populations also provides consistent results for nonlinear saturation and late-time nonlinear behavior.

V Summary

The high-order continuum kinetic methods used in this work allow for noise-free interpretation of detailed plasma dynamics in the kinetic regime. Due to their high dimensionality and significant computational expense, these methods were challenging until recently. Here, a detailed description of plasma dynamics is presented leading to the nonlinear saturation of the WI with distribution functions described well into the nonlinear phase of the instability. In agreement with previous work, the results presented here show the significance of particle trapping due to the magnetic fields. The simulation results using $v_{th}=u_{d}$ confirm magnetic trapping as the sole mechanism of the instability saturation. However, this work additionally emphasizes the role of the electrostatic potential in regimes where $v_{th}<u_{d}$. In the case of cold counter-streaming plasma beams, the electric field creates potential wells comparable to the magnetic field potential which significantly modifies the overall particle trapping.

Acknowledgements.

Authors would like to thank James Juno for many fruitful discussions. Simulations were performed at the Advanced Research Computing center at Virginia Tech (http://www.arc.vt.edu). This research was supported by the Air Force Office of Scientific Research under grant number FA9550-15-1-0193. The work of Ammar Hakim was supported by the U.S. Department of Energy under Contract No. DE-AC02-09CH11466.