© The Authors 2025

1. Introduction

There are observational indications that the solar corona and the solar wind are filled with Alfvén waves. In situ observations allow us to measure Alfvén waves directly in the solar wind (Bruno & Carbone 2013). However, in the corona they are usually only identified as spectral-line broadenings. From spectroscopic measurements, it was thought that the Alfvén waves’ energy is considerable (Banerjee et al. 1998; Hahn & Savin 2013; Pant & Van Doorsselaere 2020). Recently, a few direct measurements became available (De Pontieu et al. 2012; Shetye et al. 2021; Petrova et al. 2024). Despite these, the true energy content of Alfvén waves in the corona is not so well established.

On the other hand, high-resolution imaging observations from SDO/AIA and SolO/EUI have clearly shown that loops and plumes show ubiquitous transverse motions (Anfinogentov et al. 2013; Nakariakov et al. 2024). Since the loop is displaced from its axis, these transverse motions are often identified as kink waves (Van Doorsselaere et al. 2008; Goossens et al. 2009), in contrast to Alfvén waves, which are thought to have a torsional motion in a cylindrical structure. The transverse waves manifest as the decayless transverse waves (Tian et al. 2012; Wang et al. 2012; Nisticò et al. 2013) in coronal loops. These loops are thought to be steadily fed with energy from its convective footpoint motions (e.g. Karampelas & Van Doorsselaere 2021 and references therein) in order to counteract the strong damping that is observed for impulsively excited waves during flares (e.g. Nechaeva et al. 2019 and references therein). They have a standing character (Anfinogentov et al. 2015), even though that is not entirely clear in very short loops (Gao et al. 2022; Shrivastav et al. 2024).

In coronal plumes or long coronal loops, the transverse waves are observed as propagating waves (Tomczyk & McIntosh 2009; Thurgood et al. 2014). Also, these propagating waves are found to be damped (Tiwari et al. 2021), although the weak damping rate seems to suggest only low density contrasts (Morton et al. 2021). It is unclear how these propagating kink waves are transformed higher up in the atmosphere to Alfvén waves (as they should, as evidenced in PSP observations; see e.g. Parashar et al. 2020) when the density contrast of the plumes is vanishing.

All these transverse waves must carry a significant amount of Poynting energy flux to the higher layers of the solar atmosphere. Energy estimates from observations go from a few Watts per square meter (Tomczyk et al. 2007; Thurgood et al. 2014) to several 1000 W/m² (Petrova et al. 2023). While some of the individual observed oscillations could counteract the heating requirements of quiet-Sun or even active regions, the potential of wave heating with transverse waves was only recently studied in a statistical way. Lim et al. (2023) performed a meta-analysis of transverse wave events reported in the literature to show that their energy forms a power law of their frequency, with a supercritical slope. Thus, Lim et al. (2023) found that transverse waves can heat the corona when acting as an ensemble.

All these observational results motivated the development of 3D heating models for coronal loops and plumes (Van Doorsselaere et al. 2020a). These models rely on the development of resonant absorption (e.g. Goossens et al. 1992) kickstarting the Kelvin-Helmholtz instability (KHI) (Terradas et al. 2008; Antolin & Van Doorsselaere 2019). The KHI further develops under continued driving and encompasses the whole loop (Karampelas & Van Doorsselaere 2018). In some simulations, this mechanism was able to balance the optically thin radiative losses (Shi et al. 2021), although that seems to be strongly dependent on the specific loop conditions (De Moortel & Howson 2022). Moreover, the heating provided by propagating transverse waves was modelled in a section of the corona by Magyar & Nakariakov (2021), and similar models for coronal heating by Alfvén waves have been constructed by Suzuki & Inutsuka (2005), Shoda et al. (2019), and follow-up works. Despite the promise that these numerical models hold, it has not yet been possible to model the entire solar atmosphere heated with these wave-heating mechanisms because of the required numerical resources. So, it is necessary and instructive to parametrise the heating by transverse waves as 1D wave-heating models and apply it to the full solar atmosphere. One example of this parametrisation is shown by Verdini et al. (2010), which incorporated turbulent heating from Alfvén waves and explained the heating of plasma in an expanding coronal hole and the acceleration of the solar wind.

More recently, other models driven by Alfvén waves have been developed in 3D (AWSOM by van der Holst et al. 2014, Downs et al. 2016 and Réville et al. 2020). Such models are being used all over the world for space weather forecasting. They are excellent at predicting the overall shape of the solar corona (Riley et al. 2019). Nowadays, similar models are being used to model stellar atmospheres as well (see e.g. Alvarado-Gómez et al. 2022; Evensberget & Vidotto 2024; Cohen et al. 2024, for some recent references). However, in these models, only the wave-driving power of the Alfvén waves is taken into account. Moreover, the influence of background plasma structuring (except the structuring of the waves themselves; see L_⊥) perpendicular to the magnetic field on the wave dynamics has been ignored, despite the efforts of Evans et al. (2012). Here, we are motivated by the observations to also incorporate kink-wave driving of the solar and stellar atmospheres, following the approaches taken for the modelling of Alfvén-wave driving.

To accomplish this, we used the recent description for wave dynamics in terms of Q variables (Van Doorsselaere et al. 2024). These Q variables allowed us to track single waves in a co-propagating manner and to include non-linear effects. The latter turns out to be important for propagating kink waves, because they are known to show a turbulent evolution for a single wave pulse (Magyar et al. 2017). This phenomenon has been named ‘uniturbulence’, because contrary to Alfvén-wave turbulence in a uniform background, only one unidirectional kink wave is needed to generate small scales. This has been explained by Magyar et al. (2019), Van Doorsselaere et al. (2020b) by showing that kink waves naturally consist of both Elsässer variables simultaneously, which are constantly interacting. That leads to damping of surface Alfvén waves in numerical simulations and analytical theory (Ismayilli et al. 2022).

Van Doorsselaere et al. (2024) provided the mathematical framework of the Q variables and rewrote the MHD equations in that formalism. We show in that paper that the Q variables allow a general wave perturbation to be split in its constituent wave modes. This would thus allow us to complement the AWSOM idea (i.e. the MHD equations for the large-scale evolution of the corona, complemented by Alfvén-wave-evolution equations) with an additional equation for the kink-wave evolution.

In this paper, we first generalise the kink-wave description to a flowing plasma as in the solar wind in Sect. 2. Then we recapitulate the MHD equations in terms of Q variables (as found in Van Doorsselaere et al. 2024) in Sect. 3. From the combination of these sections, we derive a general energy evolution equation in the Q formalism in Sect. 4.2. Then, we apply this equation to the propagation and damping of kink waves (Sect. 4.2.3). We formulate the set of equations for MHD plus Alfvén- and kink-evolution equations in Sect. 5, which we call the UAWSOM (uniturbulence and Alfvén wave solar model) equations. Finally, we present the initial results of a proof-of-concept implementation in a Python package (Sect. 6).

2. Kink-wave variables and terminology

For the UAWSOM model, we considered the additional driving of kink waves in flux tubes, next to the Alfvén waves. We considered the solar atmosphere filled with flux tubes, aligned with the magnetic field. The flux tubes we considered should model the coronal plumes that carry the main magnetic wave energy out into the solar wind. As such, in this section we first need to extend the model of Van Doorsselaere et al. (2020b) to also include the background flow. However, since many of the calculations are identical to the equations and steps in Van Doorsselaere et al. (2020b), we only keep the bare minimum needed to obtain the required expressions for the calculations in Appendix A.

To model kink waves in coronal plumes with background flow, we take an equilibrium configuration of a straight cylinder with homogeneous magnetic field B₀ = B₀e_z, constant and uniform background flow of V = Ve_z, without gravity and gas pressure p = 0. In the first instance, we do not consider these equilibrium quantities (or only weakly) to vary with height. We relax this assumption in later sections in a WKB way. Our current configuration is most easily described by the cylindrical coordinate system (r, φ, z). We take a radial step function in density:

$${\rho }_0(r)=\{\begin{array}{cc}{\rho }_{\mathrm{i}}\hfill & \text{for}r\le R,\hfill \\ {\rho }_{\mathrm{e}}\hfill & \text{for}r>R.\hfill \end{array}$$

The subscripts i, e are associated with the interior and exterior regions, respectively. R is the radius of the coronal plume.

As in Van Doorsselaere et al. (2020b), we take the total pressure perturbation $P\prime =\frac{B_0\delta B_z}{\mu }$ as the basic variable. We only consider propagating waves in this paper, which we describe with

$$\begin{array}{c}\hfill P\prime (r,\phi ,z,t)=\mathcal{R}(r)cos\phi cos(k_{z}z-\omega t).\end{array}$$(1)

The function ℛ(r) has the radial dependence of the pressure variation, for which we use the thin-tube limit δ = k_zR ≪ 1 and write

$$\begin{array}{c}\hfill \mathcal{R}(r)=\{\begin{array}{cc}A\frac{r}{R}\hfill & \text{for}r\le R\hfill \\ A\frac{R}{r}\hfill & \text{for}r>R\hfill \end{array},\end{array}$$(2)

with wave amplitude A. According to Goossens et al. (1992), in the thin-tube limit for a uniform flow, the wave frequency equals the kink frequency:

$$\begin{array}{c}\hfill \omega =k_{z}V\pm {\omega }_{\mathrm{k}},\text{with}{\omega }_{\mathrm{k}}=\sqrt{\frac{2k_z^2{B}_0^2}{\mu ({\rho }_{\mathrm{i}}+{\rho }_{\mathrm{e}})}}={\omega }_{\mathrm{Ae}}\sqrt{\frac{2}{\zeta +1}}.\end{array}$$(3)

We introduced the Alfvén frequency ω_A = k_zV_A (in connection with the Alfvén speed ${\mathit{V}}_{\mathrm{A}}={\mathit{B}}_0/\sqrt{\mu \rho }$, with magnetic permeability μ) and the density contrast of the coronal plume as ζ = ρ_i/ρ_e > 1, which we consider to be always greater than one.

We compute the velocity components (v_r, v_φ, v_z) using Eqs. (10)–(12) of Goossens et al. (1992):

$$\begin{array}{c}\hfill \delta V_r=\frac{1}{{\rho }_0({\mathrm{\Omega }}^2-{\omega }_{\mathrm{A}}^2)}\frac{\partial }{\partial r}\frac{\mathrm{d}P\prime }{\mathrm{d}t},\delta V_{\phi }=\frac{1}{{\rho }_0({\mathrm{\Omega }}^2-{\omega }_{\mathrm{A}}^2)}\frac{1}{r}\frac{\partial }{\partial \phi }\frac{\mathrm{d}P\prime }{\mathrm{d}t},\delta V_z=0,\end{array}$$(4)

where the latter statement is true because of the cold plasma limit. In these equations, the total derivative d/dt = ∂/∂t + V∂/∂z is the Lagrangian derivative. Moreover, Ω = ω − k_zV is the Doppler-shifted frequency. The expression for the total derivative is also the reason why we only consider propagating waves here (in contrast to Van Doorsselaere et al. 2020b); for standing waves, the time derivative and z derivative result in a change of z − t behaviour in both terms for standing waves.

With these notations, we can also relate the amplitude A of the total pressure variations with the amplitude Υ of the velocity perturbations, generalising what we obtain in Van Doorsselaere et al. (2020b):

$$\begin{array}{c}\hfill \mathrm{{\rm Y}}=\frac{A}{R}\frac{\mathrm{\Omega }}{{\rho }_0({\mathrm{\Omega }}^2-{\omega }_{\mathrm{A}}^2)},\end{array}$$(5)

where we now include the effect of background flow.

In what follows, we also need the expressions for the magnetic-field perturbation δB. In the remainder of this paper, we consider δB_z = 0, because Van Doorsselaere et al. (2020b) showed that its contributions are always of a higher order in k_zR (even though we took P′=B₀δB_z/μ ≠ 0 before). For the other magnetic field components, we linearised the induction equation (∂B₀/∂t = ∇ × (V × B₀)) for a flowing plasma. The result for the components perpendicular to the magnetic field is

$$\begin{array}{c}\hfill \frac{\partial {\mathit{\delta }\mathit{B}}_{\perp }}{\partial t}=B_0\frac{\partial }{\partial z}{\mathit{\delta }\mathit{V}}_{\perp }-V\frac{\partial }{\partial z}{\mathit{\delta }\mathit{B}}_{\perp }-{\mathit{\delta }\mathit{B}}_{\perp }\mathrm{\nabla }·\mathit{V}.\end{array}$$(6)

In the last term, we left the dependence on the divergence of V. In this cylindrical system, it is of course zero, but in the full model later on, it is non-zero. Still, we leave this term out here because it is smaller (in a WKB sense) than the wave variations in the first two terms on the RHS.

In the end, we generalised the expressions from Van Doorsselaere et al. (2020b) for the velocity and magnetic-field perturbations, and we find that the generalisation is limited to replacing ω by the Doppler-shifted frequency Ω:

$$\begin{array}{cc}\hfill \delta V_r& =\frac{\partial \mathcal{R}}{\partial r}\frac{\mathrm{\Omega }}{{\rho }_0({\mathrm{\Omega }}^2-{\omega }_{\mathrm{A}}^2)}cos\phi sin(k_{z}z-\omega t),\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill \delta V_{\phi }& =-\frac{\mathcal{R}}{r}\frac{\mathrm{\Omega }}{{\rho }_0({\mathrm{\Omega }}^2-{\omega }_{\mathrm{A}}^2)}sin\phi sin(k_{z}z-\omega t),\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill \delta B_r& =-\frac{\partial \mathcal{R}}{\partial r}\frac{k_z{B}_0}{{\rho }_0({\mathrm{\Omega }}^2-{\omega }_{\mathrm{A}}^2)}cos\phi sin(k_{z}z-\omega t),\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill \delta B_{\phi }& =\frac{\mathcal{R}}{r}\frac{k_z{B}_0}{{\rho }_0({\mathrm{\Omega }}^2-{\omega }_{\mathrm{A}}^2)}sin\phi sin(k_{z}z-\omega t),\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill \delta B_z& =0.\hfill \end{array}$$(11)

In Van Doorsselaere et al. (2020b); corresponding to the expression in Goossens et al. (2013), we computed with these variables (in a setting without flow) that the wave energy density is given as

$$\begin{array}{c}\hfill ⟨w⟩=\pi R^2\frac{{\rho }_{\mathrm{i}}+{\rho }_{\mathrm{e}}}{2}{\mathrm{{\rm Y}}}^2,\end{array}$$(12)

which is averaged over the cross-section and the wavelength. Here, the wave energy w is given as

$$\begin{array}{c}\hfill w^{\pm }=\frac{{\rho }_0{({\mathit{\delta }\mathit{Z}}^{\pm })}^2}{4},\end{array}$$(13)

and

$$\begin{array}{c}\hfill {\mathit{\delta }\mathit{Z}}^{\pm }=\mathit{\delta }\mathit{V}\pm \frac{\mathit{\delta }\mathit{B}}{\sqrt{\mu {\rho }_0}}\end{array}$$(14)

is the linearised version of the classical Elsässer variable Z = V ± V_A. We then computed w by

$$\begin{array}{c}\hfill w=w^{+}+w^{-}.\end{array}$$(15)

We summed both of these to cover the entirety of the wave energy, because both Elsässer components of the kink wave are non-zero (Magyar et al. 2019; Van Doorsselaere et al. 2020b).

For the wave pressure P_k, Van Doorsselaere et al. (2020b) found

$$\begin{array}{cc}\hfill P_{\mathrm{ki}}& =\frac{{\omega }_{\mathrm{Ai}}^2}{2{\rho }_{\mathrm{i}}{({\mathrm{\Omega }}^2-{\omega }_{\mathrm{Ai}}^2)}^2}{\left(\frac{A}{R}\right)}^2{sin}^2(k_{z}z-\omega t),\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill P_{\mathrm{ke}}& =\frac{{\omega }_{\mathrm{Ae}}^2}{2{\rho }_{\mathrm{e}}{({\mathrm{\Omega }}^2-{\omega }_{\mathrm{Ae}}^2)}^2}\left[\frac{A^2{R}^2}{r^4}\right]{sin}^2(k_{z}z-\omega t),\hfill \end{array}$$(17)

which we generalised immediately to have the Doppler-shifted frequency Ω instead. From these generalisations of the kink-wave functions of Van Doorsselaere et al. (2020b) to plumes that now include background flow, it seems that the influence of the flow is limited to replacing the wave frequency ω with the Doppler-shifted frequency Ω in all expressions.

3. Q variables and terminology

A second ingredient that is needed to formulate the UAWSOM model is the Q variables. These were introduced in Van Doorsselaere et al. (2024) as variables co-propagating with the waves’ phase speeds. They are defined as

$$\begin{array}{c}\hfill {\mathit{Q}}^{\pm }=\mathit{V}\pm \alpha \mathit{B},\end{array}$$(18)

where α is a parameter related to the phase speed of the wave. In Van Doorsselaere et al. (2024), we show that these Q variables are well suited to tracking any wave. Thus, we used them to construct extra equations for propagating kink waves in plume regions, on top of the existing Alfvén-wave equations in the AWSOM model. Constructing these extra equations is possible, because for Q variables, waves propagating in one direction only have non-zero components in the co-propagating Q component, generalising the convenience of the Elsässer variables for Alfvén waves to all other waves. In this section, we briefly list the necessary notation and equations that form the basis for constructing the energy evolution equation for kink waves in coronal plumes.

In our previous work, Van Doorsselaere et al. (2024), we show that the MHD equations are rewritten in terms of Q variables, with the following form of the momentum and induction equation:

$$\begin{array}{cc}\hfill \frac{D^{\mp }}{\mathit{Dt}}{\mathit{Q}}^{\pm }\mp \left(\frac{{\mathit{Q}}^{+}-{\mathit{Q}}^{-}}{4}\right)\frac{D^{\mp }}{\mathit{Dt}}ln\rho {\alpha }^2& =-v_{\mathrm{s}}^2\mathrm{\nabla }ln\rho -\frac{1}{8}(1-\frac{\mathrm{\Delta }{\alpha }^2}{{\alpha }^2})\mathrm{\nabla }{({\mathit{Q}}^{+}-{\mathit{Q}}^{-})}^2+\frac{1}{4}(1-\frac{\mathrm{\Delta }{\alpha }^2}{{\alpha }^2}){({\mathit{Q}}^{+}-{\mathit{Q}}^{-})}^2\mathrm{\nabla }ln\alpha \hfill \\ \hfill & -\frac{1}{4}\frac{\mathrm{\Delta }{\alpha }^2}{{\alpha }^2}({\mathit{Q}}^{+}-{\mathit{Q}}^{-})·\mathrm{\nabla }({\mathit{Q}}^{+}-{\mathit{Q}}^{-})+\frac{1}{4}\frac{\mathrm{\Delta }{\alpha }^2}{{\alpha }^2}({\mathit{Q}}^{+}-{\mathit{Q}}^{-})\mathrm{\nabla }·({\mathit{Q}}^{+}-{\mathit{Q}}^{-})\hfill \\ \hfill & \mp \left(\frac{{\mathit{Q}}^{+}-{\mathit{Q}}^{-}}{8}\right)\mathrm{\nabla }·(3{\mathit{Q}}^{\pm }-{\mathit{Q}}^{\mp })+\left(\frac{{\mathit{Q}}^{+}-{\mathit{Q}}^{-}}{4}\right)(\left(\frac{{\mathit{Q}}^{+}-{\mathit{Q}}^{-}}{2}\right)·\mathrm{\nabla }ln\rho {\alpha }^2).\hfill \end{array}$$(19)

In this equation, we employed the notation for the co-moving derivative $\frac{D^{\pm }}{\mathit{Dt}}=\frac{\partial }{\partial t}+{\mathit{Q}}^{\pm }·\mathrm{\nabla }$; the parameter $\mathrm{\Delta }{\alpha }^2={\alpha }^2-\frac{1}{\mu \rho }$ measures the deviation of the wave from an Alfvén wave; and ${\mathit{v}}_{\mathrm{s}}^2=\frac{\gamma p}{\rho }$ is the square of the sound speed for a plasma with pressure p. Using the expression for the continuity relation

$$\begin{array}{c}\hfill \frac{D^{\pm }}{\mathit{Dt}}(ln\rho )=-\frac{1}{2}\mathrm{\nabla }·(3{\mathit{Q}}^{\pm }-{\mathit{Q}}^{\mp })\pm \frac{{\mathit{Q}}^{+}-{\mathit{Q}}^{-}}{2}·\mathrm{\nabla }ln\rho {\alpha }^2,\end{array}$$(20)

here we use Eq. (19) in the form

$$\begin{array}{cc}\hfill \frac{D^{\mp }}{\mathit{Dt}}{\mathit{Q}}^{\pm }\mp \left(\frac{{\mathit{Q}}^{+}-{\mathit{Q}}^{-}}{4}\right)\frac{D^{\mp }}{\mathit{Dt}}ln\rho {\alpha }^2& =-v_{\mathrm{s}}^2\mathrm{\nabla }ln\rho -\frac{1}{8}(1-\frac{\mathrm{\Delta }{\alpha }^2}{{\alpha }^2})\mathrm{\nabla }{({\mathit{Q}}^{+}-{\mathit{Q}}^{-})}^2+\frac{1}{4}(1-\frac{\mathrm{\Delta }{\alpha }^2}{{\alpha }^2}){({\mathit{Q}}^{+}-{\mathit{Q}}^{-})}^2\mathrm{\nabla }ln\alpha \hfill \\ \hfill & -\frac{1}{4}\frac{\mathrm{\Delta }{\alpha }^2}{{\alpha }^2}({\mathit{Q}}^{+}-{\mathit{Q}}^{-})·\mathrm{\nabla }({\mathit{Q}}^{+}-{\mathit{Q}}^{-})+\frac{1}{4}\frac{\mathrm{\Delta }{\alpha }^2}{{\alpha }^2}({\mathit{Q}}^{+}-{\mathit{Q}}^{-})\mathrm{\nabla }·({\mathit{Q}}^{+}-{\mathit{Q}}^{-})\hfill \\ \hfill & \mp \left(\frac{{\mathit{Q}}^{+}-{\mathit{Q}}^{-}}{8}\right)(-2\frac{D^{\pm }}{\mathit{Dt}}(ln\rho )\pm ({\mathit{Q}}^{+}-{\mathit{Q}}^{-})·\mathrm{\nabla }ln\rho {\alpha }^2)+\left(\frac{{\mathit{Q}}^{+}-{\mathit{Q}}^{-}}{4}\right)(\left(\frac{{\mathit{Q}}^{+}-{\mathit{Q}}^{-}}{2}\right)·\mathrm{\nabla }ln\rho {\alpha }^2).\hfill \end{array}$$(21)

For what follows, it should be noted that

$$\begin{array}{c}\hfill {\mathit{Q}}^{+}-{\mathit{Q}}^{-}=2\alpha \mathit{B}\end{array}$$(22)

is in the direction of the magnetic field.

In our previous paper, Van Doorsselaere et al. (2024), we used these equations to describe surface Alfvén waves. We found that

$$\begin{array}{c}\hfill \alpha ={\alpha }_{\mathrm{L}}={\alpha }_{\mathrm{R}}=\sqrt{\frac{2}{\mu ({\rho }_{\mathrm{R}}+{\rho }_{\mathrm{L}})}},\end{array}$$(23)

which also helps to recover the kink frequency ω_k for such surface Alfvén waves. Given the correspondence between surface Alfvén waves and kink waves (especially in the long-wavelength limit, Goossens et al. 2012), we utilised these α values for kink waves in cylinders (as in Sect. 2) as well. We previously showed that such a choice of α makes it possible to decouple and filter out the waves of interest, which in our case are the kink waves.

4. Wave-energy equations

4.1. Linearised Q variable equations

Now we can construct an energy-evolution equation that will form the basis for the kink-wave equation in the UAWSOM model. To that end, we start from the linearised Q variable equations (Eq. (21)) to create an energy-evolution equation.

We consider how the wave energy evolves in a system that has the magnetic field B₀(z), the density ρ₀(z), and the wind velocity V(z) as functions of z and slowly evolving in time. The waves only have small perturbations on the background, and a classical linearisation of the system is suitable. This is not true higher up in the solar wind (e.g. above ∼3 − 5 solar radii), but (following van der Holst et al. 2014) we take this approximation nevertheless, because higher order (non-linear) terms are included in the system anyway. Thus, we consider

$$\begin{array}{c}\hfill {\mathit{Q}}^{\pm }={\mathit{Q}}_0^{\pm }+{\mathit{\delta }\mathit{Q}}^{\pm }\text{and}ln\rho \approx ln{\rho }_0+\delta R,ln\rho {\alpha }^2\approx \delta R+ln{\rho }_0{\alpha }^2,\end{array}$$(24)

where we introduce the notation δR = δρ/ρ₀. As in Van Doorsselaere et al. (2024), the phase-speed variable α is not linearised in this process. With the above linearisation, we obtain, for the perpendicular components of δQ_⊥^±,

$$\begin{array}{cc}\hfill \frac{D^{\mp }}{\mathit{Dt}}{\mathit{\delta }\mathit{Q}}_{\perp }^{\pm }+{\mathit{\delta }\mathit{Q}}_{\perp }^{\mp }·\mathrm{\nabla }{\mathit{Q}}_0^{\pm }& =-v_{\mathrm{s}}^2\mathrm{\nabla }\delta R\pm \frac{\alpha }{2}{\mathit{\delta }\mathit{B}}_{\perp }\frac{D^{\mp }}{\mathit{Dt}}ln{\rho }_0{\alpha }^2\pm \frac{\alpha }{2}{\mathit{\delta }\mathit{B}}_{\perp }\frac{D^{\mp }}{\mathit{Dt}}\delta R-\mathrm{\Delta }{\alpha }^2{\mathit{B}}_0·\mathrm{\nabla }{\mathit{\delta }\mathit{B}}_{\perp }-\mathrm{\Delta }{\alpha }^2{\mathit{\delta }\mathit{B}}_{\perp }·\mathrm{\nabla }{\mathit{\delta }\mathit{B}}_{\perp }\hfill \\ \hfill & -\frac{1}{{\rho }_0}\mathrm{\nabla }{P}_{\mathrm{w}}\mp \left(\frac{{\mathit{\delta }\mathit{Q}}_{\perp }^{+}-{\mathit{\delta }\mathit{Q}}_{\perp }^{-}}{8}\right)[-2\frac{D^{\mp }}{\mathit{Dt}}ln{\rho }_0-2\frac{D^{\mp }}{\mathit{Dt}}\delta R\pm ({\mathit{\delta }\mathit{Q}}_{\perp }^{+}-{\mathit{\delta }\mathit{Q}}_{\perp }^{-})·\mathrm{\nabla }ln{\rho }_0{\alpha }^2\mp ({\mathit{\delta }\mathit{Q}}_{\perp }^{+}-{\mathit{\delta }\mathit{Q}}_{\perp }^{-})·\mathrm{\nabla }\delta R\hfill \\ \hfill & \pm ({\mathit{Q}}_0^{+}-{\mathit{Q}}_0^{-})·\mathrm{\nabla }ln{\rho }_0{\alpha }^2]+\left(\frac{{\mathit{\delta }\mathit{Q}}_{\perp }^{+}-{\mathit{\delta }\mathit{Q}}_{\perp }^{-}}{8}\right)[({\mathit{Q}}_0^{+}-{\mathit{Q}}_0^{-})·\mathrm{\nabla }ln{\rho }_0{\alpha }^2\hfill \\ \hfill & +({\mathit{\delta }\mathit{Q}}_{\perp }^{+}-{\mathit{\delta }\mathit{Q}}_{\perp }^{-})·\mathrm{\nabla }ln{\rho }_0{\alpha }^2+({\mathit{\delta }\mathit{Q}}_{\perp }^{+}-{\mathit{\delta }\mathit{Q}}_{\perp }^{-})·\mathrm{\nabla }\delta R],\hfill \end{array}$$(25)

where we conveniently left some terms in δB_⊥ to shorten the notation. $P_{\mathrm{w}}=\frac{\delta B^2}{2\mu }$ is the wave pressure in the equation. It will become either the kink-wave pressure P_k or the Alfvén-wave pressure P_A.

Because the kink wave is nearly incompressible, we take all terms with δR to be zero, since the kink wave’s compression scales with k_z²R². Moreover, we assume that the perpendicular variations of equilibrium quantities (e.g. δQ_⊥^∓ ⋅ ∇Q₀^±) are also zero. This is in principle incorrect, because in the kink-wave model of Sect. 2 large density variations are present across the magnetic field. However, we only account for the effect of those gradients through the use of α and averaging across the cross-section of the coronal plume. The resulting equation is thus

$$\begin{array}{cc}\hfill \frac{D^{\mp }}{\mathit{Dt}}{\mathit{\delta }\mathit{Q}}_{\perp }^{\pm }& =\pm \frac{\alpha }{2}{\mathit{\delta }\mathit{B}}_{\perp }\frac{D^{\mp }}{\mathit{Dt}}ln{\rho }_0{\alpha }^2-\mathrm{\Delta }{\alpha }^2{\mathit{B}}_0·\mathrm{\nabla }{\mathit{\delta }\mathit{B}}_{\perp }-\mathrm{\Delta }{\alpha }^2{\mathit{\delta }\mathit{B}}_{\perp }·\mathrm{\nabla }{\mathit{\delta }\mathit{B}}_{\perp }\hfill \\ \hfill & -\frac{1}{{\rho }_0}\mathrm{\nabla }{P}_{\mathrm{w}}\mp \left(\frac{{\mathit{\delta }\mathit{Q}}_{\perp }^{+}-{\mathit{\delta }\mathit{Q}}_{\perp }^{-}}{8}\right)[-2\frac{D^{\mp }}{\mathit{Dt}}ln{\rho }_0\pm ({\mathit{Q}}_0^{+}-{\mathit{Q}}_0^{-})·\mathrm{\nabla }ln{\rho }_0{\alpha }^2]+\left(\frac{{\mathit{\delta }\mathit{Q}}_{\perp }^{+}-{\mathit{\delta }\mathit{Q}}_{\perp }^{-}}{8}\right)[({\mathit{Q}}_0^{+}-{\mathit{Q}}_0^{-})·\mathrm{\nabla }ln{\rho }_0{\alpha }^2].\hfill \end{array}$$(26)

It is sometimes convenient to take the limit of ${\alpha }^2\to \frac{1}{\mu \rho },\mathrm{\Delta }{\alpha }^2\to 0,{\mathit{\delta }\mathit{Q}}_{\perp }^{\pm }\to {\mathit{\delta }\mathit{Z}}^{\pm }$ to check the results. Here, we can use this limit to see that the previous equation coincides in that limit with Eq. (20) of van der Holst et al. (2014). Since we neglected perpendicular gradients of equilibrium quantities, we do not have their terms with gradients of u and B, but we did keep the term with the wave pressure.

The linear portion of the left term in Eq. (26) is

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}{\mathit{\delta }\mathit{Q}}_{\perp }^{\pm }+{\mathit{Q}}_0^{\mp }·\mathrm{\nabla }{\mathit{\delta }\mathit{Q}}_{\perp }^{\pm },\end{array}$$(27)

which shows that the wave propagates the speed Q₀^∓ = V ∓ αB₀. This indicates that αB₀ is the phase speed of the wave and that δQ_⊥⁺ belongs to the downward-propagating wave (if the magnetic field is pointing upward).

4.2. Energy evolution

Following the procedure of van der Holst et al. (2014), we now multiply Eq. (26) by ρ₀δQ_⊥^±/2. Then, we obtain an evolution equation for the Q density W^±:

$$\begin{array}{c}\hfill W^{\pm }=\frac{{\rho }_0}{4}{\left({\mathit{\delta }\mathit{Q}}_{\perp }^{\pm }\right)}^2.\end{array}$$(28)

This expression is seemingly different than the wave-energy density w (Eq. (13)), but it is only an apparent difference. In Appendix A, we show that

$$\begin{array}{c}\hfill W^{\pm }={\rho }_0\delta V_{\perp }^2,\end{array}$$(29)

and thus that the Q density is equal to the wave-energy density w if the wave energy is in equipartition.

We now used this Q density to construct an energy-evolution equation for the kink waves. As stated before, we multiply Eq. (26) by ρ₀δQ_⊥^±/2. We then obtain

$$\begin{array}{c}\hfill \frac{D^{\mp }}{\mathit{Dt}}{W}^{\pm }-W^{\pm }\frac{D^{\mp }}{\mathit{Dt}}ln{\rho }_0-\frac{1}{2}(W^{\pm }-CC)\frac{D^{\pm }}{\mathit{Dt}}ln{\rho }_0=-\mathrm{\Delta }{\alpha }^2\frac{{\rho }_0}{2}{\mathit{\delta }\mathit{Q}}_{\perp }^{\pm }·({\mathit{B}}_0·\mathrm{\nabla }{\mathit{\delta }\mathit{B}}_{\perp })-\mathrm{\Delta }{\alpha }^2\frac{{\rho }_0}{2}{\mathit{\delta }\mathit{Q}}_{\perp }^{\pm }·({\mathit{\delta }\mathit{B}}_{\perp }·\mathrm{\nabla }{\mathit{\delta }\mathit{B}}_{\perp })+\frac{1}{2}(W^{\pm }-CC)\frac{D^{\mp }}{\mathit{Dt}}ln{\rho }_0{\alpha }^2,\end{array}$$(30)

where we introduced the shorthand $CC={\rho }_0\frac{{\mathit{\delta }\mathit{Q}}_{\perp }^{+}·{\mathit{\delta }\mathit{Q}}_{\perp }^{-}}{4}$ for the cross-correlation between upward- and downward-propagating waves. We grouped the terms not involving Δα² or lnρ₀α² (basically, the effects of Alfvén waves, as was already incorporated in AWSOM) to the left. The term with the wave pressure was dropped at this stage because (1) its net contribution is in the z direction, and (2) the time and spatially averaged term is added to the momentum equation of the background (Eq. (59)). In any case, the radial contribution of the pressure term has previously been averaged out.

We rewrite the left side of the previous equation following the AWSOM lead of van der Holst et al. (2014), using the equation for the conservation of mass of the background plasma,

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}ln{\rho }_0+\mathit{V}·\mathrm{\nabla }ln{\rho }_0+\mathrm{\nabla }·\mathit{V}=0,\end{array}$$(31)

to obtain

$$\begin{array}{cc}\hfill \frac{D^{\mp }}{\mathit{Dt}}{W}^{\pm }-W^{\pm }\frac{D^{\mp }}{\mathit{Dt}}& ln{\rho }_0-\frac{1}{2}(W^{\pm }-CC)\frac{D^{\pm }}{\mathit{Dt}}ln{\rho }_0=\hfill \\ \hfill & \frac{\partial W^{\pm }}{\partial t}+\mathrm{\nabla }·(({\mathit{Q}}_0^{\mp }+{\mathit{\delta }\mathit{Q}}_{\perp }^{\mp })W^{\pm })+\frac{1}{2}CC\frac{D^{\pm }ln{\rho }_0}{\mathit{Dt}}+W^{\pm }(\frac{1}{2}\mathrm{\nabla }·\mathit{V}\pm \frac{1}{2}\alpha {\mathit{B}}_0·\mathrm{\nabla }ln{\rho }_0{\alpha }^2),\hfill \end{array}$$(32)

where we once again consider the incompressibility of the wave mode and perpendicular gradients of background quantities to be zero. Finally, we obtain the following as the evolution equation for the Q density:

$$\begin{array}{cc}\hfill \frac{\partial W^{\pm }}{\partial t}& +\mathrm{\nabla }·(({\mathit{Q}}_0^{\mp }+{\mathit{\delta }\mathit{Q}}_{\perp }^{\mp })W^{\pm })+\frac{1}{2}CC\frac{D^{\pm }ln{\rho }_0}{\mathit{Dt}}+\frac{W^{\pm }}{2}\mathrm{\nabla }·\mathit{V}=\frac{W^{\pm }}{2}(\frac{D^{\mp }}{\mathit{Dt}}\mp \alpha {\mathit{B}}_0·\mathrm{\nabla })ln{\rho }_0{\alpha }^2-\frac{\mathit{CC}}{2}\frac{D^{\mp }}{\mathit{Dt}}ln{\rho }_0{\alpha }^2\hfill \\ \hfill & -\mathrm{\Delta }{\alpha }^2\frac{{\rho }_0}{2}{\mathit{\delta }\mathit{Q}}_{\perp }^{\pm }·({\mathit{B}}_0·\mathrm{\nabla }{\mathit{\delta }\mathit{B}}_{\perp })-\mathrm{\Delta }{\alpha }^2\frac{{\rho }_0}{2}{\mathit{\delta }\mathit{Q}}_{\perp }^{\pm }·({\mathit{\delta }\mathit{B}}_{\perp }·\mathrm{\nabla }{\mathit{\delta }\mathit{B}}_{\perp }).\hfill \end{array}$$(33)

The second term on the right side of this equation in principle groups with the third term on the left. However, this way of writing is more elegant, because the terms on the right describe how the wave-energy-density evolution is modified if $\alpha \ne 1/\sqrt{\mu \rho }$; that is, a wave other than an Alfvén wave is considered.

The rightmost term in this equation describes the energy dissipation by uniturbulence. There is self-interaction of the wave (any wave) if Δα² ≠ 0, leading to a cascade of energy in a turbulent way (Ismayilli et al. 2024).

4.2.1. Characterisation of the uniturbulence term

We further specialise these equations for single kink waves; that is, only taking into account the effect of the kink waves on themselves. The terms with lnρ₀α² will cancel out, because the α for the kink wave (Eq. (23)) shows that it is still proportional to the density (interior or exterior does not matter here). The only variation in that term would be due to spatial and temporal variations of the density contrast ζ, but we will take these to be zero in this paper, assuming that the density contrast and filling factor remain constant in the entire domain. In Sect. 6, the parameter ζ will in fact change with height, but its change is very small indeed, so the approximation here is fine.

For working out the uniturbulence terms, we considered that there is only one of δQ_⊥^±, in order to have only the self-interaction (Magyar et al. 2019). Then, we used Eq. (22) to write

$$\begin{array}{c}\hfill {\mathit{\delta }\mathit{Q}}_{\perp }^{\pm }=\pm 2\alpha {\mathit{\delta }\mathit{B}}_{\perp }.\end{array}$$(34)

With this equation, we rewrote

$$\begin{array}{c}\hfill -\mathrm{\Delta }{\alpha }^2\frac{{\rho }_0}{2}{\mathit{\delta }\mathit{Q}}_{\perp }^{\pm }({\mathit{\delta }\mathit{B}}_{\perp }·\mathrm{\nabla }{\mathit{\delta }\mathit{B}}_{\perp })=\mp \mathrm{\Delta }{\alpha }^2\mu {\rho }_0\alpha {\mathit{\delta }\mathit{B}}_{\perp }·\mathrm{\nabla }{P}_{\mathrm{k}}.\end{array}$$(35)

In the exterior region of the plume, the latter expression reduced to

$$\begin{array}{c}\hfill \mp \frac{1-\zeta }{1+\zeta }\alpha {\mathit{\delta }\mathit{B}}_{\perp }·\mathrm{\nabla }{P}_{\mathrm{k}},\end{array}$$(36)

where we used the expression for α from Eq. (23). Using Eqs. (7)–(11) (or Eqs. (40) & (44) in Van Doorsselaere et al. 2020b), we obtained that δB_⊥i ⋅ ∇P_ki = 0 up to leading orders of k_zR. The remaining contribution for uniturbulent damping was

$$\begin{array}{c}\hfill {\mathit{\delta }\mathit{B}}_{\perp \mathrm{e}}·\mathrm{\nabla }{P}_{\mathrm{ke}}=-2\frac{A^3{R}^3}{r^7}\frac{k_z{B}_0{\omega }_{\mathrm{Ae}}^2}{{\rho }_{\mathrm{e}}^2{({\mathrm{\Omega }}^2-{\omega }_{\mathrm{Ae}}^2)}^3}cos\phi {sin}^3(k_{z}z-\omega t).\end{array}$$(37)

For the radial averaging, we integrate from R to γR as for the energy (as we also did in Goossens et al. 2013; Van Doorsselaere et al. 2014, 2020b, see also the procedure in Appendix A for an example). For that integral, we have

$$\begin{array}{c}\hfill {\displaystyle {\int }_R^{\gamma R}}\frac{1}{r^7}r\mathrm{d}r=\frac{1}{5R^5}\left(\frac{{\gamma }^5-1}{{\gamma }^5}\right)=\frac{1}{5R^5}(1-f^{5/2}).\end{array}$$(38)

As discussed in Van Doorsselaere et al. (2020b), the expression Eq. (37) is an odd-powered periodic function in φ and time, and its average would result in zero. Thus, it makes sense to take the RMS average of this term to quantify the uniturbulent damping. The physical interpretation is that all small scales are cascaded to even smaller scales before they experience the inverse cascade in the second half of the period. From Van Doorsselaere et al. (2020b), we know that the RMS average will yield a $\sqrt{\pi }$ from the integral over φ, and the RMS average over t yields a factor of $\sqrt{5/16}$. We obtained the following as result:

$$\begin{array}{c}\hfill -\frac{\zeta -1}{\zeta +1}\sqrt{\pi }\sqrt{\frac{5}{16}}\frac{2A^3}{5R^2}\frac{\alpha k_z{B}_0{\omega }_{\mathrm{Ae}}^2}{{\rho }_{\mathrm{e}}^2{({\mathrm{\Omega }}^2-{\omega }_{\mathrm{Ae}}^2)}^3}(1-f^{5/2}),\end{array}$$(39)

where we only kept the minus sign to extract energy from the equation. Further substituting the relevant expressions for Υ (Eq. (5)), we then write

$$\begin{array}{c}\hfill -\frac{\zeta -1}{\zeta +1}\sqrt{\frac{\pi }{20}}\alpha k_z{B}_{0}R\frac{{\rho }_{\mathrm{e}}{\omega }_{\mathrm{Ae}}^2{\mathrm{{\rm Y}}}^3}{{\mathrm{\Omega }}^3}(1-f^{5/2}),\end{array}$$(40)

which was further simplified with the expression for Ω = αk_zB₀ to

$$\begin{array}{c}\hfill -\frac{\zeta -1}{\zeta +1}\sqrt{\frac{\pi }{20}}R\frac{{\rho }_{\mathrm{e}}{\omega }_{\mathrm{Ae}}^2{\mathrm{{\rm Y}}}^3}{{\mathrm{\Omega }}^2}(1-f^{5/2}).\end{array}$$(41)

Here, we once again used $\mathrm{\Omega }=\pm \alpha k_z{B}_0=\pm {\omega }_{\mathrm{Ae}}\sqrt{\frac{2}{\zeta +1}}$ to further simplify to

$$\begin{array}{c}\hfill -\frac{\zeta -1}{4}\sqrt{\frac{\pi }{5}}R{\rho }_{\mathrm{e}}{\mathrm{{\rm Y}}}^3(1-f^{5/2}).\end{array}$$(42)

This expression was equal to the damping rate found in Van Doorsselaere et al. (2020b), confirming the earlier result with the use of the new Q variables. For implementation in a numerical code, it is most convenient to rewrite this term as an expression of the cross-section and wavelength-averaged energy density ⟨W^±⟩ (see Van Doorsselaere et al. 2014, 2020b, and Appendix A for more details):

$$\begin{array}{c}\hfill ⟨W^{\pm }⟩=\frac{k_z}{2\pi }{\displaystyle {\int }_z}\mathrm{d}z{\displaystyle {\int }_{\phi }}\mathrm{d}\phi {\displaystyle {\int }_0^{\gamma R}}{W}^{\pm }r\mathrm{d}r=\pi R^2{\mathrm{{\rm Y}}}^2\frac{{\rho }_{\mathrm{i}}+{\rho }_{\mathrm{e}}(1-f)}{2},\end{array}$$(43)

to

$$\begin{array}{c}\hfill -{⟨W^{\pm }⟩}^{3/2}\frac{1}{\pi R^2}\frac{1}{\sqrt{{\rho }_{\mathrm{e}}}}\sqrt{\frac{2}{5}}\frac{\zeta -1}{2}\frac{1-f^{5/2}}{{(\zeta +1-f)}^{3/2}}.\end{array}$$(44)

4.2.2. Characterisation of the B₀ ⋅ ∇δB_⊥ term

In the previous subsection, we simplify the energy cascade term for the uniturbulence damping mechanism. We now simplified the remaining, unknown terms in Eq. (33). To that end, we use Eq. (34) to simplify the term in Eq. (33) containing B₀ ⋅ ∇δB_⊥:

$$\begin{array}{c}\hfill -\mathrm{\Delta }{\alpha }^2\frac{{\rho }_0}{2}{\mathit{\delta }\mathit{Q}}_{\perp }^{\pm }·({\mathit{B}}_0·\mathrm{\nabla }{\mathit{\delta }\mathit{B}}_{\perp })=\mp \mathrm{\Delta }{\alpha }^2\frac{{\rho }_0}{2}2\alpha {\mathit{\delta }\mathit{B}}_{\perp }·({\mathit{B}}_0·\mathrm{\nabla }{\mathit{\delta }\mathit{B}}_{\perp })=\mp \mathrm{\Delta }{\alpha }^2\mu {\rho }_0\alpha {\mathit{B}}_0·\mathrm{\nabla }\left(\frac{{\mathit{\delta }\mathit{B}}_{\perp }^2}{2\mu }\right)=\mp \mathrm{\Delta }{\alpha }^2\mu {\rho }_0\alpha {\mathit{B}}_0·\mathrm{\nabla }{P}_{\mathrm{k}}.\end{array}$$(45)

We previously (Eq. (36)) used that term in the region exterior to the plume, which we have as value of $\mu {\rho }_e\mathrm{\Delta }{\alpha }_{\mathrm{e}}^2=\frac{1-\zeta }{1+\zeta }$. However, to quantify the term in question, we also need the expression in the interior region:

$$\begin{array}{c}\hfill \mu {\rho }_{\mathrm{i}}\mathrm{\Delta }{\alpha }_{\mathrm{i}}^2=\frac{\zeta -1}{\zeta +1};\end{array}$$(46)

that is, it is the opposite of the exterior value.

Next, we need to average the contributing term over the entire cross-section, as we did for the uniturbulence term in Sect. 4.2.1. We must evaluate

$$\begin{array}{c}\hfill \mp \{{\displaystyle {\int }_0^R}r\mathrm{d}r\mu {\rho }_{\mathrm{i}}\mathrm{\Delta }{\alpha }_{\mathrm{i}}^2\alpha {\mathit{B}}_0·\mathrm{\nabla }{P}_{\mathrm{ki}}+{\displaystyle {\int }_R^{\gamma R}}r\mathrm{d}r\mu {\rho }_{\mathrm{e}}\mathrm{\Delta }{\alpha }_{\mathrm{e}}^2\alpha {\mathit{B}}_0·\mathrm{\nabla }{P}_{\mathrm{ke}}\}.\end{array}$$(47)

Using the results of Eqs. (36) and (46), we then have

$$\begin{array}{c}\hfill \mp \left(\frac{\zeta -1}{\zeta +1}\right)\{{\displaystyle {\int }_0^R}r\mathrm{d}r\alpha {\mathit{B}}_0·\mathrm{\nabla }{P}_{\mathrm{ki}}-{\displaystyle {\int }_R^{\gamma R}}r\mathrm{d}r\alpha {\mathit{B}}_0·\mathrm{\nabla }{P}_{\mathrm{ke}}\}.\end{array}$$(48)

Now we use the fact that the field-aligned gradient will be equal for the internal and external pressure in the WKB approximation, so it can be moved out of the integrals. We must then calculate

$$\begin{array}{c}\hfill \mp \left(\frac{\zeta -1}{\zeta +1}\right)\alpha {\mathit{B}}_0·\mathrm{\nabla }\{{\displaystyle {\int }_0^R}{P}_{\mathrm{ki}}r\mathrm{d}r-{\displaystyle {\int }_R^{\gamma R}}{P}_{\mathrm{ke}}r\mathrm{d}r\}.\end{array}$$(49)

With the aid of the Expressions (16) and (17) and Eqs. (43) and (46) in Van Doorsselaere et al. (2020b), we obtain

$$\begin{array}{c}\hfill \mp \left(\frac{\zeta -1}{\zeta +1}\right)\alpha {\mathit{B}}_0·\mathrm{\nabla }\left\{\frac{{\mathrm{{\rm Y}}}^2\pi R^2}{2{\omega }^2}({\rho }_{\mathrm{i}}{\omega }_{\mathrm{Ai}}^2-{\rho }_{\mathrm{e}}{\omega }_{\mathrm{Ae}}^2(1-f))\right\}.\end{array}$$(50)

When we take the filling factor to be small (i.e. f ≪ 1), thus, the term in the curly brackets represents exactly zero.

4.2.3. Evolution equation for cross-sectionally averaged wave energy

After simplifying all terms in the energy-evolution equation (Eq. (33)), we can gather all the previous terms and their simplified expression to achieve our goal of obtaining an evolution equation for the kink-wave-energy density. With Expressions (44) and (50), we can write the final equation for the kink-wave-energy evolution as

$$\begin{array}{c}\hfill \frac{\partial ⟨W^{\pm }⟩}{\partial t}+\mathrm{\nabla }·({\mathit{Q}}_0^{\mp }⟨W^{\pm }⟩)+\frac{⟨W^{\pm }⟩}{2}\mathrm{\nabla }·\mathit{V}=-{⟨W^{\pm }⟩}^{3/2}\frac{1}{\pi R^2}\frac{1}{\sqrt{{\rho }_{\mathrm{e}}}}\sqrt{\frac{2}{5}}\frac{\zeta -1}{2}\frac{1-f^{5/2}}{{(\zeta +1-f)}^{3/2}}.\end{array}$$(51)

This equation has inconvenient units to add up with the traditional expression for the Alfvén-wave energy density, since ⟨W^±⟩ has been integrated over the cross-section of the magnetic flux tube. Thus, it is more convenient to once again normalise it to the relevant cross-sectional area πγ²R². We therefore followed the procedure of Goossens et al. (2013) to compute the kink-wave-energy density. We define

$$\begin{array}{c}\hfill W_{\mathrm{k}}^{\pm }=\frac{⟨W^{\pm }⟩}{\pi {\gamma }^2{R}^2}\end{array}$$(52)

as the cross-sectionally averaged energy density of the kink waves to obtain the following as the evolution equation for the kink-wave-energy density:

$$\begin{array}{c}\hfill \frac{\partial W_{\mathrm{k}}^{\pm }}{\partial t}+\mathrm{\nabla }·({\mathit{Q}}_0^{\mp }{W}_{\mathrm{k}}^{\pm })+\frac{W_{\mathrm{k}}^{\pm }}{2}\mathrm{\nabla }·\mathit{V}=-\frac{1}{L_{\perp ,\mathrm{VD}}}\frac{1}{\sqrt{{\rho }_{\mathrm{e}}}}{(W_{\mathrm{k}}^{\pm })}^{3/2},\end{array}$$(53)

where we define the appropriate length for the damping of the waves as

$$\begin{array}{c}\hfill L_{\perp ,\mathrm{VD}}(R,\zeta ,f)={\left(\frac{1}{R\sqrt{f\pi }}\sqrt{\frac{2}{5}}\frac{\zeta -1}{2}\frac{1-f^{5/2}}{{(\zeta +1-f)}^{3/2}}\right)}^{-1}.\end{array}$$(54)

One may say that the introduction of such a length is no better than the infamous L_⊥ parameter for the turbulent damping of Alfvén waves. The difference is that the Alfvén turbulence parameter is solely based on a phenomenological approach, but here we offer a full description of it for uniturbulence heating. Given that R, f, and ζ all vary with distance along a magnetic-field line, it is clear that L_⊥, VD will also vary with height in the atmosphere. A reasonable prescription for this parameter can be obtained from initial atmospheric models (e.g. Sishtla et al. 2022). In the atmospheric model we use in Sect. 6, the value ranges from 6.3 Mm–7.3 Mm, which is similar to values for Alfvén waves (e.g. Sharma & Morton 2023 and references therein).

As an alternative, one may consider the theory of Hillier et al. (2019) to quantify the energy cascade rate in the saturated regime. In that case, a new expression would emerge, but it would come down to replacing L_{⊥, Hillier} in the previous equation.

An ad hoc reflection term should be added to Eq. (33). We name these terms ℛ_A,k for the respective waves. One could take the same reflection as in van der Holst et al. (2014) for the Alfvén and kink waves alike. Or, we may consider the reflection expression from Réville et al. (2020). One could also opt to take a different reflection of the Alfvén and kink waves.

5. Energy equations for the background

It was necessary to incorporate the kink-wave-energy-evolution equation (Eq. (53)) into the larger system of AWSOM equations. For this, we needed to model the slow background variation through regular MHD, but we supplemented the momentum equation with the forces exerted by the wave pressures P_k and P_A. To express the conservation of energy, we needed to compute the equation for the conservation of wave energy. To that end, we summed Eq. (33) for the upward and downward Alfvén wave and Eq. (53) for the upward and downward kink wave. Once again following van der Holst et al. (2014), we also took CC = 0. We obtain

$$\begin{array}{cc}& \frac{\partial }{\partial t}(W_{\mathrm{A}}^{+}+W_{\mathrm{A}}^{-}+W_{\mathrm{k}}^{+}+W_{\mathrm{k}}^{-})+\mathrm{\nabla }·({\mathit{Z}}_0^{-}{W}_{\mathrm{A}}^{+}+{\mathit{Z}}_0^{+}{W}_{\mathrm{A}}^{-}+{\mathit{Q}}_0^{-}{W}_{\mathrm{k}}^{+}+{\mathit{Q}}_0^{+}{W}_{\mathrm{k}}^{-})+\frac{1}{2}(W_{\mathrm{A}}^{+}+W_{\mathrm{A}}^{-}+W_{\mathrm{k}}^{+}+W_{\mathrm{k}}^{-})\mathrm{\nabla }·{\mathit{V}}_0\hfill \\ \hfill & =-{\mathrm{\Gamma }}^{+}{W}_{\mathrm{A}}^{+}-{\mathrm{\Gamma }}^{-}{W}_{\mathrm{A}}^{-}-\frac{1}{L_{\perp ,\mathrm{VD}}}\frac{1}{\sqrt{{\rho }_{\mathrm{e}}}}({(W_{\mathrm{k}}^{+})}^{3/2}+{(W_{\mathrm{k}}^{-})}^{3/2}),\hfill \end{array}$$(55)

where Z₀^± = V₀ ± V_A for the Alfvén waves and Q₀^± is the kink-phase speed in the moving frame. The term containing the uniturbulent damping is redefined to include the cross-sectionally averaged wave energy density. The terms with Γ^± are the terms of the Alfvén-wave turbulence (van der Holst et al. 2014) containing the infamous L_⊥ parameter. These terms originate in the ∇ ⋅ δQ_⊥^∓W^± terms of Eq. (33) that express the cascade due to interaction with counter-propagating waves. Thus far, we have assumed that counter-propagating kink waves do not interact (and thus that the corresponding Γ terms are 0). Indeed, it is likely that the cascade due to their interaction is much smaller than the uniturbulent self-cascade (Eq. (53)). The last term was combined with the terms for the work of the wave-pressure forces using the expression in van der Holst et al. (2014):

$$\begin{array}{c}\hfill P_{\mathrm{A}}=\frac{W_{\mathrm{A}}^{+}+W_{\mathrm{A}}^{-}}{2},\text{and}{P}_{\mathrm{k}}=\frac{W_{\mathrm{k}}^{+}+W_{\mathrm{k}}^{-}}{2\mu {\rho }_0{\alpha }^2}.\end{array}$$(56)

The terms on the right side of Eq. (55) should not be taken into account in the integrated-energy equation of MHD and the waves. Indeed, the cascade energy (via Alfvén-wave turbulence or uniturbulence) was added to the system as heat. Thus, those damping terms would be cancelled out with the heating terms.

We also took all cross-interactions to be zero; for instance, terms such as δZ⁺ ⋅ ∇W_k⁺. We could make some Fourier-based arguments as to why these do not contribute, but Guo et al. (2019) already showed that the presence of Alfvén waves does in fact facilitate the cascade of (admittedly standing) kink waves.

We now formulate the conservation equations for the entire system. We start with the conservation of mass and the induction equation:

$$\begin{array}{cc}\hfill \frac{\partial {\rho }_0}{\partial t}+\mathrm{\nabla }·({\rho }_0\mathit{V})& =0,\hfill \end{array}$$(57)

$$\begin{array}{cc}\hfill \frac{\partial {\mathit{B}}_0}{\partial t}-\mathrm{\nabla }\times (\mathit{V}\times {\mathit{B}}_0)& =0,\hfill \end{array}$$(58)

which are unchanged from straight MHD. Then, we formulate the conservation of momentum:

$$\begin{array}{c}\hfill \frac{\partial {\rho }_0\mathit{V}}{\partial t}+\mathrm{\nabla }·({\rho }_0\mathit{V}\mathit{V}-\frac{1}{\mu }{\mathit{B}}_0{\mathit{B}}_0)+\mathrm{\nabla }(p+\frac{B_0^2}{2\mu }+P_{\mathrm{A}}+P_{\mathrm{k}})=-{\rho }_0\mathit{g}(\mathit{r}),\end{array}$$(59)

where in addition to the traditional terms (including the gravity g, depending on the radial distance from the Sun), we included the wave pressures P_A and P_k.

As the equation for internal energy, we took

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\frac{p}{\gamma -1}+\mathrm{\nabla }·\frac{p\mathit{V}}{\gamma -1}+p\mathrm{\nabla }·\mathit{V}=-\mathcal{L}+{\mathcal{H}}_{\mathrm{A}+\mathrm{k}},\end{array}$$(60)

where we incorporated radiative losses ℒ and the heating ${\mathcal{H}}_{\mathrm{A}+\mathrm{k}}={\mathrm{\Gamma }}^{+}{W}_{\mathrm{A}}^{+}+{\mathrm{\Gamma }}^{-}{W}_{\mathrm{A}}^{-}+\frac{1}{L_{\perp ,\mathrm{VD}}(R,\zeta ,f)}\frac{1}{\sqrt{{\rho }_{\mathrm{e}}}}({(W_{\mathrm{k}}^{+})}^{3/2}+{(W_{\mathrm{k}}^{-})}^{3/2})$ by the Alfvén and kink waves. To obtain the final energy equation, we added Eq. (55), Eq. (60), V⋅Eq. (59), and B⋅Eq. (58), resulting in

$$\begin{array}{cc}& \frac{\partial }{\partial t}({\rho }_0\frac{V^2}{2}+\frac{p}{\gamma -1}+\frac{B_0^2}{2\mu }+{\displaystyle \sum W_{\mathrm{A},\mathrm{k}}^{\pm }})+\mathrm{\nabla }·([{\rho }_0\frac{V^2}{2}+\frac{p}{\gamma -1}+\frac{B_0^2}{2\mu }]\mathit{V}-{\mathit{B}}_0\frac{\mathit{V}·{\mathit{B}}_{\mathbf{0}}}{\mu })\hfill \\ \hfill & +\mathrm{\nabla }·({\mathit{Z}}_0^{-}{W}_{\mathrm{A}}^{+}+{\mathit{Z}}_0^{+}{W}_{\mathrm{A}}^{-}+{\mathit{Q}}_0^{-}{W}_{\mathrm{k}}^{+}+{\mathit{Q}}_0^{+}{W}_{\mathrm{k}}^{-})+\mathrm{\nabla }·\left([P_{\mathrm{A}}+P_{\mathrm{k}}]\mathit{V}\right)+(\mu {\rho }_0{\alpha }^2-1)P_{\mathrm{k}}\mathrm{\nabla }·\mathit{V}=-\mathcal{L}-{\rho }_0\mathit{V}·\mathit{g}(\mathit{r}),\hfill \end{array}$$(61)

where we used the shorthand ∑W_A, k^± = W_A⁺ + W_A⁻ + W_k⁺ + W_k⁻. Moreover, we used the expression for P_k in terms of the Q densities:

$$\begin{array}{c}\hfill P_{\mathrm{k}}=\frac{\delta B^2}{2\mu }=\frac{1}{\mu {\rho }_0{\alpha }^2}\frac{W_{\mathrm{k}}^{+}+W_{\mathrm{k}}^{-}}{2}=\frac{1+\zeta }{4}(W_{\mathrm{k}}^{+}+W_{\mathrm{k}}^{-}),\end{array}$$(62)

in the case where CC = 0. The factor (μρ₀α² − 1) is present in this equation because of the non-coincidence between the wave pressure P_k and the average of the wave energies. It reduced to $\frac{1-\zeta }{1+\zeta }$ for kink waves. For an Alfvén wave, this term reduced to zero, because the factor (μρ₀α² − 1) is zero when α² = 1/μρ₀; that is, for an Alfvén wave. The resulting extra term contained the density contrast, which varies with height in the solar atmosphere. Thus, a model parameter also appeared in the equations that need to be fixed ad hoc, just as L_⊥, VD.

Finally, the full set of equations can be written as

$$\begin{array}{cc}\hfill \frac{\partial {\rho }_0}{\partial t}+\mathrm{\nabla }·({\rho }_0\mathit{V})& =0,\hfill \end{array}$$(63)

$$\begin{array}{cc}\hfill \frac{\partial {\mathit{B}}_0}{\partial t}-\mathrm{\nabla }\times (\mathit{V}\times {\mathit{B}}_0)& =0,\hfill \end{array}$$(64)

$$\begin{array}{cc}& \frac{\partial {\rho }_0\mathit{V}}{\partial t}+\mathrm{\nabla }·({\rho }_0\mathit{V}\mathit{V}-\frac{1}{\mu }{\mathit{B}}_0{\mathit{B}}_0)+\mathrm{\nabla }(p+\frac{B_0^2}{2\mu }+P_{\mathrm{A}}+P_{\mathrm{k}})=-{\rho }_0\mathit{g}(\mathit{r}),\hfill \end{array}$$(65)

$$\begin{array}{cc}& \frac{\partial }{\partial t}({\rho }_0\frac{V^2}{2}+\frac{p}{\gamma -1}+\frac{B_0^2}{2\mu }+{\displaystyle \sum W_{\mathrm{A},\mathrm{k}}^{\pm }})+\mathrm{\nabla }·([{\rho }_0\frac{V^2}{2}+\frac{p}{\gamma -1}+\frac{B_0^2}{2\mu }]\mathit{V}-{\mathit{B}}_0\frac{\mathit{V}·{\mathit{B}}_{\mathbf{0}}}{\mu })\hfill \\ \hfill & +\mathrm{\nabla }·({\mathit{Z}}_0^{-}{W}_{\mathrm{A}}^{+}+{\mathit{Z}}_0^{+}{W}_{\mathrm{A}}^{-}+{\mathit{Q}}_0^{-}{W}_{\mathrm{k}}^{+}+{\mathit{Q}}_0^{+}{W}_{\mathrm{k}}^{-})+\mathrm{\nabla }·\left([P_{\mathrm{A}}+P_{\mathrm{k}}]\mathit{V}\right)=-\mathcal{L}-{\rho }_0\mathit{V}·\mathit{g}(\mathit{r})+\frac{\zeta -1}{\zeta +1}{P}_{\mathrm{k}}\mathrm{\nabla }·\mathit{V},\hfill \end{array}$$(66)

$$\begin{array}{cc}& \frac{\partial W_{\mathrm{A}}^{\pm }}{\partial t}+\mathrm{\nabla }·({\mathit{Z}}_0^{\mp }{W}_{\mathrm{A}}^{\pm })+\frac{W_{\mathrm{A}}^{\pm }}{2}\mathrm{\nabla }·\mathit{V}=-{\mathrm{\Gamma }}^{\mp }{W}_{\mathrm{A}}^{\pm }\mp {\mathcal{R}}_{\mathrm{A}},\hfill \end{array}$$(67)

$$\begin{array}{cc}& \frac{\partial W_{\mathrm{k}}^{\pm }}{\partial t}+\mathrm{\nabla }·({\mathit{Q}}_0^{\mp }{W}_{\mathrm{k}}^{\pm })+\frac{W_{\mathrm{k}}^{\pm }}{2}\mathrm{\nabla }·\mathit{V}=-\frac{1}{L_{\perp ,\mathrm{VD}}}\frac{1}{\sqrt{{\rho }_{\mathrm{e}}}}{(W_{\mathrm{k}}^{\pm })}^{3/2}\mp {\mathcal{R}}_{\mathrm{k}}.\hfill \end{array}$$(68)

In all of these equations, ρ₀ is actually ρ_e. When averaging the density over a cross-section, we have

$$\frac{1}{\pi {\gamma }^2{R}^2}{\displaystyle \int \mathrm{d}}\varphi {\displaystyle {\int }_0^{\gamma R}}\rho r\mathrm{d}r=\frac{1}{\pi {\gamma }^2{R}^2}(\pi R^2{\rho }_{\mathrm{i}}+\pi ({\gamma }^2{R}^2-R^2){\rho }_{\mathrm{e}})=f{\rho }_{\mathrm{i}}+(1-f){\rho }_{\mathrm{e}}={\rho }_{\mathrm{e}}+f({\rho }_{\mathrm{i}}-{\rho }_{\mathrm{e}}).$$

This averaging over the cross-section of the density would occur in the continuity equation and the non-wave pressure terms of the momentum equation. This understanding that the modelled density ρ₀ is really the external density ρ_e has no impact on the modelling that we performed (see the next section). However, appropriate values must be taken. Moreover, additional care is needed when implementing radiative losses or doing forward modelling, both of which are non-linearly weighed with the density.

In Eqs. (67)–(68) we included ad hoc reflection terms. It is important that these reflection terms adhere to the basic principles of conservation of energy; that is, what is taken out of the +-equation must be injected as a source term in the −-equation, and vice versa. Several expressions for the reflection of Alfvén waves exist in the literature (e.g. van der Holst et al. 2014; Downs et al. 2016; Réville et al. 2020). It is less clear what the reflection term for kink waves might be. Surely, it must be proportional to the strength of the longitudinal-kink-speed gradient. Inspiration for this aspect may also be found in the existence of a cut-off period for kink waves (Pelouze et al. 2023).

After a lot of calculations, it seems that the evolution equation for Alfvén and kink waves (Eqs. (67) and (68), respectively) are incredibly similar. However, there is a key difference. On the one hand, the damping term for the Alfvén waves is dependent on the amplitude of the counter-propagating wave (as has been known for a long time, at least in the homogeneous or WKB picture; conversely, in a stratified plasma the Alfvén waves also acquire a mixed character, enabling cascades; Velli et al. 1989; Verdini et al. 2009). Thus, upward-propagating Alfvén waves in plumes barely damp in the low solar atmosphere, because insufficient reflected waves have been generated. On the other hand, the damping term for the kink wave due to uniturbulence only depends on the amplitude of the kink wave itself (and not on its counter-propagating wave). Thus, the damping is at work, even if no counter-propagating waves are generated yet. This is of particular importance in the low part of the solar atmosphere. There, the kink waves damp immediately and significantly (damping lengths are between 0.25 and 10× their wavelengths; i.e. 25 Mm–4000 Mm, Van Doorsselaere et al. 2020b), providing a crucial ingredient to kick-starting the solar wind and coronal heating low down. We thus expect a major contribution of the kink-wave heating in the low corona, whereas the traditional Alfvén-wave heating takes over higher up (e.g. around 1 R_⊙).

A limitation to these wave equations (for the Alfvén wave as well as the kink wave) is that it has been derived for in the WKB limit. We (and previous authors: van der Holst et al. 2014; Réville et al. 2020) thus used the implicit assumption that the background variables only vary slowly along the magnetic field. This allows the waves to be described by a single Elsässer or Q variable and decouples the wave modes. However, if the stratification is significant, the waves are modified and will have a mixed character (Velli et al. 1989; Verdini et al. 2009). Thus, the full, true evolution equations will be much more complex and will retain their 3D nature (Zank et al. 2012).

6. Numerical implementation

To show a first proof of concept, we implemented the UAWSOM equations (Eqs. (63)–(68)) in 1D in Python using the fipy package¹ (Guyer et al. 2009). However, for this proof of concept, we neglected the Alfvén-wave equation (Eq. (67), corresponding to taking W_A^± ≡ 0 in the entire simulation). For ease of implementation, we used the equation for internal energy (Eq. (60)) instead of Eq. (66). In 1D, the system of equations (Eqs. (63)–(68)) is reduced to a coupled set of five differential equations, since the velocity only has a component along z and the induction equation is reduced to B(z) constant (as a function of time). The resulting Python implementation serves as a reference implementation for the UAWSOM concept and is available on Gitlab². In this paper, we used revision 95111d80. In order to check if the wave-propagation equations are correctly implemented, we also performed a simplified test on a stationary background. The Python script for that test is also accessible through the Gitlab repository.

In the equations (Eqs. (63)–(68)), we incorporated the radiative losses using the implementation of Hermans & Keppens (2021) of the coronal cooling curve of Dere et al. (2009). We interpolated the cooling curve’s tabulated points with a cubic spline to the temperatures needed in the simulation. On top of the terms included in the energy equation above (Eq. (60)), we also incorporated a thermal conduction term as part of its RHS:

$$\begin{array}{c}\hfill \text{Thermal}\text{conduction}=-\mathrm{\nabla }·(\kappa \mathrm{\nabla }T),\end{array}$$(69)

with κ = κ₀T^5/2 and κ₀ = 8 10⁻⁷ erg/cm/s/K^7/2. In this first implementation, we neglected the reflection terms ℛ_k = 0 in the equation for the kink-wave energy (Eq. (68)). The correct values are not known (Pelouze et al. 2023) and should be investigated in a separate modelling paper. Gravity is fixed at g(z) = 274m/s²1_z.

The Python program has an option to read in a numpy save file for its initial condition. However, here we show some results from the start of a pre-implemented plasma configuration. In that pre-implemented configuration, we used the following as equilibrium:

$$\begin{array}{cc}\hfill \rho (z)& ={\rho }_{0}exp(-z/H),\hfill \end{array}$$(70)

$$\begin{array}{cc}\hfill V(z)& =0,\hfill \end{array}$$(71)

$$\begin{array}{cc}\hfill B(z)& =B_0\frac{R_{\odot }^2}{{(z+R_{\odot })}^2},\hfill \end{array}$$(72)

$$\begin{array}{cc}\hfill R(z)& =R_0\frac{B_0^2}{B{(z)}^2},\hfill \end{array}$$(73)

$$\begin{array}{cc}\hfill \zeta & =({\zeta }_0-1)exp(-\frac{z}{5R_{\odot }})+1,\hfill \end{array}$$(74)

with the initial values ρ₀ is taken equivalent to a number density of 10⁹ cm⁻³, H = 50 Mm, B₀ = 20 G, ζ₀ = 5, and R₀ = 1 Mm. We fixed the initial temperature in the entire domain at T₀ = 1 MK.

The initial conditions for W_k^± were taken as

$$\begin{array}{ccc}\hfill W_{\mathrm{k}}^{-}& =W_0(\frac{1}{2}cos\left(\frac{\pi z}{\delta }\right)+\frac{1}{2}),\hfill & \hfill \text{for}z<\delta ,\end{array}$$(75)

$$\begin{array}{ccc}\hfill W_{\mathrm{k}}^{-}& =0,\hfill & \hfill \text{elsewhere},\end{array}$$(76)

$$\begin{array}{cc}\hfill W_{\mathrm{k}}^{+}& =0,\hfill \end{array}$$(77)

with δ = 30 Mm and W₀ = 4 10⁻⁴ J/m³. This smooth initial condition for W_k⁻ was chosen so that an initial abrupt start from a non-zero boundary condition would not lead to difficulties in this regard.

We took the following at the bottom boundary as boundary conditions:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\rho }{\mathrm{d}z}& =0,\hfill \end{array}$$(78)

$$\begin{array}{cc}\hfill V& =0,\hfill \end{array}$$(79)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{E}_{\mathrm{th}}}{\mathrm{d}z}=\frac{1}{\gamma -1}\frac{\mathrm{d}p}{\mathrm{d}z}& =-\frac{g\rho }{\gamma -1},\hfill \end{array}$$(80)

$$\begin{array}{cc}\hfill W_{\mathrm{k}}^{-}& =W_0,\hfill \end{array}$$(81)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{W}_{\mathrm{k}}^{+}}{\mathrm{d}z}& =0,\hfill \end{array}$$(82)

where the condition for the thermal energy E_th expresses a continuous isothermal stratification beyond the computational domain. At the top boundary, we set the following boundary conditions:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}\rho }{\mathrm{d}z}& =0,\hfill \end{array}$$(83)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}V}{\mathrm{d}z}& =0,\hfill \end{array}$$(84)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{E}_{\mathrm{th}}}{\mathrm{d}z}=\frac{1}{\gamma -1}\frac{\mathrm{d}p}{\mathrm{d}z}& =-\frac{g\rho }{\gamma -1},\hfill \end{array}$$(85)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{W}_{\mathrm{k}}^{-}}{\mathrm{d}z}& =0,\hfill \end{array}$$(86)

$$\begin{array}{cc}\hfill W_{\mathrm{k}}^{+}& =\frac{W_0}{200}.\hfill \end{array}$$(87)

Our aim was to have reasonably open boundary conditions at the top. However, Fig. 3 shows that there are still slow waves reflected in the simulation, as can be seen from bouncing flows between the boundaries.

We display the first 55 s of runtime (1000 iteration steps with fipy) in Fig. 1. The six top panels in the figure show the system variables (ρ, V, E_th, W_k^±) and the temperature T. The six bottom panels show the simulation, but with a much lower energy level W_k⁻ in the kink wave (by setting the left boundary condition for this variable to W₀/200). As is visible in the bottom right panel (of the six top panels), the kink wave is propagating into the system with the appropriate kink speed. The resulting energy dissipation and deposition is leading to heating, as can be seen by comparing the temperature scale of the T-panel of the six top and bottom panels. The temperature T is gradually decreasing in both because of the radiative losses. However, the top panel with the kink-wave heating has a slightly higher temperature than the non-heated simulation. This is more clearly shown in Fig. 2, where the final temperature profile as a function of height is shown for the top and bottom simulations of Fig. 1, corresponding to the case of ‘heating’ and ‘no heating’. Additionally, we display the initial temperature of 1MK and a simulation with kink-wave heating, but where the radiative losses have been switched off. The graph in Fig. 2 clearly shows that the wave heating increases the temperature, either above the initial temperature in the case of no radiative losses, or at least above the temperature in a simulation with just radiative cooling.

Fig. 1. First 1000 time steps for large simulation. Six top panels: Simulation variables for simulation with wave heating. As in the following figures, the six panels consist of density ρ, velocity V, thermal energy Eth, temperature T, and wave energies Wk±. Six bottom panels: Same, but with a much lower influx value of Wk− = W0/200 at the left boundary, resulting in almost no kink-wave heating.

Fig. 2. Temperature profile as function of height at end of short simulations shown in Fig. 1 for three different simulations: one with kink-wave heating and no radiative losses, one with kink-wave heating and radiative losses, and one without kink-wave heating and radiative losses.

The effect of the kink-wave pressure P_k is also visible in the panels of Fig. 1 showing the velocity V. The top right panel with kink wave pressure shows a gentle upflow, while the bottom panels show a gradual downflow. The latter happens because of the restructuring of the plasma due to the radiative cooling. Also visible in that panel (V of the bottom 6 panels) is the effect of the top boundary, which has a propagating signal associated with the downward kink wave, but also a signal propagating with the sound speed.

The effect of the numerical diffusivity is shown in the bottom left panel W_k⁺. The front of the W_k⁺ energy should theoretically remain a sharp interface, but it is clear that it is smoothed out during its propagation.

In Fig. 3, we show the longer duration run of the simulation with the kink-wave heating for 3100 s. In the panel with velocity V, it is shown that there is a gradual upflow generated by the kink-wave driving, of the order of 5 − 10 km/s. However, at times later than 2000 s and height above 70 Mm, numerical instabilities can be seen, which eventually crash the simulation. Some additional viscosity terms might damp out the instabilities. Still, the panel clearly shows that upflows are driven by the kink-wave pressure, because these upflows are absent in the non-driven simulation.