© The Authors 2025

1. Introduction

One of the major outstanding challenges in solar physics is solving the so-called coronal heating problem: explaining which mechanism(s) cause(s) the plasma of the solar corona to reach temperatures millions of degrees kelvin higher than the solar photosphere and chromosphere. This problem implies that one or more processes occurring on or around our closest star are poorly understood. In practice, this also means that we cannot accurately model the conditions in the solar corona, which further constrains our ability to model the solar dynamics and accurately forecast space weather. Current state-of-the-art global coronal models aimed at space weather forecasting rely on highly approximate, empirical models to add the necessary momentum and energy budgets to create a corona with realistic temperatures that drives a bi-modal solar wind (see e.g. Downs et al. 2010; Lionello et al. 2009; Baratashvili et al. 2024).

A variety of mechanisms have been proposed to explain coronal heating. An overview is provided in Klimchuk (2006). Generally, the mechanisms are divided into two classes: the processes related to magnetic reconnection and wave phenomena. Details regarding constraints on the heating due to magnetic reconnection and magnetohydrodynamic (MHD) waves can be found in the works of Aschwanden (2001) and Van Doorsselaere et al. (2020), among many others.

However, despite the substantial body of literature available on this topic, to date, there is no consensus in the community regarding the source of coronal heating. Many argue that trying to find a single process that causes coronal heating may be impossible, as it may instead be caused by several different mechanisms acting on a variety of spatial and temporal scales and interacting with each other (De Moortel & Browning 2015). Waves and instabilities can induce magnetic reconnection, and vice versa. Magnetic reconnection can lead to the formation of waves and instabilities. For example, Pueschel et al. (2014) showed that turbulence driven by magnetic reconnection on small scales gives rise to particle acceleration and a kinetic Alfvén wave (KAW) cascade, while van Ballegooijen et al. (2011) argued that Alfvénic turbulence is a strong candidate to explain coronal heating. Many phenomena observed in the solar corona, including the different types of waves and magnetic reconnection and their mutual interactions at various scales, likely contribute to the energy and momentum budgets of the coronal plasma simultaneously, albeit to a different extent.

When studying some of the possible heating processes discussed in the publications cited above, MHD and multi-fluid models may suffice. However, higher-fidelity models may be required for other processes, particularly those occurring on scales incompatible with the assumptions of the MHD and fluid models. One such process, suggested by Vranjes & Poedts (2009a,c, 2010a,b,c, 2014) as a possible mechanism through which coronal plasma could be effectively heated, are drift waves (DWs), which occur at the scale of the ion Larmor radius. As determined from thermonuclear fusion research, DWs are generated spontaneously by DW instabilities in regions where density and/or temperature gradients exist perpendicular to the background magnetic field. Many types of DWs exist, including the fusion-plasma-relevant ion- and electron-temperature-gradient-driven instabilities (Lee & Tang 1988), trapped-electron modes (Coppi & Rewoldt 1974), and microtearing modes (Hazeltine et al. 1975). The extreme ultraviolet and X-ray corona are highly filamentary and inhomogeneous, with abundant density gradients at a variety of scales (Woo 2007; Mackay et al. 2010; DeForest et al. 2018; Chitta et al. 2022; Antolin et al. 2023). Doppler measurements indicate that structures with short density-gradient scale lengths exist down to scales of 10 km (Woo 2006). Thus, in this study, we focused on density-gradient-driven DWs.

Since DWs rely solely on density gradients for their drive, it is more interesting to study them compared to mechanisms that require a specific excitation or footpoint driving (at the correct frequencies) to generate heating (see e.g. Gudiksen 2009); for DWs to be driven to instability, all that is required is a density gradient perpendicular to a sheared background magnetic field. It may be instructive to study DWs beyond the scope of coronal heating as they may be ubiquitous in the corona, just like in thermonuclear fusion plasmas. They could thus generate (drift-Alfvén) turbulence and drive secondary reconnection, which can then trigger or interfere with other coronal processes. Furthermore, since DWs produce drift motion perpendicular to the background magnetic field, their study can help us better understand and estimate cross-field transport in the corona, which is why they are studied so intensively in thermonuclear fusion research (as discussed by e.g. Voitenko & Goossens 2004).

However, unlike some of the other mechanisms mentioned above, DWs have not yet been directly identified in solar coronal observations, possibly due to the small scales of the fluctuations they excite. Drift-wave-type instabilities do, however, play an important role in other applications, such as in fusion devices, where they are the dominant source of DW turbulent transport (Terry & Diamond 1985; Rogister 1996) and have thus been widely researched in this context.

The instability mechanism generating the density-gradient-driven DWs studied in this paper is detailed in Brchnelova et al. (2024). In brief, a density perturbation ($\tilde{n}$) to the equilibrium density profile causes a differential response from the parallel motion of ions and electrons due to their disparate masses. This results in the generation of an electrostatic potential perturbation ($\tilde{\mathrm{\Phi }}$) and thus an electric field ($\tilde{\text{E}}=-\nabla \tilde{\mathrm{\Phi }}$) is generated. The resulting E×B drift then causes advection perpendicular to B and to the background density gradient (hence the term ‘drift waves’). If the response of the electrons is non-adiabatic, for example due to electron inertia, a phase difference is created between $\tilde{n}$ and $\tilde{\mathrm{\Phi }}$, and the DWs become linearly unstable.

Fluctuations produced by unstable DWs can lead to various heating mechanisms, which in turn lead to different parallel and perpendicular temperatures, as observed in the solar atmosphere. In the perpendicular direction, stochastic heating is predicted to occur (Vranjes & Poedts 2009a,c, 2010a,b,c, 2014) due to particle motion that becomes chaotic in the presence of electrostatic field fluctuations (Cerri et al. 2021). This heating has been studied both experimentally and numerically, for example in nuclear fusion experiments (McChesney et al. 1991), in the Earth's bow shock (Stasiewicz & Eliasson 2020, 2021), and in the ionosphere (Stasiewicz et al. 2000). In the ionosphere and magnetosphere, the DWs responsible for this stochastic energisation correspond to the lower hybrid drift instability between the cyclotron and lower hybrid frequencies. However, the properties of the solar corona are fairly different from those of these environments, and thus, the resulting DW frequencies can also differ significantly. Stochastic heating further requires a sufficient gradient of the electric field. Measurements of strong and dynamic electric fields in the solar atmosphere have indeed been reported (Davis 1977; Zhang & Smartt 1986), and Vranjes & Poedts (2009b) further argue that they should be omnipresent in the solar corona and sufficient to drive stochastic heating. Chandran et al. (2013) have shown that stochastic heating can explain various solar wind features, such as the temperature ratios of different plasma species. In the parallel direction, kinetic particle acceleration can take place due to a parallel electric field (E_∥) and an aligned electric current (j_∥; see Pueschel et al. 2014).

However, in the DW publications cited above, the properties of DWs in the solar corona have been estimated through approximate dispersion relations that rely on several simplifying assumptions and thus only allow qualitative conclusions regarding the potential of DWs for the coronal heating problem. By contrast, in Brchnelova et al. (2024), we conducted a linear gyrokinetic (Brizard & Hahm 2007) simulation study of DWs with the gyrokinetic code GENE (Jenko et al. 2000). We found that – given that linear mode signatures such as frequencies and cross-phases, especially at the lower wavenumbers that dominate turbulent spectra, tend to carry over to the non-linear system (Dannert & Jenko 2005) – DWs are expected to be unstable in the solar corona at frequencies in the range 0.1 mHz to 1 Hz, the exact frequency depending primarily on the considered length scale, temperature, and magnetic field strength

Linear simulations, however, cannot predict the amplitude of the turbulent fluctuations because saturation will occur due to the background variation. In addition, according to the quasi-linear scalings in Dannert & Jenko (2005) and Staebler et al. (2016), turbulent fluctuation amplitude spectra typically peak at much lower wavenumbers than linear growth rate spectra. We thus expanded on Brchnelova et al. (2024) by studying the non-linear system, which enabled more quantitative insights into fluctuation amplitudes, parallel and perpendicular heating, and fluctuation spectra; this in turn provides a better understanding of the effects that DWs can have on the surrounding coronal plasma.

In Sect. 2 we introduce the gyrokinetic solver GENE and the setup used for the non-linear simulations. In Sect. 3 we present analyses of non-linear simulations, including the corresponding spectral fluxes, frequency power spectra, fluctuation spectra, kinetic particle heating, stochastic timescales, and structure erosion due to particle fluxes. In Sect. 4 the findings from Sect. 3 are interpreted in the context of the solar corona. We also extrapolate on the findings to estimate the action of DWs at higher-density gradients. The limitations of the current study are discussed in Sect. 5, and our conclusions are presented in Sect. 6.

2. Gyrokinetic simulation setup

For this study we made use of the gyrokinetic Vlasov-Maxwell solver GENE (Jenko et al. 2000), which had previously been used for studies of astrophysical and space physics applications (Pueschel 2009; Pueschel et al. 2011, 2014, 2015; Told et al. 2015, 2016), in addition to many fusion applications. Gyrokinetics is a framework suitable to study kinetic processes in environments with strong guide fields. It assumes that any background magnetic field variation occurs on much larger scales than the Larmor radius and that any processes of dynamical relevance have frequencies much lower than the cyclotron frequency (∼1 to 100 kHz in the solar corona). Gyro-averaging eliminates the timescale of the Larmor motion and one velocity coordinate; as a result, gyrokinetic simulations typically exhibit speedups of several orders of magnitude when compared to fully kinetic simulations. Discussions about the validity regime of gyrokinetic simulations in astrophysical scenarios are provided in TenBarge et al. (2014), Schekochihin et al. (2009), and Howes et al. (2006). In addition, compared to gyrokinetic particle-in-cell simulations of DW turbulence, where noise can accumulate (Nevins et al. 2006), gyrokinetic Vlasov DW simulations do not suffer from this noise accumulation, even when run for long times (see, for instance, Fig. 5a of Peeters et al. 2016).

The magnetic geometry in which DWs are simulated here is a sheared slab that represents the entire length of the coronal loop or a similar gradient structure but only a small volume around a field line, with the macroscopic reference length L_z (where z is the coordinate parallel to the magnetic field) corresponding to the length of a structure such as a coronal loop. Ions and electrons, at hydrogen mass ratio, are assumed to have equal background temperatures T₀, setting the reference ion sound speed c_s, in a magnetic guide field B₀=B_z0, defining the ion sound Larmor radius ρ_s. The density gradient that drives the formation of DWs exists in the x-direction, giving rise to DWs drifting in the y-direction. The corresponding wavenumber k_y is normalised to 1/ρ_s, as is k_x, corresponding to a highly anisotropic process when comparing parallel and perpendicular scales.

As in Brchnelova et al. (2024), to explore the behaviour of DWs under various coronal conditions, several parameters are varied, namely the strength of the density gradient ω_n=L_z/L_n (with the density gradient scale length L_n equal for ions and electrons), electron β and magnetic shear $\hat{s}$. The magnetic shear $\hat{s}$ can be generally defined as $\hat{s}=L_z/L_s$, where L_s is the shear length, such that B=B(∇z−(x/L_s)∇y) in a slab geometry. In Brchnelova et al. (2024), we used the standard definition of $\hat{s}$ for toroidal systems, $\hat{s}=\frac{r_0}{q_0}\frac{\mathit{dq}(x)}{\mathit{dx}},$ with q(x) being the safety factor profile near a flux surface located at a distance r₀, corresponding to q₀, from the magnetic axis. The limits of $\hat{s}$ in the simulations are based on the field-line-twist profiles of coronal loops from Aschwanden (2019).

The baseline scenario is again defined by ω_n = 100, $\hat{s}=0.25$ and β = 0, with no (equal ion and electron) temperature gradient ω_T = 0. For this study, assuming the corona to be collisionless is appropriate.

The gyrokinetic distribution function is treated in Fourier space in the x-direction of the density gradient and the y-direction of wave propagation. At the same time, the z-direction along the guide field is computed using finite differencing. Velocity space covers the parallel velocity v_∥ and the magnetic moment μ. Domain sizes, mode, and grid-point resolutions were chosen such that doubling any setting would not substantially alter physical observables in the turbulent state. At $(\beta =0,{\omega }_n=100,\hat{s}=0.25)$, the number of grid points in the x direction and z direction are thus N_x×N_z = 128 × 80, the number of grid points in v_∥ and μ and N_v∥ and N_μ are 64 and 32, respectively, the converged velocity space domains in these coordinates are, respectively, L_v∥ = 3 and L_μ = 9, and the number of the resolved k_y modes is N_ky = 64, with the lowest finite $k_y^{\min }=0.0625$. The parallel boundary condition constrains the box size L_x and depends on the shear $\hat{s}$ as (Beer et al. 1995)

$$L_x=\frac{1}{\hat{s}{k}_y^{\min }}\hspace{0.5em}.$$(1)

Finally, the coefficient of the fourth-order hyper-diffusion (Pueschel et al. 2010) in the x-direction is set to D_x = 0.02, and in the z-direction, the corresponding D_z = 8; both values were determined to produce converged results.

3. Properties of drift-wave turbulence

With the code setup introduced, we now discuss the physical scenarios simulated and the corresponding results. In this section we limit ourselves to a non-dimensional discussion. More specific interpretations of these results assuming coronal conditions will be presented in Sect. 4.

3.1. Simulated scenarios

The solar corona contains a wide variety of environments, and the behaviour of DWs and the associated turbulence generally depends on all physical input parameters. Thus, in addition to the baseline scenario described in Sect. 2, simulations with varying 0≤β≤3.5·10⁻⁴ (where β = 0 is referred to as the electrostatic limit), $0.25\le \hat{s}\le 1$ and 50≤ω_n≤250 were conducted, as listed in Table 1. Note that the upper limit of the β range was chosen to facilitate a comparison with the solar-corona reconnection study in Pueschel et al. (2014).

All of the cases were evolved well into the quasi-stationary state of the DW turbulence (which is independent of the initial condition), where observables are time-averaged, typically for a period (measured in ion transit times) of 5−50 L_z/c_s.

3.2. Fluctuation properties

In the following, observables in the form of different moments of the distribution function taken at arbitrary times during the quasi-stationary, turbulent state are visualised: the electrostatic potential, Φ, the electron number density, n, the parallel flow velocity, u_∥, and the parallel and perpendicular temperatures, T_∥ and T_⊥. The corresponding normalisations read

$$\frac{\mathrm{\Phi }}{(T_0/e){\rho }^{*}},\hspace{0.5em}\frac{n}{n_0{\rho }^{*}},\hspace{0.5em}\frac{u_{\parallel }}{c_{\mathrm{s}}{\rho }^{*}},\hspace{0.5em}\frac{A_{\parallel }}{B_0{\rho }_{\mathrm{s}}{\rho }^{*}},\hspace{0.5em}\frac{T_{\parallel }}{T_0{\rho }^{*}},\hspace{0.5em}\frac{T_{\perp }}{T_0{\rho }^{*}}\hspace{0.5em}.$$(2)

Here, ρ^*=ρ_s/L_z is the ratio of the ion sound gyroradius to the loop length.

For case 1, which uses the electrostatic limit and thus does not produce magnetic fluctuations A_∥, Fig. 1 shows these fluctuations in the x-y plane, taken at the centre in the z domain along the field line. Except for Φ, these fluctuations are comparatively isotropic in x-y (recall that the z coordinate is normalised to L_z≫ρ_s, owing to an extreme anisotropy between the perpendicular and parallel directions), and given the small ρ^*, highlight the validity of the δf approach. Typical x-y structure sizes are of the order of 1−10ρ_s.

Fig. 1. Snapshot of fluctuation data in the quasi-stationary state for baseline case 1, with ωn = 100, $\hat{s}=0.25,$ and in the electrostatic limit. Shown are the electrostatic potential (Φ; where the dominant vertical structures correspond to a zonal flow), ion number density (n; the ion and electron density fluctuations are identical due to the quasi-neutrality constraint), and the parallel and perpendicular ion temperatures (T∥ and T⊥).

While a proper analysis of the DW saturation mechanism is beyond the scope of this paper, the prominent zonal flow (Diamond et al. 2005) – a vertical structure in Φ at k_y = 0 – in Fig. 1 suggests that zonal-flow-based energy transfer is a key process in saturating the linear instability, compare Makwana et al. (2014), Terry et al. (2021). By implication, strong zonal flows can lead to lower fluctuation amplitudes and vice versa for a given linear growth rate. Note that all cases studied, irrespective of driving gradient ω_n, display similarly strong zonal flows in relation to the non-zonal (k_y>0) turbulent fluctuation amplitudes. While the domain sizes L_y considered here are still on the scale of tens of gyroradii, zonal flows can in principle span much larger distances, making them a candidate for future observation efforts. Zonal flows have also been detected experimentally in fusion plasmas (Nishizawa et al. 2019), and are routinely seen in simulations of turbulence in fusion plasmas; recent work on stellarator turbulence driven by the universal instability, which is essentially the same mechanism as that studied in the present paper, reports similarly substantial amplitudes of the zonal flow (Costello et al. 2023).

For the cases with β>0, a similar zonal structure forms in the fluctuations in the magnetic vector potential A_∥, as seen in Fig. 2 for case 6; given its low amplitude, however, it remains to be determined whether this feature can exert any influence on the turbulent saturation mechanism of the DWs. More importantly, however, even at low β, the appearance of magnetic fluctuations enables interactions among and between plasmoids and current sheets, allowing them to produce secondary magnetic reconnection (see Pueschel et al. 2014).

Fig. 2. Turbulent fluctuations for simulation case 6, where β = 10−4, allowing the system to produce magnetic perturbations A∥ (right panel). Current filaments (left panel) can thus interact via magnetic reconnection.

3.3. Spectral fluxes

The fluctuation data shown in Figs. 1 and 2 highlight an inherent property of turbulence: the coexistence and interaction of multiple spatial scales. Typical structure sizes substantially exceed that of the fastest-growing linear mode (at k_y = 4) due to mixing-length effects, allowing larger structures to substantially influence observables such as fluxes. Figure 3 shows the k_y spectrum of the particle flux, normalised to n₀c_s(ρ^*)², for different parameter cases; the contribution of the higher wavenumbers is limited, with fluxes peaking in the range 0.1<k_y<0.4. Increased $\hat{s}$ or decreased ω_n, and thus reduced DW drive, shift flux peaks to slightly higher k_y. The features in the heat-flux spectra are very similar, with heat fluxes (normalised to n₀T₀c_s(ρ^*)²) scaling approximately as Q∼Γ.

Fig. 3. Particle flux, Γ, as a function of ky for different parameter cases. Due to ambipolarity, electrons and ions experience equal particle fluxes. Electromagnetic contributions to the flux are negligible for all cases.

3.4. Turbulent frequency signatures

In the absence of spectrally resolved frequency measurements, the spectral distribution of the fluctuation amplitudes – as, for instance, expressed in the fluxes – informs which frequencies are expected to be most prominent in the system overall. Power spectra for the two cases with the strongest and weakest DW drives are shown in Fig. 4. Note that in our convention, positive frequencies correspond to the drift in the ion diamagnetic direction. There, the range of frequencies widens at higher perpendicular wavenumbers, but as was seen in Fig. 3, only the lower wavenumbers, generally below k_y = 0.4, are excited to large amplitudes; on small scales, energy input due to linear instabilities tends to be weak compared to non-linear transfer, resulting in highly broadened frequency signatures of low intensity. For observational purposes, the focus, therefore, lies on low-k_y frequencies.

Fig. 4. Frequency power spectra for cases 2 (top) and 7 (bottom). Intensities are normalised separately at every ky. Narrow spectra at lower ky correspond to linear or coherent non-linear signatures, while at higher ky≳1, a strong-turbulence regime is attained.

Linear frequencies of the dominant eigenmodes (Brchnelova et al. 2024) tend to lie in the range of −3≲ω_lin≲−1 for the present parameters at wavenumbers where the non-linear frequencies are strongly coherent. However, the existence of many subdominant modes, which are substantially more difficult to probe, complicates the link between linear and non-linear physics. While future work will have to investigate this matter in more depth, at present there exists no indication that the frequency bands originating from ω≈−20 at the lowest k_y and trending downwards until k_y≈0.4, as well as the frequency band in the higher-shear case at positive ω for 0.3≲k_y≲0.8, are corresponding to (sub-dominantly unstable) linear eigenmodes. Instead, a parallel can be drawn with heuristically similar frequency signatures (Faber et al. 2015; Pueschel et al. 2016) observed in $\hat{s}\ll 1$ stellarator DW turbulence, the existence of which is not well understood.

The behaviour of the non-linear frequency ranges with changing ω_n and β, however, agrees well with the linear results (with the frequencies increasing with an increasing ω_n and decreasing for an increasing β). Preliminary simulations with a curved coronal loop geometry (which may be more representative of the solar corona) have revealed that with these three-dimensional effects included, the linear and non-linear frequencies are closer to each other, and they are lower than what is shown in Fig. 4, albeit only by an order-unity factor. For now, however, in this work, we are interested in an order-of-magnitude analysis, as the specific results are highly dependent on the exact values of $\hat{s}$, β, ω_n and B_z, L_z and T₀, which vary significantly in the corona. Thus, while we consider the non-linear frequency ranges reported above in the subsequent analysis, one should note that more realistic dominant frequencies may be lower by a factor of a few.

3.5. Kinetic Alfvén waves

An additional observable that can be compared to in situ measurements is the turbulence spectrum. The spectral slope of the magnetic energy in the quasi-inertial range (k_y≳0.7), as can be seen in Fig. 5 for case 8, is consistent with predictions in Howes et al. (2008), where KAW turbulence spectral slopes were derived for $E_{B\perp }\propto k_{\perp }^3{A}_{\parallel }^2$: using the A_∥ slopes measured from the present data, $E_{B\perp }\propto k_y^{-2.7}$ and $E_{B\perp }\propto k_x^{-2.5}$, close to the predicted −7/3. Notably, these values not only match those obtained from gyrokinetic simulations of reconnection turbulence (Pueschel et al. 2014) but are similarly consistent with observed spectra of the solar wind (Alexandrova et al. 2009; Chen et al. 2012).

Fig. 5. ky (blue) and kx (red) spectra of the perturbed magnetic vector potential, A∥, along with the respective fits of $k_y^{-5.8}$ and $k_x^{-5.4}$ in the quasi-inertial range (beyond the dotted line) for case 8.

3.6. Parallel kinetic acceleration and stochastic heating

From the perspective of solar coronal physics, a particularly interesting aspect of DWs is their ability to cause plasma heating. As discussed in the introduction, DWs may produce parallel heating via kinetic particle acceleration and perpendicular heating via triggering stochastic motion due to the electrostatic field.

To assess parallel kinetic particle acceleration, one can evaluate the DW-generated perturbed parallel electric field E_∥, which enables the determination of the average volumetric heating, j_∥E_∥, where j_∥ is the perturbed current density. This analysis was carried out for all simulation cases, with examples of the time evolution of E_∥ and j_∥E_∥ for the baseline case 1 (electrostatic) and for case 7 (β = 3.5·10⁻⁴) shown in Fig. 6. There, E_∥ is split into the static, flutter and inductive components (Pueschel et al. 2014) as

$$E_{\parallel }=E_{\parallel }^{\mathrm{es}}+E_{\parallel }^{\mathrm{fl}}+E_{\parallel }^{\mathrm{ind}}=\frac{d\mathrm{\Phi }}{\mathit{dz}}+\frac{{\text{B}}_{\perp }·{\text{E}}_{\perp }}{B_0}-\frac{{\mathit{dA}}_{\parallel }}{\mathit{dt}},$$(3)

where only the first term is finite in the electrostatic limit β = 0. For the simulations with a non-zero β (allowing the generation of perturbations of the magnetic vector potential), the perturbed magnetic field projects perpendicular electric field fluctuations along the unperturbed magnetic field (see Rechester & Rosenbluth 1978), forming the flutter component, and the time variation of the perpendicular magnetic field B_x,y∼A_∥ produces the inductive component. Even in the higher-β cases, however, the static component is dominant.

Fig. 6. Root-mean-square parallel electric field, E∥, normalised to $T_0{\rho }_{\mathrm{s}}/(e{L}_z^2$), and the heating rate, j∥E∥, normalised to $n_0{c}_{\mathrm{s}}{T}_0{\rho }_{\mathrm{s}}^2/L_z^3$, for the baseline case 1 (β = 0) and case 7 (β = 3.5·10−4). Heating primarily arises from the electrostatic component of the electric field.

It is important to acknowledge that the timescale over which this heating is able to alter the equilibrium distribution substantially and, thus, the background temperature is very long and ordered separately from the simulation time (t). The same applies to particle and heat fluxes, a common assumption for gyrokinetic flux-tube simulations.

The time-averaged root-mean-square j_∥E_∥ is listed in Table 2 in code units. The values are rounded to two significant digits due to the oscillations around the mean. As would be expected based on the linear growth rates, the volumetric heating rates most strongly scale with ω_n, spanning two orders of magnitude between ω = 50 and ω = 200. Volumetric heating is also observed to decrease with $\hat{s}$ and with the electron β.

Next, stochastic heating in the perpendicular direction is considered. Since Fig. 5 revealed that the spectrum resembles KAW turbulence, we can follow the work of Chandran (2010) on stochastic heating in Alfvén-wave (AW) turbulence. For AW and KAW fluctuations, we can express the amplitude threshold ϵ_i for strong stochastic heating, and from it the timescale t_h on which stochastic heating doubles an ion's kinetic energy, as

$${ϵ}_{\mathrm{i}}\sim \frac{q_{\mathrm{i}}\delta \mathrm{\Phi }}{m_{\mathrm{i}}{v}_{\perp \mathrm{i}}^2},t_{\mathrm{h}}\approx {ϵ}_{\mathrm{i}}^{-3}{\mathrm{\Omega }}_{\mathrm{i}}^{-1}exp\left(\frac{0.15}{{ϵ}_{\mathrm{i}}}\right),$$(4)

where v_⊥i=(2k_BT_⊥i/m_i)^1/2, Ω_i=q_iB_z/(m_ic), m_i is the ion mass, q_i the ion charge, c the speed of sound and k_B the Boltzmann constant, all in cgs units. Stochastic heating is thus relevant in environments where t_h is sufficiently small (i.e. heating occurs on short enough timescales) to be relevant for coronal heating. For this reason, the key observable from simulations in this context is δΦ, defined here as the absolute value of the peak electrostatic fluctuation amplitude. Note that since we do not have a background electrostatic potential, δΦ=Φ. This δΦ in code units is listed in Table 2. As we have seen for the fluxes, δΦ scales strongly with ω_n, resulting in stronger stochastic heating, interpreted physically later in Sect. 4 for coronal conditions.

3.7. Density structure erosion and lifetime

Finally, an essential aspect of coronal heating is the time over which this heating can be sustained – if the heating source are DWs driven by density gradients, this time is constrained by the lifetime of the density gradients. The timescale over which the gradient would erode can be estimated from the particle flux Γ (recall that due to quasi-neutrality and ambipolarity, hydrogen ions and electrons have equal number densities, equal density gradients, and equal particle fluxes) obtained from the simulations, with lifetime ${\tau }_{\mathrm{\Gamma }}\sim \frac{n_0{L}_n}{\mathrm{\Gamma }}$, in which L_n is the density gradient scale length. For each simulated case, Γ is listed for ions in Table 2 in code units. Similarly to other observables, we see strong scaling with ω_n, implying that, as expected, in the scenarios with higher density gradients, those gradients erode more quickly. Section 4 discusses these results as applied to the solar corona. As with parallel heating, we ordered out any changes to the background distributions throughout the simulation, as the simulation timescale was assumed to be much shorter than the timescale on which the structures erode. As the gradient is eroded over longer times due to the particle flux, linear DW drive will be diminished, leading to lower particle fluxes and heating rates.

In this context, however, it is helpful to recall that complete structure erosion can only occur when no restoring force exists (such a force could have been responsible for the existence of the gradient in the first place). Whether such a force persists in a given case depends on the specific scenario, and thus τ_Γ cannot be generalised to be relevant for all structures we see in the corona.

4. Implications for the solar corona

The results from Sect. 3 can now be interpreted physically by inserting normalisations appropriate for the solar corona. The solar corona contains various conditions and structures; however, as shown in this section, the simulated physical volumetric heating rates, timescales and frequencies, as per Sect. 3, depend heavily on the specific conditions. For this reason, we do not limit this discussion to just one set of reference parameters and instead perform evaluations for a wide range of conditions.

Since focused on the base of the corona, we considered T₀ temperatures in the range 0.1 MK to 10 MK, with the smaller temperatures corresponding to, for example, prominences (Brughmans et al. 2022) and the hotter ones corresponding to super-hot coronal loops above active regions (Lu et al. 2024). Magnetic fields are taken as B_z = 0.1 G to 10 G (Alissandrakis & Gary 2021) and structure lengths L_z lie between 10⁵ m and 10⁸ m (Reale 2010), with the background number density determined from the electron β for β>0, or from β∼10⁻⁵ for the electrostatic cases.

4.1. Frequencies, amplitudes, and wavelengths

The magnitudes of the frequencies at the non-linearly most strongly excited wavenumbers (k_y≲0.4) were in Sect. 3.4 found to vary between a few tens of c_s/L_z for cases 3 and 7 up to 80 c_s/L_z for case 2 (see Table 2). As c_s is temperature-dependent, this quantity determines the temperature dependence of the predicted observable frequency. The results for a range of temperatures and loop lengths are shown in the top row of Fig. 7 for cases 2 and 7, with these selected as they lie on the extreme ends of the range of cases. For the set of reference coronal conditions used in Brchnelova et al. (2024), we obtain ω = 0.145 mHz and thus frequencies in the range of 0.1 to 10 mHz. Figure 7, however, reveals that hertz-range frequencies are possible for hotter and smaller structures.

Fig. 7. For cases 2 (left) and 1 (right): the (absolute values of the) physical frequencies, parallel volumetric heating rates, timescales of stochastic heating, and structure lifetimes as functions of the length scale (Lz; on the x-axes) and either the background magnetic field (Bz) or the background temperature (T0).

The obtained range agrees with the frequencies seen in the linear simulations and corresponds to the frequencies of the oscillations observed and described by Kepko et al. (2002), Kepko & Spence (2003), and Viall et al. (2009, 2021). Future study will need to reveal how DW fluctuations generated in the solar corona propagate into the solar wind, where these periodic mesoscale structures were measured. However, if these periodic structures are associated with fluctuations in the electrostatic field, they may indeed be connected to the DW mechanism.

We can also determine whether the fluctuations resolved in our simulations are directly observable in the solar corona. The strongest fluctuations were observed in case 2 and have an amplitude of a few tens of thousands. The density fluctuations are normalised to n₀ρ^*, meaning that if we assume B₀∼10 G and L_z∼10⁸ m, these fluctuations would correspond to ∼0.5% to 5·10⁻⁵% oscillations relative to the background density. These scale linearly with the inverse of L_z, implying that shorter structures would be more straightforwardly observable. In Fig. 1 we also see that comparable or slightly smaller relative fluctuation amplitudes are produced in T_∥, T_⊥ and u_∥, in the case shown between 10⁻³% and 10⁻⁴% of the background values.

However, as evident from Figs. 1 and 2, DW oscillations occur on wavelengths of a few to a few tens of Larmor radii, which, for the investigated range of B₀, corresponds to only ∼1−100 m. Resolving these scales in observations is unfeasible with current technology. Detection of DW fluctuations may become possible when DWs propagate into regions with lower magnetic fields and farther from the surface of the Sun, enabling an expansion of the associated length scales, even if DWs are no longer unstable in those regions. However, predicting the specifics of DW signatures under such circumstances is a complex task and will have to be the subject of a separate effort.

Brchnelova et al. (2024) showed that DWs have a much larger extent in the direction parallel to the background magnetic field (typically corresponding to low-order harmonics along the coronal loop). Thus, if one could conduct a measurement aligned very precisely with the magnetic field, DW fluctuations could become discernible along the coronal loop length.

4.2. Parallel and perpendicular heating

Similarly, we can quantify to what extent DWs heat the surrounding plasma in the present simulations (see Sect. 3.6). We started with parallel heating. The baseline volumetric heating rate, j_∥E_∥, can be obtained by applying to the value in code units the normalisation of $T_0{n}_0{c}_{\mathrm{s}}{\rho }_{\mathrm{s}}^2/L_z^3$. Because, for a given β, the guide field $B_0^2$ in n₀ cancels that in ${\rho }_{\mathrm{s}}^2$, the heating rate only depends on L_z and T₀. Thus, we varied T₀ and L_z across the ranges mentioned above, with the resulting j_∥E_∥ shown in the second row of Fig. 7.

The volumetric coronal heating that global coronal models require to heat the corona sufficiently and accelerate the solar wind is in the range of ∼10⁻⁶−10⁻³ erg/(cm³s), depending on the scenario considered (see e.g. Baratashvili et al. 2024). The values displayed in Fig. 7 are well below this range. However, heating increases significantly with increasing ω_n. The maximum density gradient modelled in the present paper is ω_n = 250 due to computational-resource limitations. However, density gradients across coronal loops of thousands or more may well be realistic (see, for instance, the discussion in Brchnelova et al. 2024 about the work of Cargill et al. 2016). We thus performed an extrapolation from the existing data to determine at which ω_n parallel heating via j_∥E_∥ can become significant (see Sect. 4.4).

In the perpendicular direction, we evaluated whether stochastic heating can reach relevant levels based on the work of Chandran (2010) and the fluctuations of the electrostatic potential, δΦ, as outlined in Sect. 3.6. The maximum electrostatic potential fluctuations shown in our simulations range between 500 and 8·10⁴ in code units, normalised to (T_0,eV/e)ρ^* (see Table 2). For example, for a 10 G field with a loop length of 10⁸ m, we obtain voltage fluctuations between 5·10⁻⁴ V (for case 7) and 8·10⁻² V (for case 2).

We then determined τ_h using Eq. (4). As would be expected, τ_h is highly sensitive to the constituents of ϵ_i as a result of the exponential term, and the value of ϵ_i is primarily dependent on B₀ and L_z. For this reason, we held the plasma at 1 MK and visualise the stochastic timescale τ_h in the third row of Fig. 7. The exponential term quickly drives τ_h to very high values with increasing L_z and B₀, so we only plot those regions in which τ_h is realistic from the perspective of coronal heating (to a maximum of a few hours). Especially for the cases with stronger density gradients, almost the entire range shown results in a substantial τ_h, implying that stochastic heating is very efficient in most investigated cases. As mentioned for parallel heating, these results will be extrapolated to higher density gradients in Sect. 4.4.

4.3. Erosion timescales

Finally, we quantified the lifetime, τ_Γ, of the gradient structures based on the particle fluxes listed in Table 2. τ_Γ is again most sensitive to L_z (which for a given ω_n determines L_n) and to B₀, affect Γ through (ρ^*)². Thus, we plot τ_Γ for the aforementioned ranges of L_z and B₀ in the bottom row of Fig. 7. The results show that most of the structures can survive, even without continuous driving, on the timescales that correspond to the typical lifetimes of regular coronal loops (tens of minutes to hours and days; see e.g. Klimchuk et al. 2010 and Mulu-Moore et al. 2011) and longer.

4.4. Extrapolation to higher ω_n

The simulations conducted here do not explore the full range of normalised density gradients expected in the corona, which may reach thousands, as discussed in Brchnelova et al. (2024). As mentioned above, simulating higher-density gradients requires numerical resolutions that exceed currently available computational resources. However, from the data listed in Table 2, we can determine extrapolation power laws that allow us to estimate ω, j_∥E_∥, τ_h (via δΦ) and τ_Γ (via Γ_i) for higher ω_n. This can be determined separately for the $\hat{s}=0.25$ electrostatic and β = 3.5·10⁻⁴ sets of simulations. These are plotted in Fig. 8 along with power-law fits.

Fig. 8. For the electrostatic cases (blue squares) and β = 10−4 cases (red diamonds): scalings of the DW frequency (top left), parallel volumetric heating (top right), peak electrostatic potential (bottom right), and particle flux (bottom left) with the density gradient.

Using these extrapolation laws, we estimate that at ω_n = 1000, with L_z = 10⁷ m, T₀ = 1 MK and B₀ = 10 G, frequencies of up to ∼700 mHz are expected, along with a relatively low parallel heating of 10⁻⁸ erg/cm³/s, a fast perpendicular stochastic heating with τ_h∼7·10⁻⁶ s, and long lifetimes of τ_Γ∼5·10⁹ s. The drift velocity, v_y, for this higher gradient would reach 1 to 10 m/s. According to the same approximations, at somewhat higher gradients ω_n∼5000, the parallel heating starts to be significant (>10⁻⁵ erg/cm³/s). However, such ω_n values are far above those used to create these scalings, and higher density-gradient simulations will need to be conducted in the future to confirm these findings.

However, one factor remains to be considered before concluding the efficacy of heating: the present simulations use a sheared-slab geometry, whereas a coronal loop inherently exhibits field-line curvature and even twist (Aschwanden 2019). Preliminary gyrokinetic simulations with a three-dimensional loop geometry inspired by Tsao (2023) suggest that for a given β and ω_n, DW-driven turbulent amplitudes can be substantially boosted, with an order of magnitude higher Γ_i and two orders of magnitude higher j_∥E_∥ straightforwardly achievable. Further pursuit of DW turbulence in three-dimensional loop geometry will thus have a high priority going forwards, particularly as this boost makes parallel heating dynamically relevant for the solar corona at gradient scale lengths more common in this region.

5. Limitations of the modelling approach

While this work has demonstrated that DW turbulence constitutes a powerful mechanism to explain a range of observed features in the solar corona, further efforts are necessary to improve our confidence regarding the practical applications and interpretation of these results.

Firstly, this paper only considers slab geometry and does not consider curvature and twist effects. Preliminary simulations relaxing this constraint show stronger DW drive and heating.

Secondly, the present simulations technically require that magnetic field lines lie on closed flux surfaces and that fluctuations are thus (quasi-)periodic along the field line. More realistic parallel boundary conditions will need to be studied for both closed loops and regions with open field lines. It has already been suggested that line-tying of the (both background and perturbed) magnetic field does not substantially affect gyrokinetic turbulence (Tsao 2023).

Thirdly, the scenario considered here assumes no other equilibrium gradients are present. Linearly, the presence of additional temperature gradients can moderately stabilise density-gradient-driven DWs. Temperature gradients are also known to exist in the solar corona, especially in structures such as prominences, for instance, the measurements of Mouradian et al. (1979) or the models of Anzer & Heinzel (2008) and Brughmans et al. (2022). For example, from the simulations of Brughmans et al. (2022), one can estimate ${\omega }_T=-\frac{L_z}{T_0}\frac{{\mathit{dT}}_0}{\mathit{dx}}\lesssim 10$ in prominences, and in Brchnelova et al. (2024), we show that such ω_T may have a limited influence on DW growth rates. The effect of such additional inhomogeneities could be included in future work.

6. Conclusions

We used the gyrokinetic turbulence code GENE to determine the expected properties of density-gradient-driven DW turbulence with full non-linear simulations. The geometry considered was a sheared slab representing the immediate vicinity of a magnetic field line in a coronal loop. A prescribed density gradient perpendicular to the background magnetic field leads to the spontaneous generation of DWs via the DW instability. Multiple parameter cases were computed, in which we varied the density gradient scale length, magnetic shear, and electron β to characterise the behaviour of DWs in various coronal environments. For these cases, we have the reported flux spectra, frequency power spectra, fluctuation spectra, parallel volumetric heating rates, stochastic heating timescales, and erosion timescales that indicate the possible structure lifetime without any background driving mechanism. Thereafter, for a range of background magnetic field strengths, loop lengths, and temperatures, we computed the physical values of the expected DW frequencies, parallel volumetric heating rates, stochastic heating timescales, and structure lifetimes. The main takeaways from the analysis are:

1. The self-induced DWs in the solar corona generate zonal flows, current sheets, and plasmoids, and may thus also generate an environment suitable for secondary magnetic reconnection.

2. The spectrum of the perturbed magnetic potential at the high-k_x, k_y tail corresponds to the slope of KAW turbulence.

3. The turbulence has a spectral peak in the range 0.1≲k_yρ_s≲0.4, with stronger density gradients producing peaks at lower k_y.

4. Depending on the loop length and temperature of the plasma, the expected DW frequencies in the solar corona range from ∼0.1 mHz to a few hertz.

5. The parallel electric fields and current densities generated by DWs may give rise to non-negligible kinetic particle acceleration in the parallel direction; simulations in more realistic geometry and/or at higher gradients will be required to confirm whether the corresponding heating rates are dynamically relevant in the solar corona.

6. The amplitude of the fluctuations of the electrostatic potential indicates that in almost all of the investigated conditions, DWs lead to fast stochastic heating in the perpendicular direction, which is consistent with the observation that T_⊥>T_∥ in the solar corona. Moreover, stochastic DW heating is known (from thermonuclear fusion research and experiments) to be more efficient on ions, consistent with the observation that T_i>T_e in the solar corona.

7. The particle flux implies that most structures can exist despite erosion that lasts for minutes to days, even without a background driving mechanism.

8. Compared to other types of waves in the solar plasma, DWs can be distinguished by their largely electrostatic fluctuation signature.

9. However, since the wavelengths of DWs are only a few to a few tens of Larmor radii, our current instruments can only observe them if they are generated in or propagated into regions of very weak magnetic fields. Propagation of DWs from the corona into the solar wind should thus be further studied so that missions, such as the Solar Orbiter and Parker Solar Probe, can identify them.

As discussed in Sect. 5, future work will need to focus on DW turbulence in more realistic coronal loop geometry and determine the impact of realistic boundary conditions. Such studies could also include a moderate temperature gradient. In addition, the impact of possible non-thermal components of the background electron distribution should be assessed, as they may affect DW frequencies and growth rates.

Our results suggest that DWs can trigger both parallel and perpendicular heating and produce frequencies and turbulence spectra consistent with observations in the solar wind. This highlights the considerable potential of DWs in solar physics while simultaneously showing that gyrokinetic turbulence simulations can be a powerful tool even when applied to physical scenarios very different (Pueschel et al. 2017) from terrestrial plasma confinement experiments.

Acknowledgments

This research was funded by projects C16/24/010 (C1 project Internal Funds KU Leuven), G0B5823N and G002523N (FWO-Vlaanderen), and 4000145223 (SIDC Data Exploitation (SIDEX2). The resources and services used in this work were provided by the VSC (Flemish Supercomputer Centre), funded by the Research Foundation–Flanders (FWO) and the Flemish Government, and by the Tier 1 VSC grant project 2023/019. SP is also funded by the European Union. Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or ERCEA. Neither the European Union nor the granting authority can be held responsible for them. His project (Open SESAME) has received funding under the Horizon Europe programme (ERC-AdG agreement No 101141362).

References

All Tables

All Figures

Fig. 1. Snapshot of fluctuation data in the quasi-stationary state for baseline case 1, with ωn = 100, $\hat{s}=0.25,$ and in the electrostatic limit. Shown are the electrostatic potential (Φ; where the dominant vertical structures correspond to a zonal flow), ion number density (n; the ion and electron density fluctuations are identical due to the quasi-neutrality constraint), and the parallel and perpendicular ion temperatures (T∥ and T⊥).

In the text

	Fig. 2. Turbulent fluctuations for simulation case 6, where β = 10−4, allowing the system to produce magnetic perturbations A∥ (right panel). Current filaments (left panel) can thus interact via magnetic reconnection.
In the text

	Fig. 3. Particle flux, Γ, as a function of ky for different parameter cases. Due to ambipolarity, electrons and ions experience equal particle fluxes. Electromagnetic contributions to the flux are negligible for all cases.
In the text

	Fig. 4. Frequency power spectra for cases 2 (top) and 7 (bottom). Intensities are normalised separately at every ky. Narrow spectra at lower ky correspond to linear or coherent non-linear signatures, while at higher ky≳1, a strong-turbulence regime is attained.
In the text

	Fig. 5. ky (blue) and kx (red) spectra of the perturbed magnetic vector potential, A∥, along with the respective fits of $k_y^{-5.8}$ and $k_x^{-5.4}$ in the quasi-inertial range (beyond the dotted line) for case 8.
In the text

	Fig. 6. Root-mean-square parallel electric field, E∥, normalised to $T_0{\rho }_{\mathrm{s}}/(e{L}_z^2$), and the heating rate, j∥E∥, normalised to $n_0{c}_{\mathrm{s}}{T}_0{\rho }_{\mathrm{s}}^2/L_z^3$, for the baseline case 1 (β = 0) and case 7 (β = 3.5·10−4). Heating primarily arises from the electrostatic component of the electric field.
In the text

	Fig. 7. For cases 2 (left) and 1 (right): the (absolute values of the) physical frequencies, parallel volumetric heating rates, timescales of stochastic heating, and structure lifetimes as functions of the length scale (Lz; on the x-axes) and either the background magnetic field (Bz) or the background temperature (T0).
In the text

	Fig. 8. For the electrostatic cases (blue squares) and β = 10−4 cases (red diamonds): scalings of the DW frequency (top left), parallel volumetric heating (top right), peak electrostatic potential (bottom right), and particle flux (bottom left) with the density gradient.
In the text