Hybrid simulation of Alfvén wave parametric decay instability in a laboratory relevant plasma

Predicted by Hannes Alfvén Alfvén (1942) in 1942, the Alfvén wave (AW) is widely conceived to be the fundamental mode of a magnetized plasma Wilcox, Boley, and De Silva (1960); Hasegawa and Uberoi (1982); Gekelman (1999). Because of little geometrical attenuation, shear AWs are prevalent in space and astrophysical plasmas as a carrier of magnetic and flow energy over a large range of scales De Pontieu et al. (2007). In laboratory plasmas, AWs also exist in tokamaks and linear devices, and their interaction with energetic particles is crucial to the performance of fusion reactions Chen and Zonca (2016). Large-amplitude AWs play an important role in several nonlinear processes, such as energy cascade Kraichnan (1965) and particle acceleration Kletzing (1994), which are essential in plasmas mentioned above. Parametric instabilities are an important class of such nonlinear interactions, arising from parametric excitation of collective plasma modes. These instabilities could potentially contribute to coronal heating De Pontieu et al. (2007); Velli and Pruneti (1997), the observed spectrum and cross-helicity of solar wind turbulence Inhester (1990); Del Zanna, Velli, and Londrillo (2001); Yoon and Fang (2008), and damping of fast magnetosonic waves in fusion plasmas Lee et al. (1998); Oosako et al. (2009).

Three major types of parametric instabilities driven by a parallel (perpendicular wave vector $k_{\perp}=0$) AW have been found theoretically Sagdeev and Galeev (1969); Hasegawa and Chen (1976); Derby Jr (1978); Goldstein (1978); Wong and Goldstein (1986); Longtin and Sonnerup (1986); Hollweg, Esser, and Jayanti (1993); Hollweg (1994): modulational, beat, and decay instabilities. The modulational instability drives forward upper and lower Alfvénic sidebands, as well as a non-resonant acoustic mode at the sideband separation frequency. To allow the forward AWs to interact, the pump wave must be dispersive, i.e., requiring the finite frequency effect through inclusion of the Hall term Hollweg (1994). The beat instability drives a forward upper Alfvénic sideband and a backward lower Alfvénic sideband and is generally negligible in the low-beta plasmas to be studied here. In contrast, the parametric decay instability (PDI), involving the decay of a forward pump AW into a backward daughter AW and a forward ion acoustic wave, is best known for its prominence in low-beta plasmas.

The PDI is of special interest to the space plasma community. Firstly, the backward AW generated during the PDI can interact with the forward AW, which is an important ingredient for developing magnetohydrodynamics (MHD) turbulence Van Ballegooijen et al. (2011). Secondly, the PDI-produced ion acoustic wave can strongly heat ions through wave damping Fu et al. (2020). Recently, Fu et al. Fu et al. (2018) showed that the PDI and associated kinetic heating can be robust even in a turbulent environment, which is often the case in space plasmas. The PDI has a low-beta favored growth rate, $\gamma_{g}\propto\beta^{-1/4}$, according to a linear analysis Sagdeev and Galeev (1969). A vast region of the solar surface extending from the chromosphere to corona is dictated by low-beta plasmas Gary (2001), and several studies have suggested that the PDI play important roles in a stratified chromosphere Matsumoto and Suzuki (2014), in an expanding solar wind Shoda, Yokoyama, and Suzuki (2018); Réville, Tenerani, and Velli (2018), and in the solar wind near 1 AU Shi et al. (2017); Bowen et al. (2018), contributing to turbulence development and plasma heating.

The PDI of a shear AW was first predicted by Sagdeev and Galeev Sagdeev and Galeev (1969) in 1969 and analyzed in the low-amplitude low-beta limit for a linearly polarized AW. Derby Derby Jr (1978) and Goldstein Goldstein (1978) later independently extended the analyses to finite wave amplitudes for circular polarization in the single-fluid MHD framework. The dispersive effects arising from ion cyclotron resonance were further considered with a two-fluid framework Wong and Goldstein (1986); Hollweg (1994). Despite the firm theoretical bases, observational evidence of PDI in space has been scarce. The satellite observations in the upstream of the bow shock of Earth’s magnetosphere may have found a number of cases with possible AW decay signatures Spangler, Leckband, and Cairns (1997); Narita et al. (2007); Dorfman et al. (2017). A statistical study Bowen et al. (2018) has shown correlations of density or magnetic fluctuations with PDI over a broad parameter space. However, unambiguous observation of individual PDI events in space requires overcoming several challenges Narita et al. (2007), such as the turbulent environment, limited sampling locations, broad pump bandwidths, and various kinetic effects Spangler, Leckband, and Cairns (1997). Therefore, there is a need for laboratory experiments to validate theoretical predictions. Particularly, the active control of laboratory wave and plasma conditions may help to elucidate the underlying physical processes in a manner not possible with space measurements.

The long wavelengths of AWs makes it difficult to fit the waves in a laboratory device while maintaining a low enough plasma density to allow for the use of physical probes throughout the plasma region. Recently, with the Large Plasma Device (LAPD) at University of California, Los Angeles (UCLA), several controlled AW experiments have been conducted Gekelman et al. (1997); Gekelman (1999); Maggs and Morales (2003); Vincena and Gekelman (2006); Auerbach et al. (2010, 2011); Howes et al. (2012); Dorfman and Carter (2013, 2016). Particularly, Dorfman and Carter Dorfman and Carter (2013) excited an ion acoustic mode through beating of two counter-propagating AWs. While this setup is similar to PDI excitation in the seeded mode (i.e., seeding the PDI with an artificial daughter AW), the instability growth was not clearly identified. More recently, Dorfman and Carter Dorfman and Carter (2016) also recorded a perpendicular modulational instability of kinetic AWs; when a single finite-frequency, finite-$k_{\perp}$, AW was launched above a certain wave amplitude threshold, three daughter waves were detected: two sideband AWs co-propagating with the pump wave and a low frequency non-resonant mode at the sideband separation frequency. However, the theoretical growth rate for zero-$k_{\perp}$ modulational instability was too small to explain the observation. More critically, the PDI was missing despite having a much larger theoretical growth rate than the zero-$k_{\perp}$ modulational instability under the conditions investigated. So far, no quantitative evidence of PDI has been found in laboratory experiments. In planning future experiments toward demonstration of this important instability, it should be of considerable interest to find optimal parameters with first-principle numerical simulations.

In this paper, we report the development of novel hybrid simulation capabilities aiming at direct modeling of laboratory AW dynamics. In particular, we simulate the PDI of nonlinear AWs with physical setups and wave-plasma parameters that closely resemble LAPD conditions. As an important step, we consider LAPD-like injection of a plane AW in a bounded (along the background magnetic field direction) plasma with ion kinetics retained (i.e., collisionless Landau damping). These simulations allow us to gain insight into the threshold AW amplitudes and frequencies needed for triggering PDI in a collisionless bounded plasma, and the results can be interpreted with simple theoretical estimates. We also briefly discuss complicating effects such as ion-neutral collisions and finite AW source sizes not included in the current model. It is hoped that the present work will further developments in hybrid simulation of realistic laboratory plasmas and provide guidance for future demonstration of PDI in LAPD.

We start with introducing our hybrid simulation model. The PDI has been extensively studied via MHD simulations Del Zanna, Velli, and Londrillo (2001); Shi et al. (2017), hybrid simulations Terasawa et al. (1986); Vasquez (1995); Araneda, Marsch, and Vinas (2007); Matteini et al. (2010); Fu et al. (2018); González et al. (2020), and even full particle-in-cell simulations Nariyuki, Matsukiyo, and Hada (2008). In the context of Earth’s foreshock and solar wind, the PDI has also been invoked in simulations to explain, for example, the generation of density fluctuations Tanaka, Ogino, and Umeda (2007), the origin of low-frequency Alfvénic spectrum Réville, Tenerani, and Velli (2018), the effects of wind acceleration or expansion Shoda, Yokoyama, and Suzuki (2018), and the evolution of magnetic-field-line switchbacks Tenerani et al. (2020). However, the above-mentioned simulations have mostly considered a periodic system Shi et al. (2017); Nariyuki, Matsukiyo, and Hada (2008); Tenerani et al. (2020); González et al. (2020), limiting their applicability to laboratory plasmas where the interactions are finite in space-time and nonperiodic. The periodic boundary constraint was lifted in some MHD simulations Del Zanna, Velli, and Londrillo (2001), but they did not capture kinetic effects which are important in laboratory PDI. To address that, we have recently developed novel hybrid simulation capabilities toward direct modeling of laboratory AW dynamics.

The new developments are based on the three-dimensional hybrid code H3D  Karimabadi et al. (2006), which treats electrons as a massless fluid and ions as individual macroparticles. Therefore, the simulation captures ion kinetic effects and also saves computation by ignoring electron dynamics. The code was modified to simulate turbulent plasmas with periodic boundary conditions (Fu et al., 2018, 2020). Recently, we have relaxed the periodic boundary constraint and studied the PDI in a large system ($\sim 100\lambda_{0}$ with $\lambda_{0}$ being the reference AW wavelength) of absorption boundaries for the AWs Li, Fu, and Dorfman (2022). Compared to usual periodic interaction, several new PDI dynamics were found with the nonperiodic interaction, including reduced energy transfer, and localized density fluctuations and ion heating, as well as the formation of a stable residual wave packet that is not affected by PDI. Though the parameters (e.g., $\delta B/B_{0}\sim 0.1,\beta\sim 0.01$) used were still away from typical LAPD conditions (e.g., $\delta B/B_{0}\sim 10^{-3},\beta\sim 10^{-4}-10^{-3}$), the study uncovered several limitations of usual periodic boundary condition in addressing AW dynamics in a laboratory setting.

In the present work, we further implement wave injection at prescribed locations, similar to the wave injection from an antenna in LAPD experiments, and investigate the resultant PDI dynamics. The basic simulation setup is illustrated in Fig. 1. A typical run involves one grid cell for the $x$ direction and four cells for the $y$ direction. The $z$ direction (i.e., the AW propagation direction) has a size of $9\lambda_{0}$, corresponding to $90d_{i}$, where $d_{i}=c/\omega_{pi}$ is the ion inertial length, $c$ is the light speed in a vacuum, and $\omega_{pi}$ is the ion plasma frequency. The cell sizes are $\Delta x=1d_{i}$, $\Delta y=0.1d_{i}$, $\Delta z=0.5d_{i}$, and the time step size is $\Delta t=0.01\Omega_{ci}^{-1}$, where $\Omega_{ci}$ is the ion cyclotron frequency and satisfies $\omega_{pi}/\Omega_{ci}=c/v_{A}=300$ with $v_{A}$ being the Alfvén speed. The ions are sampled by 1000 macroparticles per cell. The ion-to-electron temperature ratio, $T_{i}/T_{e}$, can be varied to approximate different LAPD discharge conditions. A finite resistivity of $\eta=2\times 10^{-4}$ is used to model the spatial wave damping as found in LAPD Gigliotti et al. (2009). The wave injection can be prescribed at any locations within the domain and typically involves an up-ramp of $5T_{0}$ in temporal profile followed by a long plateau (100s $T_{0}$), where $T_{0}=\lambda_{0}/v_{A}$.

Our simulations essentially model the injection of a plane-wave with $k_{\perp}=0$. In this quasi-1D setup, we do not consider the finite transverse scale and large $k_{\perp}$ as found in actual experiments Gigliotti et al. (2009); Karavaev et al. (2011). Moreover, the ion-neutral collisions present in a partially ionized LAPD plasma Dorfman and Carter (2013) are not modelled in the current code. Despite the simplifications, our simulations and associated analyses should be of practical interest as we adopt realistic wave and plasma parameters (e.g., wave amplitude, polarization, and plasma beta, temperature) and geometry (e.g., plasma size, wave injection, and boundary conditions). Implementation of these features is an important step toward fully modeling LAPD experiments. The wave absorption at boundaries is achieved using field masks as reported in our previous work Li, Fu, and Dorfman (2022). Notice that the periodic boundary is still used for ion particles. To reduce the effects of particle circulation in the simulation domain, we use the field masks as a buffer area for particles, as illustrated by the shaded areas in Fig. 1(b). For typical runs, the buffer area is as large as the central unmasked region, such that the particle circulation effect on wave-plasma interactions in the central region is negligible.

Figure 1 presents a case where strong PDI is observed. A left-hand polarized AW of normalized frequency $\omega_{0}/\Omega_{ci}=0.63$ is injected at $z=3.5\lambda_{0}$ (green dashed line) by prescribing its magnetic field perturbations as

$$\frac{\delta\vec{B}}{B_{0}}=\left[\frac{\delta B_{x}}{B_{0}}\sin(\omega_{0}t)\hat{x}+\frac{\delta B_{y}}{B_{0}}\cos(\omega_{0}t)\hat{y}\right]\chi(t),$$

(1)

where $\frac{\delta B_{x}}{B_{0}}=\frac{\delta B_{y}}{B_{0}}=\frac{\delta B}{B_{0}}$ and $\chi(t)$ is the temporal profile. The velocity field perturbations, $\delta v$, are automatically generated by the field solver. The finite-$\frac{\omega_{0}}{\Omega_{ci}}$, zero-$k_{\perp}$ planar AW is verified to follow the dispersion Gekelman et al. (1997); Stasiewicz et al. (2000) $\frac{\omega_{0}}{k_{\parallel}}=v_{A}\sqrt{1-(\frac{\omega_{0}}{\Omega_{ci}})^{2}}$ and has unequal magnetic/velocity perturbations with the ratio Hollweg (1999) $R_{bv}\equiv\frac{\delta B/B_{0}}{\delta v/v_{A}}=\frac{\omega_{0}}{k_{\parallel}v_{A}}=\sqrt{1-(\frac{\omega_{0}}{\Omega_{ci}})^{2}}$. The wave amplitude of the present case, $\delta B/B_{0}=0.01$, is higher than normally achievable with LAPD, but may be realized through, for example, tuning the background magnetic field $B_{0}$ or AW maser Maggs and Morales (2003). The wave pattern is shown in Fig. 1(a-c) at $t=200T_{0}$, when the PDI has sufficiently developed; the line plots are averaged over the $y$ axis, with which the results have little dependence because of the plane-wave injection. While the wave specified by Eq. (1) propagates along the $+z$ direction, the field solver also generates a mirror wave of the same polarization and amplitude that propagates in the $-z$ direction. These waves start to be damped when they reach the mask regions on both sides. As such, the effective plasma length that the forward wave travels through is $L=2.5\lambda_{0}=25d_{i}$ as marked in Fig. 1(b); this length is typical for LAPD experiments(e.g. Dorfman and Carter, 2016).

One signature of PDI is the excitation of an ion acoustic wave (e.g., displayed as density fluctuations) and associated ion heating, as seen in Fig. 1(d). To better illustrate the whole dynamics, Figure 1(f) shows the space-time evolution of $\delta\rho/\rho_{0}$, where the density fluctuations emerge near the wave injection point in just 10s $T_{0}$ and then gradually propagate outward. The density modulation deepens with time, an indication of instability growth. A deep narrow density trough is gradually formed at $z=3.5\lambda_{0}$ due to continuous wave injection and localized plasma heating [e.g., the orange curve in Fig. 1(d)]. The density fluctuations propagates forward in the region $z>3.5\lambda_{0}$, consistent with the acoustic wave direction of PDI driven by a forward pump wave. At later times ($t\gtrapprox 300T_{0}$), the acoustic wave even forms nonlinear structures propagating at a reduced speed, due to ion trapping (not shown). In PDI, the large density fluctuations further beat with the pump wave, generating a backward daughter AW; the latter overlaps with the pump wave, resulting in the rippling modulation as found in Fig. 1(e).

To gain more insight into the daughter AW, we probe the AW signatures versus time at a fixed location in front of the injection point and further separate the forward/backward waves using the modified Elsasser variables

$$Z^{\pm}=\frac{1}{2v_{A}}(\delta v\cdot R_{bv}\mp\frac{\delta B}{\sqrt{\mu_{0}\rho}}),$$

(2)

where $Z^{+}$ and $Z^{-}$ refer to the forward wave and the backward wave, respectively. For the present case, $R_{bv}\simeq 0.78$ and the difference between the magnetic and velocity field components can be clearly seen by comparing Figs. 1(b) and (c). The separated waves and corresponding Fourier spectra are presented in Figs. 1(g) and (h). Strong bursts of backward daughter AW (red color) due to the PDI are seen at later times ($t\gtrapprox 100T_{0}$). Its spectrum shows a prominent peak around $0.96\omega_{0}$, consistent with the anticipated central frequency of a lower sideband at $\omega_{-}/\omega_{0}=1-2\sqrt{\beta}\simeq 0.96$. A broad bandwidth also occurs, inherent to the finite instability bandwidth at high wave amplitudes Li, Fu, and Dorfman (2022) and also likely enhanced by a spread in $\beta$ due to nonlinear modification of the plasma density and temperature.

The above example shows, with a simplified plane-wave setup, what the PDI may look like if it can indeed be excited in LAPD-like plasmas. Next, using the same setup, we extend the study to a large parameter space. We group the simulations by two parameters, $\beta$ and $\delta B/B_{0}$, which are most relevant in characterizing the PDI. The results are summarized in Fig. 2(a), where typical LAPD conditions correspond to the bottom-left corner (outlined by the dashed square), i.e., $\beta\sim 10^{-4}-10^{-3}$ and $\delta B/B_{0}\sim 10^{-3}-10^{-2}$. Three ion-to-electron temperature ratios, $T_{i}/T_{e}=0.1,0.25,1$, are considered and marked in the plot as different colors. Each simulation is represented by a dot in the figure, where open dots refer to simulations with no PDI and filled dots indicate those with PDI. The criteria for PDI onset are based on both density and AW signatures, i.e., whether the peak amplitude of the lower Alfvénic sideband is greater than 5% of the main peak at the fundamental frequency, and the averaged density fluctuation $\langle\delta\rho/\rho_{0}\rangle$ is greater than 5%. We have also checked that, for those open-dot cases, no PDI is found even if the AW is injected over $1000T_{0}$.

It is seen from Fig. 2(a) that the AW amplitude must exceed a threshold to trigger PDI in the system and the threshold amplitudes increase with $\beta$ and $T_{i}/T_{e}$. In general, parametric instabilities occur when the pump amplitude exceeds a threshold determined by damping effects, and both the AW waves and acoustic waves may experience damping. In LAPD experiments, AW spatial damping is due to electron-ion collisions and/or AW Landau damping; both may be varied by scanning parameters such as the electron temperature Thuecks et al. (2009). In the present simulations, the spatial damping of pump AWs due to electron-ion collisions is characterized by the finite $\eta$. While the resultant AW spatial damping is found to cause slightly weaker PDI development (compared to cases without spatial damping, i.e., $\eta\to 0$), its effect is negligible. Meanwhile, the acoustic modes in LAPD may be damped by ion-neutral collisions or collisionless Landau damping; only the latter is modelled in the current code. Therefore, the threshold AW amplitudes in the present investigation should be determined by the balance between the PDI growth and Landau damping.

The PDI growth rate may be approximated in the low-beta, low-amplitude limit Sagdeev and Galeev (1969) as

$$\gamma_{g}/\omega_{0}\simeq 0.5(\delta B/B_{0})\beta^{-1/4}.$$

(3)

It should be emphasized that this normalized growth rate was obtained in the low-frequency single-fluid limit ($\omega_{0}/\Omega_{ci}\to 0$), but it is found to approximate our finite-frequency ($\omega_{0}/\Omega_{ci}\sim 0.6$) cases well; we have checked against two-fluid growth rates Hollweg (1994) and found that noticeable differences start to appear only when $\delta B/B_{0}>10^{-2}$ for typical LAPD beta values. On the other hand, we take the Landau damping rate, $\gamma_{d}/\omega\simeq\sqrt{T_{i}/T_{e}}$, for the Cauchy ion distribution Gurnett and Bhattacharjee (2005), where $\omega$ is the acoustic wave frequency. Using the Cauchy distribution greatly simplifies the calculation of Landau damping and the Cauchy distribution reasonably approximates the actual Gaussian distribution in the ion velocity range of $|v_{i}/v_{th,i}|<1.5$, where $v_{th,i}=\sqrt{k_{B}T_{i}/M_{i}}$ is the ion thermal speed with $M_{i}$ being the ion mass and $k_{B}$ the Boltzmann constant (though the Cauchy distribution has a heavier tail at larger $v_{i}$). These simplifications are justified as we aim at providing a physical interpretation for the threshold-amplitude behavior, instead of an exact numerical prediction; the latter would require considering effects neglected by the simplifying assumptions made in Sec. II. In the low-beta regime of PDI, the acoustic wave frequency satisfies $\omega/\omega_{0}\simeq 2\sqrt{\beta}$. Cast in units of $\omega_{0}$, the Landau damping rate takes the form

$$\gamma_{d}/\omega_{0}\simeq 2\sqrt{\beta}\sqrt{T_{i}/T_{e}}.$$

(4)

In order to trigger PDI, one needs to satisfy $\gamma_{g}>\gamma_{d}$, which is translated into a threshold for the AW amplitude in a collisionless plasma

$$\delta B/B_{0}>4\beta^{3/4}\sqrt{T_{i}/T_{e}}.$$

(5)

This formula is appended as color lines in Fig. 2(a) for different $T_{i}/T_{e}$, which show reasonable agreement with the simulation results, especially in the LAPD parameter regime. The agreement confirms that the threshold is determined by the competition between PDI growth and Landau damping in our simulations. Importantly, these curves cross the bottom-left corner, an indication that the excitation of PDI might be within reach at LAPD if effects neglected by our model (see discussion in Sec. V) do not significantly modify the balance between growth and damping.

A closer comparison between Eq. (5) and the simulation data in Fig. 2 shows a few prominent differences (notice the log-scale axes). Firstly, the simulations show lower threshold amplitudes, especially when $T_{i}/T_{e}$ is small. Physically, $T_{e}$ determines the phase speed of acoustic wave and $T_{i}$ determines the ion thermal speeds. When $T_{i}/T_{e}\ll 1$, the phase speed is much larger than the thermal speed, and the Cauchy distribution used for Eq. (4) tends to overestimate the damping rate because of a heavier tail in the velocity distribution. Therefore, the theoretical curves tend to overestimate the required AW amplitude when $T_{i}/T_{e}\ll 1$. Secondly, the simulations also show threshold amplitudes below the theoretical curves in the higher-beta, higher-amplitude regime, i.e., the upper-right corner of Fig. 2(a). In this regime, despite that Eq. (3) tends to slightly overestimate the growth rate, the much stronger acoustic mode damping causes significant plasma heating which may greatly reduce Landau damping rate. To see this clearly, we show the plasma heating results of two example cases in Figs. 2(b) and (c), which are marked in Fig. 2(a) as orange/green crosses, respectively. The orange-cross case ($\beta=5\times 10^{-4},\delta B/B_{0}=1\times 10^{-2},T_{i}/T_{e}=0.25$) sits above threshold curve but falls in the low-beta, low-amplitude regime. In this case, the ion phase phase-space [Fig. 2(b1)] and ion parallel temperature distributions [Fig. 2(b3)] are largely regular, and the modification to the ion distribution function [Fig. 2(b2)] is negligible, meaning Eq. 4 is still valid. In contrast, the green-cross case ($\beta=3\times 10^{-3},\delta B/B_{0}=5\times 10^{-2},T_{i}/T_{e}=1$) has significantly stronger heating [Fig. 2(c3)]; see also the distortions found in the phase-space [Fig. 2(c1)] and velocity distributions [Fig. 2(c2)].

Aside from the requirement on AW amplitude, another important parameter for consideration in experiments is what wave frequencies one should use. This is particularly relevant to the LAPD plasma which is bounded and normally contains a few AW wavelengths Gigliotti et al. (2009). It was argued Pesme, Laval, and Pellat (1973) that the absolute instability develops in a bounded unstable region if the size of the system, $L$, satisfies

$$\gamma_{g}L/|V_{1}V_{2}|^{1/2}>\pi/2,$$

(6)

where $V_{1},V_{2}$ are the group velocity of the two daughter waves, respectively. Physically, this means the system must be large enough for the instability to emerge before the daughter waves convect away. In terms of the PDI investigated here, one has $V_{1}=\frac{d\omega_{0}}{dk_{\parallel}}$ and $V_{2}=c_{s}$, where $c_{s}$ is the acoustic wave speed. Making use of the AW dispersion relation $k_{\parallel}d_{i}=\frac{\omega_{0}/\Omega_{ci}}{\sqrt{1-(\omega_{0}/\Omega_{ci})^{2}}}$, one has $\frac{d\omega_{0}}{dk_{\parallel}}=\frac{\sqrt{1-(\omega_{0}/\Omega_{ci})^{2}}}{1+k_{\parallel}^{2}d_{i}^{2}}v_{A}=[1-(\omega_{0}/\Omega_{ci})^{2}]^{3/2}v_{A}$. The acoustic wave speed takes the form $c_{s}=\sqrt{\frac{\gamma_{e}k_{B}T_{e}+\gamma_{i}k_{B}T_{i}}{M_{i}}}=\sqrt{5\beta/6}v_{A}$ where $\gamma_{e}=\gamma_{i}=5/3$ are the electron/ion polytropic indices and we have used $\beta=\beta_{e}+\beta_{i}=(1+T_{e}/T_{i})\beta_{i}$ and $\beta_{i}=2(v_{th,i}/v_{A})^{2}$. Notice that the form of $c_{s}$ is independent of the temperature ratio $T_{i}/T_{e}$ for a fixed $\beta$. Substituting the expressions for $V_{1},V_{2}$ into Eq. (6), the requirement on the system size is translated to a threshold for the wave frequency

$$F(\omega_{0}/\Omega_{ci})\equiv\frac{\omega_{0}/\Omega_{ci}}{[1-(\omega_{0}/\Omega_{ci})^{2}]^{3/4}}\gtrapprox\frac{\pi(5/6)^{1/4}}{L/d_{i}}\frac{\sqrt{\beta}}{\delta B/B_{0}},$$

(7)

where we have used the identity $c/v_{A}=\omega_{pi}/\Omega_{ci}$ and $F(\omega_{0}/\Omega_{ci})$ is a monotonic increasing function of $\omega_{0}/\Omega_{ci}$. The result is independent of $T_{i}/T_{e}$ because the growth rate has no $T_{i}/T_{e}$ dependence [i.e., Eq. (3)], which is different from the requirement on wave amplitude. Using $\omega_{0}=2\pi v_{A}\sqrt{1-(\omega_{0}/\Omega_{ci})^{2}}/\lambda_{\parallel}$ (with $\lambda_{\parallel}=2\pi/k_{\parallel}$), the frequency threshold can be recast as

$$\frac{L}{\lambda_{\parallel}}\gtrapprox 0.48\frac{\sqrt{\beta}}{\delta B/B_{0}}[1-(\omega_{0}/\Omega_{ci})^{2}]^{1/4},$$

(8)

which demands sufficient number of wavelengths to be present in the plasma region.

To verify the threshold frequency requirement, we have performed a series of simulations with $\beta=5\times 10^{-4},\delta B/B_{0}=1\times 10^{-2},T_{i}/T_{e}=0.25$, the same as the case of Fig. 1, but having different $\omega_{0}/\Omega_{ci}$. The case of Fig. 1(f) corresponds to $\omega_{0}/\Omega_{ci}\sim 0.63$ and we also select two other cases of $\omega_{0}/\Omega_{ci}\sim 0.42,0.21$ and present their space-time density fluctuations in Figs. 3(a) and (b), respectively. It indeed shows that, as the pump frequency drops, the density fluctuations become weakened and eventually disappear for the case $\omega_{0}/\Omega_{ci}\sim 0.21$. Note that, for the present wave and plasma parameters, the threshold frequency predicted from Eq. (7) is $\omega_{0}/\Omega_{ci}>0.26$, consistent with the above simulations. The agreement is better illustrated by Fig. 3(d), which plots the maximum strength of acoustic wave (in terms of longitudinal $E_{z}$ field) for a series of $\omega_{0}/\Omega_{ci}$. A clear threshold behavior around $\omega_{0}/\Omega_{ci}\sim 0.3$ is seen, below which the wave strength is at the noise level and above which a sharp increase in the wave strength is found.

For the case with lowest frequency $\omega_{0}/\Omega_{ci}\sim 0.21$, we find in simulation $L/\lambda_{\parallel}\sim 0.83$, meaning the plasma length is shorter than one Alfvén wavelength. Eq. (8) demands $L/\lambda_{\parallel}>1.1$ for the parameters used, in order to see PDI. To verify Eq. (8), we keep the same low frequency $\omega_{0}/\Omega_{ci}\sim 0.21$ but make $L$ three times bigger, i.e., $L/\lambda_{\parallel}\sim 2.5$. As shown in Fig. 3(c), the density fluctuations indeed reappear in the larger system.

Using Eq. (7) and $L=25d_{i}$, the dependence of threshold frequency on the broad parameter space spanned by $\beta$ and $\delta B/B_{0}$ is shown in Fig. 3(e). The contours show the dependencies for a few constant $\omega_{0}/\Omega_{ci}$. Generally, the onset of PDI demands higher threshold frequencies towards the lower-right corner of the parameter space and vice versa. Because the threshold frequency in Eq. (7) is inversely proportional to the wave amplitude, the threshold frequency increases quickly with $\beta$ at small wave amplitudes. One should avoid approaching $\omega_{0}/\Omega_{ci}=1$ (bottom-right corner) where strong ion cyclotron resonances for a left-hand wave will be invoked. This plot would be useful for guiding the choice of wave frequency in future experiments.

We have so far considered two critical parameters, AW amplitude $\delta B/B_{0}$ and frequency $\omega_{0}/\Omega_{ci}$, and studied their dependencies on plasma beta, temperature, and size to excite PDI in a bounded plasma with LAPD-like wave injection. However, there are several other factors that need to be considered for optimizing PDI experiments.

Firstly, in presenting the threshold $\delta B/B_{0}$ versus $\beta$ (i.e., Fig. 2), the background field $B_{0}$ is involved in both $\delta B/B_{0}$ and $\beta=n_{0}k_{B}T/(B_{0}^{2}/2\mu_{0})$, rendering $B_{0}$ a new degree of freedom. How to choose $B_{0}$ is of practical importance in LAPD experiments. To clarify this point, we introduce the parameter $\xi$ which is defined as the ratio between PDI growth and Landau damping, i.e.,

$$\xi=\frac{\gamma_{g}}{\gamma_{d}}=\frac{1}{4}\frac{\delta B}{B_{0}}\frac{(B_{0}^{2}/2\mu_{0})^{3/4}}{(n_{0}KT_{e}+n_{0}KT_{i})^{3/4}}\sqrt{\frac{T_{e}}{T_{e}}}\propto\sqrt{B_{0}}.$$

(9)

The $\sqrt{B_{0}}$ dependence suggests that PDI favors higher-$B_{0}$ regime, which corresponds to the lower-beta, lower-normalized-amplitude regime.

Secondly, the ion-to-electron temperature ratio has been shown as an important parameter [Fig. 2(a)]. Recently, an upgrade was made in LAPD including the installation of a new $\rm LaB_{6}$ cathode, which increased the accessible temperature ranges for both electrons (1 to 15 eV) and ions (sub-eV to 10 eV) and may also indicate more flexibility in the temperature ratios. More importantly, the high electron temperatures will allow us to achieve higher pump AW amplitudes by minimizing both electron-ion collisional damping and AW Landau damping Thuecks et al. (2009).

Thirdly, actual LAPD experiments usually involve a finite wave scale in the transverse direction, which is determined by the antenna size to be used Gigliotti et al. (2009). The source size is typically much smaller than $\lambda_{\parallel}$, but the wave remains collimated over the relatively small plasma dimension due to insignificant geometric attenuation of AWs Karavaev et al. (2011), which may justify the use of the quasi-1D setup in our work. However, the finite $k_{\perp}\rho_{s}\sim 0.2$ (with $\rho_{s}$ being the ion sound gyroradius) arising from the source size could be an important factor that may potentially complicate the PDI physics Dorfman and Carter (2016). A detailed theoretical understanding of the finite $k_{\perp}$ effect on PDI has not yet been developed.

Finally, the ion-neutral collisions present in the LAPD may enhance the damping of the acoustic wave, raising the threshold curves in Fig. 2(a), i.e., requiring higher AW amplitudes. We are currently implementing 3D wave injection of finite source sizes and ion-neutral collisions in the H3D code. Studies of PDI with these additional features will be presented in a separate work.

In summary, we have presented novel hybrid simulation capabilities to model continuous wave injection in a bounded plasma with absorption wave boundaries and realistic wave-plasma parameters. Such simulations could be a powerful tool to investigate a variety of LAPD AW experiments such as beat wave excitation Dorfman and Carter (2013) and turbulence generation Howes et al. (2012). Here, we have focused on the PDI of nonlinear AWs, a fundamental process in a magnetized plasma. Using parameters that resemble LAPD conditions, our simulations and accompanied theoretical analyses have provided important insight into the requirements for exciting PDI in a LAPD-like plasma, including the threshold values for AW amplitude and frequency. Other factors not included in the present simplified model (plane wave injection, collisionless plasma) are also briefly discussed. The current work is a critical step toward direct numerical modeling of LAPD plasmas and it will guide future laboratory demonstration of the fundamental PDI process which could have significant implications in space and astrophysical plasmas.

This work is supported by a National Science Foundation and Department of Energy Partnership in Basic Plasma Science and Engineering program under the grant DE-SC0021237. X.F. is also supported by the Los Alamos National Laboratory/Laboratory Directed Research and Development program and DOE/Office of Fusion Energy Sciences. S.D. is also supported by the National Aeronautics and Space Administration (NASA) grant 80NSSC18K1235. F.L. and X.F. acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin and the National Energy Research Scientific Computing Center (NERSC) for providing HPC and visualization resources. NERSC is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

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