Introduction

Turbulent fluctuations are a likely mechanism of heating of the solar corona. However, the mode composition of the turbulence is not fully understood. It is usually thought to be Alfvén waves or Alfvénic structures that have some polarization properties of the Alfvén waves but do not obey any dispersion relations. For the turbulent fluctuations to be a source of the coronal heating, there has to be a kinetic mechanism of their dissipation. Furthermore, this mechanism has to channel their energy to ions and the ions have to be energized preferentially perpendicularly to the magnetic filed. The most straightforward damping channel to achieve that would be the cyclotron resonance. However, it appears that the Alfvénic turbulence consists mostly of fluctuations whose frequency is low compared to the cyclotron frequency. A recent review of the Alfvénic turbulence theories and their observational predictions can be found in Zank et al. (2021). In the absence of the cyclotron resonance, a suitable dissipation channel is based on demagnetization of the ions by kineticscale magnetic structures produced by the turbulence (Dmitruk et al. 2004;Dmitruk & Matthaeus 2006;Parashar et al. 2009;Markovskii & Vasquez 2011;Vasquez & Markovskii 2012;Vasquez 2015).

Alternatively, there is both theoretical and observational evidence that a fast-mode component may be present in the solar wind and coronal turbulence spectra, along with the Alfvénic fluctuations. If the fast modes propagate at a nonzero angle to the magnetic field, they can interact with the ions via the cyclotron resonance (Li & Habbal 2001;Hollweg & Markovskii 2002). This makes them an efficient means of coronal heating (Markovskii et al. 2010). The Hall magnetohydrodynamic (MHD) simulations of Ghosh et al. (1996) showed that the fast modes can be generated as a result of a turbulent cascade initiated at scales much larger than the ion kinetic scales. The fast waves can be also produced in the corona in the single-fluid MHD regime. For example, coronal density fluctuations can convert Alfvén waves into fast waves, and vice versa, through resonant three-wave interactions (Kuznetsov 2001;Chandran et al. 2009).

One of the ways to discriminate the fast/whistler and Alfvén waves in the observational data is to determine their magnetic compressibility at the ion kinetic scales. The magnetic compressibility is a measure of the magnetic fluctuation in the component parallel to the background magnetic field relative to the total fluctuation. If the electron plasma beta is small, the magnetic compressibility associated with highly oblique whistlers is much larger than that of kinetic Alfvén waves. Using the data reported by Hamilton et al. (2008) at 1 au, Gary & Smith (2009) have shown that, at least during solar wind intervals with certain properties, the observed values of the magnetic compressibility cannot be generated by the kinetic Alfvén waves and must correspond to the highly oblique whistlers. Chaston et al. (2020), Zhu et al. (2020), andZhao et al. (2021) decomposed the solar wind fluctuations close to the Sun, around several tens of solar radii, into Alfvén, fast, and slow magnetosonic modes. Chaston et al. (2020) and Zhu et al. (2020) have shown that the fast modes constitute a considerable fraction of the power spectral density, which depends on the scale but is roughly 20%. Zhao et al. (2021) also found that, at the MHD scales, the fast waves provide the second-largest contribution to the spectral energy density. Zhu et al. (2020) have further demonstrated that the fast-mode fraction increases with decreasing heliocentric distance. Therefore, one can hypothesize that it can be even higher in the inner corona, at several solar radii.

The numerical simulations of Markovskii & Vasquez (2010) have shown that the protons are heated preferentially perpendicularly to the magnetic field when the fast-mode turbulence is imbalanced. This means that the power spectral densities of the waves propagating in the direction of the magnetic field and in the opposite direction are unequal. The problem of the mode identification in the observational data is not fully resolved, and it was recently revisited by Zank et al. (2023). In particular, there is more than one way to infer the mode propagation direction. However, Chaston et al. (2020), Zhu et al. (2020), andZhao et al. (2021) all agree that the observed fast-mode fluctuations are imbalanced.

The proton heating by the fast-mode turbulence has been studied with the help of hybrid numerical simulations (particle ions and fluid electrons) in a two-dimensional (2D) configuration (Svidzinski et al. 2009;Markovskii & Vasquez 2010;Markovskii et al. 2010). In a quasi-steady state, the heating rate is approximately equal to the rate of the wave energy cascade. This latter rate is not necessarily the same in two and three dimensions (e.g., Zakharov et al. 1992). The proton heating in the three-dimensional (3D) configuration has been investigated using particle-in-cell simulations along with the electron heating (Hughes et al. 2014;Gary et al. 2016). However, due to the high computational cost of these simulations, they have not covered the entire wavenumber range where the damping of the fast modes on the protons is significant.

In this paper, we perform 3D hybrid simulations of the ion heating by the fast-mode turbulence in a plasma typical of the solar corona. We investigate if the heating is efficient enough to produce the solar wind acceleration, taking into account the turbulence imbalance. We also include alpha particles in our simulation. Despite being minor ions, they can still play a significant role in the overall energy partition.

Numerical Setup

We perform 3D hybrid simulations with particle-in-cell protons and a quasi-neutralizing electron fluid. The numerical code is described by Terasawa et al. (1986) and Vasquez (1995Vasquez ( , 2015)). The code solves the following equations:

where n e = n p + 2n α . The quantities x p,α and v p,α are positions and velocities of individual protons and alpha particles, q p,α is the ion charge, E is the electric field, c is the speed of light, e is the elementary charge, n e is the electron number density, and V e is the electron fluid velocity. The proton and alpha number density n p,α and bulk velocity V p,α for each spatial cell are calculated as moments of the distribution. The mean magnetic field B 0 is in the positive x direction. The electron temperature is set to zero, and the electron fluid is massless. The simulation box sizes are L = 100 in all directions in units of the proton inertial length W - V .

A p

1 The number of cells in the simulation grid N = 128 in all directions. The time step is 0.01 in units of the inverse proton gyrofrequency W -. p 1 Here, Ω p and the Alfvén speed V A , and thereby the spatial and temporal scales, are defined with the initial mean values of the magnetic field B 0 and the proton number density n 0p . The boundary conditions are periodic.

We initiate the turbulence in the following way. The magnetic field fluctuations ΔB(t, x) at t = 0 are given by the formula

where x is the Cartesian spatial position vector, k is the wavevector, and f(k) is a random phase. The seed spectrum is confined to the modes k 0.38 Ω p /V A and taken to be flat for simplicity. We consider the cases where the waves propagate in just the positive direction along B 0 and in both positive and negative directions. For consistency, we exclude the modes propagating exactly across B 0 because their wavevector projections on B 0 are neither positive nor negative. The fluctuations in the solar corona (and in the distant solar wind) are small at the proton kinetic scales. Since the energy cascade and plasma heating depend on the fluctuation amplitude, it is important to keep the amplitude within the observational constraints to make the simulations relevant. As a reference point, we use the interplanetary scintillation measurements at 5 solar radii (Coles & Harmon 1989;Coles et al. 1995). We can estimate that the total rms of the number density fluctuations in the spectral region under consideration should be a few percent (Markovskii & Hollweg 2002). We take the total initial rms of the magnetic field fluctuation ΔB tot (0) = 0.05 B 0 , where

This corresponds to the total rms of the proton number density fluctuation Δn tot p (0) = 0.037 n 0p in the 3D configuration. The proton and alpha density and bulk velocity fluctuations Δn p,α and ΔV p,α are defined in the same way as ΔB in Equation (6). The components of the vectors δB(0, k), δn(0, k), and δV p,α (0, k) obey the polarization relations of linear fast waves in the cold MHD limit and in the approximation n α = n p . The ion distribution functions are loaded as drifting bi-Maxwellians,

is the ion thermal speed, k B is the Boltzmann constant, and T is the total ion temperature. The initial ion number densities and temperatures are spatially uniform. The numbers of particles per cell are 200 and 100 for the protons and alphas, respectively. We set the proton beta β p to a reasonable value of 0.02 at 5 solar radii. This corresponds to the proton temperature T p = 3 • 10 6 K and the Alfvén speed V A = 1560 km s -1 . The alpha-particle temperature T α = 4T p , and the density n 0α = 0.05 n 0p . The electron temperature at this distance is difficult to measure directly, but there is evidence suggesting that it is much smaller than the proton temperature (e.g., Esser & Edgar 2000). Therefore, it is neglected in Equations (1)-( 5). After the initial conditions of the simulation are set, the system is allowed to evolve freely in time.

Simulation Results

A single fast wave would cascade the energy to higher wavenumbers by steepening, which generates higher harmonics of the initial mode. In this case, the energy transport would proceed in the "radial" direction, i.e., along the line connecting the origin of the coordinates and the initial mode location in the wavenumber space. The presence of multiple modes in the seed spectrum results in wave interactions and a nonradial and anisotropic transport preferentially across the mean magnetic field (Markovskii et al. 2010). This remains true in the 3D configuration, as demonstrated in Figure 1, which shows a 2D magnetic wavenumber spectrum S obtained by summing the full 3D spectrum over the azimuthal angle in the unidirectional propagation case.

The rate of the energy cascade is difficult to calculate directly. Therefore, we use the fluctuation energy dissipation rate as a proxy, assuming that in a quasi-steady state they balance each other. Figure 2 plots the total wave energy density, and otherwise the same parameters (Figure 2). The quantity ΔE w /(E w (0)Δt) estimated from the blue curve (c), in the 3D setup, is 3.6 • 10 -4 Ω p . Therefore, the dissipation rate scales as

This is the same scaling as obtained by Markovskii et al. (2010) in the 2D setup.

The damping of the fluctuations leads to the plasma heating. The time evolution of the ion temperatures is displayed in Figure 3. The heating of the protons and alphas is preferentially perpendicular to the background magnetic field and the alpha heating is stronger, which is consistent with the observed ion temperatures in the solar corona (Dodero et al. 1998;Kohl et al. 1998). The most likely mechanism of the perpendicular heating is a cyclotron damping of the oblique fast modes that can happen either directly or through coupling to ion Bernstein modes (Li & Habbal 2001;Hollweg & Markovskii 2002). The linear dispersion relations of these waves have a complex structure and depend, in particular, on the propagation angle θ (Li & Habbal 2001). In Figure 4, we plot the frequency ω and damping rate γ as functions of the wavenumber at q = ( ) arctan 2 , which is roughly the direction of the strongest cascade in Figure 2. The red branch starts at low wavenumbers as a fast mode and transforms to a Bernstein one through a (green line).

Figure 3. Time evolution of the parallel (black line) and perpendicular (red line) proton and parallel (green line) and perpendicular (blue line) alpha normalized to their initial values. In the case of the perpendicular alpha heating (blue curve), the quantity T ⊥α /T ⊥α (0) -1 is multiplied by 0.1 to fit the scale. linear mode coupling. Conversely, the blue branch starts as a Bernstein mode and transforms to a fast one.

Due to the low amplitude of the fluctuations in the simulation, it is reasonable to assume that they roughly follow the dispersion relations of the linear modes and damp at approximately the same rate. As the frequency increases from values ω = Ω p , the waves first encounter the alpha-cyclotron resonance at k ≈ 0.45 Ω p /V A . This can explain the stronger alpha heating. In Figure 3, the perpendicular alpha heating (blue curve) starts soon after the beginning of the simulation at around =W - t 40 .

p 1 The resonance at the proton-cyclotron frequency occurs at larger wavenumbers, k ≈ 0.85 Ω p /V A , and the perpendicular proton heating (red curve) is delayed further in time. By contrast, the parallel heating is most likely caused by the Landau resonance and has a more gradual onset.

To verify what modes are actually present in the simulation, we have derived a frequency-wavenumber magnetic power spectrum from the Fourier transform in space and time. The temporal transform was performed using 400 evenly spaced samples over an entire simulation period of 400Ω p in the 2D run. Figure 5 shows a 2D spectrum P as a function of ω and |k| at q = ( ) arctan 2 . The dashed lines are dispersion relations of the linear waves. The two top curves are the same branches as the ones shown in Figure 4, and the bottom curve is the Alfvén mode. As can be seen from the figure, almost all of the power resides in the fast branch. Traces of power are also seen in the Bernstein and Alfvén branches. However, due the small intensity of these modes, they are unlikely to account for the fluctuation energy drop in Figure 1 (red curve (b)).

The dispersion curves in Figure 4 are plotted for T e = 0.1 T p in qualitative agreement with the observational data of Esser & Edgar (2000). Note that the fast waves can also damp on electrons and heat them. The damping rate away from the ion cyclotron resonances is due entirely to the electron and ion Landau resonances, and it is much smaller than the cyclotron damping rate. The electron damping and heating are not taken into account by the hybrid simulations, which treat the electrons as a fluid. In a low-beta plasma, particle-in-cell simulations with particle electrons have been performed by Hughes et al. (2014). They found comparable electron and proton heating rates. However, in their case of relatively high wavenumbers, compared to our setup, the linear proton damping was weaker, while the electron damping was stronger.

From here, we conclude that neglecting the electron kinetics does not strongly affect our conclusions.

Our simulations show that the alpha particles produce a significant effect on the system evolution, despite being minor ions. The spectrum S from a simulation run without the alphas is displayed in Figure 6. It extends to larger wavenumbers compared to Figure 1, in the sense that the power is higher at a given |k|. This is because there is no alpha-cyclotron resonance in this case, and the damping happens at the proton-cyclotron resonance. In the absence of the alpha-particle Bernstein modes, the linear damping rate is described by the highwavenumber peak of the green curve (c) in Figure 4. The perpendicular heating of the protons without the alphas is about twice as high as shown in Figure 3.

The analyses of the observational data of Chaston et al. (2020), Zhu et al. (2020), andZhao et al. (2021) show that the spectral power of the fast modes propagating in one direction along the magnetic field is lower than in the other direction. However, it is not necessarily negligible. As discussed by Markovskii et al. (2010), the oppositely propagating fast waves Figure 5. Two-dimensional magnetic power spectrum P as a function of ω and |k| at q = ( ) arctan 2 derived from the simulations. The spectrum is normalized to B . change the partition of the heating. Figure 7 displays a bidirectional case where the initial seed spectrum of the turbulence is isotropic, apart from the excluded waves propagating exactly perpendicular to B 0 , as described in Section 2. We can see that the perpendicular proton heating is reduced by a factor of about 2.5 (during the last 100 gyrocycles of the simulation). It is still significant but it no longer dominates the parallel heating. By contrast, the perpendicular alpha heating remains dominant even though it is also reduced.

Let us now estimate the proton heating in the solar corona. Its characteristic time t heat = Δt T 0p /ΔT ⊥p can be very short under the observational constraints on the fluctuation amplitude discussed above. Here, ΔT ⊥p is the temperature increase over the time interval Δt between W - 300 p 1 and W - 400 p 1 in Figure 3. Using the figure and the dimensional parameters T p = 3 • 10 6 K, n p = 6000 cm -3 , and β p = 0.02, we obtain t heat = 24.5 s. This is much smaller than the characteristic solar wind expansion time R/V SW = 5800 s at 5 solar radii. Here, R is the heliocentric distance and V SW = 600 km s -1 is the solar wind speed. The corresponding heating rate, which can be estimated as

is almost 3 orders of magnitude higher than the rate successfully producing the solar wind acceleration in a fluid model (see, e.g., Esser et al. 1997).

The estimate of the fluctuation amplitude used in Section 2 is based on the assumption that the only contribution to the observed spectrum comes from the fast modes. The extent to which this assumption is valid at 5 solar radii is not known. We can take the fraction of the power spectral density of the fast modes to be 20%, as reported by Chaston et al. (2020) and Zhu et al. (2020) at several tens of solar radii. It will be an estimate from above, because this fraction decreases with the heliocentric distance. Even then, the heating rate will be sufficiently high, given the scaling of the dissipation rate with the fluctuation amplitude found in Section 3.

Conclusions

We have performed 3D hybrid simulations with particle protons and alphas and a quasi-neutralizing electron fluid. Realistic parameters of the background plasma and fluctuations at 5 solar radii were used. We have considered the waves propagating in just the positive direction along the background magnetic field and in both positive and negative directions. We have shown that the rate of the energy transfer by the cascade of the fast modes in a 3D case is almost the same as it is in a 2D case. The scaling of the rate with the fluctuation amplitude is also the same.

The damping of the waves at the kinetic scales results in the ion heating. Despite the fact that the relative concentration of the alpha particles is small, they play a significant role in the system evolution. As the energy cascade proceeds to higher wavenumbers, the fluctuations encounter the alpha particle before the proton-cyclotron resonance. The alphas are heated more strongly than the protons and always preferentially across the magnetic field. The presence of the alpha particles reduces the proton heating rate by about a factor of two.

The proton heating is preferentially perpendicular when all the waves propagate in the same direction along the magnetic field. Counterpropagating waves reduce the perpendicular heating rate by up to a factor of 2.5, where it would be no longer dominant but still significant. The rate in the corona will be further reduced compared to the simulation due to the fact that the fast modes contribute only a fraction of the total fluctuation power, along with other modes. Nevertheless, we find that, with all these factors taken into account, the ion heating by the fast magnetosonic waves is sufficient to produce the acceleration of the solar wind.