3.1. Overview

Figures 1 and 2 provide an overview of the time evolution of the parametric decay in the simulations listed in table 1. We show, from top to bottom, the variation of the normalized root mean square (r.m.s.) of the magnetic and total density fluctuations, the correlation between velocity and magnetic field fluctuations ($\rho _{vB} = \boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {B} / \lVert \boldsymbol {v} \rVert \lVert \boldsymbol {B} \rVert$), the r.m.s. of the field-aligned electric field $e_x$, and the variation of the mean proton and electron temperatures normalized to their initial value. The growth rate of the PDI is reported in the legend and corresponds to the slope of the best fit indicated by the red solid line in the second panel down of each figure.

Figure 1. Temporal evolution of the variation of the r.m.s. of magnetic field fluctuations (top panel), r.m.s. of total density fluctuations (second panel), the correlation between magnetic and velocity fluctuations $\rho _{VB}$ (third panel) and r.m.s. of the field-aligned electric field (fourth panel). The fifth and bottom panels show the variation of the proton and electron mean parallel and perpendicular temperatures.The vertical black dashed lines represent the end of the linear stage of PDI for Ecsim-TiTe1.

Figure 2. Temporal evolution of the variation of the r.m.s. of magnetic field fluctuations (top panel), r.m.s. of total density fluctuations (second panel), the correlation between magnetic and velocity fluctuations $\rho _{VB}$ (third panel) and r.m.s. of the field-aligned electric field (fourth panel). The fifth and bottom panels show the variation of the proton and electron mean parallel and perpendicular temperatures. Results for runs Ecsim-TiTe4 (dot-dashed line), Ecsim-TiTe1 (dashed line) and Ecsim-TiTe1B (dotted line). The vertical black dashed lines represent the end of the linear stage of PDI for Ecsim-TiTe4 and Ecsim-TiTe1B. The grey dashed lines mark the times when distribution functions are shown in figure 3.

In figure 1 we compare results from the hybrid simulation and from the explicit (VPIC) and the semi-implicit (ECsim) PIC simulations, for $T_i=T_e$ and $\beta _e=0.125$. The decay of the pump wave, marked by the exponential growth of density fluctuations (second panel), proceeds similarly in the hybrid and PIC simulations. In both models the decay starts at the same time, around $t\varOmega _{ci}\simeq 150$, with the density fluctuations growing at the same rate. The saturation of the instability occurs when the backward-propagating Alfvén wave and forward-propagating ion-acoustic wave are well developed and density fluctuations reach a constant value of $\delta n_{rms}/\langle n\rangle \simeq 0.2$ at around $t\simeq 225 \varOmega _{ci}^{-1}$ (vertical dashed line). A complete reflection of the pump wave then takes place, as can be seen from the third panel, when the sign of $\rho _{VB}$ changes from $-1$ to $+1$. During this stage, a portion of the initial magnetic energy is converted into kinetic and internal energy of both protons and electrons. In particular, protons undergo strong parallel heating along the local magnetic field with significant perpendicular heating (fifth and bottom panels), as reported in previous studies (González et al. Reference González, Tenerani, Velli and Hellinger2020). The PIC simulations show that electrons are not isothermal during the decay process and that they are also energized. The heating of electrons is on average isotropic and proceeds at nearly a constant rate. This can be seen from the bottom panels of figure 1, showing that $T_{e\parallel }\simeq T_{e\perp }$ and that temperatures keep increasing during the nonlinear stage (specifically after $t\varOmega _{ci}\simeq 500$). We noticed a small difference in the electron temperature for the explicit and semi-implicit simulations with a slightly larger final temperature in the VPIC simulation. This might be due to numerical heating associated with explicit codes, although the relative error of the total energy is less than $1\,\%$ in both simulations.

Figure 2 shows results for simulations performed with ECsim with different $\beta$ and proton-to-electron temperature ratio $T_i/T_e$. As can be seen from these plots, the total plasma $\beta$ determines the growth rate and the overall decay of the pump wave in agreement with fluid theory (Wong & Goldstein Reference Wong and Goldstein1986; Hollweg Reference Hollweg1994), which predicts that higher-beta plasmas have a lower growth rate. This dependence can be seen by comparing runs ECSim-TiTe4 (dot-dashed) and TiTe1B (dotted), which have different electron temperature but the same total plasma beta, and by comparing ECSim-TiTe1B and TiTe1, corresponding to different values of the plasma beta. We observe that the higher the initial $T_e$, the higher the electric field, in agreement with the dependence of the electric field on $\beta _e$ from the generalized Ohm's law.

The most significant difference between the simulations with the same beta but different $T_i/T_e$ is registered in the r.m.s. of the electric field (second column, fourth panel: dotted versus dot-dashed line). We observe lower r.m.s. values of $e_x$ when $T_i>T_e$. This is consistent with Landau damping of the ion-acoustic mode, which is expected to be stronger when $T_i>T_e$ than in the case where $T_i=T_e$. Linear Landau damping, however, is not strong enough to affect the growth of the instability.

To conclude our overview of the decay process, in figure 3 we show the contour plots of the velocity distribution functions (VDFs) in velocity space, $f(v_\parallel,|v_\perp |)$, integrated in $x$. We plot the VDFs of protons (top panels) and electrons (bottom panels) when the system has achieved the nonlinear stage for runs ECSim-TiTe4, ECSim-TiTe1 and ECSim-TiTe1B. This time is marked by the vertical grey dashed lines in figure 2, at $t=400 \varOmega _{ci}^{-1}$, $t=560 \varOmega _{ci}^{-1}$ and $t=400 \varOmega _{ci}^{-1}$ for each simulation, respectively. As is discussed in the next sections, protons develop an anisotropic distribution function with a clear field-aligned beam propagating backward, and an energized field-aligned population at the (positive) ion-acoustic speed (see figure 6 for more details). Electrons instead remain on average isotropic, and their heating is much less pronounced than proton heating.

Figure 3. Reduced particle VDFs for protons (top) and electrons (bottom) integrated over the entire domain. Shown, from left to right, are the VDFs for Ecsim-TiTe4, Ecsim-TiTe1 and Ecsim-TiTe1B, at $t=400 \varOmega _{ci}^{-1}$, $t=560 \varOmega _{ci}^{-1}$ and $t=400 \varOmega _{ci}^{-1}$, respectively.

3.2. Energy conversion and particle heating

The energy lost by the pump wave powers the formation of daughter waves and results in thermal and parallel bulk energy gain for protons and electrons. The PIC simulations show that, after the decay of the pump wave, most of the particle internal energy gain goes into proton internal energy. In figure 4, we analyse the pump wave energy conversion and partition by breaking down the total energy (per unit volume) into its contributions from Alfvén wave energy (top panel), parallel kinetic energy (thus including the contribution from the ion-acoustic wave; second panel) and internal energy of protons and electrons (third and bottom panels). All quantities are expressed as variations from their initial value, and they are normalized to the initial pump wave energy $E_{w0}=E_{B_\bot }(0)+E_{v_\bot }(0)$, where

(3.1)\begin{equation} E_{B_\bot}(t)=1/2\langle B_y(t)^2+B_z(t)^2\rangle \end{equation}

is the magnetic energy density and

(3.2)\begin{equation} E_{v_\bot}(t)= \sum_\alpha 1/2 \langle m_\alpha n_\alpha (u_{y_\alpha}(t)^2+u_{z_\alpha}(t)^2) \rangle \end{equation}

is the bulk kinetic energy and the brackets denote spatial average. The $x$ component of the kinetic energy is defined as

(3.3)\begin{equation} E_{v_\parallel}(t)= \sum_\alpha 1/2\langle m_\alpha n_\alpha u_{x_\alpha}(t)^2\rangle \end{equation}

and the internal energy of species $\alpha$ as

(3.4)\begin{equation} E_{th,\alpha}=1/2\langle Tr(P_{ij,\alpha})\rangle. \end{equation}

Figure 4. Energy conversion and partition. Variation of the wave energy (top panel), of the field-aligned bulk kinetic energy (second panel) and the internal energy variations for protons and electrons (third and bottom panels, respectively). Each panel shows results for runs Ecsim-TiTe4 (dashed), Ecsim-TiTe1 (solid) and Ecsim-TiTe1B (dotted). The vertical dashed line indicates the end of the linear stage for the case $T_i/T_e=4$.

Energy balance is quantitatively similar in all our simulations. The top panel of figure 4 shows that the energy contained in the Alfvén waves decreases with $|\Delta E_w|/E_{w0} \simeq 60\,\%$. Notice that this does not correspond to the total pump wave energy decrease. The pump wave fully decays in the simulations considered here (see the correlation $\rho _{vb}$ in figure 2). After the decay ($\rho _{vb}=+1$ after $t\varOmega _{ci}\gtrsim 500$), $40\,\%$ of the pump wave energy goes into a reflected Alfvén wave, and the remaining $60\,\%$ goes in an ion-acoustic mode and other forms of internal energy. The second panel, as well as the trend of $\rho _{rms}$ reported in figure 2, shows that at the end of the linear stage (indicated by the vertical dashed line for the case $T_i/T_e=4$), $\Delta E_\parallel$ has increased due to the growth of the ion-acoustic mode. The subsequent peak of $E_\parallel$, however, is due to the ongoing particle trapping by the ion-acoustic mode that generates a population of protons moving at approximately the ion-sound speed, as discussed in § 3.3 and, for example, in Matteini et al. (Reference Matteini, Landi, Velli and Hellinger2010b). Particle trapping and subsequent phase-space mixing ultimately lead to the dissipation of the ion-acoustic mode and particle heating. After the complete decay of the pump wave, most of the Alfvén wave energy has gone into proton internal energy. In particular, we find that approximately $50\,\%$ of the initial pump wave energy goes into proton internal energy and $10\,\%$ to electron internal energy. Notice that, unlike protons, electrons keep heating at a constant rate when the instability has saturated and the pump wave has decayed. At the end of the simulation, the electron internal energy is a factor of 2 larger than the energy at the end of the complete decay ($t\varOmega _{ci}\gtrsim 500$). Still, the good energy conservation in the simulations and similar results between explicit and semi-implicit simulations suggest that the electron heating is indeed physical and not a numerical artefact. Such an energy balance appears to be unaffected by the total plasma beta and only slightly by $T_i/T_e$ and $m_i/m_e$ (not shown).

3.3. Phase-space dynamics

Compressibility plays an important role in proton and electron energization and in providing paths to particle heating. In the low-beta simulations presented in this work, we were able to identify and distinguish two types of compressible fluctuations generated during parametric decay: an ion-acoustic mode and a compressible steepened front of the fast-mode type propagating at the Alfvén speed. Such signatures were difficult to discern in higher-beta simulations (González et al. Reference González, Tenerani, Velli and Hellinger2020). Both types of compressible fluctuations contribute to the increase of internal energy. However, they lead to fundamentally different phase-space dynamics and to different kinetic features. The ion-acoustic mode contributes nonlinearly to phase-space mixing via particle trapping. On the other hand, the development of the steepened front with non-constant $|{\boldsymbol B}|$ leads to the formation of strong parallel heat fluxes and ultimately causes the acceleration of a population of protons into a well-defined beam of particles and the enhancement of proton perpendicular heating. In the following, we analyse the signatures of these processes in the fields and particle VDFs.

In figure 5 (left column) we show contour plots in the $(t,x)$ plane of $B_y$ (figure 5a), of the magnetic field magnitude $|{\boldsymbol B}|$ (figure 5b) and of the electric field $e_x$ (figure 5c). We also plot the scalar pressure agyrotropy (Aunai, Hesse & Kuznetsova Reference Aunai, Hesse and Kuznetsova2013; Cazzola et al. Reference Cazzola, Innocenti, Goldman, Newman, Markidis and Lapenta2016),

(3.5)\begin{equation} D_{ng}=\sqrt{8\left( P_{xy}^2 + P_{xz}^2 + P_{yz}^2 \right)}/(P_\parallel{+} 2 P_\perp), \end{equation}

and the parallel heat flux ($q_\parallel$) for protons (figure 5d,e) and electrons (figure 5f,g), respectively. The thermal energy flux is obtained by calculating the total energy flux for each species ($\boldsymbol {Q}_s = \tfrac {1}{2}m_s \int f_s(\boldsymbol {x},\boldsymbol {v},t) v^2 \boldsymbol {v} \,{\rm d}^3\boldsymbol {v}$) and subtracting the energy flux components in the Eulerian frame ($\boldsymbol {q}_s = \boldsymbol {Q}_s - \boldsymbol {H}_s - \boldsymbol {K}_s$), namely the contribution from the bulk energy flux ($\boldsymbol {K}_s = \boldsymbol {u}_s ({\rho u_s^2}/{2})$) and the enthalpy flux ($\boldsymbol {H}_s = P_s \boldsymbol {\cdot } \boldsymbol {u}_s + \boldsymbol {u}_s U_{th}$). This requires the calculation of the lower moments of the velocity distribution (see details in Lapenta et al. (Reference Lapenta, El Alaoui, Berchem and Walker2020)). We show results from ECSim-TiTe4 as a reference, since ECSim-TiTe1 and ECSim-TiTe1B are similar. In the right column of figure 5 we show the r.m.s. of the divergence of the parallel heat flux (normalized to its initial value) and the standard deviation of the scalar agyrotropy for both species and for all simulations.

Figure 5. Left panels: contour plot in the $(t,x)$ plane of the $b_y$ component (a), the magnitude of the magnetic field (b), the field-aligned component of electric field (c), and the scalar agyrotropy (d) and heat flux (e) for protons and the scalar agyrotropy (f) and heat flux (g) for electrons. Results are shown for ECSim-TiTe4. Right panels: r.m.s. of the divergence of the parallel heat flux for protons (top panel) and electrons (second panel), and variance of the agyrotropy for protons (third panel) and electrons (bottom panel) for all simulations.

From figure 5 (left column) one can see that a dominant backward-propagating mode forms at $t\varOmega _{ci}\simeq 200$, when the PDI has reached saturation. The signature of the forwards ion-acoustic mode can be seen in the plots of $e_x$ in the time interval $t\varOmega _{ci}=200\unicode{x2013}400$, before it dissipates, corresponding to a wave vector $k_s=7 \times 2{\rm \pi} /L_x$. At about $t\varOmega _{ci}\simeq 400$, a steepened front propagating backwards at the Alfvén speed also develops. This steepened front corresponds to a localized increase of $|\boldsymbol {B}|$, positively correlated with a density enhancement. As discussed in our previous work (González et al. Reference González, Tenerani, Velli and Hellinger2020, Reference González, Tenerani, Matteini, Hellinger and Velli2021), compressive fluctuations in $|{\boldsymbol B}|$ are associated with a field-aligned electric field $e_x$. Enhancements of $e_x$ at the steepened front can also be seen in figure 5(c). The strongest non-thermal features (agyrotropy and heat fluxes) of particles develop at such a steepened front and within the ion-acoustic mode, as shown in figure 5(c–g).

The cumulative contribution of non-Larmor particle motion produces non-zero values of the off-diagonal terms of the pressure tensor, quantified by $D_{ng}$, where local values of up to $|D_{ng,e}|\simeq 0.3$ and $|D_{ng,i}|\simeq 0.65$ are reached. The variance of $D_{ng}$ for both species is reported in the third and fourth panels on the right-hand side of figure 5, and the r.m.s. of the divergence of the parallel heat flux in the first two panels. Agyrotropy in particular is considered one of the VDF-based dissipation measures (Pezzi et al. Reference Pezzi, Liang, Juno, Cassak, Vásconez, Sorriso-Valvo, Perrone, Servidio, Roytershteyn and TenBarge2021). From figure 5 (right panels) one can see that agyrotropy first builds up during the linear stage of the PDI for both species in correspondence with the ion-acoustic mode, before decaying during the nonlinear stage as the ion-acoustic mode is dissipated and, correspondingly, the plasma is heated. However, protons, unlike electrons, undergo a two-step heating process. Besides particle trapping, protons subsequently interact with the steepened front, where they develop a second peak of agyrotropy. Interaction of protons with a field-aligned discontinuity and electric field causes both proton scattering in phase space, contributing to heating, and the acceleration of protons into a field-aligned beam (González et al. Reference González, Tenerani, Matteini, Hellinger and Velli2021). This is consistent with the observed time evolution of proton agyrotropy and parallel heat flux (see top and third panels on the right-hand side of figure 5). Electrons, which have a smaller gyroradius than protons, display on average $D_{ng,e}\ll D_{ng,i}$ (bottom right panel) and, accordingly, they do not gain much internal energy.

Kinetic features in phase space corresponding to the dynamical evolution described above can be seen in figure 6, where we show the proton and electron VDFs in the $(x,v_\parallel )$ space (top) and $(x,v_\perp )$ space (bottom) at $t\varOmega _{ci}=400$. We show results for EcsimTiTe4 only since similar features develop in the other simulations. Protons in phase space are shown on the left and electrons on the right. For reference, we plot $B_y(x)$ (green colour) and $|\boldsymbol {B}(x)|$ (red colour) at the same time to show that the particle beam forms at the steepened front. The contour plot of $f(x,v_\parallel )$ shows the two sources of the proton internal energy increase discussed above. First, saturation of the growth of the ion-acoustic mode is determined by nonlinear particle trapping. The latter causes the formation of a second proton population that propagates forwards at the ion-acoustic speed, and develops well-known signatures of phase-space mixing (vortexes). Second, a proton beam is generated locally at the backward-propagating steepened Alfvén front. The latter can be identified in figure 6 at around $x \simeq 85 \, d_i$, and has a size of a few $d_i$. Both particle trapping and beam formation cause the large increase of the parallel proton temperature and parallel heat flux. One can see in the bottom left-hand panel that also perpendicular heating takes place at the steepened edge, where particles with high perpendicular velocities are observed. The perpendicular energization is due to particle scattering at those discontinuities (González et al. Reference González, Tenerani, Matteini, Hellinger and Velli2021). The electrons are also characterized by a local deviation from Maxwellian distribution at the steepened front. The inset in the top right-hand panel shows the electron VDF integrated in $x$ over an interval $\Delta x=3 \, d_i$, at different positions in the simulation domain: centred at the location of discontinuity $x=85 \, d_i$ (blue line) and outside the discontinuity, at $x=60 \, d_i$ and $x=95 \, d_i$ (orange and green lines), where the VDF is Maxwellian. The electrons at the discontinuity show a flat-top distribution with enhanced number of electrons in that region as to maintain the quasi-neutrality of the system (we do not find evidence of charge–space separation in the simulations). We do not observe any electron beam in the simulations, and the electrons undergo on average isotropic heating.

Figure 6. The phase space ($x,v_\parallel /v_\bot$) for protons (left) and electrons (right). The phase space is shown at $t\varOmega _{ci}=400$ for Ecsim-TiTe4. The magnetic field's magnitude (red line) and $B_y$ (green line) are also presented. The inset at top right represents the electron distribution function integrated in space around the discontinuity (blue) and outside the discontinuity (orange and green).