© The Authors 2023

1. Introduction

Recent advances in observations and modelling have shown that magnetohydrodynamic (MHD) waves could significantly contribute to the heating of the solar corona (see review by Van Doorsselaere et al. 2020). In particular, transverse waves are ubiquitously observed, and several kinds have been identified. The type that was first discovered are the transverse waves that are impulsively excited after a flare (Nakariakov et al. 1999). However, these transverse waves are only sporadically excited and do not play an important role in the energy budget of the solar corona (Terradas & Arregui 2018). Later on it was discovered that the corona is filled by small-amplitude transverse waves (Tomczyk et al. 2007; Tomczyk & McIntosh 2009; McIntosh et al. 2011; Tian et al. 2012). These were observed in coronal loops as propagating (Tiwari et al. 2019) or standing waves (Anfinogentov et al. 2015). These low-amplitude transverse waves were also observed as propagating waves in open-field regions (Thurgood et al. 2014; Morton et al. 2015). These low-amplitude waves show little to no decay (Morton et al. 2021), and are thus named ‘decayless’.

Because the flare-excited standing waves rapidly decay (Goddard et al. 2016; Nechaeva et al. 2019) due to resonant absorption (Goossens et al. 2002) and non-linear Kelvin–Helmholtz instability (KHI) damping (Terradas et al. 2008; Antolin et al. 2014; Van Doorsselaere et al. 2021; Arregui 2021), it is generally thought that the decayless waves must be continuously supplied with energy to counteract the strong damping. Several mechanisms for excitation have been proposed: slip-stick driving with steady flows (Nakariakov et al. 2016; Karampelas & Van Doorsselaere 2020), vortex shedding (Nakariakov et al. 2009; Karampelas & Van Doorsselaere 2021), and footpoint driving (Nisticò et al. 2013; Karampelas et al. 2017) through p-modes (Morton et al. 2019) or convective shuffling. The option of footpoint driving has had some success in generating standing mode decayless waves (Afanasyev et al. 2020), which counterbalance the non-linear damping through the KHI (Guo et al. 2019) and heat the loops (Shi et al. 2021).

However, for the driving of decayless waves through their footpoints, it is not well understood how the transverse waves propagate through the complicated structure of the chromosphere and transition region. The simulations of transverse-wave-induced KHI heating (e.g., Karampelas et al. 2019) only take into account the coronal part of the loop; in other words, a driver is imposed at the top of the transition region. To properly model the whole loop evolution due to the wave heating, it is essential to also model the wave driver in the photosphere and accurately capture its influence on the coronal loop dynamics.

In plane-parallel atmospheres the propagation of fast and slow waves has been well studied. It was found that these modes couple efficiently to Alfvén waves through resonant absorption (Hansen & Cally 2009; Cally & Andries 2010; Khomenko & Cally 2012). Currently, investigations into what happens if the cross-field structuring is included in the wave propagation model are ongoing (Cally & Khomenko 2019; Riedl et al. 2019, 2021). Another crucial ingredient is the wave’s behaviour in strong (i.e., non-Wentzel-Kramers-Brillouin, non-WKB) stratification. It is well known that slow waves experience a cutoff while propagating through a stratified medium (Bel & Leroy 1977). This has been verified observationally (Jess et al. 2013) and numerically (Felipe et al. 2018), although it is still unknown whether a similar cutoff exists for transverse waves in structured media. For the driving of the observed decayless waves in the corona, this is a crucial property to understand.

Several analytical works predict that transverse waves are cut off in the transition below a given frequency. The first formula was derived by Spruit (1981):

$$\begin{array}{c}\hfill {\omega }_{\mathrm{Sp}81}^2=\frac{g}{8H}\frac{1}{2\beta +1}.\end{array}$$(1)

Here g is the gravity projected along the loop, H the pressure scale height, and β the ratio of the gas to the magnetic pressures. For a typical isothermal atmosphere, this corresponds to a cutoff period of 700 s (Spruit 1981). However, Lopin et al. (2014) showed that this classical cutoff is suppressed when the radial component of the magnetic field is taken into account. Lopin & Nagorny (2017) later showed that transverse waves can still be cut off, provided a non-isothermal atmosphere. They predict a cutoff frequency of

$$\begin{array}{c}\hfill {\omega }_{\mathrm{LN}17}^2=\frac{c_{k0}^2}{4H_{0}H(z)}({\delta }_B^2\frac{\mathrm{d}H(z)}{\mathrm{d}z}+\frac{H^2(z)}{z^2}),\end{array}$$(2)

where z is the altitude, c_k0 is the kink speed at the base of atmosphere (z = z₀), H is the pressure scale height, H₀ = H(z₀), and ${\delta }_B^2=\left(B_{\text{0i}}^2-B_{\text{0e}}^2\right)/\left(B_{\text{0i}}^2+B_{\text{0e}}^2\right)$ is the relative difference between the magnetic field inside (B_0, i) and outside (B_0, e) the flux tube at z = z₀. Finally, an alternative formula was derived by Snow et al. (2017),

$$\begin{array}{c}\hfill {\omega }_{\mathrm{Sn}17}^2=\frac{v_{\mathrm{A}}^2(z)}{4z^2},\end{array}$$(3)

where z is the altitude and v_A is the Alfvén speed.

In this article we modelled the propagation of kink waves in an open magnetic flux tube, embedded in a non-isothermal atmosphere. The atmosphere extends from the chromosphere to the corona, and includes gravitational stratification and thermal conduction (Sect. 2). We drove kink waves at different periods, and determined whether they experienced a cutoff (Sect. 3). We compare these results to the three analytical formulas given above in Sect. 4, and summarize our conclusions in Sect. 5.

2. Numerical model: Magnetic flux tube through the transition region

We modelled a vertical magnetic flux tube of radius R = 1 Mm embedded in a stratified atmosphere, starting in the chromosphere (altitude z = 0 Mm) and extending through the transition region (z ≈ 4 Mm) into the corona. Kink waves were excited in the flux tube by applying a monoperiodic driver at the bottom of the domain (z = 0 Mm). In the upper half of the domain (z > 50 Mm), we implemented a ‘velocity rewrite layer’ to absorb the kink waves. The driver and the velocity rewrite layer are described in Sect. 2.1. A sketch of the domain is shown in Fig. 1. We solved the 3D MHD evolution of this tube using the PLUTO code (Mignone et al. 2007), version 4.3. This code solves the conservative MHD equations (mass continuity, momentum conservation, energy conservation, and induction equation). We used the corner transport upwind finite volume scheme, where characteristic tracing is used for the time stepping and a linear spatial reconstruction with a monotonized central difference limiter is performed. The magnetic field divergence was kept small using the extended divergence cleaning method (generalized Lagrange multiplier, or GLM), and the flux was computed with the linearized Roe Riemann solver. We did not include explicit viscosity, resistivity, or cooling. However, numerical dissipation results in higher effective viscosity and resistivity than expected for the solar corona, as discussed by Karampelas et al. (2019). We included a modified thermal conduction, as described below.

Fig. 1. Sketch of the simulation domain, showing the magnetic flux tube, the location of the kink wave driver (bottom boundary), chromosphere, transition region, corona, and velocity rewrite layer.

The transition region between the chromosphere and the corona is characterized by a very sharp temperature gradient. Resolving this gradient requires a very high resolution along the tube (∼1 km in the transition region). In order to keep computational costs reasonable, we artificially broadened the transition region (thus reducing the temperature gradient). To that end, we modified the thermal conductivity using the method developed by Linker et al. (2001), Lionello et al. (2009), and Mikić et al. (2013). Below the cutoff temperature T_c = 2.5 × 10⁵ K, the parallel thermal conductivity was set to ${\kappa }_{\Vert }=C_0{T}_{\mathrm{c}}^{5/2}$ with C₀ = 9 × 10⁻¹² W m⁻¹ K^−7/2; above T_c, it was set to κ_∥ = C₀T^5/2. This allowed us to use a resolution of 98 km along the tube. This grid allows us to fully resolve the broadened transition region, which has a minimum temperature scale length of 1.6 Mm (see Johnston & Bradshaw 2019). The dimensions of the domain were (L_x, L_y, L_z) = (16, 6, 100) Mm. We used a uniform grid of 400 × 150 × 1024 cells, with a size of 40 km in the x and y directions, and 98 km in the z direction. Furthermore, we verified that the results did not change significantly when using a resolution of 40 km in the z direction. To that end, we ran a separate simulation and verified that the resulting cutoff altitude and comparison to the analytical formulas (see Sect. 4) were not strongly modified. We note that this resolution is too costly in terms of running time to be used for all simulations in this work.

The strong stratification in the transition region makes it challenging to obtain a relaxed initial state for the model. We first initialized the domain with a field-aligned hydrostatic equilibrium (Sect. 2.2). We then let the simulation relax in 2D for 47 ks (Sect. 2.3). Finally, we filled the 3D domain with this relaxed state through cylindrical symmetry, where we drove kink waves of different periods for a duration up to 2.7 ks (Sect. 2.4).

2.1. Boundary conditions and driver

We first describe the boundary conditions used for the relaxation (2D) and kink wave (3D) simulations.

Bottom boundary. At the bottom boundary (base of the chromosphere, z = 0), the density and pressure were extrapolated using the hydrostatic equilibrium equation. The magnetic field was extrapolated using the zero normal-gradient condition described by Karampelas et al. (2019, Sect. 2.4). For v_z, we either imposed a reflective boundary condition (2D relaxation, see Sect. 2.3), or imposed v_z = 0 (in 3D, see Sect. 2.4). We verified that both boundary conditions give the same results in 3D simulations. The parallel velocity components v_x and v_y were set to obey either a zero-gradient boundary condition (2D relaxation), or to follow a driver that excites kink waves (in 3D). We used a monoperiodic, dipole-like, driver developed by Pascoe et al. (2010) and updated by Karampelas et al. (2017). Inside the tube the driver imposes

$$\begin{array}{c}\hfill \{v_x(x,y,t),v_y(x,y,t)\}=\{v(t),0\},\end{array}$$(4)

where v(t) = v₀cos(2πt/P₀), with v₀ the driver amplitude set to 2 km s⁻¹. The driver period, P₀, was set to different values in order to test the cut off of kink waves. Outside the tube the driver imposes

$$\begin{array}{c}\hfill \{v_x(x,y,t),v_y(x,y,t)\}=v(t)R^2\frac{\{{(x-x_0(t))}^2-y^2,2(x-x_0(t))y\}}{{({(x-x_0(t))}^2+y^2)}^2},\end{array}$$(5)

where x₀(t) = v₀P₀/(2π)⋅sin(2πt/P₀) is the centre of the tube’s footpoint at time t. This driver generates a kink wave polarized in the x direction.

Upper boundary. At the upper boundary (top of the corona, z = 100 Mm), the magnetic field was kept symmetric. All other variables obeyed a reflective boundary condition. In order to absorb the upward waves excited by the driver, we artificially modified the velocity in the upper half of the domain (z > 50 Mm). At each time step, after solving the MHD equations, we decreased each component of the velocity v_i by multiplying it by a quantity α_v ≲ 1:

$$\begin{array}{c}\hfill v_{\mathrm{i}}^{\prime }={\alpha }_v(t,z)v_{\mathrm{i}}.\end{array}$$(6)

In the driven 3D simulations α_v was kept constant in time and varied linearly along the loop, from 1 at z = z_v = 50 Mm to α_v,min = 0.9995 at z = L = 100 Mm:

$$\begin{array}{c}\hfill {\alpha }_{v,3\mathrm{D}}(z)=\{\begin{array}{cc}1\hfill & \text{if}z\le z_v,\hfill \\ 1-(1-{\alpha }_{v,\min })\left(\frac{z-z_v}{L-z_v}\right)\hfill & \text{else.}\hfill \end{array}\end{array}$$(7)

In the 2D relaxation run, the first third of the simulation (t_1/3 = 15.7 ks) was run without modifying the velocity (i.e., α_v = 1). During the second third α_v was linearly ramped down in time to match the profile α_v, 3D(z) described above. Finally, the last third of the simulation was run with the constant α_v, 3D(z):

$$\begin{array}{c}\hfill {\alpha }_{v,2\mathrm{D}}(z,t)=\{\begin{array}{cc}1\hfill & \text{if}t\le t_{1/3},\hfill \\ 1-(1-{\alpha }_{v,3\mathrm{D}}(z))\left(\frac{t-t_{1/3}}{t_{1/3}}\right)\hfill & \text{if}{t}_{1/3}<t\le 2t_{1/3},\hfill \\ {\alpha }_{v,3\mathrm{D}}(z)\hfill & \text{else.}\hfill \end{array}\end{array}$$(8)

The evolution of α_v is shown in Fig. 2. This velocity rewrite layer can successfully absorb the kink waves that are excited by the driver at the bottom of the chromosphere. As a result, these waves are not reflected at the upper boundary, and do not propagate downwards back into the domain. We note that the solution obtained inside the velocity rewrite layer (i.e., above z = 50 Mm) is not physical, and that this layer should be considered a part of the upper boundary.

Fig. 2. Velocity-rewrite coefficient αv, applied to the velocity above 50 Mm so that upper-propagating waves are not reflected back into the domain. αv is shown for different times of the 2D relaxation run. The last profile (t ≥ 31.3 ks) is also applied in the 3D driven simulations.

Side boundaries. At the side boundaries (x and y axes), all variables obeyed a zero-gradient boundary condition. In the 2D relaxation run, we only simulated half of the tube radius (x > 0). For these simulations, we imposed a reflective boundary condition on all variables at the centre of the tube (x = 0).

2.2. Initial conditions: Field-aligned hydrostatic equilibrium

The simulation was initialized with a uniform vertical magnetic field of magnitude B₀ = 42 G. Along the tube, we imposed the temperature profile derived from Aschwanden & Schrijver (2002)

$$\begin{array}{c}\hfill T(x,y,z)=\{\begin{array}{cc}{T}_{\mathrm{ch}}\hfill & \text{if}z\le {\mathrm{\Delta }}_{\mathrm{ch}},\hfill \\ T_{\mathrm{ch}}+(T_{\mathrm{cor}}(x,y)-T_{\mathrm{ch}}){(1-{\left(\frac{L-z}{L-{\mathrm{\Delta }}_{\mathrm{ch}}}\right)}^2)}^{0.3}\hfill & \text{else,}\hfill \end{array}\end{array}$$(9)

where z is the altitude, L is the height of the computational domain, Δ_ch = 4 Mm is thickness of the chromosphere, and T_ch = 20 000 K is the temperature in the chromosphere. We defined the transverse temperature profile at the top of the domain, T_cor(x, y), as

$$\begin{array}{c}\hfill T_{\mathrm{cor}}(x,y)=T_{\mathrm{cor},\mathrm{ext}}+(T_{\mathrm{cor},\mathrm{int}}-T_{\mathrm{cor},\mathrm{ext}})\zeta (x,y),\end{array}$$(10)

where T_cor, int = 1.2 MK is the temperature inside the tube, and T_cor, ext = 3.6 MK is the temperature outside the tube. The shape of the profile was set by ζ(x, y),

$$\begin{array}{c}\hfill \zeta (x,y)=\frac{1}{2}[1-tanh\left((\sqrt{x^2+y^2}/R-1)b\right)],\end{array}$$(11)

where R = 1 Mm is the tube radius and b = 5 is a dimensionless number setting the width of the inhomogeneous layer between the interior and exterior of the tube (l ≈ 6R/b). The value of ζ(x, y) is close to 1 inside the tube, and to 0 outside.

We also set the density at the bottom of the chromosphere (z = 0) to

$$\begin{array}{c}\hfill {\rho }_{\mathrm{ch}}(x,y,z=0)={\rho }_{\mathrm{ch},\mathrm{ext}}+({\rho }_{\mathrm{ch},\mathrm{int}}-{\rho }_{\mathrm{ch},\mathrm{ext}})\zeta (x,y),\end{array}$$(12)

where ρ_ch, int = 3.51 × 10⁻⁸ kg m⁻³ is the density inside the tube and ρ_ch, ext = 1.17 × 10⁻⁸ kg m⁻³ is the density outside. We then integrated the field-aligned hydrostatic equilibrium equation numerically using a Crank–Nicholson scheme. The profiles of the imposed temperature and of the density resulting from the integration are shown in Fig. 3a. The temperature contrast (interior temperature divided by exterior temperature) is 1 in the chromosphere, and decreases to 1/3 in the corona. The density contrast is 3 in the chromosphere, increases to around 7 in the transition region, and decreases again to about 4 in the upper corona. The pressure contrast is 3 in the chromosphere, and slowly decreases to reach 1.2 in the upper corona.

Fig. 3. Temperature (black), density (red), and magnetic field magnitude (blue) profiles inside (r = 0 Mm; solid lines) and outside (r = 8 Mm; dashed lines) the flux tube. (a) After solving the field-aligned hydrostatic equilibrium. (b) After the 2D magnetohydrodynamic relaxation.

However, this initial state is not in magnetohydrostatic (MHS) equilibrium because the pressure varies across the flux tube, while the magnetic field does not. To fix this, we let the tube relax by running a 2D magnetohydrodynamic simulation (Sect. 2.3). We then used this relaxed state to initialize the 3D simulation of kink waves (Sect. 2.4).

2.3. Flux tube relaxation (2D)

In order to obtain a flux tube in MHS equilibrium, we first ran a 2D simulation, initialized with the initial state described in Sect. 2.2. The MHD equations were solved in a longitudinal plane at y = 0 (see Fig. 1), with x ∈ [0, 8.56] Mm and z ∈ [0, 100] Mm. We used a uniform grid of 64 × 2048 cells with a size of 134 km × 49 km. The resolution along z was higher than in the 3D runs in order to resolve the sharper gradients in the transition region (see Fig. 3). We verified that a resolution of 40 km in the x direction yielded the same results by running a separate 2D simulation followed by a 3D driven simulation (P₀ = 200 s), and verifying that the cutoff altitude and comparison to the analytical formulas (Sect. 4) were not significantly modified.

We let the system evolve for 47 ks, during which the velocity rewrite parameter α_v varied as described in Eq. (8). As a result of the relaxation, periodic longitudinal flows with a velocity of about 15 km s⁻¹ develop along the tube. They are damped during the later stages of the simulation, as the velocity rewrite layer is gradually introduced. At the end of the relaxation run, residual velocities are lower than 0.5 km s⁻¹ everywhere in the domain. The resulting temperature, density, and magnetic field profiles are shown in Fig. 3b. Compared to the initial state (Fig. 3a), the transition region is significantly broadened, with a thickness of about 7 Mm. This is the direct result of the modified thermal conductivity used in this set-up, and allows for a coarser resolution along the loop in the 3D simulations. In addition, the temperature and density decrease, both inside and outside the tube. Overall, the density contrast (ρ_int/ρ_ext) decreases: it reaches 1 in the chromosphere, 1.2 in the transition region, and 1.8 in the corona. The temperature contrast also changes to about 1.3 in the transition and about 0.8 in the corona. Finally, the magnetic field amplitude contrast remains very close to 1 everywhere in the domain (0.97 in the chromosphere and 1 in the corona), with a magnitude of about 11 G everywhere in the domain. Compared to the initial uniform magnetic field, the magnitude is divided by about four, while the contrast remains close to 1. The final temperature and density profile significantly differ from the initial conditions of the 2D relaxation run. However, this is not an issue as the goal of this study is to investigate how the analytical formulas we consider (Spruit 1981; Lopin & Nagorny 2017; Snow et al. 2017) predict the cutoff frequency for a given temperature and density profile. By using the relaxed profiles as input to these analytical formulas, we obtained predictions for the relaxed system.

This relaxed 2D simulation was then mapped onto the 3D domain through cylindrical symmetry. We used a rotation about the line x = 0 (i.e., the centre of the loop), and a trilinear interpolation to project onto the 3D Cartesian grid.

2.4. Kink wave propagation (3D)

In order to simulate the propagation of kink waves from the chromosphere to the corona, we drove the 3D simulations with the monoperiodic, dipole-like, driver described in Eqs. (4) and (5). We ran four simulations, with different driver periods P₀: 200 s, 335 s, 700 s, and 2000 s. The propagating kink waves generated by the driver are absorbed by the velocity rewrite layer at the top of the domain, and are thus not reflected downwards. The first three simulations were run for a duration of 5P₀. The last simulation was run for 1.75P₀. At the beginning of the simulations, the system goes through an initial transitory phase before the propagating kink wave is fully established (i.e. its amplitude does not change with time). We waited for 2P₀ (0.42P₀ for P₀ = 2000 s) for the kink wave to enter a stable sinusoidal regime. After this duration we saved high-cadence snapshots at the centre of the loop (line x = y = 0). For all further analyses we used the snapshots saved after the transitory phase. The transverse velocity v_x at the loop centre is shown in Fig. 4. As can be seen in this figure, the amplitude of the kink wave decreases as the period increases. For the two longer driver periods (700 and 2000 s), the amplitude of the kink wave is small enough for some perturbations to become visible. They travel at the Alfvén speed, and appear to be triggered by the flows remaining after the relaxation (see Sect. 2.3). These perturbations have amplitudes smaller than 0.2 km s⁻¹, and should thus have no effect on the wave.

Fig. 4. Kink wave transverse velocity (vx) at the loop centre (x = y = 0) as a function of altitude and time. The velocity is shown for four 3D simulations with different driver periods P0, after an initial settling time of 2P0 (for P0 = 200 s, 335 s, and 700 s), or 0.42P0 (for P0 = 2000 s). The dashed black lines represent a propagation at the kink speed (see Eq. (13)), and are independent of the driver period.

3. Results: Cut off and tunnelling of transverse waves

In order to determine whether the kink waves driven in the 3D simulations are experiencing a cut off, we looked at the evolution of the velocity amplitude (Sect. 3.1), as well as the phase speed (Sect. 3.2) as a function of altitude. The analysis of these profiles allows us to establish that the transverse waves are subject to a low-frequency cutoff in the transition region.

3.1. Wave amplitude increases with frequency

In order to compute the velocity amplitude of the kink wave, we fitted the function A_x(z)sin(ω(z)t+ϕ(z)) to the transverse velocity v_x(z, t), at each altitude (z). Here A_x(z) is the velocity amplitude, ω(z) is the kink wave frequency, and ϕ(z) is the phase. The frequency varies by less than 1% with altitude, confirming the theoretical understanding. The velocity amplitude is shown in Fig. 5. In all simulations, the wave amplitude increases with altitude because of the density decreases with altitude and energy conservation. Across the simulations the amplitude at a given altitude increases with the frequency of the wave. This means that kink waves with higher frequencies propagate better from the chromosphere to the corona. This would be consistent with the low-frequency cutoff predicted by analytical models (see Sect. 1).

Fig. 5. Velocity amplitude of kink waves as a function of altitude. The velocity is shown for four different driver periods (P0). The inset has the same axes as the main figure, with a zoom-in on the vertical axis.

3.2. Evanescent waves in the transition region

To determine the altitude at which the waves are cut off, we compared their phase speed v_p(z) to the kink speed of the flux tube c_k(z). The inverse phase speed is equivalent to the phase difference Δϕ(z) between two altitudes separated by Δz: 1/v_p(z) = Δϕ(z)/(ωΔz). The phase difference has been successfully used to determine the cutoff frequency of acoustic and slow-magnetosonic waves in observations (Centeno et al. 2006; Felipe et al. 2010; Krishna Prasad et al. 2017; Felipe et al. 2018), and in simulations (Felipe & Sangeetha 2020). In these articles the authors determine the phase speed for a wide range of frequencies, but at a limited number of altitude positions. In the present study, however, we could only examine four frequencies because of the high computational cost of a simulation. However, we computed the phase difference at all altitudes of the simulation domain. This allows us to determine the altitude at which the wave is cut off.

The phase speed at a given altitude z was computed from the transverse velocity in the cells above and below: v_x(t, z + Δz/2) and v_x(t, z − Δz/2), where Δz = 98 km is the cell size. We apodized these velocity time series with a Hann window, and computed the cross-correlation: C(τ, z) = v_x(t, z + Δz/2)⋆v_x(t, z − Δz/2). We then determined the time delay Δτ(z) by finding the maximum of C(τ, z). To that end, we fitted the function A + Bcos(ω(τ−Δτ)/δ) to C(τ, z) with τ ∈ [ − P₀/4, +P₀/4]. Finally, the phase difference was given by Δϕ(z) = ωΔτ(z), and the inverse phase speed by 1/v_p(z) = Δτ(z)/Δz. The inverse phase speed is shown in Fig. 6 along with the inverse kink speed for the simulated flux tube. The kink speed c_k is calculated using

$$\begin{array}{c}\hfill c_{\mathrm{k}}^2(z)=\frac{{\rho }_i(z)v_{\mathrm{A}i}^2(z)+{\rho }_e(z)v_{\mathrm{A}e}^2(z)}{{\rho }_i(z)+{\rho }_e(z)},\end{array}$$(13)

Fig. 6. Inverse phase speed of the kink wave (1/vp) and inverse kink speed of the flux tube (1/ck) as a function of altitude. The phase speed is given for four different driver periods (P0).

where ρ(z) is the density, $v_{\mathrm{A}}(z)=B(z)/\sqrt{{\mu }_0\rho (z)}$ is the Alfvén speed, B(z) is the magnetic field amplitude, and μ₀ is the magnetic permittivity of vacuum. The indices i and e correspond, respectively, to internal and external quantities relative to the flux tube and are taken at x = 0 and x = 8 Mm.

In simulations with short driver periods, the inverse phase speed is somewhat smaller than the inverse kink speed in the chromosphere and transition region (v_p/c_k ≈ 2 for P₀ = 200 s, and 5 for P₀ = 335 s), and equals the inverse kink speed in the corona. On the other hand, in simulations with longer periods, the inverse phase speeds are much lower than the inverse kink speed below a given altitude. For P₀ = 700 s, 1/v_p is about 250 times smaller than 1/c_k below z = 1 Mm. For P₀ = 2000 s a similar drop occurs below z = 20 Mm.

For a propagating kink wave the inverse phase speed is expected to be equal to the inverse kink speed. Conversely, standing and evanescent (i.e., cutoff) waves have inverse phase speeds lower than the inverse kink speed. Thus, the decreased inverse phase speed for higher periods indicates that the waves are cut off in at least some regions.

To distinguish between the standing and evanescent cases, we also looked at the wave amplitude (Fig. 5). In the absence of vertical stratification, the amplitude of evanescent waves decreases with altitude. However, in a stratified atmosphere (our case), the amplitude increases with altitude because of the density decrease, even for evanescent waves. In Fig. 5, the amplitude of waves with longer periods (for which 1/v_p ≪ 1/c_k) increases less with altitude compared to waves with shorter periods (for which 1/v_p ≲ 1/c_k). We thus conclude that the waves with longer periods are evanescent in parts of the low atmosphere where their inverse phase speed is much lower than the inverse kink speed. This means that these long-period waves are cut off in the transition region.

3.3. Wave tunnelling at higher frequencies

Waves with shorter periods (P₀ = 200 and 335 s) also show signs of cutoff at low altitudes. Below z = 3 Mm, the inverse phase speed 1/v_p is lower than the inverse kink speed 1/c_k (Fig. 6), and the amplitude increase with altitude is smaller for P₀ = 335 s than for P₀ = 200 s (Fig. 5). However, this cutoff is significantly weaker than in the long-period case. This is explained by the fact that the cutoff region (where 1/v_p < 1/c_k) is narrower for short periods (∼1 Mm) than for long periods (∼10 Mm). As a result, short-period waves can tunnel through the cutoff region, and propagate into the corona. Furthermore, the weak attenuation in the cutoff region results (1/v_p ≲ 1/c_k) further reduces the effect of the cutoff.

4. Discussion: Comparison to analytical formulas

In order to compare our simulations to the analytical models, we quantified the cutoff frequency as a function of altitude. We define z_c, the altitude at which c_k/v_p goes above a given threshold t_r. This corresponds to the altitude where the wave leaves the cutoff regime and enters the propagating regime (i.e., the cutoff altitude). We computed z_c for four values of t_r between 0.2 and 0.5. Considering the four simulations with different driver frequencies, ω, we obtained the cutoff altitude as a function of the frequency, z_c(ω). We compare this to the cutoff frequency as a function of altitude, ω_c(z), predicted by the analytical models presented in Sect. 1.

In Fig. 7 we show the cutoff frequency and altitude computed in our simulations, for different values of t_r (black points). In the same figure we show the predictions of the analytical formulas of Spruit (1981, Eq. (1)), Lopin & Nagorny (2017, Eq. (2)), and Snow et al. (2017, Eq. (3)) (coloured lines) computed for the temperature and density profiles used in our simulations. We implement the formula of Lopin & Nagorny (2017) for different values of z₀, defined by the authors as ‘the base of the atmosphere’, with no further details. Because this quantity is not accurately defined, we used four values of z₀ in the range of 24 km (bottom cell of our simulation domain) to 1978 km. This loosely defined parameter broadens the range for the cutoff frequencies predicted by this formula. While the match is rather loose, the cutoff altitude z_c(ω) measured in our simulations matches the overall variation the cutoff frequency ω_c(z) predicted by the Lopin & Nagorny (2017) formula. In particular, the shape of the profiles are in good agreement. On the contrary, the Snow et al. (2017) model correctly predicts the cutoff frequency only in the lower transition region, but fails to do so in the upper transition region and corona. In particular, their model predicts a slower decrease in the cutoff frequency above 20 Mm, while the simulations and the Lopin & Nagorny (2017) formula show a continued decrease. Finally, the Spruit (1981) predictions are off by almost an order of magnitude at all altitudes. Thus, the formula of Lopin & Nagorny (2017) best predicts the cutoff frequency of transverse waves at different altitudes.

Fig. 7. Kink wave cutoff frequency as a function of altitude, from analytical models (see legend, left column) and from our numerical simulations (right column). We show the analytical predictions of Spruit (1981, SP81), Snow et al. (2017, Sn17), and Lopin & Nagorny (2017, LN17) (coloured lines). For the last model, we computed the cutoff frequency for different values of z0, the ‘base of the atmosphere’. We show the cutoff altitude (zc) for the four simulations that we ran with different driver frequencies (black symbols). The cutoff altitudes are computed with different thresholds tr, indicated in the legend and described in the text.

While the broadened transition region in our simulations could affect the altitude-dependence of the cutoff frequency, this should have little impact on the validation of the analytical formulas. These formulas include the atmospheric stratification through altitude-dependent profiles of either the pressure scale height or the Alfvén speed (see Sect. 1). Because they make no hypothesis on these profiles, they should be valid regardless of the atmosphere considered. As such, the agreement with the simulations should not depend on the broadening of the transition region, provided the appropriate profile is fed into the formulas. After validating the Lopin & Nagorny (2017) formula by comparing it to our simulations, it should be applicable to other stratification profiles.

We note that while analytical formulas can predict the kink cutoff frequency, this is not sufficient to know whether a kink wave with a given frequency will propagate into the corona. To that end, the thickness of the cutoff region and the strength of the attenuation have to be taken into account. As shown by our simulations, kink waves with higher frequencies (≥3 mHz) can propagate into the corona by tunnelling through a region where they are cut off (Sect. 3.3). Furthermore, these waves only experience a weak attenuation because their frequency is close to the cutoff frequency. The cutoff frequency does not constitute a clear-cut boundary between oscillatory and non-oscillatory solutions. This was also reported for sound waves by Felipe & Sangeetha (2020). Although the question of whether a solution is oscillating is well defined mathematically, this is not straightforward to translate into a single cutoff frequency (Schmitz & Fleck 1998). For this reason there are several canonical definitions for cutoff frequencies, set within the continuous variation between the oscillating and non-oscillating regimes (see e.g., Schmitz & Fleck 1998 for sound waves in the solar atmosphere). As a result, cutoff frequencies are bound to be mere indications, rather than strong constraints on the physical behaviour of a wave (Chae & Litvinenko 2018).

5. Conclusions

Transverse waves are a candidate mechanism for heating the solar corona. However, several analytical models predicted that they are cut off in the transition region. In order to assess whether transverse waves can indeed heat the corona, it is thus crucial to determine whether they can propagate through the transition region. To that end, we have simulated the propagation of transverse kink waves in an open magnetic flux tube, embedded in an atmosphere extending from the chromosphere to the corona. We find that transverse waves are indeed cut off in the lower solar atmosphere. However, only waves with low frequencies (ν ≲ 2 mHz) are significantly affected. At higher frequencies, the cutoff occurs in a very thin layer (∼1 Mm), and results in a weak attenuation. In this case waves can tunnel through the cutoff layer, experiencing little to no amplitude attenuation. This means that transverse waves with high frequencies are able to transport energy from the chromosphere to the corona, where it can be dissipated and result in heating.

Furthermore, we compared our simulations to several analytical models that predict the cutoff frequency of transverse waves. We conclude that the formula proposed by Lopin & Nagorny (2017) gives the best prediction. While our simulations use a broadened transition, we expect it to have little impact on the validation of analytical formulas. As such, the formula by Lopin & Nagorny (2017) should be able to predict the cutoff frequency for any atmospheric stratification profile. We note that while the cutoff frequency is a good first indicator of whether a wave can propagate into the corona, it alone cannot predict the whole behaviour of the wave. In particular, waves with frequencies just below the cutoff frequency (which should thus be cut off) can still reach the corona thanks to a combination of tunnelling and weak attenuation.

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 724326). GP was supported by a CNES postdoctoral allocation. TVD was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 724326) and the C1 grant TRACEspace of Internal Funds KU Leuven. K.K. recognises support from a postdoctoral mandate from KU Leuven Internal Funds (PDM/2019), from a UK Science and Technology Facilities Council (STFC) grant ST/T000384/1, and from a FWO (Fonds voor Wetenschappelijk Onderzoek – Vlaanderen) postdoctoral fellowship (1273221N). The results received support from the FWO senior research project with number G088021N. Software: Astropy (Astropy Collaboration 2013, 2018).

References

All Figures

	Fig. 1. Sketch of the simulation domain, showing the magnetic flux tube, the location of the kink wave driver (bottom boundary), chromosphere, transition region, corona, and velocity rewrite layer.
In the text

	Fig. 2. Velocity-rewrite coefficient αv, applied to the velocity above 50 Mm so that upper-propagating waves are not reflected back into the domain. αv is shown for different times of the 2D relaxation run. The last profile (t ≥ 31.3 ks) is also applied in the 3D driven simulations.
In the text

	Fig. 3. Temperature (black), density (red), and magnetic field magnitude (blue) profiles inside (r = 0 Mm; solid lines) and outside (r = 8 Mm; dashed lines) the flux tube. (a) After solving the field-aligned hydrostatic equilibrium. (b) After the 2D magnetohydrodynamic relaxation.
In the text

Fig. 4. Kink wave transverse velocity (vx) at the loop centre (x = y = 0) as a function of altitude and time. The velocity is shown for four 3D simulations with different driver periods P0, after an initial settling time of 2P0 (for P0 = 200 s, 335 s, and 700 s), or 0.42P0 (for P0 = 2000 s). The dashed black lines represent a propagation at the kink speed (see Eq. (13)), and are independent of the driver period.

In the text

	Fig. 5. Velocity amplitude of kink waves as a function of altitude. The velocity is shown for four different driver periods (P0). The inset has the same axes as the main figure, with a zoom-in on the vertical axis.
In the text

	Fig. 6. Inverse phase speed of the kink wave (1/vp) and inverse kink speed of the flux tube (1/ck) as a function of altitude. The phase speed is given for four different driver periods (P0).
In the text

Fig. 7. Kink wave cutoff frequency as a function of altitude, from analytical models (see legend, left column) and from our numerical simulations (right column). We show the analytical predictions of Spruit (1981, SP81), Snow et al. (2017, Sn17), and Lopin & Nagorny (2017, LN17) (coloured lines). For the last model, we computed the cutoff frequency for different values of z0, the ‘base of the atmosphere’. We show the cutoff altitude (zc) for the four simulations that we ran with different driver frequencies (black symbols). The cutoff altitudes are computed with different thresholds tr, indicated in the legend and described in the text.

In the text