# Parametric decay of Alfvénic wave packets in nonperiodic low-beta plasmas

Feiyu Li New Mexico Consortium, Los Alamos 87544, USA  Xiangrong Fu New Mexico
Consortium, Los Alamos 87544, USA  Seth Dorfman Space Science Institute,
Boulder, Colorado 80301, USA  University of California Los Angeles, Los
Angeles, California 90095, USA

###### Abstract

The parametric decay of finite-size Alfvén waves in nonperiodic low-beta
plasmas is investigated using one-dimensional hybrid simulations. Compared
with the usual small periodic system, a wave packet in a large system under
the absorption boundary condition shows different decay dynamics, including
reduced energy transfer and localized density cavitation and ion heating. The
resulting Alfvén wave dynamics are influenced by several factors of the
instability including the growth rate, central wave frequency, and unstable
bandwidth. A final steady state of the wave packet may be achieved when the
instability does not have enough time to develop within the residual packet,
and the packet size shows well-defined scaling dependencies on the growth
rate, wave amplitude, and plasma beta. Under the proper conditions enhanced
secondary decay can also be excited in the form of a narrow, amplified
wavepacket. These results may help interpret laboratory and spacecraft
observations of Alfvén waves, and refine our understanding of associated
energy transport and ion heating.

Alfvén waves — instabilities — hybrid simulation

##  1 Introduction

The Alfvén waves (AWs) are of fundamental importance in magnetized plasmas.
Due to their little geometrical attenuation, the waves may efficiently carry
fluctuating magnetic energy over long distances, playing crucial roles in the
energy transport found in space plasmas (De Pontieu et al., 2007),
astrophysical plasmas (Jatenco-Pereira, 1995), and fusion plasmas (Chen &
Zonca, 2016). However, at large amplitudes, the AWs may dissipate energy via
nonlinear wave-wave and wave-particle interactions, which might be relevant to
several important problems in space plasmas such as corona heating (Heyvaerts
& Priest, 1983), solar wind acceleration (Suzuki & Inutsuka, 2005), and minor
ion heating (Fu et al., 2020). Parametric instabilities are a class of such
nonlinear interactions that have gained considerable interest for the past
half century (Sagdeev & Galeev, 1969). Dependent on the plasma beta
($\beta=c_{s}^{2}/v_{A}^{2}$ where $c_{s}$ is the sound wave speed and $v_{A}$
the Alfvén speed), several instability types have been identified, including
the decay, modulational, and beat instabilities (Hollweg, 1994). In low-beta
plasmas, the parametric decay instability (PDI) predominates over the other
two for its considerably higher growth rates. In the PDI, a forward pump AW
decays to a backward AW and a forward sound wave. The generation of strong
backward AWs in the PDI may trigger MHD (magnetohydrodynamic) turbulence via
the interaction of counter-propagating AWs (Goldreich & Sridhar, 1995). The
PDI may also strongly heat plasma through damping of the compressible sound
wave (e.g., Fu et al., 2018; González et al., 2020).

The basic physics of PDI of AWs have been theoretically studied under various
conditions. Sagdeev & Galeev (1969) first analyzed the instability for
linearly polarized AWs in the low-amplitude low-$\beta$ limit. Derby (1978)
and Goldstein (1978) extended the theory to finite wave amplitudes for
circular polarization under the single-fluid MHD framework. Dispersive effects
from ion cyclotron motion were later included under the two-fluid framework,
where the dispersion gives rise also to the modulational instability by
enabling the interaction of forward propagating AW sidebands (Sakai &
Sonnerup, 1983; Longtin & Sonnerup, 1986; Wong & Goldstein, 1986; Hollweg,
1994). The PDI process has also been studied via either MHD (e.g., Del Zanna
et al., 2001), hybrid (e.g., González et al., 2020), or particle-in-cell
simulations (e.g., Nariyuki et al., 2008), under different dimensions (Ghosh
et al., 1993; Del Zanna et al., 2001), driven by either monochromatic or
broadband waves (Malara & Velli, 1996; Shoda et al., 2018), coherent or
incoherent waves (Umeki & Terasawa, 1992), with either parallel or oblique
wave propagation (Verscharen et al., 2012; Gao et al., 2013a), single or
multiple ion species (Gao et al., 2013b), and in quiescent or turbulent
background (Shi et al., 2017; Fu et al., 2018). In the context of the Earth’s
foreshock and the solar wind, studies also show observational evidence of PDI
(e.g., Spangler et al., 1997; Narita et al., 2007; Dorfman et al., 2017; Bowen
et al., 2018), which has been simulated to illustrate, for example, the
generation of density fluctuations (Tanaka et al., 2007), the origin of low-
frequency AWs (Réville et al., 2018), the effect of wind
acceleration/expansion (Shoda et al., 2018), and the evolution of switchbacks
(Tenerani et al., 2020).

Despite these extensive studies, most of the theoretical treatment and
simulations have essentially considered a periodic infinite system. This is
achieved using the periodic boundary condition where the waves repeatedly
interact with a localized plasma (e.g., Del Zanna et al., 2001; Nariyuki et
al., 2008; Gao et al., 2013a; Shi et al., 2017; Tenerani et al., 2020;
González et al., 2020). While this setup is useful for addressing the temporal
dynamics of the instability, e.g., confirming the temporal growth rates (Shi
et al., 2017), actual AWs observed in space and laboratory are usually finite
in time and space, and the associated interactions nonperiodic (Gigliotti et
al., 2009). These discrepancies may complicate the interpretation of
laboratory or spacecraft observation of AW dynamics; for example, a
modulational-like instability was recently observed in a laboratory plasma,
while the PDI was mysteriously missing despite its significant growth rates in
the low-beta environments (Dorfman & Carter, 2016). Pruneti & Velli (1997)
first lifted the periodic-boundary constraint by considering wave injection
into an open system using 1D/2D MHD simulations. Del Zanna et al. (2001)
further extended it to 3D. However, these simulations involved only a few
Alfvén wavelengths along the wave direction, thus were limited in capturing
any large-scale dynamics pertinent to the nonperiodic interactions. Notice
that the open boundary condition has been used in some solar-wind simulations
to account for the global wind structure, where the MHD framework is adopted
to handle the vast space-time scales at the expense of ignoring kinetic
effects (e.g., Tanaka et al., 2007; Shoda et al., 2018; Réville et al., 2018).

The objective of this paper is to address how a finite-size AW packet decays
in a large-scale nonperiodic system of uniform plasmas. This basic setup has
fundamental interest and can also be realized in laboratory experiments with
well-controlled plasma and wave conditions (e.g., Gigliotti et al., 2009). To
that end, we consider a wave packet of tens of Alfvén wavelengths in size and
examine its nonperiodic interaction with a low-beta plasma over a large time
scale. This study is carried out using hybrid simulations with an absorption
boundary condition for the wave fields. The large simulation domain considered
here removes the constraints on the integer wave numbers as imposed under the
periodic system. The hybrid simulation also allows us to keep the ion kinetic
effects (e.g., wave damping, ion heating) that are absent in the MHD
simulations.

The paper is organized as follows. In Section 2, we describe the simulation
model, in particular the implementation of field masking for the absorption
boundary condition. Novel PDI dynamics from the nonperiodic interactions are
then detailed in Section 3, as well as its influences on the wave propagation.
In Section 4, we discuss the results in relation to spacecraft and laboratory
observations, and briefly mention future directions.

##  2 Simulation model

The simulations are performed with H3D, a three-dimensional hybrid code that
self-consistently advances individual ion particles while treating electrons
as a massless neutralizing fluid (Karimabadi et al., 2006). Periodic boundary
conditions were assumed by the code. In order to realize the absorption
boundary condition, we apply a masking technique (Umeda et al., 2001) to the
field equations,

$E=-\frac{u_{i}\times B}{c}-\frac{B\times(\nabla\times B)}{4\pi
q_{i}n_{i}}-\frac{\nabla P_{e}}{q_{i}n_{i}}+\frac{c}{4\pi}\eta\nabla\times B,$
(1) $B^{n+1}=\overbrace{B^{n}-c\Delta t(\nabla\times E)}^{\times F_{m}},$ (2)

where $E$ is the electric field, $B$ magnetic field, $\eta$ resistivity,
$P_{e}$ electron pressure, and $u_{i},n_{i},q_{i}$ are the ion fluid velocity,
density, and charge, respectively. The superscript $n$ represents the time
step and $\Delta t$ the step size. The mask consists of a layer on either side
of the $z$ boundaries ($z$ is along the mean magnetic field), and its profile
$F_{m}$ is defined as

$F_{m}=\begin{dcases}1-\left(r\frac{z-\delta_{m}}{\delta_{m}}\right)^{2},&z\leq\delta_{m}\\\
1,&\delta_{m}\leq z<L_{z}-\delta_{m}\\\
1-\left(r\frac{z-L_{z}+\delta_{m}}{\delta_{m}}\right)^{2},&L_{z}-\delta_{m}\leq
z\leq L_{z}\\\ \end{dcases}$

where $L_{z}$ is the domain size along $z$, $\delta_{m}$ the length of either
mask layer, and $r$ a ratio factor controlling the slope of the mask. In H3D
time subcycling is used for advancing the $B$ field with Equation (2), thus
the masking is applied only at the final subcycle. It is important that the
background magnetic field $B_{z}$ is not masked. Equation (1) is also left
unmasked as the $E$ field will be automatically damped with $B$. It is found
that the masking works significantly better with smaller $r$ (i.e., smaller
gradients of $F_{m}$) but requires thicker mask layers to fully damp the
waves, as was also observed for full Maxwellian field solvers (Umeda et al.,
2001). Using $r=0.1,\delta_{m}=5d_{i}$ we achieved 99% absorption of the field
energy upon $\sim 3d_{i}$ propagation into the mask layers, where
$d_{i}=c/\omega_{pi}$ is the ion inertial length and $\omega_{pi}$ the ion
plasma frequency. Notice that particles are still periodic making the
simulation less reliable at high wave amplitudes where strong plasma heating
may induce particle circulation in the simulation domain.

We carry out essentially 1D simulations by using one cell for the $x$
direction, four cells for $y$, and considerably more cells along the wave
direction $z$. The results have little dependence on $y$ and are therefore
averaged over this coordinate. We take $z_{\rm max}=4400d_{i}\equiv
110\lambda_{0}$ for simulations with the field masking, where
$\lambda_{0}=2\pi v_{A}/\omega_{0}=40d_{i}$ is the Alfvén wavelength. In all
simulations, the cell size is $\Delta x=\Delta y=\Delta z=1d_{i}$ and time
step $\Delta t=0.01\Omega_{ci}^{-1}=4\omega_{pi}^{-1}$ where $\Omega_{ci}$ is
the ion cyclotron frequency satisfying $\omega_{pi}/\Omega_{ci}=400$. This
implies weakly dispersive AWs with $\omega_{0}/\Omega_{ci}=0.157$. A forward
propagating (along $+z$) circularly-polarized AW packet of magnetic/velocity
field

$\displaystyle\begin{aligned} \delta\vec{B}/B_{0}&=(\delta
B_{x}/B_{0})\chi(z)\cos(k_{0}z-\omega_{0}t)\hat{x}\\\ &+(\delta
B_{y}/B_{0})\chi(z)\sin(k_{0}z-\omega_{0}t)\hat{y},\end{aligned}$ (3)
$\displaystyle\delta\vec{u}=-v_{A}\delta\vec{B}/B_{0},$

is initialized inside the domain next to the left mask layer (i.e.,
$z\geq\delta_{m}$), where $\delta B_{x}=\delta B_{y}=\delta B$ and $B_{0}=1$.
The wave $E$ field is self-consistently determined from the field equations.
The wave envelope $\chi(z)$ consists of sine-square-shaped ramps on sides and
a plateau in between, i.e.,

$\chi=\begin{dcases}\sin^{2}\left(\frac{\pi(z-\delta_{m})}{2l_{r}}\right),&0\leq
z-\delta_{m}<l_{r}\\\ 1,&0\leq z-\delta_{m}-l_{r}<l_{p}\\\
\cos^{2}\left(\frac{\pi(z-z^{\prime})}{2l_{r}}\right),&0\leq
z-z^{\prime}<l_{r}\\\ \end{dcases}$

where $z^{\prime}=\delta_{m}+l_{r}+l_{p}$, and $l_{r},l_{p}$ are the length of
the ramp and plateau segment, respectively. The full packet size is
$S_{z}=l_{p}+2l_{r}$. The plasma ions (protons) are sampled by 1000
macroparticles per cell. The initial density noise, $\delta\rho/\rho_{0}\sim
0.5\%$, is suitable for triggering the PDI over a reasonable time scale. We
have confirmed very similar results using 64 particles per cell where the PDI
develops slightly earlier due to larger initial noises. Electrons have the
same initial temperature as ions and follow the adiabatic equation of state
$T_{e}/n_{e}^{\Gamma-1}=\rm const$ where $\Gamma=5/3$.

##  3 Results

Figure 1: Space-time evolution of (a, d) density fluctuation
$\delta\rho/\rho_{0}$ and (b, e) transverse $x$-component of the wave magnetic
field $\delta B_{x}/B_{0}$. (c, f) Time evolution of forward wave energy
$E^{+}$ (black solid), backward wave energy $E^{-}$ (red dashed), and cross
helicity $\sigma$ (green dotted) where the wave energies are normalized to the
initial value of $E^{+}$. The upper panels (a-c) are obtained for the
absorption boundary condition and the bottom panels (d-f) for a periodic
system. Detailed parameters are provided in the main text and notice the
different axis ranges of the space-time plots.

Figure 1 first presents an overview of two contrast simulations, highlighting
the difference between the more realistic absorption boundary condition (upper
panels a-c) and the usual periodic system (lower panels d-f). In the former
case, a square wave packet of full size $S_{z}=40\lambda_{0}$ and small ramps
$l_{r}=1\lambda_{0}$ is launched at $5\lambda_{0}\leq z\leq 45\lambda_{0}$.
The latter periodic case has a smaller simulation domain, $z_{\rm
max}=5\lambda_{0}$, which is fully occupied by a wave of constant amplitude
initially. Both cases have the same wave amplitude $\delta B/B_{0}=0.1$ and
plasma beta $\beta=0.01$. Significant differences are found in the space-time
evolution of density fluctuation (a, d) and wave magnetic field (b, e). Under
the absorption boundary condition the wave packet undergoes a strong decay at
the early stage ($t<25T_{0}$ where $T_{0}=\lambda_{0}/v_{A}$) with the packet
head being little distorted. Correspondingly, the density fluctuation is
largely localized and constrained to $z<50\lambda_{0}$. In contrast, the
periodic case shows mainly temporal dynamics where the density fluctuation
grows or decays almost uniformly across the domain. The wave field is
dominated by either forward- or backward-going waves at a time.

To measure the strength of PDI, we separate the waves into forward- and
backward-going components using the normalized Elsasser variables
$z^{\pm}=\frac{1}{2v_{A}}(\delta u\mp\delta B/\sqrt{4\pi\rho})$, and define
the associated wave energies and cross helicity as
$E^{\pm}=\frac{1}{2}\sum|z^{\pm}|^{2},\sigma=\frac{E^{+}-E^{-}}{E^{+}+E^{-}}$,
where the summation integrates over space and, in the absorption boundary
case, includes also the absorbed waves in order to ensure energy conservation.
By comparing between Figure 1(c) and (f), it is seen that about 60% of the
pump energy ($E^{+}$) is decayed during the early stage in the absorption
boundary case, while virtually the entire $E^{+}$ is depleted for the periodic
system. This dramatic difference is also revealed by the evolution of the
cross helicity which flips exactly in the latter case, indicating the
dominance of backward waves as also seen in Figure 1(e) during
$15T_{0}<t<30T_{0}$. Despite this difference, 80% of the energy lost by the
pump is transferred to the backward wave in both cases. To interpret that, we
cast the frequency matching condition for the PDI in the quantum form
$\hbar\omega_{0}=\hbar\omega+\hbar\omega_{-}$. It illustrates the energy
apportioning between the sound wave ($\omega$) and backward lower-sideband AW
($\omega_{-}$). For $\beta=0.01$ in the weakly nonlinear regime, the mode of
the maximum growth rate corresponds to $\omega/\omega_{0}=18\%$ and
$\omega_{-}/\omega_{0}=82\%$ (Derby Jr, 1978). The consistency in energy
apportioning thus suggests that both cases contain the same PDI nature despite
of very different dynamics. The two cases also share similar growth rates of
the density wave (not shown), but their maximum density fluctuation amplitudes
differ as the two start with different pump energies.

The absorption boundary case quickly evolves towards a steady state
($\sigma=0.1$ in this case) after the initial strong decay. However, the
periodic system could undergo several flips of $\sigma$ until the wave
energies in both directions exhaust. This latter phenomenon is apparently due
to the repeated wave interaction with the local plasma forced by the periodic
boundary condition. Thus, investigations using a periodic system will
significantly overestimate the energy conversion from pump to daughter waves
in the PDI, and are therefore of limited use in unraveling the PDI dynamics of
finite-size AWs in more realistic nonperiodic systems.

Figure 2: Snapshots of (a) the forward wave envelope $|z^{+}|$ (black),
backward wave envelope $|z^{-}|$ (red), the initial wave envelope (translated
to the current time step at the Alfvén speed, dashed line), density
fluctuation $\delta\rho/\rho_{0}$ (blue), and (b) parallel electric field
$E_{z}$ (black) and parallel temperature $T_{\parallel}$ (orange) at
$t=5,20,40,60T_{0}$ (from top to bottom) for the absorption boundary case
presented in Figure 1. The magenta curves in (a) are smoothed density
fluctuations, indicating any large-scale modulations.

Let us now focus on the absorption boundary case and identify how the
instability emerges, cascades, and saturates. Figure 2(a) shows progressively
the wave separation and associated density fluctuation. Initially, the wave is
dominated by the square pump sitting close to the left boundary. The maximum
growth rate for the present case is $\gamma_{\rm max}\sim 0.13\omega_{0}$. As
a result, the PDI grows exponentially following $\exp(\gamma_{\rm max}t)$ in
the linear stage. At $t=20T_{0}$, it is seen that the instability is already
strongly developed: the back portion of the wave is much depleted, a
significant backward wave appears, and a large-amplitude density fluctuation
excited. Despite that, the pump front remains largely intact, but is gradually
tailored to have a low-amplitude precursor and a smooth down-ramp profile
towards the back. We term this entire front wave portion as the residual
packet, because it stays rather stable against further propagation as shown at
$t=40,60T_{0}$. The formation of this steady state may be qualitatively
understood as follows. In the frame comoving with the AW packet, its front
portion is always fed with unperturbed plasmas. During the transit of the
plasma through the wave, the density perturbation grows in a convective
manner, becoming stronger at later times or towards the back of the wave. The
density perturbation develops to sufficient levels once it reaches the
position corresponding to the back end of the residual packet. This leads to
the decay of the back portion of the pump wave, leaving only the residual
packet intact. Notice that the residual packet propagates at a speed smaller
than the nominal Alfvén speed $v_{A}=B_{0}/\sqrt{4\pi\rho_{0}}$ as seen by the
delay relative to the dashed envelopes. About half of the delay is found to be
attributed to the finite-frequency effect, i.e., we obtain the group velocity
$v_{g}=d\omega/dk_{\parallel}=(\omega/k_{\parallel})/(1+k_{\parallel}^{2}d_{i}^{2})\simeq
0.96v_{A}$ from the dispersion relation
$\omega/k_{\parallel}=v_{A}\sqrt{1-(\omega/\Omega_{ci})^{2}}$ where
$k_{\parallel}$ is the parallel wave vector along the background magnetic
field. The wave dispersion of the finite packet is among the other effects
that would cause the additional delay.

Another remarkable feature is the emergence of a wave spike, trailing behind
the pump (i.e., outside the assumed unperturbed envelope outlined by the
dashed line). It has a large amplitude, twice that of the pump wave. This is
the secondary decay, or inverse cascade as termed in some studies with the
periodic system (e.g., Del Zanna et al., 2001). It was generated by a
secondary PDI process where the backward AW serves as the pump. This secondary
decay mostly occurs when the initial decay was strongly developed around
$t=20T_{0}$. While secondary decays are also found in the periodic system,
they never form such isolated shapes, involving both an amplification in
amplitude and a narrow duration/size. Thus, it may imply new potentials for
secondary decays in real nonperiodic environments. The detailed mechanism of
their formation is worthy of a separate study elsewhere.

Once the final steady state of the pump is reached, the PDI does not have
enough time to develop as evidenced by the negligible density fluctuations
($\delta\rho/\rho_{0}<5\%$) excited within the modulated wave envelope. On the
other hand, the large-amplitude density waves produced during the initial
decay remain localized because of the small wave speed under the low-beta
condition ($c_{s}=\sqrt{\beta}v_{A}=0.1v_{A}$); see also the slope of the wave
trace in Figure 1(a). However, at later times, the density fluctuations
gradually evolve into large-scale cavitation as indicated by the magenta
curves in Figure 2(a). This is caused by localized ion heating. To see that,
Figure 2(b) presents how the parallel electric field (black curves) develops
with the compressible sound wave. Eventually the parallel wave field vanishes
due to Landau damping of the sound wave. As a consequence, the ions are heated
up to a parallel temperature $T_{\parallel}$ (orange) several times of the
initial value $T_{\parallel}^{0}$, whereas the perpendicular temperature
remains low throughout. The locally heated ions spread out longitudinally due
to thermal expansion, causing the density cavitation. This long-term dynamics
is also a reflection of the pressure balance, $P\propto nT$, where a higher-
temperature region tends to evolve into lower densities.

The spatio-temporal (instead of just temporal) dynamics, along with the
reduced energy transfer and localized density fluctuation/cavitation and ion
heating, constitute our central new insights into the PDI process under the
more realistic nonperiodic conditions. The spatial separation of Alfvénic
fluctuations from density fluctuations and our new estimates of the efficiency
of energy conversion and new insights on ion heating may have important
implications for laboratory or spacecraft observations. For these
applications, it may be imperative to understand better the Alfvénic
propagation, especially how the final steady state depends on the plasma and
wave conditions.

Figure 3: Snapshots of $|z^{+}|$ (black), $|z^{-}|$ (red),
$\delta\rho/\rho_{0}$ (blue), and the initial wave envelope (translated to the
current time step at the Alfvén speed, dashed line) for different cases (from
top to bottom) at $t=60T_{0}$, where $S_{z}$ is the full packet size and
$S_{z}^{\rm eff}$ the effective packet size measured at FWHM.

As a first step, we investigate the effects of the wave envelope profile and
packet size on the Alfvénic propagation. Figure 3 displays the final steady
states of several cases with different wave profiles or packet sizes (e.g.,
illustrated by the dashed lines). All of them are obtained under $\beta=0.01$
and $\delta B/B_{0}=0.1$. In the case of sine-square waves (i.e., no wave
plateau region or $l_{p}=0$) the wave amplitude refers to the peak value. Case
(a) is the baseline that we have discussed above, and it has an effective
packet size (measured at full width at half maximum, FWHM) of $S_{z}^{\rm
eff}=39\lambda_{0}$. Compared with that, case (b) is only differed by a
smaller packet size ($S_{z}^{\rm eff}=19\lambda_{0}$). We see virtually the
same residual packet size for the two cases, $S_{r}^{\rm eff}\simeq
12\lambda_{0}$. This illustrates the irrelevance of initial packet size on the
steady state so along as the initial packet is larger than its final residual
packet, because the latter corresponds to only the front portion of the wave.
However, larger initial packet sizes means more energies to be converted due
to the PDI. This is evidenced by the stronger density fluctuations and more
intense cascade spike found in case (a) as compared with (b). For an initial
packet smaller than the final state size, no PDI should develop and the wave
is stable throughout because of the short transit time. This is confirmed by
case (c) where $S_{z}^{\rm eff}=7.3\lambda_{0}<12\lambda_{0}$. To clarify that
the result is not due to the change of the initial wave envelope (i.e., from
square to sine-square), we double the packet size [case(d)] and find that the
sine-square wave becomes unstable again because its effective size $S_{z}^{\rm
eff}=14.6\lambda_{0}>12\lambda_{0}$. These additional comparisons shows the
importance of packet size (or transit time), rather than the envelope, on the
PDI excitation. The envelope may still play a role through the ponderomotive
force, but its effect is generally small here due to the large packet size (or
small envelope gradients) used. Finally, Figure 3(e) presents a periodic case
having a wave size comparable to the simulation domain of the absorption
boundary cases. Again, the instability develops uniformly across the domain,
and no final steady state or enhanced wave cascades can be identified.

Figure 4: (a-c) Snapshots of $|z^{+}|$, $|z^{-}|$, and
$\langle\delta\rho/\rho_{0}\rangle$ (smoothed over $1\lambda_{0}$ along $z$
axis) for cases of increasing $\delta B/B_{0}$ (from top to bottom) and of
increasing $\beta$ (from left to right), shown at $t=60T_{0}$. The cases in
the second row correspond to the same maximum growth rate (marked as dashed
lines in d-f) despite their different $\delta B/B_{0}$ and $\beta$. (d)
Maximum growth rates versus $\delta B/B_{0}$, solved from single-fluid MHD
theory under different plasma beta $\beta$. (e) The wave number $k$ and (f)
frequency $\omega$ of the compressible ion acoustic wave corresponding to the
maximum growth rates versus $\delta B/B_{0}$. The shaded areas indicate the
range within which the growth rates are non-negative.

With the above insights, we next focus on the simple square wave profile and
extend the investigation to a broad range of wave amplitude and plasma beta.
The simulation results are summarized in Figure 4(a-c) in terms of the final
states and their corresponding density fluctuations. Furthermore, following
the single-fluid theory (Derby Jr, 1978; Goldstein, 1978) which is justified
here due to the small pump dispersion, the parameters $\delta B/B_{0}$ and
$\beta$ adequately define the dispersion relation of the PDI,
$\omega=\omega(k)$. The latter can be solved to find the growth rate, central
frequency, and unstable bandwidth of the instability measured by the sound
wave. These calculations, central to understanding the PDI and the resulting
Alfvénic propagation, are presented in Figure 4(d-f) in terms of the maximum
growth rate $\gamma_{\rm max}$ and the bandwidth corresponding to non-negative
$\gamma$. Of the most importance is the growth rate which defines the
timescale relevant for the PDI development. In general, more significant
decays happen with increasing growth rates, corresponding to larger $\delta
B/B_{0}$ and lower $\beta$ following Figure 4(d). This is verified by the
overall pattern of the simulation results orderly arranged in Figure 4(a-c).
Closer to the top-right having smaller $\delta B/B_{0}$ and higher $\beta$,
the wave is less prone to PDI. The trend reverses towards to bottom-left
(having larger $\delta B/B_{0}$ and smaller $\beta$) where strong decays are
found even at the early stage of the interactions.

Despite its importance, the growth rate does not exclusively define the PDI-
driven Alfvénic evolution. As seen in Figure 4(e-f), the central $\omega$ and
$k$ of the sound wave (solid curves) corresponding to the maximum growth rates
also show strong dependence with $\beta$. For weakly nonlinear regimes (i.e.,
not too large $\delta B/B_{0}$), one has $\omega/\omega_{0}<1$ and
$k/k_{0}\sim 2/(1+\sqrt{\beta})$, where $\omega/\omega_{0}$ increases and
$k/k_{0}$ decreases with $\beta$, respectively. Therefore, following the
energy apportioning relation $N\hbar\omega_{-}=N\hbar\omega_{0}-N\hbar\omega$,
more energy goes into the backward AWs at smaller $\beta$ for a given pump
depletion (i.e., fixed $N$ units of energy quanta). In addition, following
$|k_{-}|=k-k_{0}$, one has
$|k_{-}/k_{0}|\sim\frac{1-\sqrt{\beta}}{1+\sqrt{\beta}}<1$ which means
$|k_{-}/k_{0}|$ keeps cascading into smaller numbers upon each secondary
decay. The more limited $k$-space at higher $\beta$ (which corresponds to
smaller initial $|k_{-}/k_{0}|$) may make secondary decays more difficult.
This phenomenon was identified in periodic systems (Del Zanna et al., 2001)
where the minimum wavenumber is limited by the box size. It may also hold for
the present nonperiodic condition when the packet size is large enough to
allow for multiple decays. From the above analysis, both the frequency and
wavenumber matching conditions suggest that the secondary wave decays favor
smaller $\beta$ where more energy goes into the backward wave having larger
$k_{-}/k_{0}$. This is supported by the simulation results where the secondary
decays are stronger from column (c) to column (a) with decreasing $\beta$; for
example, notice the prominent spike observed for $\beta=0.01,\delta
B/B_{0}=0.1$. We have checked that the wave power in the low-frequency range
(corresponding to the secondary decays) indeed increases with decreasing
$\beta$ even in cases without a clearly visible cascade spike.

A final contributing factor to the Alfvénic dynamics is the unstable mode
bandwidth, i.e., the shaded areas in Figure 4(e-f). The theory predicts
rapidly increasing bandwidths with $\delta B/B_{0}$, becoming even larger than
the central frequency/wavenumber as $\delta B/B_{0}\to 1$. The presence of
multiple modes may cause strong competition among each other and make the wave
envelopes more spiky. These features can be seen by comparing top to bottom
panels in Figure 4(a-c).

The above analyses may point out the criteria for strong secondary decays with
good wave coherency (e.g., narrow bandwidths and smooth envelopes). These
phenomena would favor low-beta plasmas where the central frequency is small
and the wavenumber is large, such that more energy goes into the backward
waves with sufficiently large wavenumbers. Besides, the wave amplitude should
not be too large, involving much mode competition and making the interactions
less coherent. It should neither be too small requiring a large timescale for
the PDI to develop.

Figure 5: (a) The effective residual packet size, $S_{r}^{\rm eff}$, versus
the maximum growth rate of the PDI under different $\beta$. All simulation
results (point markers) are obtained with the square wave of an initial size
$S_{z}=60\lambda_{0}$. The wave amplitudes $\delta B/B_{0}$ are chosen to have
the maximum growth rate fall between $0.07\omega_{0}$ and $0.22\omega_{0}$.
The solid curves represent fitting by Equation (4), where the coefficients
$C(\beta)$ are presented in (b) as a function of $\beta$ (vertical lines
indicate errorbars). The solid curve in (b) is a fitting of $C(\beta)$ using
the $\sqrt{\beta}$ dependence. (c) All simulation data presented in (a) are
rearranged here using the combined axis $\beta/(\delta B/B_{0})^{2}$. The
solid line represents Equation (5). The shaded areas in (b, c) come from the
uncertainties in fitting $C(\beta)$.

Finally, we have assumed in the above analyses a certain pump depletion, which
corresponds to the total energy lost to both the backward wave and sound wave.
It can be defined as $\mathcal{D}=(S_{z}^{\rm eff}-S_{r}^{\rm eff})/S_{z}^{\rm
eff}$. The quantity $1-D=S_{r}^{\rm eff}/S_{z}^{\rm eff}$ determines the
portion of energy that can be stably transported over large distances and
should be useful for estimating the energy transport efficiency and the limit
of PDI-related plasma heating. As identified above, several factors may
contribute to the Alfvénic dynamics, making analytical treatment of $D$ or
$S_{r}^{\rm eff}$ a challenge. Here, we resort to simulations for its
dependencies on the wave amplitude and plasma beta. We scan over a range of
parameters where the maximum growth rate falls within $0.07<\gamma_{\rm
max}/\omega_{0}<0.22$. Higher growth rates typically involve large $\delta
B/B_{0}$ and broad mode competition, making the final state less
distinguishable. Even lower growth rates may require a long time for PDI to
develop sufficiently, a challenge to the simulations. Our results are
summarized as point markers in Figure 5(a). Despite the complexity of the
interactions, a simple inverse-square dependence of the residual packet size
versus the maximum growth rate is found as

$\frac{S_{r}^{\rm eff}}{\lambda_{0}}=C(\beta)\left(\frac{\gamma_{\rm
max}}{\omega_{0}}\right)^{-2},$ (4)

where the coefficient $C(\beta)$ has a square-root dependence on $\beta$,
i.e., $C(\beta)\simeq 1.8\sqrt{\beta}$, as displayed in Figure 5(b). To more
vividly illustrate the $\beta$-dependence, the second row of Figure 4(a-c)
shows three cases of different $\delta B/B_{0}$ and $\beta$ designed to have
the same growth rate [marked as the dashed lines in Figure 4(d-f)]. Despite
this, the final evolved states of the three cases are quite different, with
significantly more pump depletion at smaller $\beta$.

In the low-amplitude low-beta limit, an approximated form of the growth rate
was suggested to be $\gamma_{\rm max}/\omega_{0}\simeq\frac{1}{2}\frac{\delta
B}{B_{0}}\frac{1}{\beta^{1/4}}$ (Sagdeev & Galeev, 1969). By substituting it
into Equation (4), we arrive at

$\frac{S_{r}^{\rm eff}}{\lambda_{0}}=8C\frac{\sqrt{\beta}}{(\delta
B/B_{0})^{2}}\simeq 14.4\frac{\beta}{(\delta B/B_{0})^{2}},$ (5)

which explicitly associates the residual packet size with the plasma beta and
wave amplitude. Notice a difference of $\sqrt{2}$ in wave amplitude when
applying Sagdeev & Galeev’s scaling who studied linearly polarized AWs. We
have verified Equation (5) by arranging all simulation data in the combined
axis $\beta/(\delta B/B_{0})^{2}$ as displayed in Figure 5(c). It is seen that
they follow the linear scaling reasonably well, especially when the
uncertainty in the coefficient of Equation (5) is considered (the shaded
area). The $\beta=0.1$ cases are slightly off, confirming the applicability of
Equation (5) to only the low-beta limit. In addition, high-amplitude cases
typically involve broad mode bandwidths and spiky modulations of the wave
envelope, making identification of the residual state more difficult; their
detailed effects on long-term AW dynamics need more refined studies in the
future.

##  4 Discussion

In this study, we have relaxed the usual periodic boundary condition and
investigated the parametric decay dynamics of finite-size AWs in a large-scale
low-beta plasma under the more realistic nonperiodic boundary condition. By
simulating the dynamics over large space and time scales, our 1D hybrid
simulations unravel novel dynamics that were not appreciated in the periodic
system before, including strongly space-time dependent evolution, reduced pump
energy transfer, and localized density fluctuation/cavitation and ion heating.
These results may have important implications for spacecraft and laboratory
observations.

The robustness and importance of the PDI have stimulated the search of their
existence in space over the past decades, for example, in the ultra-low-
frequency wave field upstream of the Earth’s bow shock (Spangler et al., 1988;
Narita et al., 2007). However, observational signature in the magnetic
spectrum has been limited due to the dominance of the pump wave and the close
separation of the pump and sidebands in frequency space (Tsurutani et al.,
1987; Agim et al., 1995). The density spectral feature, on the other hand, may
provide another independent signature of PDI because it is clearly separated
from the magnetic signature in both wavenumber and frequency in low-beta
plasmas (Spangler et al., 1997). As our simulations have shown, however,
simultaneous observation of both magnetic and density signatures may not
always be possible because the density fluctuation can be separated in space
and time from the Alfvénic fluctuation. This may explain the missing “decay
line” signature reported in Spangler et al. (1997) and help improve the search
by establishing proper detection criteria. Furthermore, the conjunction
observation from multiple spacecraft in the heliosphere (Velli et al., 2020)
may identify missing signatures for PDI by correlating magnetic and density
signatures at different locations along the same field line.

In addition, our simulations have shown that the front portion of the wave
packet may evolve into a steady state under suitable wave and plasma
parameters. This is because the instability has not enough time to grow in the
wave packet. It also indicates efficient energy transport of a short wave
train over long distances even in an environment where the PDI has a
considerable growth rate. By combining simplified theories and numerical
scans, we have identified well-behaved scaling dependencies of the residual
packet size with both the growth rate, wave amplitude, and plasma beta. These
scaling laws are valid for relatively low-beta low-amplitude regimes and
should be useful for further analytical treatment on the energy transport. The
AW propagation under the absorption boundary condition in higher-beta higher-
amplitude regimes, which may apply to certain regions of the heliosphere
(Gary, 2001), is beyond the scope of the present work.

Finally, we have elucidated the major factors and corresponding parameter
space that would influence the strength of decays and cascades. These insights
should stimulate new interest in the basic PDI dynamics, especially in their
highly nonlinear regime, for example, the modification to the spectral
evolution of inverse cascade (Chandran, 2018). However, it should be pointed
out that while our study has adopted a basic setup of fundamental importance,
it is limited by the 1D geometry. High-dimensional effects such as finite wave
$k_{\perp}$ (Dorfman & Carter, 2016), oblique wave propagation (Viñas &
Goldstein, 1991; Verscharen et al., 2012), filamentation, magnetosonic decays,
and perpendicular turbulence cascades (Del Zanna et al., 2001; González et
al., 2020) are among immediate next steps to check on the features identified
in this work. We may also consider more realistic wave injection from the
boundary for a finite period of time (Del Zanna et al., 2001). The latter
would have direct correspondence with laboratory experiments (Gigliotti et
al., 2009) and may help resolve some puzzles related to the AW parametric
instabilities (e.g., Dorfman & Carter, 2016).

This work is supported by NSF/DOE Partnership in Basic Plasma Science and
Engineering under the grant DE-SC0021237. The authors 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.

## References

  * Agim et al. (1995) Agim, Y., Vinas, A., & Goldstein, M. 1995, Journal of Geophysical Research: Space Physics, 100, 17081
  * Bowen et al. (2018) Bowen, T. A., Badman, S., Hellinger, P., & Bale, S. D. 2018, The Astrophysical Journal Letters, 854, L33
  * Chandran (2018) Chandran, B. D. 2018, Journal of plasma physics, 84
  * Chen & Zonca (2016) Chen, L., & Zonca, F. 2016, Reviews of Modern Physics, 88, 015008
  * De Pontieu et al. (2007) De Pontieu, B., McIntosh, S., Carlsson, M., et al. 2007, Science, 318, 1574
  * Del Zanna et al. (2001) Del Zanna, L., Velli, M., & Londrillo, P. 2001, Astronomy & Astrophysics, 367, 705
  * Derby Jr (1978) Derby Jr, N. 1978, The Astrophysical Journal, 224, 1013
  * Dorfman & Carter (2016) Dorfman, S., & Carter, T. 2016, Physical review letters, 116, 195002
  * Dorfman et al. (2017) Dorfman, S., Hietala, H., Astfalk, P., & Angelopoulos, V. 2017, Geophysical Research Letters, 44, 2120
  * Fu et al. (2020) Fu, X., Guo, F., Li, H., & Li, X. 2020, The Astrophysical Journal, 890, 161
  * Fu et al. (2018) Fu, X., Li, H., Guo, F., Li, X., & Roytershteyn, V. 2018, The Astrophysical Journal, 855, 139
  * Gao et al. (2013a) Gao, X., Lu, Q., Li, X., Shan, L., & Wang, S. 2013a, Physics of Plasmas, 20, 072902
  * Gao et al. (2013b) Gao, X., Lu, Q., Tao, X., Hao, Y., & Wang, S. 2013b, Physics of Plasmas, 20, 092106
  * Gary (2001) Gary, G. A. 2001, Solar Physics, 203, 71
  * Ghosh et al. (1993) Ghosh, S., Viñas, A., & Goldstein, M. 1993, Journal of Geophysical Research: Space Physics, 98, 15561
  * Gigliotti et al. (2009) Gigliotti, A., Gekelman, W., Pribyl, P., et al. 2009, Physics of Plasmas, 16, 092106
  * Goldreich & Sridhar (1995) Goldreich, P., & Sridhar, S. 1995, The Astrophysical Journal, 438, 763
  * Goldstein (1978) Goldstein, M. L. 1978, The Astrophysical Journal, 219, 700
  * González et al. (2020) González, C., Tenerani, A., Velli, M., & Hellinger, P. 2020, The Astrophysical Journal, 904, 81
  * Heyvaerts & Priest (1983) Heyvaerts, J., & Priest, E. R. 1983, Astronomy and Astrophysics, 117, 220
  * Hollweg (1994) Hollweg, J. V. 1994, Journal of Geophysical Research: Space Physics, 99, 23431
  * Jatenco-Pereira (1995) Jatenco-Pereira, V. 1995, Physica Scripta, 1995, 113
  * Karimabadi et al. (2006) Karimabadi, H., Vu, H., Krauss-Varban, D., & Omelchenko, Y. 2006, in Numerical Modeling of Space Plasma Flows, Vol. 359, 257
  * Longtin & Sonnerup (1986) Longtin, M., & Sonnerup, B. Ö. 1986, Journal of Geophysical Research: Space Physics, 91, 6816
  * Malara & Velli (1996) Malara, F., & Velli, M. 1996, Physics of Plasmas, 3, 4427
  * Narita et al. (2007) Narita, Y., Glassmeier, K. H., Fränz, M., Nariyuki, Y., & Hada, T. 2007, Nonlinear Processes in Geophysics, 14, 361, doi: 10.5194/npg-14-361-2007
  * Nariyuki et al. (2008) Nariyuki, Y., Matsukiyo, S., & Hada, T. 2008, New Journal of Physics, 10, 083004
  * Pruneti & Velli (1997) Pruneti, F., & Velli, M. 1997, in Fifth SOHO Workshop: The Corona and Solar Wind Near Minimum Activity, Vol. 404, 623
  * Réville et al. (2018) Réville, V., Tenerani, A., & Velli, M. 2018, The Astrophysical Journal, 866, 38
  * Sagdeev & Galeev (1969) Sagdeev, R. Z., & Galeev, A. A. 1969, Nonlinear Plasma Theory
  * Sakai & Sonnerup (1983) Sakai, J.-I., & Sonnerup, B. Ö. 1983, Journal of Geophysical Research: Space Physics, 88, 9069
  * Shi et al. (2017) Shi, M., Li, H., Xiao, C., & Wang, X. 2017, The Astrophysical Journal, 842, 63
  * Shoda et al. (2018) Shoda, M., Yokoyama, T., & Suzuki, T. K. 2018, The Astrophysical Journal, 860, 17
  * Spangler et al. (1988) Spangler, S., Fuselier, S., Fey, A., & Anderson, G. 1988, Journal of Geophysical Research: Space Physics, 93, 845
  * Spangler et al. (1997) Spangler, S. R., Leckband, J. A., & Cairns, I. H. 1997, Physics of Plasmas, 4, 846
  * Suzuki & Inutsuka (2005) Suzuki, T. K., & Inutsuka, S.-i. 2005, The Astrophysical Journal Letters, 632, L49
  * Tanaka et al. (2007) Tanaka, S., Ogino, T., & Umeda, T. 2007, Journal of Geophysical Research: Space Physics, 112
  * Tenerani et al. (2020) Tenerani, A., Velli, M., Matteini, L., et al. 2020, The Astrophysical Journal Supplement Series, 246, 32
  * Tsurutani et al. (1987) Tsurutani, B. T., Thorne, R. M., Smith, E. J., Gosling, J., & Matsumoto, H. 1987, Journal of Geophysical Research: Space Physics, 92, 11074
  * Umeda et al. (2001) Umeda, T., Omura, Y., & Matsumoto, H. 2001, Computer Physics Communications, 137, 286
  * Umeki & Terasawa (1992) Umeki, H., & Terasawa, T. 1992, Journal of Geophysical Research: Space Physics, 97, 3113
  * Velli et al. (2020) Velli, M., Harra, L. K., Vourlidas, A., et al. 2020, A&A, 642, A4, doi: 10.1051/0004-6361/202038245
  * Verscharen et al. (2012) Verscharen, D., Marsch, E., Motschmann, U., & Müller, J. 2012, Physical Review E, 86, 027401
  * Viñas & Goldstein (1991) Viñas, A. F., & Goldstein, M. L. 1991, Journal of plasma physics, 46, 129
  * Wong & Goldstein (1986) Wong, H., & Goldstein, M. 1986, Journal of Geophysical Research: Space Physics, 91, 5617