Cut-off of transverse waves through the solar transition region

Context. Transverse oscillations are ubiquitously observed in the solar corona, both in coronal loops and open magnetic flux tubes. Numerical simulations suggest that their dissipation could heat coronal loops, counterbalancing radiative losses. These models rely on a continuous driver at the footpoint of the loops. However, analytical works predict that transverse waves are subject to a cut-off in the transition region. It is thus unclear whether they can reach the corona, and indeed heat coronal loops. Aims. Our aims are to determine how the cut-off of kink waves affects their propagation into the corona, and to characterize the variation of the cut-off frequency with altitude. Methods. Using 3D magnetohydrodynamic simulations, we modelled the propagation of kink waves in a magnetic flux tube, embedded in a realistic atmosphere with thermal conduction, that starts in the chromosphere and extends into the corona. We drove kink waves at four different frequencies, and determined whether they experienced a cut-off. We then calculated the altitude at which the waves were cut-off, and compared it to the prediction of several analytical models. Results. We show that kink waves indeed experience a cut-off in the transition region, and we identified the analytical model that gives the best predictions. In addition, we show that waves with periods shorter than approximately 500 s can still reach the corona by tunnelling through the transition region, with little to no attenuation of their amplitude. This means that such waves can still propagate from the footpoints of loop, and result in heating in the corona.


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 openfield 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'.
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 transversewave-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. 2019Riedl et al. , 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): 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 where z is the altitude, c k0 is the kink speed at the base of atmosphere (z = z 0 ), H is the pressure scale height, H 0 = H(z 0 ), and is the relative difference between the magnetic field inside (B 0,i ) and outside (B 0,e ) the flux tube at z = z 0 .Finally, an alternative formula was derived by Snow et al. (2017), 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.

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.
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), andMikić et al. (2013).Below the cutoff temperature T c = 2.5 × 10 5 K, the parallel thermal conductivity was set to κ = C 0 T 5/2 c with C 0 = 9 × 10 −12 W m −1 K −7/2 ; above T c , it was set to κ = C 0 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).

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 A105, page 2 of 8 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 where v(t) = v 0 cos (2πt/P 0 ), with v 0 the driver amplitude set to 2 km s −1 .The driver period, P 0 , was set to different values in order to test the cut off of kink waves.Outside the tube the driver imposes where x 0 (t) = v 0 P 0 /(2π) • sin (2πt/P 0 ) 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: 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: 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): else. (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.
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).

Initial conditions: Field-aligned hydrostatic equilibrium
The simulation was initialized with a uniform vertical magnetic field of magnitude B 0 = 42 G.Along the tube, we imposed the temperature profile derived from Aschwanden & Schrijver ( 2002) 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 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), where R = 1 Mm is the tube radius and b = 5 is a dimensionless number setting the width of the inhomogeneous layer between A105, page 3 of 8 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 where ρ ch,int = 3.51 × 10 −8 kg m −3 is the density inside the tube and ρ ch,ext = 1.17 × 10 −8 kg m −3 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.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).

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 0 = 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 −1 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 −1 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.

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 0 : 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 A105, page 4 of 8   13)), and are independent of the driver period.three simulations were run for a duration of 5P 0 .The last simulation was run for 1.75P 0 .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 (0.42P 0 for P 0 = 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 −1 , and should thus have no effect on the wave.

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.

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).

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 A105, page 5 of 8 then determined the time delay ∆τ(z) by finding the maximum of C(τ, z).To that end, we fitted the function A+B cos (ω(τ − ∆τ)/δ) to C(τ, z) with τ ∈ [−P 0 /4, +P 0 /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 where ρ(z) is the density, v A (z) = B(z)/ µ 0 ρ(z) is the Alfvén speed, B(z) is the magnetic field amplitude, and µ 0 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 0 = 200 s, and 5 for P 0 = 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 0 = 700 s, 1/v p is about 250 times smaller than 1/c k below z = 1 Mm.For P 0 = 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.

Wave tunnelling at higher frequencies
Waves with shorter periods (P 0 = 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 0 = 335 s than for P 0 = 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.

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 0 , 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 0 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.
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 A105, page 6 of 8  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 z 0 , the 'base of the atmosphere'.We show the cutoff altitude (z c ) for the four simulations that we ran with different driver frequencies (black symbols).The cutoff altitudes are computed with different thresholds t r , indicated in the legend and described in the text.
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).

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.

Fig. 1 .
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.

Fig. 4 .
Fig. 4. Kink wave transverse velocity (v x ) 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 P 0 , after an initial settling time of 2P 0 (for P 0 = 200 s, 335 s, and 700 s), or 0.42P 0 (for P 0 = 2000 s).The dashed black lines represent a propagation at the kink speed (see Eq. (13)), and are independent of the driver period.

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

Fig. 6 .
Fig.6.Inverse phase speed of the kink wave (1/v p ) and inverse kink speed of the flux tube (1/c k ) as a function of altitude.The phase speed is given for four different driver periods (P 0 ).
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 ofSpruit  (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 z 0 , the 'base of the atmosphere'.We show the cutoff altitude (z c ) for the four simulations that we ran with different driver frequencies (black symbols).The cutoff altitudes are computed with different thresholds t r , indicated in the legend and described in the text.