Kinetic Alfvén wave cascade in sub-ion range plasma turbulence

Sharma, Johan; Makwana, Kirit D.

doi:10.3389/fspas.2024.1423642

ORIGINAL RESEARCH article

Front. Astron. Space Sci., 29 October 2024

Sec. Space Physics

Volume 11 - 2024 | https://doi.org/10.3389/fspas.2024.1423642

Johan Sharma*

Kirit D. Makwana

Department of Physics, Indian Institute of Technology Hyderabad, Sangareddy, Telangana, India

Kinetic Alfvén waves (KAWs) are simulated with a 3D particle-in-cell (PIC) code by using the eigenvector relations of density, velocity, electric, and magnetic field fluctuations derived from a two-fluid KAW model. Similar simulations are also performed with a whistler waves setup. The 2D two-fluid eigenvector relations are converted into 3D by using rotation of the reference frame. The initial condition for the simulations is a superposition of several waves at scales slightly larger than the ion skin depth. The nonlinear interactions produce a transfer of energy to smaller scales. The magnetic field perturbation ratios, velocity perturbation, and density perturbation ratios are calculated from the simulation at higher wavenumbers and compared with the analytically expected ratios for KAWs and whistler waves. We find that in both types of simulations, initialized either with an ensemble of KAWs or with whistlers, the observed polarization relations at later times match better with the KAW relations compared to whistlers. This indicates a preference for excitation of KAW fluctuations at smaller scales. The power spectrum in the perpendicular direction is calculated, and it shows similar indices as measured in the solar wind power spectrum in the transition (sub-ion) region. The power law extends to smaller scales when a higher ion-to-electron mass ratio is taken. The 2D magnetic power spectrum in magnetic field parallel and perpendicular directions shows typical anisotropy where the power spreads more in the perpendicular direction than in the parallel direction. This study shows that KAWs can explain features of the sub-ion range plasma turbulence in the solar wind.

1 Introduction

The solar wind is an excellent laboratory for studying plasma turbulence (Bruno and Carbone, 2013). Observations of frequency spectra of velocity and magnetic fluctuations in the solar wind in an intermediate range of frequencies ( $\sim 1 0^{- 4} - 1 0^{- 1}$ Hz) (Coleman, 1968; Podesta et al., 2007; Roberts, 2007; Bamert et al., 2008) have revealed scaling similar to that of hydrodynamic turbulence of Kolmogorov $(- 5 / 3)$ (Kolmogorov, 1941) and/or Iroshnikov-Kraichnan $(- 3 / 2)$ scaling (Kraichnan, 1965; Iroshnikov, 1963). Many attempts have been made to understand these observations through theory and simulation models. There are a wide variety of models ranging from the observationally driven large-scale models of magnetohydrodynamics (MHD) to the small-scale models of kinetic scales. MHD models provide the overall picture of the solar system plasma, while kinetic models promise to explain the physical processes that result in the heating and acceleration of solar wind (Ofman, 2010).

The shear Alfvén waves (Alfvén, 1942) are often found in the solar wind (Belcher and Davis, 1971). The incompressible MHD model is often used to describe the turbulence generated by Alfvén waves (Goldreich and Sridhar, 1995; Galtier et al., 2000; Boldyrev, 2006; Matthaeus et al., 1983; Cho and Vishniac, 2000; Bhattacharjee and Ng, 2001; Schekochihin, 2022). There could also be compressible slow and fast magnetosonic modes present in MHD turbulence (Makwana and Yan, 2020). As the energy cascades to smaller scales anisotropically (Shebalin et al., 1983), the perpendicular scales become comparable to the ion gyro-radius. As a result, ions cannot follow magnetic field lines, whereas electrons can continue to do so, leading to finite Larmor radius effects and kinetic Alfven waves (KAWs) (Hasegawa and Uberoi, 1982). KAWs are thought to play a role in the heating and acceleration of the solar wind (Ofman, 2010; Voitenko and Goossens, 2006; Narita et al., 2020). Recent observations also show a turbulent power law spectrum below the ion-scales (Alexandrova et al., 2009). The existence of a KAW cascade is posited to explain this kinetic range turbulence (Shaikh and Zank, 2009; Voitenko and Goossens, 2004; Howes et al., 2008; Sahraoui et al., 2009; Schekochihin et al., 2009; Sahraoui et al., 2010; Howes et al., 2011; Boldyrev and Perez, 2012; TenBarge et al., 2012; Told et al., 2015; Chen, 2016; Chen and Boldyrev, 2017; Cerri et al., 2019; Passot and Sulem, 2019). At the same time whistler waves have also been proposed to play a role in solar wind heating (Chang et al., 2011; Saito and Peter Gary, 2012) and in the kinetic range cascade (Gary et al., 2012; Gary and Smith, 2009). Salem et al. (2012) compared the kinetic Alfvén wave and whistler waves model predictions with magnetic and electric field measurements in the solar wind. They found that the turbulent fluctuations at small scales are consistent with the KAW spectrum. Similar results were obtained by Chen et al. (2013), Sahraoui et al. (2009, 2010) in their studies.

Kobayashi et al. (2017) and others have shown that there is a steeper spectral law in the transition range (near the ion characteristic scale) than at smaller scales [which show a slope of around $- 2.8$ (Alexandrova et al., 2009; Sahraoui et al., 2010; Sahraoui et al., 2013; Matteini et al., 2017; David and Galtier, 2019)]. The physics behind this steeper slope in the transition range is not well understood. Some studies give possible explanations for this transition range like the weak or imbalanced KAW turbulence (Voitenko and De Keyser, 2011; Passot and Sulem, 2019), ion-cyclotron waves (Bruno and Trenchi, 2014; Bruno and Telloni, 2015; Squire et al., 2022; Bowen et al., 2024), weak dispersive range (WDR) formed by co-collisions between co-propagating waves along mean magnetic field $B_{0}$ (Voitenko and De Keyser, 2016; Passot et al., 2022; Passot et al., 2024) and helicity-barrier hypothesis (Meyrand et al., 2021). Huang et al. (2020) showed the coexistence of KAWs and Alfvén ion cyclotron waves (AICW) in the slow solar wind in the inner heliosphere using NASA’s Parker Solar Probe. They measured the magnetic power spectra of KAW in the ion-transition range and found the slope of around $k_{⊥}^{- 3.73}$ . Huang et al. (2021) used the high-resolution data from the Parker Solar Probe and found that the power spectra in the sub-ion range do not follow a single turbulence scaling and found the spectral index varying in the range −3 to −5.7. These results are similar to the reported results by Bruno and Trenchi (2014), Bruno et al. (2014), Kobayashi et al. (2017), Huang et al. (2020).

First-principle numerical simulations can help in disentangling the different modes present in this turbulence. Earlier simulation studies on sub-ion range turbulence have involved initial conditions which are quite ideal and nonlinear (Olshevsky et al., 2018), or uses Alfven modes for initialization (Franci et al., 2018; Grošelj et al., 2018) or there is random forcing (Kobayashi et al., 2017), or with forcing containing compressive terms (Cerri et al., 2016) or using OU forcing (Arzamasskiy et al., 2019) or Alfven forcing (Grošelj et al., 2019) or a general initial perturbation (Camporeale and Burgess, 2011). There are other simulations that are restricted to 2D (Matteini et al., 2020; Makwana et al., 2023; Parashar et al., 2010; Parashar et al., 2011; Valentini et al., 2014; Franci et al., 2015; Servidio et al., 2015; Franci et al., 2016; Cerri and Califano, 2017; Cerri et al., 2017; Grošelj et al., 2017; Perrone et al., 2018; Pezzi et al., 2018; Matthaeus et al., 2020; Franci et al., 2022), or discusses the effect of collisions (Pezzi et al., 2017) or discusses the generation of sub-ion scale turbulence by means of reconnection (Franci et al., 2017), or only looks at electron scales (Califano et al., 2020; Muñoz et al., 2023; Gary et al., 2012). In some studies, they also obtain polarization ratios but either look at more general compressibility ratios (Matteini et al., 2020) or do not compare with different wave modes (Grošelj et al., 2017) or compare ratios at different $β$ values (Cerri et al., 2016). There are also 3D simulations but with the main focus on spectrum (Howes et al., 2011; Gary et al., 2012; TenBarge et al., 2013; Chang et al., 2015; Cerri et al., 2018; Cerri et al., 2016) and do not discuss the polarization ratios. Some studies look at dissipation and energy spectrum with gyrokinetic simulations (Told et al., 2015). Cerri et al. (2019) used collection of simulations described above and found similar conclusions of KAW. Other simulations analyze current sheet structures (Sisti et al., 2021) or ion cyclotron waves (Squire et al., 2022) or discuss linear and nonlinear timescales in the kinetic range turbulence (Gary et al., 2020) or look at stochastic in heating with predictions from Alfven and KAW scalings (Cerri et al., 2021).

In this study, we perform 3D fully-kinetic PIC simulations of sub-ion range turbulence. Instead of arbitrary perturbations or forcing, we carefully initialize these simulations by using a superposition of either KAWs or whistlers setup according to their eigenvector (i.e., polarization) relations derived from a two-fluid model. This allows for separately studying the nonlinear interactions of these modes. We find that these larger-scale waves interact with each other to produce a cascade of energy to smaller scales. We then compare the various polarization relations at higher wavenumbers with the expected polarization relations of whistlers and KAWs. This allows us to identify the modes produced by the forward cascade. We find a preference for the KAW cascade, even in simulations that are initialized with whistlers. We also measure the power spectrum and compare it with observations. Hence our study is complementary to all the above studies, which have also hinted at the presence of KAWs in the sub-ion range turbulence. In Section 2 we describe the setup of the numerical simulations, including the two-fluid relations for KAW and whistlers and the way they are used to provide the initial conditions for the simulations. In Section 3 we compare the analytical dispersion relation of KAW from the two-fluid model and that obtained from 2.5D and 3D simulations. Comparison of magnetic polarization ratios obtained from simulations with the analytically expected ratios of KAWs and whistlers are shown in Section 3.1. In Sections 3.2, 3.3, we compare the ion-density perturbation and velocity perturbation ratios, respectively. Section 4.1 shows the magnetic and velocity field spectrum obtained from the simulation at different time steps and its spectral index. Here, simulations with different plasma beta and ion-to-electron mass ratios are also shown. In Section 4.2 we present the 2D magnetic energy spectrum in $k_{⊥} - k_{‖}$ plane and its anisotropy. The conclusions and implications are discussed in Section 5.

2 Setup of numerical simulations

We aim to simulate the nonlinear interactions of KAW and whistler waves in 3D particle-in-cell simulations. The initial condition of the simulations is a superposition of these waves at scales close to the ion skin-depth $(k d_{i} ≲ 1)$ with the goal of seeing a cascade of energy to smaller scales. In order to set up the KAW modes, we use the Hollweg (1999) (from now on referred to as H99) relations of two-fluid KAW. The two-fluid eigenvector relations in H99 paper are derived by considering $k_{x} = 0$ , i.e., the relations are valid only for wave-vectors in the $k_{y} - k_{z}$ plane. We are using these eigenvector relations in a 3D simulation set in a reference frame $S^{'}$ . In order to do this we need to generalize the eigenvector relations to 3D. To accomplish this we take the initialization wavevectors in 3D and then rotate our frame of reference around the $z^{'}$ -axis such that the $k_{x}$ becomes zero in this new reference frame $S$ . In this rotated frame we can use the eigenvector relations from the H99 paper. Let the primed symbols $k_{x}^{'}, k_{y}^{'}, k_{z}^{'}$ and $δ E_{x}^{'}, δ E_{y}^{'}, δ E_{z}^{'}$ represent the wave-vector and the electric field perturbations respectively in the simulation frame $S^{'}$ . The unprimed symbols $k_{x}, k_{y}, k_{z}$ and $δ E_{x}, δ E_{y}, δ E_{z}$ represent the corresponding quantities in the rotated reference frame $S$ in which the eigenvectors from H99 hold. In the simulation, we take $k_{x}^{'}, k_{y}^{'}$ , and $k_{z}^{'}$ as input, and then we rotate our frame of reference along the z-axis by an angle $ϕ$ to unprimed coordinates. The $δ E$ and other quantities are obtained in this frame from H99 relations, and then rotated back to frame $S^{'}$ to use them in initializing the simulations (see Supplementary Appendix for more details). Similar process is used for whistler wave as well.

2.1 KAW polarization relations

The eigenvector relations express the ratio of the different perturbed physical quantities in a KAW mode. Their amplitude is defined up to a multiplicative constant, which we set by setting the amplitude of $δ E_{x}$ . The rest of the perturbed quantities are defined by the eigenvector relationships. These relations are used to setup the initial wave modes in the simulation, as well as to identify the wave modes produced by the cascade. In frame $S$ the KAW electric field eigenvector relations are taken from H99, Makwana et al. (2023).

\frac{δ E_{z}}{δ E_{x}} = i \frac{k_{z} (ω^{2} - k^{2} v_{A}^{2})}{k_{y} ω ω_{c, i}} \frac{(G_{e} - m_{e} ω^{2} / k_{z}^{2})}{G_{t}}, (1)

\frac{δ E_{y}}{δ E_{x}} = i \frac{(ω^{2} - k^{2} v_{A}^{2}) (Ω^{2} (m_{e} v_{A}^{2} - G_{i}) - G_{e})}{G_{t} (Ω^{2} - 1) ω ω_{c, i}} . (2)

Here the index $i$ represents ions, $e$ represents electrons, $s$ represents species and $t$ represents total. G_t, G_i, and G_e are given by Equation 3,

G_{t} = G_{i} + G_{e} = γ_{i} κ T_{0, i} + γ_{e} κ T_{0, e}, (3)

$Ω = \frac{ω}{k_{z} v_{A}}$ , ion cyclotron frequency $ω_{c, i} = \frac{e B_{0}}{m_{i} c}$ , and Alfvén speed is given by Equation 4,

v_{A}^{2} = \frac{B_{0}^{2}}{4 π n_{0} (m_{i} + m_{e})} . (4)

where, $B_{0} = B_{0} \hat{z}$ is mean magnetic field.

The electric field polarization in the simulation reference frame $(S^{'})$ is obtained using rotation matrix (see Supplementary Appendix). The same rotation relations are used to obtain the eigenvector relations for the magnetic and velocity fields in the simulation frame. Since density is a scalar, it remains unchanged $δ n_{s}^{'} = δ n_{s}$ . The dispersion relation of KAW from the 2-fluid model is

\begin{align} Ω^{6} \frac{k_{z}^{2}}{k^{2}} (1 + k_{y}^{2} L_{e}^{2}) - Ω^{4} [\frac{k_{z}^{2}}{k^{2}} (1 + k_{y}^{2} L^{2}) + (1 + β + k_{y}^{2} L_{e}^{2})] \\ + Ω^{2} (1 + 2 β + k_{y}^{2} L^{2}) - β = 0, \end{align} (5)

where, β, ω_p,e, L, and L_e are given by Equations 6–9 respectively

β = β_{t o t} = \frac{G_{t}}{m_{i} v_{A}^{2}}, (6)

ω_{p, e}^{2} = 4 π q_{e}^{2} n_{0} / m_{e}, (7)

L = \frac{1}{| ω_{c, i} |} \sqrt{\frac{G_{t}}{m_{i}}} = \frac{v_{s}}{| ω_{c, i} |}, (8)

L_{e} = \frac{c}{ω_{p, e}} = \frac{L}{\sqrt{β}} \frac{\sqrt{1 + m_{r}}}{m_{r}} . (9)

Here $m_{r} = m_{i} / m_{e}$ is the mass ratio of ion to electron. Equation 5 has three roots; the intermediate value of these three roots is the KAW branch. The magnetic perturbation ratios give some important characteristics of the waves in plasma. The magnetic field perturbations $δ B$ are obtained from the electric field perturbations using Faraday’s law as

δ B = \frac{c}{ω} (k \times δ E) . (10)

The magnetic perturbations $δ B$ in frame $S$ are obtained from Equation 10 and then rotated to obtain $δ B^{'}$ in the simulation frame. The theoretically expected magnetic perturbation ratios $δ B_{x}^{'} / δ B_{y}^{'}$ and $δ B_{x}^{'} / δ B_{z}^{'}$ are then obtained from this.

KAWs also have density perturbations. The relation for ion number density perturbation ratio for KAW is given by Equation 11,

\frac{δ n_{s}}{n_{0, s}} = \frac{i e δ E_{z}}{k_{z} (G_{e} - m_{e} ω^{2} / k_{z}^{2})} . (11)

We have kept only the terms with lowest order in $ω / ω_{c, s}$ since $ω ≪ ω_{c, s}$ . For the electron-proton plasma $n_{0} = n_{0, i} = n_{0, e}$ , and we use the quasineutrality condition for electron density perturbation, $δ n_{i} = δ n_{e}$ . The velocity perturbations for KAWs from the 2-fluid model are

δ v_{x, s} = \frac{q_{s}}{m_{s} ω_{c, s}} (δ E_{y} - i \frac{ω}{ω_{c, s}} δ E_{x}) - i \frac{G_{s} k_{y}}{m_{s} ω_{c, s}} \frac{δ n_{s}}{n_{0}}, (12)

δ v_{y, s} = - \frac{q_{s}}{m_{s} ω_{c, s}} (δ E_{x} + i \frac{ω}{ω_{c, s}} δ E_{y}) - \frac{G_{s} k_{y} ω}{m_{s} ω_{c, s}^{2}} \frac{δ n_{s}}{n_{0}}, (13)

δ v_{z, s} = i \frac{q_{s}}{m_{s} ω_{c, s}} δ E_{z} + \frac{G_{s} k_{z}}{m_{s} ω} \frac{δ n_{s}}{n_{0}} . (14)

Boldyrev et al. (2013) derive the KAW relations from kinetic theory. They obtain KAW relations from full kinetic equations, by use of certain approximations. In their paper, to obtain KAW relations they make simplifications to Equation 18 and obtain the polarization ratios given in Equations 42, 43. Our relations are similar to those, as shown below. Dividing Equations 1, 2, then in the denominator, we have $G_{i} ≫ m_{e} v_{A}^{2}$ , and also, we are analyzing in the range $Ω ≳ 2$ . Also, $G_{e} ≫ m_{e} ω^{2} / k_{z}^{2}$ since $v_{t h, e} ≫ v_{A}$ . Now if we take $Ω^{2} - 1 \sim Ω^{2}$ , then our polarization ratio reduces to

\frac{δ E_{z}}{δ E_{y}} \approx \frac{k_{z} G_{e} Ω^{2}}{- k_{y} Ω^{2} G_{i}} = - \frac{k_{z}}{k_{y}} \frac{T_{e}}{T_{i}} (15)

Since Boldyrev et al. (2013) assume $k_{y} = 0$ , whereas, in this work, $k_{x} = 0$ , a rotation of frame is required for comparison. Their frame transformed by $y \to - x$ , and $x \to y$ gives our reference frame. Then Equation 42 of Boldyrev et al. (2013) becomes Equation 15. Now Equation 43 in Boldyrev et al. (2013) is the inverse of our Equation 2. In the denominator of Equation 2, we can say $ω^{2} - k^{2} v_{A}^{2} \sim - k_{⊥}^{2} v_{A}^{2}$ and using above approximations we get Equation 43 of Boldyrev et al. (2013). So our equations reduce to Boldyrev et al. (2013) expressions except for the frequencies, which are their respective frequencies. By comparing the frequency expression of Boldyrev et al. (2013) with our frequencies, we find that in the range $1 < k_{⊥} d_{i} < 3$ these frequencies differ by less than 30%. The Boldyrev et al. (2013) frequency is obtained by assuming $k_{⊥} d_{i} ≫ 1$ , whereas for our case where we are looking at the modes excited close to the ion skin depth scale, $k_{⊥} d_{i} ≳ 1$ . In this case, the H99 expressions are more appropriate.

2.2 Whistler wave polarization relations

For obtaining the whistler wave relations we use Swanson (2003) (SW03). The mean magnetic field $B_{0}$ in SW03 is the same as in our case, but the wave vector $k$ in SW03 is in $x$ - $z$ plane, but it is in $y$ - $z$ plane in our case. To express the relations of SW03 we give a rotation by 90 $°$ so that $y \to - x$ and $x \to y$ , so the electric field components also become $δ E_{x} \to δ E_{y}$ and $δ E_{y} \to - δ E_{x}$ . Also, we use a signed definition of $ω_{c, e}$ , which is negative. Then the dispersion relation for the whistler wave at an arbitrary angle $θ$ is given by Equation 16,

ω = - \frac{k^{2} c^{2} ω_{c, e} \cos θ}{ω_{p, e}^{2}}, (16)

where $θ$ is the angle between the mean magnetic field. The wave equation for cold plasma is also modified to our frame of reference, given by Equation 17,

(\begin{matrix} i D & S - n^{2} \cos^{2} θ & n^{2} \cos θ \sin θ \\ n^{2} - S & i D & 0 \\ 0 & n^{2} \cos θ \sin θ & P - n^{2} \sin^{2} θ \end{matrix}) (\begin{matrix} δ E_{x} \\ δ E_{y} \\ δ E_{z} \end{matrix}) = 0 . (17)

As $δ E_{z} = 0$ for whistler wave, so the equation becomes

(\begin{matrix} i D & S - n^{2} \cos^{2} θ \\ n^{2} - S & i D \end{matrix}) (\begin{matrix} δ E_{x} \\ δ E_{y} \end{matrix}) = 0 . (18)

So we get the electric field perturbation ratio $δ E_{y} / δ E_{x}$ , given by Equation 19,

\frac{δ E_{y}}{δ E_{x}} = \frac{- i D}{S - n^{2} \cos^{2} θ} . (19)

Here, S and D are given by Equations 20, 21,

S = 1 + \frac{ω_{p, i}^{2}}{ω_{c, i}^{2} - ω^{2}} + \frac{ω_{p, e}^{2}}{ω_{c, e}^{2} - ω^{2}}, (20)

i D = - \frac{i ω_{c, e} ω_{p, e}^{2}}{ω (ω_{c, e}^{2} - ω^{2})} - \frac{i ω_{c, i} ω_{p, i}^{2}}{ω (ω_{c, i}^{2} - ω^{2})}, (21)

and $n = \frac{c k}{ω}$ . Also $δ E_{z} / δ E_{x} = 0$ for whistler wave.

We can obtain the magnetic field perturbation $δ B$ for whistler waves again by using Faraday’s law as in Equation 10 for KAWs. Then, we use the same rotation of the reference frame as done for 2-fluid KAW relations and obtain the magnetic perturbation ratios in full 3D. Also, the velocity perturbation for the whistler wave has been obtained from SW03. They obey the same relations as Equations 12–14, except that the number density perturbations for whistlers are zero. Using the same method of the rotation of the reference frame, these relations are then used to obtain the velocity perturbation ratios for the 3D case.

2.3 Simulation setup

We use the implicit particle in cell code iPIC3D (Lapenta, 2012; Markidis et al., 2010) for our simulations. The simulations are initialized with either an ensemble of superimposed KAWs or an ensemble of superimposed whistlers. This allows us to study the nonlinear interactions of the wave modes amongst themselves. The polarization relations shown in Sections 2.1, 2.2 are used to set up the KAWs and whistlers, respectively. Normalization in the simulation has been taken to make length in units of ion-skin depth $d_{i}$ , velocity in units of speed of light $c$ , and time in units of inverse plasma ion frequency $ω_{p, i}^{- 1}$ $(ω_{p, i} = c / d_{i})$ . The mass and charge for ions are taken as unity ( $e = 1$ and $m_{i} = 1$ ). The number density is normalized to the background density, $n_{0}$ , and due to the unit mass of ions, the normalized mass density of ions is unity. The magnetic field is normalized to $c \sqrt{4 π ρ_{0}}$ , where $ρ_{0}$ represents the total mass density. The electric field follows the same normalization. In these units, the Alfvén speed normalized to the speed of light numerically corresponds to the normalized background field strength, ${\hat{B}}_{0} = v_{A} / c = 0.01$ . The eigenvector relations described above are converted into these normalized units for utilization within the code.

We have used periodic boundary conditions in the three-dimensional simulation box. The uniform plasma is initialized with a uniform magnetic field $B_{0}$ along $z$ -axis. The information about the length of the box, resolution, and thermal velocities of electrons and ions for the simulations is listed in Table 1. The modes for simulation initialized with KAW are chosen such that $k_{⊥} ≫ k_{‖}$ for each wave. For whistler simulation, we have used a more isotropic distribution of wave modes. In both simulations, we have taken care that the wavenumber modes have small values, $k d_{i} < 1$ (because we are injecting modes at large scales). We have used 64 quasi-particles per species per cell for all simulations. A superposition of different wave vectors is used as the initial condition where each physical quantity, $ψ$ , is given by Equation 22,

ψ (r) = ψ_{0} + R [\sum_{k} δ ψ \exp (i k \cdot r + ϕ_{k})], (22)

with $ψ$ representing the components of the electric field, magnetic field, ion velocity, electron velocity, ion density, and electron density. A random phase $ϕ_{k}$ is put for each wave-vector, $ψ_{0}$ is the equilibrium field value and $δ ψ$ is the perturbation amplitude derived from the polarization relations. The electric and magnetic fields are directly initialized on the grid. The density and velocity perturbations are used in the initial Maxwellian distribution of particles. For initialization, the Maxwellian distribution for a given temperature $T_{s}$ is used where $k_{B} T_{s} = m_{s} v_{t h, s}^{2} / 3$ with thermal speed $v_{t h, s}$ . The fluid velocities are added as a drift velocity in the exponential term of the Maxwellian. The wave-vector is $k = k_{x} \hat{i} + k_{y} \hat{j} + k_{z} \hat{k}$ , where $k_{x, y, z} = n_{x, y, z} (2 π / L_{x, y, z})$ with wavenumbers $n_{x, y, z}$ . The initial wavevectors are large-scale, of the size of the simulation box, such that $k d_{i} ≲ 1$ . This allows us to look at the cascade of energy to smaller scales by the nonlinear interactions. For KAWs, we have used $| n_{z} | < \sqrt{n_{x}^{2} + n_{y}^{2}}$ so that they satisfy $k_{⊥} > k_{‖}$ . 16 wavevectors are used in the initial conditions, where S0 wave input is used for KAW simulations, and the S1 wave input is used for whistler simulations. They are listed in Table 2. We use the amplitude of electric field perturbation in $x$ -direction $(δ E_{x})$ as an input parameter and all the other perturbation amplitudes are obtained from the eigenvector relations. The angular frequency in these relations is used from the 2-fluid dispersion relations. The simulation then solves for the time evolution of the plasma with these initial conditions.

Table 1

Table 1. Simulation parameters.

Table 2

Table 2. List of simulations, the number of wavenumbers in the initial condition, amplitude of waves, initial perturbation amplitude, and the initial mode-numbers.

3 Dispersion and polarization properties

We check the linear dispersion properties of the KAW by comparing the frequency measured in KAW simulations with the frequency calculated from the two-fluid model and frequency from the NHDS (New Hampshire Dispersion Solver) which is a dispersion solver for the fully kinetic case (Verscharen and Chandran, 2018; Verscharen et al., 2013). We take the KAW root of Equation 5 and evaluate it for different values of the wavenumbers. Then simulations are run with single KAW modes corresponding to these wavenumbers. The frequencies from the simulation are obtained by using the magnetic field $B_{x}$ propagation in y-direction. The wave is observed at two different time steps, and the distance traveled by this wave in the y-direction, $Δ y$ , in this time $Δ t$ is measured. The observed frequency is then $ω = k_{y} Δ y / Δ t$ . This is done for both 2.5D and 3D simulations and the results are shown in Figure 1 where the wavenumbers used are $n_{y}$ = 2, 3, 4, 6, 8, 10, $n_{z}$ = 1, and $n_{x} = 0$ . We find that the frequency from the 2.5D simulations matches very closely with the 3D simulation results which means that the 3D simulation gives expected results. We have also run a 2.5D KAW simulation with input wavenumber $n_{x}, n_{y}, n_{z} = 0,5,1$ which gives a frequency of 0.00394 $ω_{p, i}$ and a 3D KAW simulation with $n_{x}, n_{y}, n_{z} = 3,4,1$ giving a frequency of 0.00402 $ω_{p, i}$ . Both the simulations have the same $k_{‖}$ and $k_{⊥}$ and give almost the same result for frequency. This further verifies the accuracy of the rotation and, hence, the 3D simulation.

Figure 1

Figure 1. Comparison of frequencies $(ω / ω_{p, i})$ from 2-fluid KAW, NHDS, 2.5D and 3D KAW simulation with varying $k_{y} d_{i}$ at constant $k_{z}$ . Here $n_{z} = 1$ and $n_{y}$ = 2, 3, 4, 6, 8, 10 and A0 simulation parameters has been used.

Figure 1 shows that the observed frequency follows the trend but is smaller than that expected from the two-fluid dispersion relation. This is very similar to the result obtained in Makwana et al. (2023) for similar simulation parameters. It was shown there that the observed dispersion relation matches much better with a hot plasma dispersion relation, even though we use a two-fluid description in the initialization. That showed that such an initialization is suitable for PIC simulations. Here also, we find that the NHDS dispersion relation results are very close to those observed in our 3D simulations. Hence we conclude that the 2-fluid eigenvector is a good initialization for 3D simulation also, while being easier to implement than a fully kinetic initialization. We also checked for the value of $k^{'} \cdot r^{'}$ before and after the rotation, and it remains the same, i.e., $k^{'} \cdot r^{'}$ = $k \cdot r$ , as expected. This means that the waveform factor $\cos (k_{x} x + k_{y} y + k_{z} z + ϕ_{k})$ remains unchanged as we go from frame $S$ to $S^{'}$ .

Figure 2 shows the RMS amplitude of the magnetic field perturbations, i.e., $δ B_{rms} / B_{0}$ at different time steps for the 3D KAW simulation. As this is an initial perturbation simulation, the turbulence it produces decays, and the fluctuation amplitude decreases with time. At $ω_{p, i} t = 0$ , the RMS amplitude is 0.339, and it reaches a value of around 0.2 at $ω_{p, i} t = 1700$ . The Alfvén crossing time is $250 ω_{p, i}^{- 1}$ . We see that the amplitude reduces by more than 40% in a few Alfvén crossing times. We analyze the development of turbulence over this time period.

Figure 2

Figure 2. $δ B_{rms} / B_{0}$ for 3D KAW simulation at different time steps. A2 simulation parameters with S0 input waves are used to run the simulation.

3.1 Magnetic perturbation ratios

The ratio $| δ B_{x} / δ B_{y} |$ is calculated from simulation data. The magnetic field data $B (r)$ at a given time step is Fourier transformed to $\hat{B} (k)$ and then the ratio of the Fourier amplitudes is calculated at specific wavevectors $k$ . These ratios are shown in Figure 3 for both the simulations initialized with KAWs and whistlers. The S0 wave input is used for KAW simulations, and the S1 wave input is used for whistler simulations. The results of $| δ B_{x} / δ B_{y} |$ from the analytical calculations of two-fluid KAW, and two-fluid whistler modes are also shown in this figure. Figures 3A–C shows the variation of ratio $| δ B_{x} / δ B_{y} |$ with perpendicular wavevector along $k_{x}$ (with $k_{y} = 0$ and $n_{z} = 1$ ). Figures 3D–F show this ratio along $k_{y}$ (with $k_{x} = 0$ and $n_{z} = 1$ ).

Figure 3

Figure 3. $| δ B_{x} | / | δ B_{y} |$ ratios from analytic KAW relations (black), analytic whistler (dotted blue), KAW simulation (red) and whistler simulation (green) (A) at $t = 0$ (B) at $t = 800 ω_{p, i}^{- 1}$ and (C) at $t = 1600 ω_{p, i}^{- 1}$ along $k_{y} = 0$ and (D–F) along $k_{x} = 0$ at their respective times are shown. To run the simulation A2 parameters are used with S0 wave input for KAW simulation and S1 wave input for whistler simulation.

At $t = 0$ , Figure 3D shows that the observed ratio from the KAW simulation matches with the analytical expectation at $k_{y} d_{i} = 0.5$ . Similarly, the observed ratio from whistler simulation matches with the analytical expectation for whistlers at $k_{y} d_{i} = 0.5$ . This is because the input wavevectors in both simulations include the wavenumber (0, 2, 1) which corresponds to this wavevector. On the other hand, the larger wavenumbers are not initialized with any specific wave mode. At $t = 0$ , Figure 3A shows that for these higher wavenumbers, the KAW simulation shows polarization ratios quite close to the whistler expectation, while the whistler simulation shows fluctuations falling between the KAW and whistler expectation values. However, as time progresses to $ω_{p, i} t = 800$ and $ω_{p, i} t = 1600$ , Figures 3B, C show that the polarization ratios from both the KAW as well as whistler simulation show values quite close to the KAW expectation. Similarly, in Figure 3D, we see that at $t = 0$ , the higher wavenumbers in the KAW simulation show polarization ratios closer to the KAW expectation, whereas the whistler simulation lies closer to the whistler expectation. However, as time advances, the higher wavenumbers in both KAW and whistler simulations show ratios closer to the KAW expectation (Figures 3E, F). This shows that the magnetic polarizations of modes excited non-linearly at smaller scales tend to behave more like KAWs. We have also looked at this ratio along the $k_{x} = k_{y}$ and $k_{x} = 2 k_{y}$ line of wavevectors. For $k_{x} = k_{y}$ , the KAW and whistler modes have the same analytical expectations, so they do not distinguish the modes. Along $k_{x} = 2 k_{y}$ , also we find the ratios in simulations match closely with KAW. However, the gap between KAWs and whistlers is largest for the cases shown in Figure 3.

We also look at the ratio $| δ B_{x} / δ B_{z} |$ from the simulation data. These results are obtained by the same process of Fourier transforming and then taking the ratio as done for $| δ B_{x} / δ B_{y} |$ . Figure 4 shows these ratios for both simulations as well as their analytically expected values from the two-fluid KAW and whistler models. It shows the ratios as a function of $k_{y}$ (with $k_{x} = 0$ and $n_{z} = 1$ ). Since the initial wavenumbers include $n_{x} = 0$ , $n_{y} = 2$ , and $n_{z} = 1$ which corresponds to the first point $(k_{y} = 0.5)$ , there the simulation ratios match with their analytical counterparts at $t = 0$ . Figure 4A shows that initially, at the higher wavenumbers, the ratios from KAW simulation are similar to the analytical expectation for KAWs, and the whistler simulation lies closer to the whistler analytical expectation. At higher time steps Figures 4B, C it can be seen that the $| δ B_{x} / δ B_{z} |$ both from the whistler and KAW simulations shows a closer match with analytical KAW than to analytical whistler, indicating the cascade of KAW to smaller scales in both cases. The ratio $| δ B_{x} / δ B_{z} |$ along $k_{x}$ (with $k_{y} = 0$ and $n_{z} = 1$ ), gives the same results for analytical expectations of both KAWs and whistlers since $δ B_{x} / δ B_{z} = - k_{z} / k_{x}$ for this case. Both the simulations also lie close to this analytical expectation. Hence it does not help to distinguish between KAW and whistler and is not shown.

Figure 4

Figure 4. $| δ B_{x} | / | δ B_{z} |$ ratios from analytic KAW relations (black), analytic whistler (dotted blue), KAW simulation (red) and whistler simulation (green) at (A) $t = 0$ , (B) $t = 800 ω_{p, i}^{- 1}$ and (C) $t = 1600 ω_{p, i}^{- 1}$ are shown. Here $n_{z} = 1$ , $k_{x} = 0$ and $k_{y}$ is varied. A2 simulation parameters are used.

The ratio $| δ B_{z} / δ B | = | δ {\hat{B}}_{z} (k) | / \sqrt{δ {\hat{B}}_{x}^{2} (k) + δ {\hat{B}}_{y}^{2} (k) + δ {\hat{B}}_{z}^{2} (k)}$ is also obtained from the simulation data. This is related to the compressibility, $C_{‖} = δ B_{‖}^{2} / δ B^{2}$ , which is measured in solar wind observations and thought to be related to KAW fluctutations (Kiyani et al., 2013). Figure 5 shows these ratios from the whistler and KAW simulations as well as it is analytic expectation for KAWs and whistlers obtained from their 2-fluid models. The top row shows the variation of $| δ B_{z} / δ B |$ along $k_{x}$ (with $k_{y} = 0$ and $n_{z} = 1$ ) and second row shows it along $k_{y}$ (with $k_{x} = 0$ and $n_{z} = 1$ ). As before Figure 5D shows the simulations matching with analytical expectations at $k_{y} d_{i} = 0.5$ . At $t = 0$ , along the $k_{x}$ direction, Figure 5A, the KAW simulation shows ratios closer to the whistler expectation, while the whistler simulation fluctuates between the two modes. As the simulations evolve in time, Figures 5B, C, we see that the KAW simulation lies closer to the KAW analytical expectation at higher wavenumbers. The whistler simulation is closer to the whistler analytical expectation up to $k_{x} d_{i} = 1$ , but beyond that, it also appears to go towards the KAW expectation. Along $k_{y}$ direction, Figures 5E, F, the KAW simulation again seems closer to the analytical KAW prediction. At $t = 800 ω_{p, i}^{- 1}$ the whistler simulation shows ratios matching with KAW predictions, while at $t = 1600 ω_{p, i}^{- 1}$ it shows a fluctuating behavior between the 2 modes. This could mean that some fraction of whistler modes may be present in the simulation. However, the difference between the analytical value of $δ B_{z} / δ B$ for KAW and whistler is very small (between 0.6 and 1), while for $δ B_{x} / δ B_{y}$ and $δ B_{x} / δ B_{z}$ the difference was of an order of magnitude. Hence it is difficult to distinguish the two modes using $δ B_{z} / δ B$ .

Figure 5

Figure 5. $| δ B_{z} | / | δ B |$ ratios from analytic KAW relations (black), whistler (dotted blue), KAW simulation (red) and whistler simulation (green) at (A) $t = 0$ , (B) $t = 800 ω_{p, i}^{- 1}$ and (C) $t = 1600 ω_{p, i}^{- 1}$ along k_y = 0, and (D–F) along k_x = 0 at their respective times are shown, for A2 simulation parameters.

We have also run these simulations with 125 particles per cell per species, and we get similar results. This indicates that noise is not affecting our results. From the spectra, we can see that noise becomes important at $k_{⊥} d_{i} > 4$ for the velocity fields (and further smaller scales for magnetic fields), but we are comparing polarization ratios in the region $k_{⊥} d_{i} < 3$ . We have also compared these results for different plasma $β$ mentioned in Table 1 and found similar results, i.e., in this case also we find the KAWs as the dominant mode for these beta values. However, it is possible that other modes may be important at other plasma beta. We initialize simulations with either KAWs or whistlers at large scales and allow modes to be excited at smaller scales. This can be seen from the energy spectra in Section 4.1. The ratios $| δ B_{x} / δ B_{y} |$ , $| δ B_{x} / δ B_{z} |$ and $| δ B_{z} / δ B |$ show that the modes excited at these smaller scales tend to be in agreement with the ratios expected of KAWs, even in the simulations initialized with whistler modes. This indicates a natural preference for the excitation of a KAW cascade in the sub-ion range. Typically “energy cascade” refers to the local, nonlinear energy transfer from large to small scales in the inertial range over a vast range of scales. While the kinetic range of scales in a plasma is quite limited and even more so in simulations, there also “energy cascade” is often used to describe the transfer of energy (Franci et al., 2017; Passot et al., 2022). Here, by “cascade”, we are referring to the transfer of energy from large scales to smaller scales (as seen in the energy spectra) and the excitation of wave modes at smaller scales (as seen in the polarization relations which match with KAW). This is referred to as the KAW cascade. Earlier we had shown in 2.5D simulations the locality of KAW nonlinear energy transfer (Makwana et al., 2023) by using bispectral analysis. Next we take a further look at what relations the density and velocity perturbations follow.

3.2 Number density perturbation

The ratio of density perturbation to mean number density for ions along $k_{y}$ (with $k_{x} = 0$ and $n_{z} = 1$ ) is obtained from simulation data. We take the Fourier transform of the density field and then take the ratio of the Fourier amplitude at a given $k$ , $\hat{n} (k)$ with the Fourier amplitude at $k = 0$ , since that would correspond to the mean number density. We use A2 simulation parameters to run the simulations. We also obtain the ratio from the analytical KAW 2-fluid model. The results are shown in Figure 6. There is no density perturbation in the case of the whistler waves, so analytical whistler expectations are not plotted here. Figure 6A shows that at time t = 0, the ratios from analytical KAW relations and KAW simulation are quite close at the input wavenumber $k_{y} d_{i} = 0.5$ , but they do not match at higher wavenumbers, with the higher wavenumbers showing ratios lower by more than 3 orders of magnitude. The whistler simulation shows density perturbations lower by more than five orders of magnitude compared to KAW expectations. As the KAW simulation evolves, Figures 6B, C, these density perturbations grow. In some cases, like at $ω_{p, i} t = 800$ for $k_{y} d_{i} = 0.75 - 1.75$ , the ratio has increased by a factor of around 100 compared to $t = 0$ . Although it does not reach the level of KAW expectation, it tends to get closer to it. It maybe that particle noise might affect the density perturbations and bring them to a lower level. On the other hand, we know that whistler waves have no density perturbation, but we are still getting density perturbation for the whistler simulation at higher times. We see that the density perturbation for whistler simulation increases from the order of $1 0^{- 7}$ at $ω_{p, i} t = 0$ to $1 0^{- 3}$ at $ω_{p, i} t = 1600$ which is as much as the level seen in the KAW simulation. Figures 6D–F show the ratio of ion-density perturbation to the parallel magnetic field at $t = 0$ , 800 and $1600 ω_{p, i}^{- 1}$ respectively along $k_{y} (d_{i})$ . We see that at $ω_{p, i} t = 0$ , the KAW simulation results match with analytic KAW, but the whistler simulation has a very low value. At $ω_{p, i} t = 800$ , both the simulation results are close to analytic KAW ratios, and at $ω_{p, i} t = 1600$ , both the simulation results match with analytic KAW at higher wavenumbers. We see that the whistler simulation ratio in this case rises from 0.1 at $t = 0$ to $30 - 100$ at $ω_{p, i} t = 1600$ . This shows that KAW-like density perturbations are being produced even in the whistler simulation. We also find almost the same spectra for electron density fluctuations showing that the quasineutrality condition works well.

Figure 6

Figure 6. The ratio of ion density perturbation to the mean ion density at time (A) $t = 0$ , (B) 800 and (C) $1600 ω_{p, i}^{- 1}$ with varying $k_{y}$ values. (D–F) shows the ratio of ion-density perturbation to parallel magnetic field at $t = 0$ , 800 and $1600 ω_{p, i}^{- 1}$ , respectively. Green and red curves are obtained from whistler and KAW simulation, respectively, and black from analytic KAW. A2 simulation parameters are used to run the simulation.

3.3 Velocity perturbation ratios

We now compare the velocity perturbation ratio $| δ v_{x} / δ v_{y} |$ for the analytical whistler wave, analytical KAW, and the simulations of KAWs and whistlers. The simulation data provides the current density for both species. The fluid velocity perturbation of a species is obtained by taking the ratio of its current density to its charge density and then taking the Fourier transform. The ratio of these Fourier-transformed velocities is then taken to obtain the required ratios. Figure 7 shows these different ratios at the three different time steps. The simulation results are obtained using the A2 simulation parameters. At $t = 0$ in Figure 7A along $k_{x}$ , both the KAW and whistler simulations lie in between their analytical expectations. In Figure 7D, along $k_{y}$ at $t = 0$ , the whistler simulation results are closer to analytical whistler but KAW simulation does not match with analytical KAW, except at the input wavevector $k_{y} d_{i} = 0.5$ . Along the $k_{x}$ direction, as the simulation evolves Figure 7B shows that at time $t = 800 ω_{p, i}^{- 1}$ the results from both KAW and whistler simulations seem closer to the whistler expectation at lower wavenumbers, while at higher wavenumbers they lie in between the two modes’ expectations. However, at $t = 1600 ω_{p, i}^{- 1}$ , Figure 7C shows that both the whistler and KAW simulations show ratios closer to the KAW expectation at higher wavenumbers. Along $k_{y}$ direction as well, at $t = 800 ω_{p, i}^{- 1}$ Figure 7E shows that the simulation ratios of both KAW and whistlers lie between the analytical expectations. Then, later in Figure 7F, the simulated ratios seem closer to the KAW expectation at higher wavenumbers. After looking at all these characteristic polarization ratios from simulations, it seems that the KAW modes are preferentially excited in the sub-ion range cascade.

Figure 7

Figure 7. The ratio of velocity perturbation $| δ v_{x} | / | δ v_{y} |$ at time t = 0 (A, D), t=800ω_P,i^-1 (B, E) and t=1600ω_P,i^-1 (C, F) (from left to right) with varying $k_{x}$ ( $k_{y}$ ) values in top (bottom) row from KAW simulation (red), whistler simulation (green), analytic KAW (black) and for analytic whistler (dotted green). A2 simulation parameters are used.

4 Energy spectrum

4.1 Magnetic energy and velocity spectrum

Now we look at the spectra that are produced by the KAW simulations. We have calculated the magnetic energy and velocity spectrum along the perpendicular wavevector $(k_{⊥})$ for the KAW simulation. To obtain the magnetic energy 2D power spectrum we use the discrete Fourier transform of magnetic field components, $| \hat{B} (k_{⊥}, k_{z}) |^{2} = | {\hat{B}}_{x} (k_{⊥}, k_{z}) |^{2} + | {\hat{B}}_{y} (k_{⊥}, k_{z}) |^{2} + | {\hat{B}}_{z} (k_{⊥}, k_{z}) |^{2}$ . The spectrum along the perpendicular direction is the spectrum in the $k_{x} - k_{y}$ plane. To make the spectrum along the $k_{⊥}$ direction, we take an average of energy spectrum over a ring in the $k_{x} - k_{y}$ plane for a fixed value of $k_{⊥}$ . Then the spectrum along the $k_{⊥}$ direction can be obtained by taking different values of $k_{⊥}$ and applying the same averaging operation. The perpendicular wavenumber spectrum is obtained by summing over all the parallel wavenumbers, i.e., along the wavevector $k_{z}$ . Figure 8A shows the magnetic spectrum in the perpendicular wavenumber obtained from the simulation with 16 waves injection using A2 simulation parameters. The spectrum is shown at $ω_{p, i} t = 0,300,500,900,1300,1700$ with $48 0^{3}$ resolution for S0 input wave. The dotted lines are made to see the visual comparison of the simulation slope with these reference slopes. At t = 0 the spectrum peaks at the input wavenumbers and drops steeply for $k_{⊥} d_{i} > 1$ . As the time increases, the higher wavenumber modes gain energy, indicating energy cascade to smaller scales. As can be seen in the figure at $t = 500 ω_{p, i}^{- 1}$ the higher wavenumber modes get more energy than at $t = 300 ω_{p, i}^{- 1}$ and similar behavior can be seen for higher wavenumbers. At $t = 1700 ω_{p, i}^{- 1}$ the slope becomes saturated. We find a slope of $- 3.02$ up to $k_{⊥} d_{i} \approx 3$ . Beyond that, there appears to be a break in the slope. However, for simulation parameters A3 and A4 that use a mass ratio of 50 (implying smaller electron scale) the KAW spectrum appears continuous up to $k_{⊥} d_{i} \sim 6$ (Figure 9). So, for the A2 simulation, the break could be due to some electron scale physics.

Figure 8

Figure 8. The (A) magnetic energy and (B) velocity spectrum for KAW simulation, and (C) magnetic energy and (D) velocity spectrum for whistler simulation at different times with respect to the perpendicular wavevector $k_{⊥} d_{i}$ . A2 simulation parameters are used with S0 wave input for KAW and S1 wave input for whistler simulation. Vertical dashed line is $k_{⊥} d_{e} = 1$ .

Figure 9

Figure 9. The magnetic energy spectrum $| \hat{B} (k) |^{2}$ at different times with respect to the perpendicular wavevector $k_{⊥} d_{i}$ . (A) is for A3 and (B) A4 simulation parameters.

Previously Makwana et al. (2023), we have shown perpendicular magnetic energy spectrum obtained by 2.5D PIC simulations, and the spectrum is similar to 3D case but less smooth and clear. The reason can be given to the fact that in the 3D case, for obtaining the perpendicular spectrum, we have to average over a ring on the perpendicular wavevector plane, while for the 2D case, there is no averaging. Huang et al. (2020) have also obtained the magnetic power spectrum for the solar wind with data taken from NASA’s Parker Solar Probe. They have obtained the power spectrum for both the parallel and the perpendicular direction and have got different slopes. For the perpendicular magnetic spectrum in the transition range (also called the sub-ion range), they have obtained a power law slope of −3.73. There is a large variation in the spectral slope in the sub-ion range (Sahraoui et al., 2009; Sahraoui et al., 2010; Bruno and Trenchi, 2014; Bruno et al., 2014) depending on various physical conditions of the solar wind. The spectral slopes measured in our simulations also lie within this range.

The velocity field spectrum in the perpendicular direction is also obtained and its behaviour is similar to the magnetic spectrum. Initially, at $ω_{p, i} t = 0$ , all the energy is in the initial wavenumbers and very little energy in higher wavenumbers. As the simulation evolves in time, the higher wavenumbers start to get energy, as can be seen in Figure 8B. As can be seen from the figure, the slope at higher times saturates from $ω_{p, i} t = 500$ to 1,700, giving an index of around $- 4.4$ . We see the noise in velocity spectrum becomes dominant after $k_{⊥} d_{i} \sim 4$ . However, we are analyzing the polarization ratios in the sub-ion range up to $k_{⊥} d_{i} \sim 3$ . Also, we have run 125ppc simulation, and we get similar polarization ratio graphs. So our results are not affected by the particle noise. Figure 8C shows the magnetic spectrum for the whistler simulation, which is steeper than the KAW simulation, and Figure 8D shows the velocity spectrum for the whistler simulation.

We have also run KAW simulations with different mass ratios and different ion and electron thermal velocities. The magnetic power spectra in the perpendicular direction for these simulations are shown in Figure 9. Figure 9A shows the magnetic spectrum with a ion-to-electron mass ratio of 50 and A3 parameters. This results in the decrease of electron beta $β_{e}$ shown in Table 1 in A3 parameters. The magnetic spectrum in Figure 9B is the result of simulation with a mass ratio of 50 with A4 parameters. These changes lead to a decrease in both electron and ion plasma beta, the values for which are shown in Table 1, A4 parameters. From Figure 9A, we observe a spectral index close to $- 3$ in the saturated state, which is similar to that observed in Figure 8A. However, now we get a continuous power-law spectrum down to smaller scales and do not see a break. This extension of spectra to smaller scale may be because electron scales become smaller with increased mass ratio. It is known that in the sub-ion range the spectral slope is quite variable while it settles to around −2.8 value as we go deeper towards electron scales Alexandrova et al. (2013) In Figure 9B, the steepness of the slope further decreases, and we observe a spectral index of around $- 2.8$ in the saturated state, which is the index sometimes observed in the solar wind turbulence dissipation range (Verscharen et al., 2019). This shows that the spectral index could depend on the plasma beta parameter. There could be various reasons for this decrease in steepness like higher intermittency at ion scales at higher $β$ (Parashar et al., 2018) and that could cause a steeper slope at higher beta (Boldyrev and Perez, 2012).

4.2 2D magnetic energy power spectrum

To understand these interactions in more detail we obtain the 2-dimensional magnetic field power spectrum for KAW simulation in the $k_{⊥} - k_{z}$ plane. Figure 10 shows the 2D magnetic field power spectrum at $t = 0,800$ and $1600 ω_{p, i}^{- 1}$ (from left to right, respectively) for A2 simulation parameters for KAW (top row) and for Whistler simulation (bottom row). As can be seen in Figure 10, at $t = 0$ , all the magnetic energy is in the initial wavenumbers, but as the simulation evolves in time, the energy starts to spread to higher wavenumbers. The energy is spread more in the perpendicular direction than in parallel direction, showing the typical anisotropic energy cascade. We get similar results for whistler simulation but there appears to be lesser anisotropy compared to the KAW simulation case. In Makwana et al. (2023), even stronger anisotropy was observed with the energy concentrated mainly along $k_{z} = 0$ , whereas here, it is not so prominent. The reason might be the partial information in 2.5D simulations where we had taken $k_{x} = 0$ , and the spectrum was only w.r.t. $k_{y}$ . This means that some nonlinear interactions would be suppressed, and it was not able to show the behavior of the full 3D spectrum. In the 3D simulations, all three-wave interactions are accounted for. In the 2.5D simulations, we saw the non-resonant 3-wave interactions responsible for the energy spread to higher wavenumbers. Using bispectral analysis we showed that these interactions were local. We suspect that similar 3-wave interactions are playing an important role here. Another non-linear energy transfer possibility is magnetic reconnection-induced energy transfer from reconnection scales back to larger scales, as in Franci et al. (2017). Other possibilities could be nonlocal transfers from large to small scales, by shear (Gorman and Klein, 2024) or other effects (Friedrich et al., 2024). One unknown point is the observation of modes at the high $k_{z} d_{i} \sim 2 - 3$ in KAW simulation case. This region is $k_{⊥} \sim k_{‖}$ , so these may not be KAWs. These were observed in 2.5D simulations also, and they could be ion-cyclotron waves, or some other modes. The identity of these modes and their generation mechanisms will be studied in future work.

Figure 10

Figure 10. The logarithm of 2 dimensional magnetic energy spectrum $l n | \hat{B} (k_{⊥}, k_{z}) |^{2}$ at time (A) $t = 0$ , (B) $800 ω_{p, i}^{- 1}$ and (C) $1600 ω_{p, i}^{- 1}$ (from left to right) in $k_{⊥} - k_{z}$ plane for KAW simulation. And the logarithm of 2-dimensional magnetic energy spectrum for whistler simulation at the time (D) $t = 0$ , (E) $800 ω_{p, i}^{- 1}$ and (F) $1600 ω_{p, i}^{- 1}$ . A2 simulation parameters are used to run the simulation.

5 Conclusion

In this study, we found that the two-fluid KAW eigenvector relations work well as an initial condition in a fully kinetic 3D PIC simulation. Previously, in Makwana et al. (2023), we did 2.5D PIC simulations and showed that the two-fluid eigenvector relations gave a linear dispersion relation close to the hot plasma kinetic dispersion relation. In this work we extend that study to full 3D PIC simulations. Comparison of the analytical two-fluid KAW dispersion relation, 2.5D simulations, and 3D simulations shows that the frame rotation works appropriately in order to set up initial conditions for different wave-mode simulations.

The simulations are initialized at small wavenumbers (large scales) by either KAWs or whistlers. Therefore, initially, the magnetic perturbation ratios $| δ B_{x} / δ B_{y} |$ , $| δ B_{x} / δ B_{z} |$ , and $| δ B_{z} / δ B |$ at the input wavenumbers match with the analytical expectations of those modes. The higher wavenumbers (smaller scales) are not initialized with any specific mode so the perturbation ratios for those do not show any correspondence with the KAW or whistler modes. As the simulation evolves and nonlinear interactions occur we see that energy is transferred and modes are excited at higher wavenumbers. The simulations initialized with KAWs, as well as those initialized with whistlers, show excitation of modes that show perturbation ratios closer to those expected from KAW. Similar results are seen in the velocity and density perturbation ratios. These results are more in line with the KAWs as compared to the whistler wave. In the KAW simulations, this can be explained as a cascade of KAWs in the simulation with the evolution of time. Their presence in simulations initialized with whistlers shows that we can expect KAWs to be the dominant wave-modes at these scales. For the generation of KAW modes in the whistler simulation, there could be several possibilities. There is the conversion of incident wave mode to KAW at the resonant layer in tokamaks (Hasegawa and Chen, 1976), parametric decay of whistler modes can lead to the generation of KAWs (Sharma et al., 1986) or other non-linear couplings of whistlers could also generate KAWs (AC-L, 1995). KAWs are also generated by phase mixing (Vásconez et al., 2015). What processes are responsible for the generation of KAWs in these simulations needs to be investigated.

The perpendicular magnetic energy spectrum obtained for the KAW simulation by running the simulation with $48 0^{3}$ resolution shows different behavior from time $t = 0$ to $1700 ω_{p, i}^{- 1}$ . We see that at $t = 0$ , all the energy is contained in the injection wavenumbers. As the simulation evolves in time, the higher wavenumbers get energy. Although the turbulence is decaying in these simulations, the magnetic spectrum comes close to a steady slope of index $- 3$ at later times, and the velocity spectrum reaches a steady slope close to $- 4.4$ . This implies that the non-linear interactions in the simulation produce a cascade of KAWs, resulting in this spectrum. Huang et al. (2020) and Huang et al. (2021) have obtained the energy spectrum in parallel and perpendicular directions in the transition (sub-ion) range from the data of NASA’s Parker Solar Probe. In their study, they found magnetic spectra with slopes ranging from $- 3$ to $- 5.7$ , which was also attributed to KAWs. Our results of the KAW cascade show similar behavior and are able to produce the magnetic energy perpendicular spectrum of solar wind in the ion-transition range. They support the existence of KAWs in the ion-transition range of solar wind turbulence. By increasing the ion-to-electron mass ratio, we see extension of this sub-ion range power law to smaller scales. When we decrease plasma $β$ we observe slightly shallower power-law indices. Earlier studies have also shown that the spectral indices can vary with plasma $β$ (Franci et al., 2016), although the behavior may be complicated.

A more detailed picture of the non-linear interactions of KAWs in the simulation is shown by the 2D magnetic energy spectrum in the $k_{⊥} - k_{z}$ plane. The spectrum shows that there is a spread of power to higher wavenumbers in both the parallel and perpendicular direction, but the spread is anisotropically more in the perpendicular than in the parallel direction. By comparing these results to our previously obtained 2.5D simulation results (Makwana et al., 2023), we find that the energy spread in the perpendicular direction in 3D simulations is not as strongly anisotropic as in the 2.5D case. The reason could be the absence of some of the non-linear interactions in the case of 2.5D simulations, as the waves are confined only to the 2D $k_{y} - k_{z}$ plane.

There have been earlier studies showing nonlinear generation of whistlers by pump KAWs (Dwivedi et al., 2012) and also excitation of low-frequency KAWs by whistlers (Goyal et al., 2017). There is also the possibility of the co-existence of whistler modes and KAWs (Mithaiwala et al., 2012). However, this generally happens closer to the electron scales. The cascade of whistler waves has also been simulated in PIC at close to the electron scales (Gary et al., 2008). In this study we find that at the sub-ion scales, KAWs are preferentially excited compared to whistlers. It needs to be studied whether it is simply conversion of the Alfvén cascade at magnetohydrodynamic scales into a KAW cascade (Xiang et al., 2019) and/or is it also excitation of the KAWs by some instabilities. One drawback of this study is that it simulates decaying turbulence. More general forced simulations will provide better characteristics of turbulence. Also, there could be other modes playing an important role in addition to the two modes studied here. To study further smaller scales, we would need to use kinetic dispersion properties. Furthermore, the resolution is quite limited. We require higher resolution in both space and time, with higher mass ratios, to obtain more realistic kinetic range turbulence properties.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

JS: Data curation, Formal Analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing–original draft, Writing–review and editing. KM: Conceptualization, Funding acquisition, Investigation, Supervision, Validation, Writing–review and editing, Methodology, Resources, Formal Analysis, Writing–original draft.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by Science and Engineering Research Board (SERB) grant SRG/2021/001439. We would like to thank DST – FIST (SR/FST/PSI-215/2016) for the financial support alongwith Indian Institute of Technology, Hyderabad seed grant. The support and the resources provided by the “PARAM Sanganak facility” under the National Supercomputing Mission, Government of India at the Indian Institute of Technology, Kanpur, as well as the “PARAM Seva facility” at Indian Institute of Technology, Hyderabad are gratefully acknowledged.

Acknowledgments

KM acknowledges useful discussions with Dr. M. K. Verma.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fspas.2024.1423642/full#supplementary-material

References

AC-L, C. (1995). Nonlinear alfvén wave phenomena in the planetary magnetosphere. doi:10.1088/0031-8949/1995/T60/005