© The Authors 2024

1. Introduction

Eighty years since their discovery (Alfvén 1942), Alfvén waves continue to attract much research interest from laboratory and space plasma physicists. The waves represent an interplay between plasma inertia and the tension of the magnetic field lines, similar to waves on a string instrument.

Alfvén waves have long been thought to be a possible candidate in resolving two major and longstanding problems in solar physics and astrophysics: the heating of the solar and stellar atmospheres and the acceleration of solar and stellar winds. The waves can be generated through random motions or magnetic reconnection in the photosphere and chromosphere.

The main reason why Alfvén waves have been favoured in comparison to other types of magnetohydrodynamic (MHD) waves is that incompressible Alfvénic disturbances travelling along magnetic field lines should be less susceptible to steepening and dissipation before reaching the upper atmosphere of the sun. They could therefore transport significant amounts of energy along the magnetic field that extends from the visible surface of the sun (photosphere) into the hot and tenuous atmosphere (corona) and further (Liu et al. 2019; Soler et al. 2019; Alielden & Taroyan 2022).

However, early studies by Ferraro (1954) and Ferraro & Plumpton (1958) indicated that the chromosphere and transition region between the photosphere and the corona should act as strong reflectors of wave energy due to steep temperature and density gradients. Follow-up theoretical modelling endorsed the view that Alfvén waves are attenuated by the abrupt increase in their propagation speed by more than a factor of 10 (over about 1 Mm) in passing from the chromosphere to the corona (An et al. 1989; Parker 1991).

The reflection is especially strong at long periods. On the other hand, Alfvén waves with shorter periods are efficiently damped in the chromosphere by a variety of mechanisms, for example by ion-neutral collisions (Zaqarashvili et al. 2013).

These views have been challenged in the new era of solar observations. High-resolution spectroscopic observations with Hinode/EIS, SST, IRIS, and other instruments have demonstrated the ubiquity of Alfvén waves in the solar atmosphere.

Large-amplitude magnetic twists and Alfvénic waves appear to be propagating along different magnetic structures carrying significant amounts of energy into the chromosphere and the corona (Jess et al. 2009; McIntosh et al. 2011; De Pontieu et al. 2014). Recent work has revealed and characterised ubiquitous and continuous wave activity in extreme ultraviolet image time series through the application of a novel method (Morgan & Korsós 2022).

Several studies have proposed that the reflection may lead to resonances, and thus to a very large energy flux at those frequencies (Hollweg 1978, 1984; Matsumoto & Shibata 2010). However, Leer et al. (1982) pointed out that such large energy fluxes can only be produced if the wave-driver matches the nodes of the resonator and remains phase-coherent over many wave periods. In addition, the resonator should remain fixed in time over many wave periods. These requirements impose strong restrictions on the transport of large amounts of energy from the photosphere to the corona through resonances.

Studies carried out by Wedemeyer-Böhm et al. (2012), Battaglia et al. (2021), Breu et al. (2023), and Silva et al. (2024) show the presence of ubiquitous vortical or swirling motions both in observations and numerical simulations. These events are correlated with perturbations of the magnetic field. At photospheric and chromospheric levels, they form Alfvén pulses that propagate upward and may contribute to chromospheric heating. Energy is transported following a vortical motion.

Theoretical modelling of Alfvén wave propagation has traditionally been carried out in standard magnetic field geometries. These include a uniform field or a thin flux tube whereby the tube thickness is small compared with the characteristic scale of the magnetic field variation. However, the magnetic field lines can expand in different ways as the density drops from the photosphere to the corona. In general, the shape of the field lines can be determined from a divergence-free condition in curvilinear coordinates (Taroyan et al. 2021).

Here we demonstrate that the above-mentioned standard geometries inhibit the flux of Alfvén wave energy into the upper atmosphere. Examples of structures with such geometries are a stratified atmosphere with a uniform magnetic field, a radially expanding magnetic field, or a thin flux tube. Reflection approaches 100% at low frequencies. We show that the wave reflection strongly depends on the magnetic field geometry and that it may vanish in certain cases, resulting in very efficient transport of Alfvén wave energy into the upper atmosphere.

The paper is structured as follows. Section 2 summarises some important results for thin flux tubes in a cylindrical coordinate system. In Sect. 3 we demonstrate that the problem of Alfvén wave propagation in a general field line geometry can be reduced to a set of Klein–Gordon equations in local curvilinear coordinates. Three different cases of field line expansion rates in the lower atmosphere are introduced and compared in Sect. 4; the corresponding Alfvén wave variables are derived in Sect. 5 and compared in Sect. 6. The wave energy fluxes associated with the three cases are derived and compared in Sect. 7. Section 8 summarises the presented results and draws our overall conclusions.

2. Alfvén waves on thin flux tubes

Models of Alfvén wave propagation on thin flux tubes are based on the well-known relationship (9) between the background magnetic field and the tube cross-section (Hollweg et al. 1982). Studies carried out in this approximation show that Alfvén waves are strongly reflected by density and temperature gradients especially at low frequencies (An et al. 1989; Murawski & Musielak 2010; Soler et al. 2017; Tsap & Kopylova 2021). Therefore, very little wave energy is available at coronal heights.

The relationship between the background field and the cross-sectional radius (9) permits a wide range of thin flux tubes with different variations in the tube cross-section. It also includes the case of a uniform field with vertical field lines. The propagation and reflection of Alfvén waves in an atmosphere with a uniform field has been studied since Ferraro (1954) and Ferraro & Plumpton (1958). The wave behaviour in a uniform field is similar to the results obtained for a thin flux tube.

The thin flux tube approximation considers field lines close to the tube’s vertical axis of symmetry, where the radial component, B_r, vanishes. On field lines near the tube axis we can Taylor expand B_r with respect to r,

$$\begin{array}{c}\hfill B_r=0+\frac{\partial B_r(0,z)}{\partial r}r+\frac{{\partial }^2{B}_r(0,z)}{\partial r^2}\frac{r^2}{2}+\dots ,\end{array}$$(1)

where z is the vertical coordinate. Dividing by r and ignoring terms of order r, we obtain

$$\begin{array}{c}\hfill \frac{B_r}{r}=\frac{\partial B_r(0,z)}{\partial r}.\end{array}$$(2)

The derivative of B_r can be Taylor expanded in a similar way:

$$\begin{array}{c}\hfill \frac{\partial B_r}{\partial r}=\frac{\partial B_r(0,z)}{\partial r}+\frac{{\partial }^2{B}_r(0,z)}{\partial r^2}r+\dots .\end{array}$$(3)

Ignoring terms of order r in expansion (3) and combining it with expression (2), we obtain

$$\begin{array}{c}\hfill \frac{B_r}{r}=\frac{\partial B_r(0,z)}{\partial r}=\frac{\partial B_r(r,z)}{\partial r}.\end{array}$$(4)

The divergence-free condition in cylindrical coordinates reads

$$\begin{array}{c}\hfill \mathbf{\nabla }\mathbf{·}\mathit{B}=\frac{\partial B_z}{\partial z}+\frac{\partial B_r}{\partial r}+\frac{B_r}{r}=0.\end{array}$$(5)

For field lines near the axis of symmetry we apply the obtained relationship (4), to rewrite the divergence-free condition in the following form:

$$\begin{array}{c}\hfill \frac{\partial B_z}{\partial z}+2\frac{B_r}{r}=0.\end{array}$$(6)

Along the field lines the following relationship is satisfied:

$$\begin{array}{c}\hfill \frac{\mathrm{d}r}{\mathrm{d}z}=\frac{B_r}{B_z}.\end{array}$$(7)

Using the obtained Eq. (6) for field lines close to the symmetry axis, we have

$$\begin{array}{c}\hfill \frac{\mathrm{d}r}{\mathrm{d}z}=-\frac{r}{2}\frac{\mathrm{d}{B}_z}{B_z\mathrm{d}z}\text{or}2\frac{\mathrm{d}r}{r}+\frac{\mathrm{d}{B}_z}{B_z}=0.\end{array}$$(8)

Integrating Eq. (8) along an arbitrary field line, we obtain the condition of flux conservation for a thin flux tube:

$$\begin{array}{c}\hfill B_z{r}^2=\mathrm{const}\mathrm{along}\mathrm{field}\mathrm{lines}.\end{array}$$(9)

The above condition (9) can be derived in different ways. For example, in the case of a potential field, the flux function is conserved along the field lines, and it reduces to the left-hand side of (9) for field lines close enough to the symmetry axis (Ruderman et al. 2008; Verth et al. 2010; Soler et al. 2017). Condition (9) has been extensively used in studies of Alfvén waves on thin flux tubes.

The governing equation for the Alfvén waves on thin flux tubes can be written in the form of a Bessel equation, with a general solution in terms of the Hankel functions (see e.g. Hollweg (1978)):

$$\begin{array}{c}\hfill x=a{H}_0^{(1)}\left(\frac{h\omega }{c_A}\right)+b{H}_0^{(2)}\left(\frac{h\omega }{c_A}\right).\end{array}$$(10)

Here x = v_φ/r, and c_A = c_A0exp(z/h) is the exponentially increasing Alfvén speed with a scale height h. The general solution (10) is derived assuming condition (9) is satisfied. In the case of an atmosphere with a uniform magnetic field and an exponential Alfvén speed, the governing equations are reduced to the wave equation, which, in turn, can be reduced to the Bessel equation, and the general solution is given by the same formula (10) (Hollweg 1978; Verth et al. 2010; Tsap & Kopylova 2021).

The model with an exponentially increasing Alfvén speed has been used to investigate the propagation of Alfvén waves in stellar atmospheres since Ferraro (1954). Ferraro & Plumpton (1958) expressed the solutions in terms of Bessel functions of the first and second kind of order zero. Cally (2012) pointed out that dropping the Bessel function Y₀ used by Ferraro & Plumpton (1958) and An et al. (1989), among others, leads to unphysical reflection at infinity. Therefore, solutions are usually expressed in terms of the Hankel functions (10) that represent upgoing and downgoing waves, respectively.

Fig. 1. Axisymmetric poloidal magnetic field described by cylindrical coordinates and local orthogonal curvilinear coordinates. The angle ϕ representing the angle of rotation around the vertical symmetry axis, z, is not shown. The figure shows distance, s, along a curved magnetic field; the angle, θ, between the poloidal field and the vertical axis; and incremental increases δs and δa, along and across the poloidal field, respectively.

3. Klein–Gordon equations for Alfvén waves

In cylindrical coordinates (r, ϕ, z) axisymmetric motions are characterised by the condition ∂/∂ϕ = 0. We introduced local orthogonal curvilinear coordinates (a, ϕ, s) along a given magnetic field line tailored to the geometry of the magnetic field. Figure 1 shows that s is a curvilinear coordinate because it follows the curved path of the magnetic field line rather than a straight line. The coordinate s measures distance along the poloidal field line, and a is the distance perpendicular to the poloidal field line. The angle of rotation, ϕ, around the vertical z-axis is not shown. The magnetic field may be decomposed either into cylindrical components (B_r, B_ϕ, B_z) or into toroidal and poloidal components (0, B_ϕ, B_s), where B_s denotes the poloidal field and there is no component in the transverse direction, a.

The propagation of axisymmetric Alfvén waves in general field geometries can be studied using the equations of motion and induction in local curvilinear coordinates (Hollweg et al. 1982; Taroyan et al. 2021):

$$\begin{array}{c}\hfill \frac{\partial V_{\phi }}{\partial t}=\frac{B_s}{{\mu }_{0}r\rho }\frac{\partial }{\partial s}\left(r{B}_{\phi }\right),\frac{\partial B_{\phi }}{\partial t}=r{B}_s\frac{\partial }{\partial s}\left(\frac{V_{\phi }}{r}\right).\end{array}$$(11)

Here B_s is the time-independent background magnetic field and s is the local curvilinear coordinate along a given field line.

It is well known that one-dimensional wave propagation in the presence of gravity typically gives rise to the Klein–Gordon equation and with it the introduction of a cut-off frequency (Roberts 2004; Taroyan & Erdélyi 2008). The cut-off frequency is related to the scale-height, which is a measure of the strength of gravity. The significance of the Klein–Gordon equation is that it introduces a timescale imposed by the equilibrium. For a Klein–Gordon equation with constant coefficients, wave propagation occurs only for frequencies above the cut-off frequency. Frequencies below this value result in no propagation and the motion is evanescent (Roberts 2004).

Noble et al. (2003) derived a Klein–Gordon equation for the velocity perturbation, V_φ, and Routh et al. (2020) show that the magnetic field perturbation, B_φ, is determined from the Klein–Gordon equation. These studies are carried out under the assumption of a thin flux tube approximation (i.e. the condition of flux conservation (9) is satisfied).

Here we show that the general problem of Alfvén wave propagation along field lines of arbitrary geometry can be reduced to a set of Klein–Gordon equations for the perturbations of the magnetic field, B_φ, and the velocity, V_φ. This approach allows us to explore the relationship between the waves, the corresponding energy fluxes, and the cut-off frequencies for different geometries.

Firstly, we differentiate the momentum and the induction Equation (11) with respect to t and s, respectively, to obtain

$$\begin{array}{c}\hfill \frac{{\mu }_{0}r\rho }{B_s}\frac{{\partial }^2{V}_{\phi }}{\partial t^2}=\frac{{\partial }^2}{\partial t\partial s}\left(r{B}_{\phi }\right),\frac{{\partial }^2}{\partial t\partial s}\left(r{B}_{\phi }\right)=\frac{\partial }{\partial s}\left(B_s{r}^2\frac{\partial }{\partial s}\left(\frac{V_{\phi }}{r}\right)\right).\end{array}$$(12)

Equation (12) can be combined to eliminate the variable B_φ

$$\begin{array}{c}\hfill \frac{{\partial }^2}{\partial t^2}\left(\frac{V_{\phi }}{r}\right)=\frac{c_A^2}{A^2}\frac{\partial }{\partial s}\left(A^2\frac{\partial }{\partial s}\frac{V_{\phi }}{r}\right),\end{array}$$(13)

where $c_A^2=\frac{B_s^2}{{\mu }_0\rho }$ and A² = B_sr². By introducing the variable ${\displaystyle Q_1=\frac{A}{r}{V}_{\phi }}$, we reduce the right-hand side of Eq. (13) to

$$\begin{array}{c}\hfill \frac{c_A^2}{A^2}\frac{\partial }{\partial s}\left(A^2\frac{\partial }{\partial s}\frac{V_{\phi }}{r}\right)=\frac{c_A^2}{A^2}(A\frac{{\partial }^2{Q}_1}{\partial s^2}-Q_1\frac{{\partial }^{2}A}{\partial s^2}).\end{array}$$(14)

Equation (13) is therefore reduced to a Klein–Gordon equation:

$$\begin{array}{c}\hfill \frac{{\partial }^2{Q}_1}{\partial t^2}-c_A^2\frac{{\partial }^2{Q}_1}{\partial s^2}+\frac{c_A^2}{A}\frac{{\partial }^{2}A}{\partial s^2}{Q}_1=0.\end{array}$$(15)

The second Klein–Gordon equation can be derived in a similar way. We differentiate the momentum and the induction Equation (11) with respect to s and t, respectively, to obtain

$$\begin{array}{c}\hfill \frac{{\partial }^2}{\partial t\partial s}\left(\frac{V_{\phi }}{r}\right)=\frac{\partial }{\partial s}\left(\frac{c_A^2}{A^2}\frac{\partial }{\partial s}\left(r{B}_{\phi }\right)\right),\frac{1}{A^2}\frac{{\partial }^2}{\partial t^2}\left(r{B}_{\phi }\right)=\frac{{\partial }^2}{\partial t\partial s}\left(\frac{V_{\phi }}{r}\right).\end{array}$$(16)

Equation (16) can be combined to eliminate the variable V_φ:

$$\begin{array}{c}\hfill \frac{1}{A^2}\frac{{\partial }^2}{\partial t^2}\left(r{B}_{\phi }\right)=\frac{\partial }{\partial s}\left(\frac{c_A^2}{A^2}\frac{\partial }{\partial s}\left(r{B}_{\phi }\right)\right).\end{array}$$(17)

By introducing the variable ${\displaystyle Q_2=\frac{c_A}{A}r{B}_{\phi }}$, we reduce the right-hand side of Eq. (17) to

$$\begin{array}{c}\hfill \frac{\partial }{\partial s}\left(\frac{c_A^2}{A^2}\frac{\partial }{\partial s}\left(r{B}_{\phi }\right)\right)=\frac{c_A}{A}\frac{{\partial }^2{Q}_2}{\partial s^2}-Q_2\frac{{\partial }^2}{\partial s^2}\left(\frac{c_A}{A}\right).\end{array}$$(18)

Equation (17) is therefore reduced to a Klein–Gordon equation of the form

$$\begin{array}{c}\hfill {\displaystyle \frac{{\partial }^2{Q}_2}{\partial t^2}-c_A^2\frac{{\partial }^2{Q}_2}{\partial s^2}+A{c}_A\frac{{\partial }^2}{\partial s^2}\left(\frac{c_A}{A}\right)Q_2=0.}\end{array}$$(19)

In summary, we have demonstrated that the propagation of Alfvén waves along field lines of arbitrary geometry is determined by the Klein–Gordon equation

$$\begin{array}{c}\hfill \frac{{\partial }^2{Q}_i}{\partial t^2}-c_A^2\frac{{\partial }^2{Q}_i}{\partial s^2}+{\mathrm{\Omega }}_i^2{Q}_i=0,\end{array}$$(20)

where s is the coordinate along the poloidal field, c_A is the Alfvén speed,

$$\begin{array}{c}\hfill Q_i=\{\begin{array}{cc}{\displaystyle \frac{A}{r}{V}_{\phi },}\hfill & \text{for}i=1,\hfill \\ \hfill \\ {\displaystyle \frac{r{c}_A}{A}{B}_{\phi },}\hfill & \text{for}i=2\hfill \end{array}\end{array}$$(21)

is the wave variable, A² = B_sr², and

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_i^2=\{\begin{array}{cc}{\displaystyle \frac{c_A^2}{A}\frac{{\partial }^{2}A}{\partial s^2},}\hfill & \text{for}i=1\hfill \\ \hfill \\ {\displaystyle A{c}_A\frac{{\partial }^2}{\partial s^2}\left(\frac{c_A}{A}\right),}\hfill & \text{for}i=2\hfill \end{array}\end{array}$$(22)

is the square of the cut-off frequency.

By setting Q_i = q_iexp(iωt), we can reduce the Klein–Gordon Eq. (20) to the following form:

$$\begin{array}{c}\hfill \frac{d^2{q}_i}{d{s}^2}+\frac{{\omega }^2-{\mathrm{\Omega }}_i^2}{c_A^2}{q}_i=0.\end{array}$$(23)

In the WKB approximation, eigenfrequencies below the cut-off, ω < Ω_i, correspond to growing and decaying solutions, whereas those above the cut-off, ω > Ω_i, correspond to propagating solutions.

The thin flux tube approximation represents a particular case of a general class of field lines that can be self-consistently determined using the solenoidal condition in curvilinear coordinates (Taroyan et al. 2021). Here we show that the commonly used models, including the thin flux tube, correspond to a degeneracy in the Klein–Gordon Eq. (20), which has important implications.

For field lines close to the symmetry axis we have B_r ≪ B_z so that B_s ≈ B_z. The conservation of magnetic flux for thin flux tubes (9) is then reduced to A² = B_sr² ≈ B_zr²= const along the field lines.

Therefore, the thin flux tube approximation assumes A= const along the poloidal field, in which case the cut-off frequency Ω₁ vanishes. It is important to note that this result is true for all thin flux tubes regardless of the particular field line geometry.

The condition A= const remains true for the cases of a uniform field and a radially expanding field. Therefore, the cut-off frequency vanishes (Ω₁ = 0) and the corresponding Klein–Gordon equation again becomes degenerate.

Below we demonstrate that vanishing cut-off frequencies (cases I and II) strongly restrict the wave energy flux into the upper atmosphere. A counterexample (case III) shows that more generally (i.e. in non-degenerate cases) the transport of the wave energy by Alfvén waves into the corona can be much more efficient.

Fig. 2. Equilibrium parameters in the atmosphere. The panels on the left display the variation of density, Alfvén speed and magnetic field. The panels on the right show the radial distance, r = r(s), and the corresponding parameters that determine the cut-off frequencies Ω1 and Ω2. Logarithmic scales are used.

Fig. 3. Geometry of the magnetic field in the r − z plane. The field lines exponentially expand in the lower atmosphere and are vertical in the uniform upper atmosphere.

4. The background atmosphere

Following previous studies, we split the atmosphere into two parts (e.g. Hollweg et al. 1982). The upper atmosphere is uniform and all the parameters here are constant. In the lower atmosphere, the Alfvén speed increases exponentially with height. In order to explore in more detail the behaviour of the waves and the role of the magnetic field geometry, we consider three different cases for the lower atmosphere that are summarised in Table 1.

Table 1 shows the variation of the Alfvén speed, density, magnetic field, and the radial distance in the lower atmosphere. All the parameters are treated as functions of the field line coordinate, s. The photosphere is located at s = −2 Mm, and the coronal base is located at s = 0 Mm. The subscript 0 indicates the values of the corresponding parameters at the coronal base, s = 0. Table 1 shows that all the parameters are the same for the three cases except the distance between the field line and the symmetry axis, r.

We assume that the density decreases from the photosphere to the base of the corona by eight orders of magnitude (Vernazza et al. 1981). Further, in agreement with the measurements, Table 1 suggests a magnetic field decrease by two orders of magnitude (Aschwanden 2019). The Alfvén speed should then increase by two orders of magnitude. The scale height corresponding to these variations is h = 869.5 km. The main difference between the three cases is the expansion factor of the field lines.

Observations suggest that magnetic flux tubes originating from concentrated sources in the photosphere can expand significantly, by a factor of 100 or more, as they rise to the corona (Gary 2001; Solanki et al. 2006). The cross-sectional area of these flux tubes increases due to the divergence of field lines (Aschwanden 2019). We consider three possibilities: the field lines expand by a factor of 10 in case I, a factor of 1000 in case II, and a factor of 100 in case III.

The parameter variations described above are visualised in Fig. 2. We note that logarithmic scales are used for the vertical axes. We adopt a two-layer model of the solar atmosphere with a constant scale height of less than 1 Mm in the lower layer and an infinite scale height corresponding to high temperatures in the upper layer (Leer et al. 1982; Hollweg 1984; Tsap & Kopylova 2021). The background atmosphere can be approximated more accurately. For example, Hollweg (1978) introduced a model with 16 exponential layers mimicking the Munro–Jackson model to study the behaviour of Alfvén waves from the lower atmosphere to the outer corona and their role in solar wind acceleration. However, the Klein–Gordon equations derived in the previous section do not assume any specific structure of the atmosphere. Therefore, we could argue that our main conclusions regarding the relationship between the cut-off frequencies, the field line geometry, and the energy flux carried by the Alfvén waves should remain valid for more complex structures.

Figure 2 shows that the atmosphere is split into a lower layer, s < 0, with an exponentially increasing Alfvén speed, density, and magnetic field, and an upper layer, s > 0, with a constant value for the same parameters. In agreement with previous studies, the solution in the infinite and constant upper layer is determined by an outgoing wave condition.

The panels on the right show the variations in the radial distance from the symmetry axis, r = r(s), and the parameters that determine the cut-off frequencies for the three cases. The parameter A is constant for case I (orange dashed line), and leads to a vanishing cut-off frequency Ω₁. For case II (green dotted line) it is the parameter c_A/A that remains constant, leading to vanishing cut-off frequency Ω₂. For the case III (blue line) both cut-off frequencies remain different from 0.

Figure 3 displays the geometry of the magnetic field in arbitrary units. Equations (7) and (8) from Taroyan et al. (2021) establish relationships between the field components in local curvilinear and cylindrical coordinates. These relationships are used to construct and plot the field lines in the r − z plane. The axisymmetric geometry implies no dependence on the azimuthal angle ϕ.

The uniform field, B_s0, in the upper atmosphere, s > 0, corresponds to vertical field lines that are identical for the three cases. In the linear regime, the Alfvénic motions are decoupled from the other motions. However, the dynamic nature of the transition region or the inclination of the field lines in the upper atmosphere could play an important role as the waves become non-linear. The three cases show differences in the lower atmosphere, s < 0. In each case, the rate of the field line expansion is determined by the corresponding scale height: h for case I, h/3 for case II, and h/2 for case III. These three expansion rates result in very different wave behaviour.

The expressions from Table 1 show that the three cases have the same poloidal field, B_s, from the photospehre to the corona. These expressions and the relationships between the field components in local curvilinear and cylindrical coordinates can be used to determine the horizontal, r, and the vertical, z, components of the magnetic field. Figure 3 shows that Case I has the strongest horizontal component and the weakest vertical component at the photospheric level, s = −2 Mm, whereas case II has the strongest vertical component and the weakest horizontal component at the photospheric level.

5. Solutions of the Klein–Gordon equations

Using the equilibrium parameters introduced above, we derive and examine the solutions of the Klein–Gordon equations for the three cases. We first consider the two degenerate cases I and II that correspond to vanishing cut-off frequencies Ω₁ and Ω₂, respectively, in the Klein–Gordon equations. We then examine a counterexample (case III) where neither of the cut-off frequencies vanishes. The obtained results are then used to compare the magnetic field and velocity perturbations in the next section.

The solutions to the Klein–Gordon equations are constructed based on three key principles:

• 1.

  The net wave energy flux is finite at all frequencies.

• 2.

  The net wave energy flux is zero when the upgoing flux equals the downgoing flux.

• 3.

  The upgoing flux from the photosphere is independent of the magnetic field geometry above and remains consistent across all three cases.

5.1. Case I

In the first case, which also includes thin flux tubes, we have the condition A = const, and the first cut-off frequency vanishes: Ω₁ = 0. The corresponding Klein–Gordon equation for the lower atmosphere becomes degenerate. It is reduced to the wave equation with solutions that represent propagating waves.

The condition Ω₁ = 0 holds for all thin flux tubes regardless of the particular dependence r = r(s) or c_A = c_A(s). In the lower atmosphere, s < 0, where the Alfvén speed increases exponentially with a scale height of h/2, the solutions are expressed in terms of the Hankel functions, similar to the expression (10) for cylindrical geometry:

$$\begin{array}{c}\hfill q_1=\sqrt{\omega }a{H}_0^{(1)}\left(\frac{h\omega }{2c_A}\right)+\sqrt{\omega }b{H}_0^{(2)}\left(\frac{h\omega }{2c_A}\right),s<0.\end{array}$$(24)

One difference between our solution (24) and the solution (10) used by Hollweg (1978) is the Alfvénic scale height: h/2 in (24) and h in (10). Another difference is the presence of an extra coefficient $\sqrt{\omega }$ in (24) in order to have physically acceptable solutions for all frequencies. We show in Sect. 7 that due to the coefficient $\sqrt{\omega }$ there is no net energy flux when the upgoing and downgoing wave energy fluxes become equal to each other as ω → 0. On the other hand, the net flux remains finite as ω → ∞. The solution (24) is also applicable to a uniform or a radial magnetic field, where the same condition, A = const, is satisfied.

In addition to the variables q_i, we introduce the new variables v_φ and b_φ: V_φ = v_φexp(iωt), B_φ = b_φexp(iωt). The velocity perturbation, v_φ, can easily be found from the solution (24) and the relationship (21):

$$\begin{array}{c}\hfill v_{\phi }=\frac{i\sqrt{\omega }}{\sqrt{B_s}}[a{H}_0^{(1)}\left(\frac{h\omega }{2c_A}\right)+b{H}_0^{(2)}\left(\frac{h\omega }{2c_A}\right)],s<0.\end{array}$$(25)

The magnetic field perturbation, b_φ, can be determined from the induction in Eq. (11):

$$\begin{array}{c}\hfill b_{\phi }=-\frac{ir{B}_s}{\omega }\frac{\mathrm{d}}{\mathrm{d}s}\left(\frac{q_1}{A}\right)=-\frac{i\sqrt{B_s}}{\omega }\frac{\mathrm{d}{q}_1}{\mathrm{d}s}.\end{array}$$(26)

Using the obtained solution (24), the above expression is reduced to

$$\begin{array}{c}\hfill b_{\phi }=-\frac{i\sqrt{\omega B_s}}{c_A}[a{H}_1^{(1)}\left(\frac{h\omega }{2c_A}\right)+b{H}_1^{(2)}\left(\frac{h\omega }{2c_A}\right)],s<0.\end{array}$$(27)

The cut-off frequencies vanish in the upper atmosphere, and therefore only propagating waves can exist in this open region. Consistent with previous studies, we require the solution to represent an outgoing wave in the open region, s > 0,

$$\begin{array}{c}\hfill q_1=cexp(-\frac{i\omega s}{c_{A0}}),s>0,\end{array}$$(28)

where c_A0 = c_A(s = 0) and c is the wave amplitude in the upper layer. The magnetic field and velocity perturbations, b_φ and v_φ, in the region s > 0 can be determined from expression (28):

$$\begin{array}{c}\hfill v_{\phi }=\frac{c}{\sqrt{B_{s0}}}exp(-\frac{i\omega s}{c_{A0}}),b_{\phi }=-\frac{c\sqrt{B_{s0}}}{c_{A0}}exp(-\frac{i\omega s}{c_{A0}}).\end{array}$$(29)

Here the index 0 indicates that the functions are evaluated at s = 0.

A relationship between the constant coefficients a and b can be established using the continuity of the magnetic field and the velocity perturbations across the interface at s = 0. The result is

$$\begin{array}{c}\hfill \frac{b}{a}=-\frac{i{H}_1^{(1)}-H_0^{(1)}}{i{H}_1^{(2)}-H_0^{(2)}},\end{array}$$(30)

where the argument of the Hankel functions is at hω/(2c_A0). Hollweg (1978) identified ${\left|\frac{b}{a}\right|}^2$ as the ratio of downgoing energy flux to upgoing energy flux (see also Leer et al. (1982), Tsap & Kopylova (2021)).

5.2. Case II

For the second case we have c_A/A= const along the field lines, corresponding to Ω₂ = 0. Therefore, the second Klein–Gordon equation becomes degenerate. In the lower atmosphere, s < 0, the field lines expand with a scale height h/3, and the Alfvén speed increases with a scale height h/2. The equation for q₂,

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{q}_2}{\mathrm{d}{s}^2}+\frac{{\omega }^2}{c_A^2}{q}_2=0,\end{array}$$(31)

is reduced to a Bessel equation with solutions similar to (24):

$$\begin{array}{c}\hfill q_2=\sqrt{\omega }a{H}_0^{(1)}\left(\frac{h\omega }{2c_A}\right)+\sqrt{\omega }b{H}_0^{(2)}\left(\frac{h\omega }{2c_A}\right),s<0.\end{array}$$(32)

Solutions of this type have been studied by Hollweg (1984), Cally (2012), and others in different contexts. Similarly to case I, the extra coefficient, $\sqrt{\omega }$, ensures that there is no net flux of wave energy when the upgoing and downgoing fluxes balance each other (see Sect. 7). The net flux remains finite for all frequencies. The coefficients a and b are the same as in case I to satisfy the third principle (i.e. the upgoing flux from the photospehre remains independent of the magnetic field geometry above).

The magnetic field perturbation, b_φ, in the lower atmosphere is expressed through Eq. (21):

$$\begin{array}{c}\hfill b_{\phi }=\frac{\sqrt{\omega B_s}}{c_A}[a{H}_0^{(1)}\left(\frac{h\omega }{2c_A}\right)+b{H}_0^{(2)}\left(\frac{h\omega }{2c_A}\right)].\end{array}$$(33)

The velocity perturbation in the region s < 0 can be expressed using the equation of motion (11):

$$\begin{array}{c}\hfill v_{\phi }=-\frac{i{c}_{A}r}{A\omega }\frac{\mathrm{d}{q}_2}{\mathrm{d}s}=-\frac{i\sqrt{\omega }}{\sqrt{B_s}}[a{H}_1^{(1)}\left(\frac{h\omega }{2c_A}\right)+b{H}_1^{(2)}\left(\frac{h\omega }{2c_A}\right)].\end{array}$$(34)

Here we use the well-known differentiation rule for the Hankel functions: ${\displaystyle \frac{\mathrm{d}}{\mathrm{d}z}{H}_0^{(i)}(z)=-H_1^{(i)}(z),i=1,2.}$

As in case I, we require the solution in the open and uniform region, s > 0, to represent an outgoing wave,

$$\begin{array}{c}\hfill q_2=cexp(-\frac{i\omega s}{c_{A0}}),s>0,\end{array}$$(35)

where c_A0 = c_A(s = 0), and c is the constant amplitude in the upper layer.

The magnetic field and velocity perturbations, b_φ and v_φ in the region s > 0 can be determined from expression (35),

$$\begin{array}{c}\hfill v_{\phi }=-\frac{c}{\sqrt{B_{s0}}}exp(-\frac{i\omega s}{c_{A0}}),b_{\phi }=\frac{c\sqrt{B_{s0}}}{c_{A0}}exp(-\frac{i\omega s}{c_{A0}}),\end{array}$$(36)

where the index 0 indicates that the functions are evaluated at s = 0. The magnetic field and velocity perturbations, b_φ and v_φ, are continuous across the interface at s = 0,

$$\begin{array}{c}\hfill a\sqrt{\omega }{H}_0^{(1)}+b\sqrt{\omega }{H}_0^{(2)}=c,-i\sqrt{\omega }[a{H}_1^{(1)}+b{H}_1^{(2)}]=-c,\end{array}$$(37)

where the Hankel functions are evaluated at $\frac{h\omega }{2c_{A0}}$, corresponding to the interface s = 0. The ratio b/a can be determined by solving the above set of algebraic equations (37). The result is expressed by the same formula (30). It can be used to determine the ratio of the upgoing and downgoing fluxes, ${\left|\frac{b}{a}\right|}^2$.

5.3. Case III

More generally, for non-degenerate cases the cut-off frequencies Ω₁ and Ω₂ are different from 0. Here we consider an example of a non-degenerate case with both A and c_A/A being variable, as shown in Fig. 2. In the lower atmosphere the Alfvén speed increases and the field lines expand with a scale height of h/2. The second Klein–Gordon equation is reduced to

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{q}_2}{\mathrm{d}{s}^2}+(\frac{{\omega }^2}{c_A^2}-\frac{1}{h^2})q_2=0,s<0.\end{array}$$(38)

By introducing the new variables

$$\begin{array}{c}\hfill z=exp(-\frac{2s}{h}),Z=\sqrt{z}{q}_2,\end{array}$$(39)

we transform Eq. (38) into a second-order ordinary differential equation with constant coefficients:

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^{2}Z}{\mathrm{d}{z}^2}+{\left(\frac{\omega h}{2c_{A0}}\right)}^{2}Z=0.\end{array}$$(40)

Equation (40) has solutions of the form ${\displaystyle Z=exp(\pm \frac{ih\omega }{2c_A}z)}$. The general solution of the original Eq. (38) is given by

$$\begin{array}{c}\hfill q_2=\sqrt{\frac{4\pi c_A}{h}}[aexp\left(\frac{ih\omega }{2c_A}\right)+bexp(-\frac{ih\omega }{2c_A})].\end{array}$$(41)

The wave energy flux corresponding to the above solution (41) remains finite for all real frequencies. The constants outside the brackets in Eq. (41) ensure that the same amount of upgoing energy flux from the photosphere is produced in all three cases (see Sect. 7).

The magnetic field perturbation in the lower atmosphere, b_φ, is determined from Eq. (21):

$$\begin{array}{c}\hfill b_{\phi }=\frac{\sqrt{B_s}}{c_A}\sqrt{\frac{4\pi c_A}{h}}[aexp\left(\frac{ih\omega }{2c_A}\right)+bexp(-\frac{ih\omega }{2c_A})].\end{array}$$(42)

The velocity perturbation, v_φ, is determined from the equation of motion (11):

$$\begin{array}{c}\hfill v_{\phi }=-\frac{ir{c}_A^2}{A^2\omega }\frac{\partial }{\partial s}(b_{\phi }r)=-\frac{ir{c}_A^2}{A^2\omega }\frac{\partial }{\partial s}\left(\frac{A}{c_A}{q}_2\right).\end{array}$$(43)

Using the obtained expression (41) for q₂, we find the velocity perturbation in the region s < 0:

$$\begin{array}{c}\hfill v_{\phi }=-\frac{1}{\sqrt{B_s}}\sqrt{\frac{4\pi c_A}{h}}[aexp\left(\frac{ih\omega }{2c_A}\right)-bexp(-\frac{ih\omega }{2c_A})].\end{array}$$(44)

The solutions in the upper atmosphere, s > 0, are again required to represent an outgoing wave. They are determined by the expressions (35) and (36) with a different constant coefficient $\stackrel{\sim }{c}$. Using the continuity conditions for b_φ and v_φ at the interface s = 0, we have

$$\begin{array}{c}\hfill \sqrt{\frac{4\pi c_{A0}}{h}}[aexp\left(\frac{ih\omega }{2c_{A0}}\right)+bexp(-\frac{ih\omega }{2c_{A0}})]=\stackrel{\sim }{c},\\ \hfill -\sqrt{\frac{4\pi c_{A0}}{h}}[aexp\left(\frac{ih\omega }{2c_{A0}}\right)-bexp(-\frac{ih\omega }{2c_{A0}})]=-\stackrel{\sim }{c}.\end{array}$$(45)

From the set of equations (45) we have b = 0. We show in Sect. 7 that the obtained result indicates an absence of any downgoing wave energy flux for case III.

Fig. 4. Comparison of the velocity, vϕ, and magnetic field, bϕ, perturbations for a zero frequency wave corresponding to vortex motion. Both the real and imaginary components, as well as the amplitudes are shown. The dotted vertical line indicates the interface between the lower and the upper regions.

Fig. 5. Comparison of the velocity, vϕ, and magnetic field, bϕ, parameters for a low-frequency wave at ω = 0.00175 rad s−1 corresponding to a period of 1 hour. Both the real and imaginary components, as well as the amplitudes are shown.

Fig. 6. Comparison of the velocity, vϕ, and magnetic field, bϕ, parameters for a high-frequency wave at ω = 0.209 rad s−1 corresponding to a period of 30 s. The real and imaginary components of the perturbations, as well as their amplitudes are shown.

6. Comparison of the perturbations

Based on the obtained solutions in the previous section, we compare the velocity and magnetic field perturbations for the three cases. The SciPySpecial package is used to calculate the Hankel functions. The spatial profiles of the solutions are plotted between the photosphere at s = −2 Mm and the corona at s = 1 Mm. The coronal base lies at s = 0 Mm and is indicated with a dotted vertical line.

A scale height of h = 869.5 km is used for the plots. This is to allow the equilibrium parameters between the photosphere at s = −2 Mm and the coronal base at s = 0 Mm to be consistent with those shown in Fig. 2: a density and magnetic field decrease by 8 and 2 orders of magnitude, respectively, and a corresponding Alfvén speed increase by two orders of magnitude, from 10 to 1000 km s⁻¹. The field lines expand by a factor of 10 for case I, a factor of 1000 for case II, and a factor of 100 for case III. It should be noted that the present model is restricted to a constant scale height. However, the actual scale height in the solar atmosphere is expected to vary with temperature.

In order to correctly compare the spatial behaviour of the Alfvén waves and their amplitudes, we have assumed that the amount of upward wave energy flux from the photosphere does not depend on the field line geometry in the atmosphere above. We show in the next section that |a|² is a measure of the upward energy flux in the lower region. Therefore, the solutions of the linearised equations are normalised with respect to the constant a and are displayed in arbitrary units.

A torsional Alfvén wave with zero frequency does not propagate and instead forms a twisted steady-state configuration of the magnetic field, which can resemble a vortex in the plasma. The twist in the magnetic field induces a rotational motion of the plasma around the vertical symmetry axis, similar to a vortex in traditional fluid dynamics. While a vortex in fluid dynamics typically refers to the swirling motion of fluid around a central core, a zero-frequency torsional Alfvén wave can produce a similar swirling motion in a magnetised plasma, where the magnetic field lines themselves are twisted.

Figure 4 displays the behaviour of the perturbations with zero frequency representing a vortex. No vortical motions are observed for cases I and II. However, case III displays a vortical motion that reaches the corona. For cases I and II the upgoing energy flux is balanced by the downgoing flux, whereas case III exhibits no downgoing energy flux due to the particular magnetic field geometry (see Sect. 7). Vortical motions have been detected at different atmospheric heights in numerical simulations and in solar observations (Bonet et al. 2008; Battaglia et al. 2021; Silva et al. 2024). Only the real parts of the perturbations are present. There is an increase in the velocity perturbation and a decrease in the magnetic field perturbation from the photosphere to the corona.

The imaginary parts of the perturbations appear as the frequency begins to increase. Figure 5 shows the real and imaginary parts of v_ϕ and b_ϕ for the three cases in a low-frequency oscillation regime, ω = 0.00175 rad s⁻¹. The upper panels show the real parts of v_ϕ and b_ϕ, the middle panels show their imaginary parts, and the lower panels display the corresponding amplitudes.

The plots show no oscillations in v_ϕ or b_ϕ. Instead, their spatial profiles are characterised by a continuous increase or a decrease between the photosphere and the coronal base.

The amplitudes in v_ϕ increase with distance s, whereas b_ϕ decreases in magnitude. Case III possesses the largest increase in v_ϕ, especially in the real part. The corresponding magnetic field in the lower atmosphere is comparatively high among the three cases. We show in the next section that this distinctive behaviour for case III is consistent with the wave energy flux having no downgoing component. Cases I and II both end up at the same low values at the coronal base at s = 0 Mm. In the lower region, s < 0 Mm, their values vary; the magnetic field for case II is higher.

Figure 6 displays the spatial profiles for v_ϕ and b_ϕ for a higher frequency oscillation of 0.209 rad s⁻¹ corresponding to a period of 30 s. The three cases still display an increase in |v_φ| and a decrease in |b_φ| with distance s, with case III showing the largest increase. The imaginary parts of the perturbations increase with increasing frequency. In contrast to Figs. 4 and 5 for the low-frequency cases, the spatial profiles in Fig. 6 are more similar with moderate variations from the photosphere to the corona. This trend can be explained by the gradual decrease in the downgoing flux with increasing frequency for cases I and II (see Sect. 7). Another feature of Fig. 6 corresponding to a higher frequency is the oscillatory nature of the curves. Case III does not show oscillations in the magnitude plots as the real and imaginary components are out of phase.

Figures 4–6 demonstrate that the oscillatory nature of the perturbations becomes more prominent with increasing frequency. This behaviour is consistent with the predictions of the WKB theory. In the WKB approximation, the solutions of Eq. (23) can be expressed in the following form:

$$\begin{array}{c}\hfill q_i\sim exp\pm \int \frac{\sqrt{{\mathrm{\Omega }}_i^2-{\omega }^2}}{c_A}\mathrm{d}s.\end{array}$$(46)

Here Ω_i represent the corresponding cut-off frequencies. At high frequencies, the difference Ω_i² − ω² is negative and the solutions become oscillatory, which is consistent with the profiles of the curves in Figs. 4–6. At low frequencies the difference Ω_i² − ω² is positive, and leads to non-oscillatory solutions.

In certain physical contexts, non-oscillatory solutions represent evanescent waves, which decay exponentially with distance from a localised source. In the present context, the non-oscillatory solutions contain both exponentially growing and decaying terms corresponding to the plus and minus signs in Eq. (46). These solutions contrast with the more familiar oscillatory wave solutions, which describe propagating waves.

7. The wave energy flux

In the previous section we explored the spatial profiles of the velocity and the magnetic field perturbations as well as the amplitudes in the solar atmosphere for three different cases. To further understand how the wave amplitudes in the corona are affected by the field line geometry, we derive expressions for the wave energy flux and establish relationships between the flux and the amplitude for each case.

Equations (11) can be combined to derive an equation of wave energy (Taroyan & Williams 2016; Taroyan et al. 2021).

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\left(\frac{W}{B_s}\right)+\frac{\partial }{\partial s}\left(\frac{F}{B_s}\right)=0,\end{array}$$(47)

where

$$\begin{array}{c}\hfill W=\frac{\rho V_{\phi }^2}{2}+\frac{B_{\phi }^2}{2{\mu }_0},F=-B_s\frac{V_{\phi }{B}_{\phi }}{{\mu }_0}\end{array}$$(48)

are the wave energy density and the wave energy flux, respectively.

A similar equation of energy can be derived for the real parts of the perturbations using the real parts of equations (11). The result is

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\left(\frac{W_r}{B_s}\right)+\frac{\partial }{\partial s}\left(\frac{F_r}{B_s}\right)=0,\end{array}$$(49)

where

$$\begin{array}{c}\hfill W_r=\frac{\rho \mathfrak{R}{(V_{\phi })}^2}{2}+\frac{\mathfrak{R}{(B_{\phi })}^2}{2{\mu }_0},F_r=-B_s\frac{\mathfrak{R}(V_{\phi })\mathfrak{R}(B_{\phi })}{{\mu }_0}.\end{array}$$(50)

By taking a time average over a single wave period we have

$$\begin{array}{c}\hfill \langle \frac{\partial }{\partial t}\left(\frac{W_r}{B_s}\right)\rangle =-\langle \frac{\partial }{\partial s}\left(\frac{F_r}{B_s}\right)\rangle .\end{array}$$(51)

For the left-hand side we have

$$\begin{array}{c}\hfill \langle \frac{\partial }{\partial t}\left(\frac{W_r}{B_s}\right)\rangle =\frac{1}{B_{s}P}\underset{0}{\overset{P}{\int }}\frac{\partial W_r}{\partial t}\mathrm{d}t=\frac{1}{B_{s}P}{W_r|}_{t=0}^P=0\end{array}$$(52)

because the variable W_r is sinusoidal in time. Therefore, for the time-averaged energy flux on the right-hand side of Eq. (51) we have

$$\begin{array}{c}\hfill \frac{⟨F_r⟩}{B_s}=\mathrm{const}.\mathrm{along}\mathrm{field}\mathrm{lines}.\end{array}$$(53)

In practice, it is more convenient to deal with the complex variables rather than their real parts. We therefore express the flux, F_r, through the complex perturbations, v_φ and b_φ:

$$\begin{array}{cc}\hfill \frac{⟨F_r⟩}{B_s}& =-\frac{1}{{\mu }_{0}P}\underset{0}{\overset{P}{\int }}\mathfrak{R}(V_{\phi })\mathfrak{R}(B_{\phi })\mathrm{d}t\hfill \\ \hfill & =-\frac{1}{4{\mu }_{0}P}\underset{0}{\overset{P}{\int }}[v_{\phi }exp\left(i\omega t\right)+v_{\phi }^{\ast }exp(-i\omega t)]\hfill \\ \hfill & \times [b_{\phi }exp\left(i\omega t\right)+b_{\phi }^{\ast }exp(-i\omega t)]\mathrm{d}t=-\frac{1}{2{\mu }_0}\mathfrak{R}\left(v_{\phi }^{\ast }{b}_{\phi }\right).\hfill \end{array}$$(54)

Here the asterisk (*) denotes complex conjugation, and the integral is evaluated using the relationship between period and frequency, P = 2π/ω. Combining Eqs. (53) and (54) we obtain

$$\begin{array}{c}\hfill \frac{⟨F_r⟩}{B_s}=-\frac{1}{2{\mu }_0}\mathfrak{R}\left(v_{\phi }^{\ast }{b}_{\phi }\right)=\mathrm{constant}\mathrm{along}\mathrm{field}\mathrm{lines}.\end{array}$$(55)

The obtained relationship (55) is valid for all s and can be used to derive expressions for the energy flux in each of the three cases shown in Table 1. This derivation also verifies that the ratio of the energy flux to the magnetic field strength, B_s, is conserved along the field lines.

Previous applications of the conservation law (55) mainly focused on models with a uniform magnetic field (Hollweg et al. 1982; Tsap & Kopylova 2021). Equation (55) can be used to explore the relationships between the cut-off frequencies, the wave energy flux, and the corresponding amplitudes in general magnetic field geometries.

7.1. Case I

The expressions (25) and (27) for the velocity and magnetic field perturbations, v_φ and b_φ, in the region s < 0 can be substituted in the formula (55) for the energy flux to obtain

$$\begin{array}{cc}\hfill \frac{⟨F_r⟩}{B_s}& =-\frac{\omega }{2{\mu }_0}\mathfrak{R}\left(v_{\phi }^{\ast }{b}_{\phi }\right)\hfill \\ \hfill & =\frac{\omega }{2{\mu }_0{c}_A}\mathfrak{R}\left[i(a^{\ast }{H}_0^{(1)\ast }+b^{\ast }{H}_0^{(2)\ast })(a{H}_1^{(1)}+b{H}_1^{(2)})\right]=\frac{\omega }{2{\mu }_0{c}_A}\hfill \\ \hfill & \times \mathfrak{R}\left[i(a^{\ast }{J}_0-i{a}^{\ast }{Y}_0+b^{\ast }{J}_0+i{b}^{\ast }{Y}_0)(a{J}_1+ia{Y}_1+b{J}_1-ib{Y}_1)\right]\hfill \\ \hfill & =\frac{\omega }{2{\mu }_0{c}_A}\mathfrak{R}\left[{i|a|}^2(J_0{J}_1+i{J}_0{Y}_1-i{Y}_0{J}_1+Y_0{Y}_1)\right]\hfill \\ \hfill & +\frac{\omega }{2{\mu }_0{c}_A}\mathfrak{R}\left[{i|b|}^2(J_0{J}_1-i{J}_0{Y}_1+i{Y}_0{J}_1+Y_0{Y}_1)\right]\hfill \\ \hfill & +\frac{\omega }{2{\mu }_0{c}_A}\mathfrak{R}\left[i{a}^{\ast }b(J_0{J}_1-i{J}_0{Y}_1-i{Y}_0{J}_1-Y_0{Y}_1)\right]\hfill \\ \hfill & +\frac{\omega }{2{\mu }_0{c}_A}\mathfrak{R}\left[ia{b}^{\ast }(J_0{J}_1+i{J}_0{Y}_1+i{Y}_0{J}_1-Y_0{Y}_1)\right],\hfill \end{array}$$(56)

where the arguments of the functions, $\frac{h\omega }{2c_A}$, are omitted for brevity. The sum of the last two terms in (56) is 0 because ℜ(iz) + ℜ(iz * ) = 0. Therefore,

$$\begin{array}{cc}& \frac{⟨F_r⟩}{B_s}=\frac{\omega }{2{\mu }_0{c}_A}{(|a|}^2-{|b|}^2)(J_1{Y}_0-J_0{Y}_1).\hfill \end{array}$$(57)

Using the well-known property of the Wronskian for the Bessel functions, W(J₀, Y₀) = J₁Y₀ − J₀Y₁ = 2π/z, where z represents the argument, we obtain

$$\begin{array}{c}\hfill \frac{⟨F_r⟩}{B_s}=-\frac{1}{2{\mu }_0}\mathfrak{R}\left(v_{\phi }^{\ast }{b}_{\phi }\right)=\frac{2\pi }{{\mu }_{0}h}{(|a|}^2-{|b|}^2).\end{array}$$(58)

A similar expression was presented by Hollweg (1978), and it demonstrates that the ratio of the time-averaged energy flux to the field strength, ⟨F_r⟩/B_s, is conserved along the magnetic field lines. The formula can be interpreted as the difference between the upward and downward fluxes. In the present paper the net flux vanishes when upward and downward fluxes balance each other.

On the other hand, for the upper atmosphere, s > 0, we have the expressions (29). The ratio of the flux to the background magnetic field can therefore be expressed in an alternative format:

$$\begin{array}{c}\hfill \frac{⟨F_r⟩}{B_s}=-\frac{1}{2{\mu }_0}\mathfrak{R}\left(v_{\phi }^{\ast }{b}_{\phi }\right)=\frac{{|c|}^2}{2{\mu }_0{c}_{A0}}.\end{array}$$(59)

Therefore, the net wave energy flux can be expressed either as the difference between the upward and downward fluxes in the lower atmosphere (58) or as an upward flux in the upper atmosphere (59). The latter expression establishes a direct link between the energy flux and the amplitudes in the corona.

7.2. Case II

The expressions (33) and (34) for the magnetic field and velocity perturbations, b_φ and v_φ, in the region s < 0 can be substituted in the formula for the energy flux to obtain

$$\begin{array}{cc}\hfill \frac{⟨F_r⟩}{B_s}& =-\frac{1}{2{\mu }_0}\mathfrak{R}\left(v_{\phi }^{\ast }{b}_{\phi }\right)\hfill \\ \hfill & =-\frac{\omega }{2{\mu }_0{c}_A}\mathfrak{R}\left[i(a^{\ast }{H}_1^{(1)\ast }+b^{\ast }{H}_1^{(2)\ast })(a{H}_0^{(1)}+b{H}_0^{(2)})\right]=-\frac{\omega }{2{\mu }_0{c}_A}\hfill \\ \hfill & \times \mathfrak{R}\left[i(a^{\ast }{J}_1-i{a}^{\ast }{Y}_1+b^{\ast }{J}_1+i{b}^{\ast }{Y}_1)(a{J}_0+ia{Y}_0+b{J}_0-ib{Y}_0)\right]\hfill \\ \hfill & =-\frac{\omega }{2{\mu }_0{c}_A}\mathfrak{R}\left[{i|a|}^2(J_0{J}_1+i{J}_1{Y}_0-i{Y}_1{J}_0+Y_1{Y}_0)\right]\hfill \\ \hfill & +\frac{\omega }{2{\mu }_0{c}_A}\mathfrak{R}\left[{i|b|}^2(J_0{J}_1-i{J}_1{Y}_0+i{Y}_1{J}_0+Y_0{Y}_1)\right]\hfill \\ \hfill & +\frac{\omega }{2{\mu }_0{c}_A}\mathfrak{R}\left[i{a}^{\ast }b(J_0{J}_1-i{J}_1{Y}_0-i{Y}_1{J}_0-Y_0{Y}_1)\right]\hfill \\ \hfill & +\frac{\omega }{2{\mu }_0{c}_A}\mathfrak{R}\left[ia{b}^{\ast }(J_0{J}_1+i{J}_1{Y}_0+i{Y}_1{J}_0-Y_0{Y}_1)\right],\hfill \end{array}$$(60)

where the arguments of the functions, $\frac{h\omega }{2c_A}$, are once again omitted for brevity. As in the previous case, the sum of the last two terms in (60) is 0 because ℜ(iz) + ℜ(iz * ) = 0. Combining this with the expression for the Wronskian, W(J₀, Y₀) = J₁Y₀ − J₀Y₁ = 2π/z, where z represents the argument, we obtain the same expression (58) derived for the previous case. It confirms that the ratio of the time-averaged energy flux to the field strength is conserved along the magnetic field lines, in agreement with (55).

A procedure similar to the one for case I can be carried out to establish expression (59) for case II. Therefore, expressions (58) and (59) are applicable to both degenerate cases I and II.

7.3. Case III

The magnetic field and velocity perturbations, b_φ and v_φ, in the lower atmosphere, s < 0, are given by the expressions (42) and (44). In Sect. 5, we establish that the coefficient b = 0 when the perturbation in the upper atmosphere is outgoing. Therefore, the following expression is obtained for the energy flux:

$$\begin{array}{c}\hfill \frac{⟨F_r⟩}{B_s}=-\frac{1}{2{\mu }_0}\mathfrak{R}\left(v_{\phi }^{\ast }{b}_{\phi }\right)=\frac{{2\pi |a|}^2}{{\mu }_{0}h}.\end{array}$$(61)

The above formula indicates that for the field line expansion representative of case III there is only an upward flux of energy. Using the solutions in the region s > 0 the same energy flux can be expressed in the form (59) with a different constant $\stackrel{\sim }{c}$.

Fig. 7. Ratio of the downgoing to upgoing energy fluxes as a function of frequency. The overlaying curves indicate that the ratio, R, is identical for cases I and II.

Fig. 8. Comparison of the time-averaged flux (⟨Fr⟩/Bs,) for the three cases over a range of frequencies up to a frequency of 0.209 rad s−1 representing a period of 30 s.

7.4. Comparison of fluxes

The three principles set out in Sect. 5 allow us to construct physically acceptable solutions and to compare the wave amplitudes and the corresponding energy fluxes for the three different cases. Hollweg (1978) introduced the ratio of the downgoing energy flux to the upgoing energy flux, $R={\left|\frac{b}{a}\right|}^2$. Previous theoretical studies suggested that vortices and low-frequency waves are strongly suppressed in the atmosphere (Hollweg 1978; Leer et al. 1982; Tsap & Kopylova 2021).

Figure 7 demonstrates the change in the flux ratio, R, with frequency, ω. It gradually decreases with increasing frequency as the proportion of downgoing flux decreases, and the waves transfer more energy into the upper atmosphere. Cases I and II are identical, as they are both given by Eq. (30).

On the other hand, for case III R = 0, indicating no downward flux, in agreement with the coefficient b being 0 from Equation (45). The energy flux enters the corona unhindered. This behaviour is consistent across the entire frequency spectrum.

The ratio of the time-averaged energy flux to the background field (55) is plotted in Fig. 8. The three cases are shown on log scales for various frequencies up to a frequency of 0.2 rad s⁻¹ or a corresponding period of 30 s. Most notably, it can be seen for case III that the flux is significantly higher than the fluxes for cases I and II, especially at low frequencies. Cases I and II have the same fluxes, which is in agreement with the ratio R being the same. Both vanish in the limit of ω = 0, whereas the flux for case III does not change with frequency and remains consistently high. The fluxes for the three cases merge when the frequency increases because the proportion of the downgoing flux becomes small.

8. Summary and conclusions

The role of Alfvén waves in heating solar or stellar atmospheres and in wind acceleration remains one of the most important questions in solar and stellar physics. Conclusions based on theoretical models suggest that low-frequency Alfvén waves, which are the least susceptible to damping (Zaqarashvili et al. 2013), are unable to transport sufficient amounts of energy to the upper atmosphere of the Sun: a strong downward flux of wave energy is developed due to strong spatial gradients.

These views have been challenged by recent observations, which indicate the important role of these waves in the dynamics and energetics of the upper atmosphere (De Pontieu et al. 2014). Additionally, both observations and simulations indicate that vortices and swirling motions can lead to efficient transport of energy with subsequent heating (Bonet et al. 2008; Wedemeyer-Böhm et al. 2012; Breu et al. 2023). Revised models are required to reassess established views on wave behaviour and to reconcile observations and simulations with theory.

In the present paper we addressed the role of the magnetic field line geometry. It was shown that the propagation of Alfvén waves in arbitrary magnetic field geometries is governed by the Klein–Gordon equations for the magnetic field and velocity perturbations. The obtained results demonstrate that the energy flux reaching the corona and the corresponding perturbation amplitudes strongly depend on the magnetic field line expansion rate in the lower atmosphere. The flux is suppressed when for a given expansion rate either of the governing Klein–Gordon equations becomes degenerate and the corresponding cut-off frequency vanishes. A strong downward flux is developed that balances the upgoing flux (i.e. the ratio of the downgoing to upgoing fluxes, R = 1). This may explain why low-frequency Alfvén waves and vortices do not transport sufficient energy to heat the atmosphere in commonly adopted models. Examples include models based on a uniform magnetic field, a thin flux tube, or a radially expanding field. As a counterexample, a model is constructed where none of the cut-off frequencies vanish, and the Klein–Gordon equations remain non-degenerate, and lead to large perturbation amplitudes.

The energy flux is proportional to the perturbation amplitudes. Expressions for the time-averaged flux are derived and compared for three different cases. These cases have different expansion rates for the background magnetic field in the lower atmosphere. For cases I and II the flux has both downgoing and upgoing components, whereas the intermediate case III only has an upward component.

The amplitudes and the associated energy flux are strongly suppressed when the cut-off frequencies are zero. The cut-off frequencies mark the transition from higher frequency propagating waves to lower frequency growing and decaying waves. A vanishing cut-off frequency implies that even at low frequencies the solutions to the corresponding Klein–Gordon equations have an oscillatory behaviour. These propagating waves are reflected off the interface between the lower and the upper layers. A strong downward flux associated with the reflected counterpart balances the upward flux. Therefore, the amplitudes of the vortices and low-frequency waves become suppressed, in agreement with Figs. 4 and 5. Reflection is less efficient at higher frequencies (Fig. 7), and, therefore, the wave amplitudes and the energy flux reaching the corona increase with increasing frequency for cases I and II (see Figs. 6 and 8).

The intermediate case III demonstrates that non-zero cut-off frequencies in the Klein–Gordon equations are associated with large perturbation amplitudes and an energy flux that is significantly higher compared to the degenerate cases. The perturbations do not propagate and there is no reflection at the interface. The solutions have growing and decaying spatial profiles. The magnetic field perturbation decreases with distance, whereas the velocity perturbation increases as the energy flux should be conserved.

An increase in the velocity perturbations with increasing height (Fig. 4) is consistent with the presence of swirls and torsional motions in simulations and observations of the solar atmosphere (De Pontieu et al. 2014; Battaglia et al. 2021; Silva et al. 2024). When the frequency of the perturbations increases they become more oscillatory. The net energy flux remains high across the entire frequency spectrum as there is no downgoing component.

When Alfvénic motions are confined to a symmetry axis, they are usually seen to be highly reflective. This conclusion is also valid for models with a uniform magnetic field. The counterexample presented in this study shows the important role of the magnetic field line geometry. The energy flux reaching the corona can be consistently high across the entire frequency spectrum.

We focused on a particular field line geometry with different expansion rates in the lower atmosphere and exponential spatial profiles corresponding to a constant scale height. We need to ask the question of how sensitive the obtained results are to the background atmosphere. Future studies should look at wave propagation in an atmosphere with a more complex magnetic field line geometry. The preceding discussion is based on the general properties of the Klein–Gordon equations. We could therefore argue that the link between the cut-off frequencies and the energy flux carried by the Alfvén waves should hold both in simple and in more complex structures. There should be efficient transfer of energy (even if not 100%) as long as the cut-off frequencies do not vanish and there is no degeneracy in the corresponding Klein–Gordon equations.

The present study is restricted to the linear regime. An important question is how non-linearity can affect the transfer of energy flux and the reflectivity of the Alfvén waves. As the wave amplitude grows, Alfvén waves can become distorted, leading to steepening and the formation of shock-like structures. Non-linear interactions can also lead to the coupling of Alfvén waves with other wave modes, such as fast and slow magnetoacoustic waves, which can alter the reflectivity (Hollweg et al. 1982; Matsumoto & Shibata 2010). These coupled modes may have different reflectivity characteristics than pure Alfvén waves.

The transition from linear to non-linear behaviour is expected to lead to a more complex interaction with the solar atmosphere, where more energy is dissipated or converted into other forms. Both long- and short-period waves can experience non-linear effects, but short-period waves are more likely to steepen and dissipate. Detailed studies in the non-linear regime with different field line geometries would be required for more accurate conclusions.

Acknowledgments

T. Borradaile would like to thank the STFC for their financial support (ST/S505225/1).

References

All Tables

All Figures

Fig. 1. Axisymmetric poloidal magnetic field described by cylindrical coordinates and local orthogonal curvilinear coordinates. The angle ϕ representing the angle of rotation around the vertical symmetry axis, z, is not shown. The figure shows distance, s, along a curved magnetic field; the angle, θ, between the poloidal field and the vertical axis; and incremental increases δs and δa, along and across the poloidal field, respectively.

In the text

	Fig. 2. Equilibrium parameters in the atmosphere. The panels on the left display the variation of density, Alfvén speed and magnetic field. The panels on the right show the radial distance, r = r(s), and the corresponding parameters that determine the cut-off frequencies Ω1 and Ω2. Logarithmic scales are used.
In the text

	Fig. 3. Geometry of the magnetic field in the r − z plane. The field lines exponentially expand in the lower atmosphere and are vertical in the uniform upper atmosphere.
In the text

	Fig. 4. Comparison of the velocity, vϕ, and magnetic field, bϕ, perturbations for a zero frequency wave corresponding to vortex motion. Both the real and imaginary components, as well as the amplitudes are shown. The dotted vertical line indicates the interface between the lower and the upper regions.
In the text

	Fig. 5. Comparison of the velocity, vϕ, and magnetic field, bϕ, parameters for a low-frequency wave at ω = 0.00175 rad s−1 corresponding to a period of 1 hour. Both the real and imaginary components, as well as the amplitudes are shown.
In the text

	Fig. 6. Comparison of the velocity, vϕ, and magnetic field, bϕ, parameters for a high-frequency wave at ω = 0.209 rad s−1 corresponding to a period of 30 s. The real and imaginary components of the perturbations, as well as their amplitudes are shown.
In the text

	Fig. 7. Ratio of the downgoing to upgoing energy fluxes as a function of frequency. The overlaying curves indicate that the ratio, R, is identical for cases I and II.
In the text

	Fig. 8. Comparison of the time-averaged flux (⟨Fr⟩/Bs,) for the three cases over a range of frequencies up to a frequency of 0.209 rad s−1 representing a period of 30 s.
In the text