© The Authors 2025

1. Introduction

Measuring quantities such as the mean magnetic field or the density in the solar atmosphere is a task that requires a combination of observations, theory, and numerical modelling to extract information about the local plasma parameters. Coronal seismology (e.g. Roberts et al. 1984; Nakariakov & Verwichte 2005) uses observed waves in the solar atmosphere, which are understood to be magnetohydrodynamic (MHD) waves in inhomogeneous plasma (Edwin & Roberts 1983), to perform such calculations. Space-based telescopes such as the Transition and Coronal Explorer (TRACE; Handy et al. 1999), the Solar Dynamics Observatory/Atmospheric Imaging Assembly (SDO/AIA; Lemen et al. 2012), and the Solar Orbiter/Extreme Ultraviolet Imager (EUI; Rochus et al. 2020) have revealed an omnipresence of transverse oscillations that are understood to be kink waves (Van Doorsselaere et al. 2008) and which have been extensively used for plasma diagnostics. Observations of kink waves have revealed both propagating motions in open and closed field regions (for a review, see Morton et al. 2023) and standing oscillations in structures such as coronal loops (see Nakariakov et al. 2021, for a review).

For the case of standing kink oscillations in loops, two different regimes have been discovered. The first is the decaying regime first reported in Aschwanden et al. (1999) and Nakariakov et al. (1999). It is characterised by quickly decaying oscillation amplitudes, with the oscillations being excited by external energetic phenomena (e.g. Nechaeva et al. 2019), coronal rain (Verwichte et al. 2017; Verwichte & Kohutova 2017; Kohutova & Verwichte 2017), and nanojets (Sukarmadji & Antolin 2024). In addition, the coexistence of multiple harmonics in decaying oscillations of loops has been demonstrated (e.g. De Moortel & Brady 2007; Duckenfield et al. 2019). The second regime consists of decayless oscillations (e.g. Wang et al. 2012; Tian et al. 2012), which are characterised by a near-constant oscillation amplitude over many periods (Nisticò et al. 2013). These standing oscillations are observed in the solar atmosphere in loops of different scales (e.g. Wang et al. 2012; Nisticò et al. 2013; Anfinogentov et al. 2013, 2015; Anfinogentov & Nakariakov 2019; Zhong et al. 2022a,b, 2023a; Li & Long 2023). However, their interpretation as standing waves is still a matter of debate for most observations in shorter transition region loops and coronal loops (Gao et al. 2022, 2024; Petrova et al. 2023; Shrivastav et al. 2024). It was recently shown in Lim et al. (2024a) that under-sampling due to low cadence can lead to an overestimation of the oscillation periods in shorter loops, giving rise to the observed lack of a correlation between the periods and loop lengths for shorter loops.

Simulations of decayless oscillations have hinted at an energy content sufficient to counterbalance the radiative losses in the quiet Sun corona (e.g. Karampelas et al. 2019; Shi et al. 2021; De Moortel & Howson 2022). Recent meta-analysis studies of observations (e.g. Lim et al. 2023, 2024b) have found that the high-frequency regime of decayless oscillations can heat the quiet Sun corona, but cannot support the wave heating in active regions unless unresolved oscillations with frequencies up to 0.17 Hz are considered. Decayless standing waves are also reported to exhibit linear polarisation (see e.g. Zhong et al. 2023b, for a combined observation with SDO/AIA and the High Resolution Imager HRI_EUV of Solar Orbiter/EUI), while the coexistence of multiple harmonics (first and second) was also reported in Duckenfield et al. (2018). Equally important is that decayless oscillations still have an unidentified excitation mechanism, with different hypotheses proposing it to be broadband footpoint drivers (Afanasyev et al. 2020; Ruderman & Petrukhin 2021; Ruderman et al. 2021; Howson & De Moortel 2023; Karampelas & Van Doorsselaere 2024), p-modes (Skirvin et al. 2023; Gao et al. 2023), plasma flows along the loop (Kohutova & Verwichte 2018), and external flows (Nakariakov et al. 2009, 2016; Karampelas & Van Doorsselaere 2020, 2021). Additionally, Antolin et al. (2016) suggest that apparent observations of decayless oscillations are actually observations of decaying oscillations muddled by line-of-sight effects.

As stated earlier, transverse kink oscillations have been used extensively as tools in coronal seismology. For example, by measuring the period of the fundamental kink mode and the loop length, we can calculate the kink speed of the loop and then estimate the average Alfvén speed and the average magnetic field using common seismology techniques (e.g. Edwin & Roberts 1983; Roberts et al. 1984; Aschwanden & Schrijver 2011; Anfinogentov & Nakariakov 2019; Zhong et al. 2023a; Gao et al. 2024). Longitudinal density stratification and the expansion of the magnetic field can also lead to the modification of the spatial profile and the respective frequency of the fundamental kink mode and its overtones (e.g. Andries et al. 2005a, see also Andries et al. 2009 for a review). For example, Verth et al. (2007) showed that the antinode of the second harmonic is shifted towards the loop footpoints as the density scale height decreases, for a flux tube with a straight magnetic field (see also Erdélyi & Verth 2007). If multiple kink mode harmonics are present in the same oscillating loop, then their ratios can be used to determine the density scale height, assuming again that there is no variation in the magnetic field (Andries et al. 2005b; Safari et al. 2007). The effects of the lower solar atmosphere and the transition region have started to be incorporated in 3D simulations of open flux tubes and coronal loops in studies of wave propagation (Pelouze et al. 2023), excitation of decayless oscillations (e.g. Gao et al. 2023; Karampelas & Van Doorsselaere 2024), and wave energy dissipation (e.g. Guo et al. 2023; Karampelas et al. 2024). Howson & Breu (2023) have numerically shown from the 1D Sturm-Liouville problem describing an oscillating loop (Dymova & Ruderman 2005) that the existence of the transition region and the chromosphere will greatly modify the profile of the harmonics of standing Alfvén and kink waves. It was also suggested that this deformation could be potentially misleading when trying to identify the oscillation modes in observations of loop systems.

We studied the effects of the lower solar atmosphere on the spatial profiles and frequencies of standing kink modes of coronal loops driven by noisy footpoint motions. Our aim is to identify the individual modes in our simulation results and create synthetic emission maps that can be compared with observations. We also applied the techniques of coronal seismology to understand the importance of the lower atmosphere when using kink oscillations as plasma diagnostic tools for quantities such as the kink speed and mean magnetic field. In Sect. 2 we describe the initial conditions of our 3D coronal loop model and how we excite the oscillations of interest. Section 3 describes the results of the simulations, the creation and analysis of our synthetic observations, the equivalent 1D eigenvalue problem, and the use of seismological tools. Finally, Sect. 4 includes a discussion of the most important points of this study.

2. Numerical setup

This numerical study is based on the analysis of our setup that consists of a stratified 3D straight magnetic flux tube, modelling a coronal loop with chromospheric footpoints. The numerical setup is similar to the ones in Karampelas & Van Doorsselaere (2024) and Karampelas et al. (2024). Following our past work, we initialised a 2D model in hydrostatic equilibrium along the vertical z direction, using sinusoidal gravity (g_z(z) = 274 cos(π z/200) in m s⁻²), to emulate the gravity variation along a semi-circular loop. We ran a 2.5D resistive MHD simulation to let our model reach a quasi-equilibrium state and we then interpolated into a 3D domain.

The initial profiles for the magnetic field (B), density (ρ), and temperature (T) are

$$\begin{array}{cc}& \{B_r,B_{\theta },B_z\}=\{0,0,30\}\mathrm{G}\hfill \end{array}$$(1)

$$\begin{array}{cc}& T(r,z)=T_{\mathrm{Ch}}+(T_{\mathrm{C}}(r)-T_{\mathrm{Ch}}){(1-{[(L-z)/(L-{\mathrm{\Delta }}_{\mathrm{Ch}})]}^2)}^{0.3},\hfill \end{array}$$(2)

$$\begin{array}{cc}& T_{\mathrm{C}}(r)=T_{\mathrm{C},\mathrm{e}}+(T_{\mathrm{C},\mathrm{i}}-T_{\mathrm{C},\mathrm{e}})\zeta (r),\hfill \end{array}$$(3)

$$\begin{array}{cc}& {\rho }_{\mathrm{Ch}}={\rho }_{\mathrm{Ch},\mathrm{e}}+({\rho }_{\mathrm{Ch},\mathrm{i}}-{\rho }_{\mathrm{Ch},\mathrm{e}})\zeta (r),\hfill \end{array}$$(4)

$$\begin{array}{cc}& \zeta (r)=0.5[1-tanh\left(([r/R]-1)20\right)],\hfill \end{array}$$(5)

where the subscripts C and Ch correspond to the coronal (z = 100 Mm) and chromospheric (z = 0, 200) values. The temperature in the transition region and the corona is calculated for height z ∈ [5, 195] Mm. The loop length is L = 200 Mm, while Δ_Ch = 5 Mm is the width of our chromosphere. For z < 5 Mm and z > 195 Mm, we considered a uniform chromospheric temperature of T_Ch = 0.02 MK. The temperature at the apex (T_C(r)) is T_C, i = 1 MK and T_C, e = 1.5 MK inside and outside the flux tube, respectively. The density is ρ_Ch, i = 3.51 × 10⁻⁸ kg m⁻³ inside the footpoints and ρ_Ch, e = 1.17 × 10⁻⁸ kg m⁻³ outside the footpoints. The function ζ(r) calculates the loop cross-section, with R = 1 Mm being the loop minor radius.

For the 2.5D domain (r ∈ [0, 8] Mm and z ∈ [0, 200] Mm), we took 200 × 2048 uniformly spaced grid points in cylindrical coordinates. We considered axisymmetry at r = 0 and open boundary conditions at r = 8 Mm. At the z-axis boundaries, we used zero-gradient conditions for the magnetic field components, anti-symmetric conditions for the velocity components, symmetric conditions for the density, and a constant temperature. We then let the setup reach a semi-equilibrium state for a total time of t = 3890 s, reducing the velocity values per iteration as v_{{r, θ, z}} = v_{{r, θ, z}}/1.0001 for t ∈ [778, 3112] s, similarly to Pelouze et al. (2023) and Guo et al. (2023). Figure 1 shows the post relaxation 2D temperature profile as well as 1D profiles for the density and magnetic field inside (i,  r = 0) and outside (e,  r = 8 Mm) the flux tube. After the relaxation, the magnetic field has been restructured through advection and dissipation and is no longer purely vertical, but also has a radial component (B_r) that has its a maximum value near the loop footpoints. This maximum value of the B_r is an order of magnitude lower than the vertical component (B_z) there. The vertical component gets its minimum value near r = 0 at the footpoint, and its maximum value at r ∼ 1.5 Mm near the footpoint, as was shown for a similar setup in Karampelas et al. (2024). At larger radii, as well as towards the transition region and the corona, the radial component is almost zero, and the magnetic field is near-vertical there. Similarly to Karampelas et al. (2024), the loop minor radius does not change post relaxation, but is approximately R = 1 Mm, with the exception of the footpoints (z = 0, 200 Mm) where it has increased to R = 1.2 Mm. Finally, in Fig. 1 we also plot the calculated kink speed (C_Kink) along the z-axis, with the coronal part (z ∈ [7, 193] Mm) highlighted.

Fig. 1. Post 2.5D MHD relaxation profiles for the temperature, density, Bz magnetic field, and kink speed of our stratified flux tube. Left: 2D temperature profile for a section of our model up to r = 3 Mm. Top right: 1D profiles along the z direction of the density (black lines) and Bz magnetic field (red lines) at r = 0 (solid lines) and at r = 8 Mm (dashed lines). Bottom right: Profile of the kink speed along the z-axis, with the coronal part (z ∈ [7, 193] Mm) highlighted (dashed orange line).

We interpolated the post 2.5D relaxation slice into a Cartesian grid of size x ∈ [ − 4, 4] Mm, y ∈ [0, 4] Mm and z ∈ [0, 200] Mm, with δx = δy = 0.040 Mm everywhere, δz = 0.098 Mm for z ≤ 10 Mm and z ≥ 190 Mm and δz = 0.8 Mm for 24 Mm ≤ z ≥ 176 Mm. We considered stretched grids with 40 cells each for 10 Mm ≤ z ≤ 24 Mm and 176 Mm ≤ z ≤ 190 Mm.

For the 3D cases, we took open boundary conditions at x = −4, 4 Mm and y = 4 Mm and reflective boundary conditions at the y = 0, which allowed us to simulate only half of the loop. Along the z-axis, we considered symmetric boundary conditions for both density and pressure, as well as zero-gradient conditions for the magnetic field and anti-symmetric boundary conditions for the v_x, v_y and v_z velocity components, apart from the case where a velocity driver is employed for v_x and v_y at z = 0, for t ≤ 202 s.

Two different cases are considered in the numerical part of this study. For the first case, we used a model similar to that in Karampelas & Van Doorsselaere (2024), where at z = 0 for t > 202 s, we applied a broadband, linearly polarised velocity driver:

$$\begin{array}{c}\hfill \}{v}_x,v_y\{=\}V(t)\zeta (r_d,t),0\{,\end{array}$$(6)

with the function ζ(r) from Eq. (5) calculating the driver cross-section and $r_d=\sqrt{{(x-x_0(t))}^2+y^2}$ being the location of the driver centred at x₀(t) with radius R = 2.5 Mm. The velocity signal V(t) follows a red noise spectrum (∼f^−1.66, where f is the frequency), constructed through the colorednoise 2.1.0 python package (see also Afanasyev et al. 2020), but with reduced power at the lower frequencies (< 1.5 mHz) after using a high-pass filter, similar to that used in Karampelas et al. (2024). The left panels of Fig. 2 show the velocity driver profile and its respective spectrum. The RMS velocity (V_RMS) of the signal (∼0.94 km s⁻¹) is comparable to the RMS velocities of the horizontal motions of magnetic bright points (∼1.32 km s⁻¹) derived from the Swedish 1 m Solar Telescope (SST) and Hinode observations (Chitta et al. 2012), which could be treated as loop footpoints motions.

Fig. 2. Profiles of the velocity driver and pulses used in the three 3D simulations. Left panels: Temporal profile of the vx driver at z = 0 for the first simulation (top) and its corresponding power spectrum (bottom). Right panel: Two initial velocity conditions used in the two simulations without the footpoint driver.

For the second case, we ran two simulations in which we impulsively perturbed the flux tube with an initial condition in the v_x velocity. We used a sinusoidal pulse, which is symmetric with respect to the apex and an off-centre Gaussian pulse, which was asymmetric with the respect to the apex:

$$\begin{array}{cc}& v_x=15\zeta (r)sin(\pi z/L),\hfill \end{array}$$(7)

$$\begin{array}{cc}& v_x=45\zeta (r)exp{((z-80)/30)}^2.\hfill \end{array}$$(8)

Both velocity pulses are given here in km s⁻¹, and z is in Mm. The radius of the velocity pulse cross-section, defined by ζ(r) from Eq. (5), is R = 1.2 Mm to ensure that the entire flux tube cross-section is perturbed. We note that a different amplitude was used so that the initial displacement of the flux tube would be of the same order. In the right panel of Fig. 2, we see the velocity initial condition used to perturb our flux tube in the two simulations of that setup.

We solved the compressible MHD equations for a hydrogen plasma using the PLUTO code (Mignone et al. 2007). We enforced the ∇ ⋅ B = 0 condition through a hyperbolic divergence cleaning method. We used anisotropic thermal conduction, with the thermal conduction coefficient (in J s⁻¹ K⁻¹ m⁻¹) from the Spitzer conductivity (Orlando et al. 2008), and the method for artificially broadening the transition region during the simulation, described in Linker et al. (2001), Lionello et al. (2009), and Mikić et al. (2013):

$$\begin{array}{cc}& {\kappa }_{\Vert }=5.6\times {10}^{-12}\{\begin{array}{cc}{T}^{5/2},\hfill & \mathrm{if}T\ge 2.5\times {10}^5\mathrm{K},\hfill \\ {(2.5\times {10}^5)}^{5/2},\hfill & \mathrm{if}T\le 2.5\times {10}^5\mathrm{K},\hfill \end{array}\hfill \end{array}$$(9)

$$\begin{array}{cc}& {\kappa }_{\perp }=3.3\times {10}^{-21}{n}_{\mathrm{H}}^2/(\sqrt{T}{B}^2).\hfill \end{array}$$(10)

Our setup has a different value for the parallel thermal conduction coefficients than the ones used in Karampelas & Van Doorsselaere (2024) and Karampelas et al. (2024) in order to match the values used in Guo et al. (2023) and Gao et al. (2023). In the 2.5D setup, explicit magnetic diffusivity (η = R_m⁻¹ = 10⁻⁴) was used to improve the code stability, while none was used in the 3D case. Our numerical scheme also introduces effective numerical diffusivity, with a corresponding effective magnetic Reynolds number of R_m, eff = 10⁴ − 10⁵.

3. Results

We ran three 3D MHD simulations, in which our stratified loop performed transverse oscillations. In the first case, our loop is perturbed by a footpoint broadband and a linearly polarised driver, while in the second and third case we used a sinusoidal pulse symmetric along the loop velocity and an off-centre Gaussian velocity pulse asymmetric along the loop, respectively. We limited our 3D study mostly to the coronal part of our system, and we treated the chromosphere as a mass reservoir of near-constant temperature. However, we focused on the importance of the large gradients in the transition region for our understanding the dynamics of the standing transverse oscillation modes.

3.1. 3D loop simulations and forward modelling

Our reference step is to revisit and expand upon the result in our past work, as this is the basis for the upcoming analysis. In Fig. 3 we plot the time-distance maps of the loop displacement for the three simulations, alongside their respective spectra derived from Fourier analysis. The loop displacement was calculated by tracking the position of the loop centre of mass at each xy-plane through a weighted surface averaging, with (ρ(z)−ρ_e(z))² as the weight. Subtracting the external background density (ρ_e) and squaring reduced the effects of the background plasma in the calculation of the centre of mass. The left panels of Fig. 3 show the displacement for a small selection of coronal heights, while we used the entire array when calculating the spectra in the right panels. We also note that while the loop displacement is along the x-direction, the left panels depict the magnified signal projected on the z-axis for better visualisation.

Fig. 3. Centre of mass displacement (left panels) and corresponding power spectra density profiles (right panels) for the three simulations of our 3D loop. The time-distance maps of the displacement along the x-direction show the magnified (×20) signal projected along the z-axis for visualisation purposes. Only the coronal part of the loop is depicted. Shown are the simulations with the driver (top), with the symmetric initial pulse (middle), and with the asymmetric pulse (bottom).

As was shown in Karampelas & Van Doorsselaere (2024) for a red noise driver and in Karampelas et al. (2024) for a similar driver with less power in the lower frequencies, such footpoint drivers excite persisting transverse oscillations with a non-decaying, fluctuating amplitude. This can be seen on the top left panel of Fig. 3. Examining the spectra in the top right panel, we again identify a frequency band at ∼5 mHz that spatially resembles the fundamental standing kink mode of the coronal part of the loop (and was understood as such in Karampelas & Van Doorsselaere 2024). We also identify many higher harmonics, as well as wide frequency band at ∼2 − 3.5 mHz, the nature of which was not fully understood in our past work (see Karampelas & Van Doorsselaere 2024) and which was called descriptively as ‘half harmonic’.

The middle panels of Fig. 3 show the time-distance map and spectra for the loop excited by symmetric sinusoidal velocity pulse. When considering the same sinusoidal pulse in Karampelas & Van Doorsselaere (2024), wavelet spectra for the oscillation signal at the apex revealed a strong signal at ∼5 mHz and a very weak one at ∼2.5 mHz, which was one of the primary reasons why the harmonic at 5 mHz was understood as the fundamental mode, the other reason being its spatial distribution along the loop for the driven oscillations. In this study, we repeated this experiment, and, as shown on the middle panels of Fig. 3, we calculated the spectra along the loop length and again see a very strong harmonic at ∼5 mHz that resembles the fundamental mode of the oscillating coronal part of the loop. Alongside it, we also detect a weaker signal at ∼2.5 mHz matching the previously mentioned half harmonic.

In the bottom panels of Fig. 3 we plot the same quantities as before but for the loop perturbed by an asymmetric, off-centre Gaussian pulse. This type of pulse was chosen to model in a more realistic way a transverse oscillation in a loop excited by a random pulse rather than aiming to excite a specific harmonic. As stated earlier, the amplitude of the velocity pulse was chosen so that the resulting oscillation has a similar amplitude to the case where the symmetric pulse was used. From the spectra, we again see that the main signal is coming from the mode around 5 mHz, with a much weaker signal than before seen at 2.5 mHz. Additionally, we also see higher harmonics at ∼8, 12 and 16 mHz, whose frequencies and spatial profiles make them appear as the second, third, and fourth harmonics with respect to the mode at ∼5 mHz, for a gravitationally stratified coronal loop (e.g. Andries et al. 2005b; Safari et al. 2007).

To compare our results with observations of decayless oscillations, we performed forward modelling to create synthetic observations for our 3D simulations. We used the FoMo code (Van Doorsselaere et al. 2016) to calculate the extreme ultraviolet (EUV) emission for the Fe IX 171 Å line from the optically thin coronal part of our oscillating loops. The 171 Å line is ideally suited to the temperature range of the coronal part of our models, since it has its maximum formation temperature of log T = 5.93. In our analysis we targeted the spatial resolution of current and upcoming instruments, specifically the SDO/AIA, the Solar Orbiter/EUI HRI_EUV and the Spectograph (SG) on board the Multi-slit Solar Explorer (MUSE; De Pontieu et al. 2022; Cheung et al. 2022). While SDO/AIA covers the 171 Å line and MUSE/SG is designed to do so as well, Solar Orbiter/HRI_EUV has a channel that covers instead the 174 Å line. This line, however, covers a similar temperature range as the 171 Å line, so we can use the latter as a proxy in our analysis, similarly to Petrova et al. (2023). We degraded the original spatial resolution of the synthetic images to $0\stackrel{″}{.}6$ for SDO/AIA, to $0\stackrel{″}{.}4$ for MUSE/SG, and to 200 × 200 km² for Solar Orbiter/HRI_EUV, assuming that Solar Orbiter is located at a distance of 0.52 AU from the Sun. That was the distance when the observations of high-frequency decayless oscillations discussed by Petrova et al. (2023) took place with Solar Orbiter. We note here that the time resolution of our results (∼15.52 s) was not changed. Therefore, caution is needed when directly comparing to Solar Orbiter/HRI_EUV data (4 s cadence) and to a lesser extent to SDO/AIA and MUSE/SG data (∼12 s).

The time-distance maps of the 171 Å line emission for the three simulations of the oscillating loops at the apex, for the resolution of the three instruments (SDO/AIA, Solar Orbiter/HRI_EUV, and MUSE/SG), are shown in Fig. 4. We note here that we synthesised observations at each height along the three loops but here only show the results at the apex. The top panels show the results for the loop perturbed by the footpoint driver, while the middle and bottom panels show the results for the impulsive oscillations excited by the symmetric Gaussian velocity pulse and the off-centre exponential, asymmetric velocity pulse. We used the Automatic Northumbria University Wave Tracking code (Auto-NUWT; Morton et al. 2013; Weberg et al. 2018) to track the transverse motions of our loops. The code tracks the loop position by measuring the intensity gradients. The position of each loop along the x-axis is overplotted on the panels of Fig. 4. The different bin size used for resampling the data at each resolution can affect the readings of the intensity gradient, leading to slight differences in the measured oscillation signal from each targeted instrument. While the simulation with the broadband driver clearly shows a non-decaying, fluctuating amplitude, the oscillations excited by the two different pulses show no qualitative difference regarding the amplitude and the oscillation period. This implies that the two different pulses are expected to have similar power spectra density profiles, as was shown for the displacement of the centre of mass in Fig. 3.

Fig. 4. Time-distance maps of the synthetic emission in the 171 Å line, for the three cases of oscillating loops, at the apex. Left column: SDO/AIA resolution. Middle column: Solar Orbiter/HRIEUV resolution. Right column: MUSE/SG resolution. Top panels: Loop with the broadband driver. Middle panels: Loop with the sinusoidal, symmetric pulse. Bottom panels: Asymmetric, off-centre Gaussian pulse. Overplotted in each panel are the oscillating signal of the loop centre (black asterisks) and a fitted sinusoidal function (solid red line) with the four prominent frequencies for each signal, as detected by Auto-NUWT.

Since our focus in this study is the spectra of the excited transverse oscillations, we wanted to extract the relevant information from the synthetic emission and compare it with the results in Fig. 3. Using the location of the loop centre, we then used Auto-NUWT to identify the frequencies detected for the transverse oscillations through a Fourier analysis with a 95% significance level. In each panel of Fig. 4 we show the fit of the following sum of sine functions:

$$\begin{array}{c}\hfill x(t)={\displaystyle \sum _{i=1}^4}{c}_{i}sin(\frac{2\pi t}{P_i}+\varphi )+x_0,\end{array}$$(11)

where the coefficient c_i is the displacement amplitude for each respective period (P_i) and x₀ is just the position of the loop at t = 0. As we see in Fig. 4, this fit seems to not precisely track the loop centre (black asterisks) for the model with the broadband footpoint driver (top panels). This is because only the four frequencies with the highest amplitudes (c_i) are used for the fit in Eq. (11), while additional frequencies with less power are omitted. The models with the initial velocity perturbation (middle and bottom panels) show a better agreement between the fit and the detected loop centre, because they only exhibit one dominant frequency, as was also shown in Fig. 3. Here, we also note that damping coefficients are not included in our analysis, as shown in Eq. (11).

The distribution and relative strength of the detected frequencies along each loop, for each instrument resolution is shown in Fig. 5. The panels here follow the layout of Fig. 4. Focusing on first in the middle and bottom panels, we see that the synthetic images for the oscillations excited by the symmetric and asymmetric velocity pulses respectively exhibit very similar spectra, with one dominant frequency at ∼5 mHz and with very small differences between the different instrument resolutions. This is the same frequency band detected for these two cases in the spectra in Fig. 3. Some additional frequencies with very low power are detected near footpoints, which can be attributed to noise near the bottom of the corona – top of the transition region in our setups. This is also at the tail end of the temperature range of the 171 Å line. A few points are highlighted near the apex at about 0.5 mHz, for the asymmetric pulse case, as seen in the AIA resolution. However, this can also be attributed to noise, since no respective signal has been detected in the spectra in Fig. 3. We should highlight here that no other harmonic or mode in general can be detected in the middle and bottom panels, apart from the one at ∼5 mHz.

Fig. 5. Distribution of oscillation frequencies, detected with Auto-NUWT at each height, z, from the synthesised 171 Å line emissions. The colour bars correspond to the normalised amplitude of the oscillation for each detected frequency. As in Fig. 4, the left, middle, and right columns correspond to SDO/AIA, Solar Orbiter/HRIEUV, and MUSE/SG resolution, while the top, middle, and bottom panels correspond to the simulations with the broadband driver, symmetric pulse, and asymmetric pulse, respectively.

For the simulation with the footpoint driver, the detected frequencies from the synthetic observations are shown in the top panels of Fig. 5. Unlike the other two cases, multiple frequency bands are detected along the loop, including the one at ∼5 mHz. Traces of higher frequency modes are present here, but the main signal comes from the frequency range of the previously mentioned half harmonic. We detect distinct frequency bands along z at ∼2.5, 3.5 and 4 mHz. Increasing the resolutions reduces the number of frequencies with power above the 95% significance level, again showing the effects that the resolution has on the detection algorithms. It is also interesting to note that while the frequency band at 5 mHz has a power distribution along z with the maximum at the apex and minima near the footpoints, the other frequency bands at lower frequencies exhibit additional strong power near the driven footpoint at z = 0. This implies that they are at least partially associated with the driving frequencies of the footpoint driver. A similar phenomenon was observed in simulations of short loops in Gao et al. (2023), where the frequency of an inclined p-mode driver also manifested in wavelet spectra of the transverse oscillations alongside the fundamental mode and it harmonics. Also, in Karampelas & Van Doorsselaere (2024), the addition of a different driver to the second footpoint did not modify the spectra, but seemed to reinforce these frequencies near the second footpoint.

3.2. 1D wave equation eigenvalue problem

To better understand the nature of the observed frequency bands in the spectra of our oscillating loops and conclude which are associated with the oscillating system itself and which with the broadband nature of the footpoint driver, we needed to calculate the natural frequencies of the standing transverse modes of our system. An analytical approach is not feasible for the scope of this study, due to the non-uniform nature of the kink speed of our system and its strong gradients along the loop. Instead, we employed the approach used in Dymova & Ruderman (2005), and the approach used in Howson & Breu (2023) for the case of Alfvén waves. For the thin tube approximation, the modes and frequencies of non-axisymmetric, transverse oscillations of loops were calculated using the Sturm-Liouville problem

$$\begin{array}{c}\hfill \frac{{\partial }^{2}u}{\partial z^2}+\frac{{\omega }^2}{v_{\mathrm{ph}}^2(z)}u=0,u(z=0,L)=0.\end{array}$$(12)

The phase speed along the loop v_ph is taken equal to the kink speed C_Kink that we have calculated and show in Fig. 1. We numerically solved Eq. (12) to find the eigenfunctions u that correspond to the standing kink mode harmonics and the eigenvalues ω² = (f/2π)², which give the respective frequencies f for each harmonic. We note here that the thin tube approximation used to derive Eq. (12) is only valid for the coronal part of our loop model. For the part of the loop below the corona, the scale height and thus also the vertical scale are expected to drop significantly from their coronal values, reducing the validity of the thin tube approximation. Therefore, the results of Eq. (12) are taken as an approximation of the actual transverse modes of the system.

Figure 6 shows the normalised transverse oscillation modes and their frequencies for our loop model. The left panel shows the solutions for our entire loop, or when the entire C_Kink profile from Fig. 1 was used, and the right panel shows the solutions only for the coronal part of the loop (i.e. when only the coronal part of the C_Kink profile was used, shown with the dashed orange line on the bottom right panel of Fig. 1). Focusing on the left panel, we notice that the normalised eigenfunctions are deformed, with the nodes of the higher-order harmonics being pushed towards the footpoints of our loop. This effect has been been reported in past studies of straight flux tubes and was deemed to be caused by the density stratification along the tube axis for models with a straight magnetic field along height (e.g. Erdélyi & Verth 2007; Verth et al. 2007), that is, for loops with a non-expanding magnetic field. Similarly to Howson & De Moortel (2022), we see that the existence of the transition region leads to an extreme manifestation of this effect. Considering a coronal loop rooted at the base of the corona, we get eigenfunctions with less extreme deformation with respect to the case of no density stratification, as we see in the right panel of Fig. 6. From the mode profiles for the two cases considered here, we see that the first two kink oscillation harmonics, which are the fundamental mode and the first overtone, are significantly deformed, with both of them exhibiting large amplitudes near the base of the corona (z ∼ 7 and 193 Mm). As expected, such a qualitative behaviour is absent from the coronal loop model, for which all harmonics exhibit nodes at the base of the corona. Additionally, modes like the third, sixth, ninth, and twelfth harmonics of the full loop have nodes located below, but near the base of the corona and show similar spatial profiles to the first, second, third, and fourth harmonics for the coronal loop.

Fig. 6. Normalised eigenfunctions of the 1D oscillating string system and its corresponding frequencies. Left: Solution for the entire system. Right: Solution when only the coronal part of the system is used. The cyan shaded region highlights the coronal part in the left panel.

From the same two panels of Fig. 6, we also get the frequencies of each oscillations mode, shown in the respective legends. Focusing on the case of the full loop (left panel), we see that the fundamental mode and the first overtone have frequencies f₁ = 2.69 mHz and f₂ = 2.76 mHz, respectively. These frequencies fall within the range of the half harmonic mode that is detected in the spectra shown in Figs. 3 and 5. These two modes exhibit very similar frequencies, due to the extreme density gradient in the transition region (Andries et al. 2005a,b). The third harmonic has an eigenfrequency equal to f₃ = 4.75 mHz, which is shown to be very similar to that of the fundamental standing kink mode of the coronal part of our loop (f₁ = 5, 09 mHz), as shown in the right panel of Fig. 6. It is also very interesting to note that the eigenfrequencies included higher harmonics for the full loop (f₆ = 8.34 mHz, f₉ = 12.20 mHz and f₁₂ = 16.03 mHz) are again similar to those of the second (f₂ = 8.93 mHz), third (f₃ = 13.07 mHz) and fourth harmonics (f₄ = 17.22 mHz) of the coronal loop, in the same way as their respective spatial profiles.

The results of the 1D analysis can thus explain the findings in our past work (Karampelas & Van Doorsselaere 2024; Karampelas et al. 2024) as well as in Figs. 3 and 5. What we have descriptively been referring to as half harmonic is in fact the combined result of the loop expected eigenmodes with the added effects of the continuous broadband driving, when the footpoint driver is employed. The thin frequency band around 2.5 mHz detected for the oscillating loops without the footpoint driver correspond to the fundamental and first overtone of the system. The additional frequencies between 2.5 and 3.5 mHz detected for the loop with the footpoint driver are associated with the spectrum of the broadband driver itself. This is further reinforced by the results in simulations of short loops (Gao et al. 2023), where the frequency of an inclined p-mode driver also manifested in wavelet spectra of the transverse oscillations alongside the fundamental mode and the second and third harmonics. Additionally, in simulations of a stratified but purely coronal loop perturbed with monochromatic footpoint drivers of various frequencies, the 3D flux tube would exhibit a resonant response to the monochromatic drivers, even for small deviations from the resonant frequency (Afanasyev et al. 2019). Although not focused on in that study, this could be the result of the non-linear evolution of the loop profile due to the development of the instabilities.

The frequency band at ∼5 mHz, detected in Fig. 3 and in Fig. 6 is the third harmonic of the full-loop system with frequency f₃ = 4.75 mHz. In past studies (Karampelas & Van Doorsselaere 2024; Karampelas et al. 2024) this mode was misidentified as the fundamental mode of the system, due to its spatial profile and as we see from the right panel of Fig. 6, this mode closely resembles the fundamental mode of the coronal part of our loop model. This mode gives a strong signal in the spectra from the synthetic observations of the driven model and is the dominant mode in the spectra for the oscillating loops without the footpoint driver. The prominence of this mode could therefore lead to it being misidentified with the fundamental mode, which in turn can have implications for coronal seismology, as was discussed in Howson & Breu (2023), as well as in Karampelas & Van Doorsselaere (2024) for a a system driven by a broadband footpoint driver. Our 1D analysis also helps us identify the nature of the higher harmonics detected in the panels of Fig. 3, the frequencies and spatial profile of which seems to closely match that of the harmonics expected from the 1D eigenvalue system.

3.3. Using the 1D eigenfunctions for seismology

To explore how the effects of the density stratification can influence the results derived by seismological techniques, we employed various tools to study the results from our 1D analysis. In a recent paper by Chen et al. (2023), spatial seismology methods were used to infer the density scale and magnetic field expansion factor. This was done by fitting observations of transverse oscillations of loops with the expected profiles of transverse standing oscillations modes, using a χ² goodness of fit test. Here we employed a similar approach to calculate the density scale height and magnetic field for various cases. For the role of observations, we considered randomly selected data points (η_i, measurements) from (a) the first, second, third and sixth harmonics of the full-loop system, as well as (b) from the first and second harmonic of the coronal loop system. For both cases, we used the results derived from the 1D eigenvalue problem of Eq. (12). As predictions ψ, we used the first and second harmonics from the eigenvalue problem of Eq. (12), calculated for a phase speed (v_ph(z)) equal to the kink speed of a coronal loop model with an exponential density profile, in sinusoidal gravity like in our 3D simulations, with a straight, uniform magnetic field, B:

$$\begin{array}{cc}& \mathit{B}=\}0,0,B_z\{=\}0,0,B\{\hfill \end{array}$$(13)

$$\begin{array}{cc}& {\rho }_{\mathit{ie}}={\rho }_{0,ie}exp(-Lsin(\pi z/L)/(\pi H))\hfill \end{array}$$(14)

$$\begin{array}{cc}& ϵ=({\rho }_{0,e}/{\rho }_{0,i})=0.1\hfill \end{array}$$(15)

$$\begin{array}{cc}& v_{\mathrm{ph}}^2(z)=\frac{{\rho }_i{V}_{\mathrm{A},i}^2+{\rho }_e{V}_{\mathrm{A},e}^2}{{\rho }_i+{\rho }_e}=\frac{2}{1+ϵ}\frac{B_z^2}{{\mu }_0{\rho }_i}·\hfill \end{array}$$(16)

Here, the subscripts ‘i’ and ‘e’ indicate the values at the centre of the loop cross-section and at a radial distance of r = 8 Mm from the loop, respectively, V_A is the Alfvén speed, L = 186 Mm is the length of the fitted coronal loop model and ρ_0, i = 20.63 × 10⁻¹² kg m⁻³ is the density value inside the loop at the bottom of the corona, for the considered model.

To fit the model, we employ the same technique as described in Chen et al. (2023). We first coded an iterative stochastic global search optimisation algorithm, known as a simulated annealing algorithm. Starting from an initial guess for the density scale height (H = 40 Mm) and the uniform magnetic field (B_z = 30 G), we performed at each iteration a χ² goodness of fit test:

$$\begin{array}{c}\hfill {\chi }^2(H,c)={\displaystyle \sum _{i=1}^N}\frac{{(\psi (H,z_i)-c{\eta }_i)}^2}{\psi (H,z_i)}·\end{array}$$(17)

The coefficient c is used to normalise the observations η_i and generally has dimensions the inverse of the dimensions of the observations. As an initial guess for this coefficient we take $c={\eta }_{i,\max }^{-1}$. The total number of measurements are equal to N. At each iteration, we solve Eq. (12) for randomly adjusted values of H and c, and we use the new prediction ψ to minimise the χ² statistic. We then compare the final value of the χ² against the critical value for the χ² test (χ_crit²), at 95% significance level (or α = 0.05) with N − 1 degrees of freedom. If χ² is above the critical value, then we reject the null hypothesis that the observations are part of predicted values. If the critical value is above the χ² statistic, then we fail to reject the hypothesis, and thus we consider the model to be a good fit for the observations. Once we reach a good fit, we then vary the value of the uniform magnetic field B, for the predicted eigenfrequency of the model to match the frequency of the observations.

The results for our set of successful χ² fitting tests are shown in the panels of Fig. 7. For our tests, we have a critical value χ_crit² = 14 for 7 degrees of freedom (N = 8 data points) and a 95% significance level. We took data points from (a) the first, second, third and sixth harmonics of the full-loop system and (b) from the first and second harmonics of the coronal loop system. We then fit those data points with the first and second harmonics calculated for the coronal loop model with the exponential density profile, described in Eq. (15).

Fig. 7. Successful tests of χ2 goodness of fit between the eigenfunctions of the 1D oscillating string system (blue crosses, observations) and the eigenfunctions of the 1D system using a coronal model for the density and magnetic field distribution (solid orange line, prediction). Our models have a χ2 critical value of χcrit2 = 14 at α = 0.05 (95% significance level) with seven degrees of freedom. Bottom panels: Purely coronal system. Top and middle: Full-length 1D system. Shown are observations fit with the first harmonic of our predicted 1D model (left column) and the second harmonic (right). The first, second, third, and sixth harmonic are used as observations. The predicted scale height, magnetic field, scaling parameter, and χ2 test result are included in each panel.

Starting with the modes for the full loop, when fitting the observations of the third (sixth) harmonic of the full loop with the fundamental (first overtone) of the fitted 1D model gives us H = 41.9 Mm (H = 46.2 Mm) and B = 37.4 G. The magnetic field in particular is very close to the mean value in our 3D models (B_z = 30 G), showing the strength of this seismological technique. On the other hand, considering the fundamental and first overtone of the full-loop system as observations leads to vastly different results for the scale height (H = 9.2 and 2.1 Mm) and the magnetic field (B = 5.2 and 1.2 G). This big difference in the results is due to the fact that the first two harmonics of the full-loop system are the most influenced by the sharp gradients in the transition region, rendering the usual approach in seismology of fitting simple coronal model profiles problematic.

To test the validity of our χ² goodness of fit tests, we compared our results with more traditional seismological techniques, using the formulas for the ratio of periods derived in Andries et al. (2005b) and Safari et al. (2007):

$$\begin{array}{cc}& P_1=P_{\mathrm{Kink}}{(1+L/(3{\pi }^{2}H))}^{-1},\hfill \end{array}$$(18)

$$\begin{array}{cc}& P_2=P_{\mathrm{Kink}}{(2+2L/(15{\pi }^{2}H))}^{-1},\hfill \end{array}$$(19)

$$\begin{array}{cc}& P_3=P_{\mathrm{Kink}}{(3+3L/(35{\pi }^{2}H))}^{-1},\hfill \end{array}$$(20)

where P_Kink = 2L/C_Kink is the period of the fundamental standing kink mode for a loop of the same length, mean magnetic field, and no density stratification along its axis, with density equal to its value at the footpoints (Edwin & Roberts 1983). Using the modes for the full loop, the periods of the fundamental (P₁) and the second harmonic (P₂) give P₁/P₂ = 1.03 (for an observed loop length in coronal lines equal to L = 186 Mm) and a corresponding scale height of H = 4.03 Mm. The latter is the same order of magnitude to the values calculated from the χ² tests. If we instead misidentify the third and sixth harmonic as the first and second, then we get P_1, false/P_2, false = 1.76 and H = 34.9 Mm, which again is similar to the results of the χ² fit tests.

Applying the same χ² goodness of fit tests when using the first (second) harmonics of the coronal part of the loop as measurements, we then get H = 37.7 Mm (H = 39.1 Mm) for the scale height and B = 37.6 G (B = 39.1 G) for the magnetic field, respectively. The results can be seen in the bottom two panels of Fig. 7. These values for the derived scale height and the magnetic field are very close to those we get for the third and sixth harmonics of the full-loop system. Additionally, when filling in the values from the modes of the coronal part of the loop in Eqs. (18)–(20), and considering the same loop length as in the case of the full loop, we get P₁/P₂ = 1.75 and H = 34.6 Mm. These results further hint to the possibility of observations misidentifying the order of the observed standing kink modes, if we assume that the oscillations occurred along the entire loop length and not just the coronal part. For reference, when using the modes for the coronal part of the loop and considering the same loop length we get P₁/P₂ = 1.75 and H = 34.6 Mm.

3.4. Seismological implications for transverse oscillations

In the context of coronal seismology, it is also useful to revisit the results from older observational studies and try to understand them in light of our current results. In Duckenfield et al. (2018), both the first and second harmonics were detected in an observation of a decayless transverse oscillation of a coronal loop with a period ratio equal to P₁/P₂ = 1.4. The third harmonic was detected in an observation of a decaying transverse oscillation of coronal loops in Duckenfield et al. (2019), with the period ratio of the fundamental to the third harmonic equal to P₁/P₃ = 2.61. Considering our 1D eigenvalue problem for our full-loop system, the corresponding period ratios obtained are P₁/P₂ = 1.03 and P₁/P₃ = 1.77. These values seem to significantly deviate from the observations. However, if we assume that the third, sixth, and ninth harmonics were misidentified for the first, second, and third harmonics, we obtain ratios that match the observations: P_1, false/P_2, false = 1.76 and P_1, false/P_3, false = 2.57. Unsurprisingly by this point, if we consider the first three harmonics of the coronal loop system, the same ratios give P₁/P₂ = 1.75 and P₁/P₃ = 2.57, which are again a good match to the observations. Although a larger statistical sample is needed for observations and simulations of oscillating stratified loops, the period ratios above indicate that the observed modes for standing oscillations of longer coronal loops either are purely coronal in nature, with the part of the loop below the corona not being perturbed, or are the modes that exhibit nodes near the base of the corona.

Our analysis so far has focused only on our loop setup, which models a longer coronal loop with length 200 Mm. However, recent observations of short coronal and transition region loops with length ≤50 Mm have also revealed transverse decayless oscillations. Unlike with the longer coronal loops, most studies of decayless oscillations in short loops (e.g. Gao et al. 2022, 2024; Petrova et al. 2023; Shrivastav et al. 2024) could not conclusively characterise these waves as standing or propagating, leaving their interpretation still open. Despite that, the assumption of them being standing waves has been used extensively in these studies when considering them for seismology. Lim et al. (2024a) show that under-sampling can lead to overestimation of the oscillation periods, producing values that match the observations for shorter loops. The difference with respect to the longer coronal loops is that a higher percentage of the total length of these short loops will be located below the corona. Therefore, we expect these short loops to exhibit qualitatively different spatial profiles of their kink oscillation modes. To get better insight into the oscillatory profiles of these shorter loops, we repeated our analysis with the 1D wave equation eigenvalue system for a short loop of length 30 Mm. The kink speed profile that we used as phase speed in Eq. (12) was taken from the numerical model of a short coronal loop model that was used in Gao et al. (2023). This profile of the kink speed is shown in the left panel of Fig. 8. The approximate height of the top of the transition region is this model is located at height z = 4.2 Mm above each footpoint. As a reference, the coronal part of the kink speed profile is also shown. The right panel of Fig. 8 shows the spatial profiles of the standing kink modes for this system and their respective frequencies. For reference, the derived frequencies of the first, second and third harmonics reported in Gao et al. (2023) are f₁ = 13.98 mHz, f₂ = 18.62 mHz and f₃ = 24.04 mHz, with P₁/P₂ = 1.33 and P₁/P₃ = 1.72. From this panel, we see that it is easier to distinguish between the first and third harmonics, due to the proximity of the loop apex to the base of the corona. In fact, the oscillation of such a low-lying loop would likely be visible in transition region lines, as was the case with the observed decayless oscillations in transition region loops reported in Gao et al. (2024). Therefore, in the case of shorter, low-lying loops we expect less ambiguity over the order of the observed standing kink modes, than in longer coronal loops.

Fig. 8. Profiles of the kink speed along the z-axis (left) and the corresponding normalised eigenfunctions of the 1D system of Eq. (12) (right) for the short loop model used in Gao et al. (2023). The coronal part of the kink speed profile is highlighted in the left panel (z ∈ [4.2, 25.8] Mm, dashed orange line) and the eigenfunctions for the entire loop length in the right.

The aforementioned possible ambiguity over the order of standing kink modes in coronal loops, as it translates to the observations of transverse decaying and decayless oscillations can have seismological implications. For example, the observed periods of the fundamental mode and the loop length can be used to calculate the kink speed and through that average magnetic field, when combined with estimates of the plasma density and the ratio of densities inside and outside of the oscillating loops (see for example Anfinogentov & Nakariakov 2019; Zhong et al. 2023a; Gao et al. 2024). In Fig. 9 we plot the oscillation period and corresponding kink speed with respect to the loop length, for data taken from a number of past studies (Wang et al. 2012; Nisticò et al. 2013; Anfinogentov et al. 2013, 2015; Anfinogentov & Nakariakov 2019; Zhong et al. 2022a,b; Gao et al. 2022, 2024; Li & Long 2023; Petrova et al. 2023; Shrivastav et al. 2024). The left panel shows the scatter plot of the period versus the loop length. The blue circles represent the original observations for short loops, with a length of less than 50 Mm (dashed vertical line), the value of which was arbitrary chosen. We show the original observations of long coronal loops and the best linear fit (P = a L) with coefficient a = 1.103 ± 0.043 Mm s⁻¹. Estimating the kink speed as C_Kink = 2L/P_Kink, where P_Kink is considered as the period of the fundamental standing kink mode, we get an average kink speed C_Kink = 1.81 Mm s⁻¹. The estimated kink speed with respect to the loop length for each observation is shown in the right panel of Fig. 9. Considering now the possible ambiguity between the first and the third kink modes in long loops, we corrected the values of the observed periods in Fig. 9 by a factor of 1.77, since P₁ = 1.77 P₃. Here we used the results for our system as a suggested average correction for the observations of the long loops (L > 50 Mm). The corrected values for the period and the estimated kink speed are shown in both panels. From the new linear fit we get a coefficient a = 1.953 ± 0.076 and a new average kink speed of C_Kink = 1.02 Mm s⁻¹. Using common seismology techniques (e.g. Anfinogentov & Nakariakov 2019; Zhong et al. 2023a; Gao et al. 2024), we can calculate the average magnetic field of an oscillating loop with the following formula:

Fig. 9. Scatter plots for various parameters for observations of loops undergoing decayless oscillations, as reported in past works (Wang et al. 2012; Nisticò et al. 2013; Anfinogentov et al. 2013, 2015; Anfinogentov & Nakariakov 2019; Zhong et al. 2022a,b; Gao et al. 2022, 2024; Li & Long 2023; Petrova et al. 2023; Shrivastav et al. 2024). Left: Relation between the period and the loop length, as well as the best linear fit for the long loops. Right: Relation between the kink speed and the loop length. The dashed lines at 50 Mm divide the short- and long-loop data. The blue circles correspond to the short loops, the black circles to the long loops (original data), and the red circles to the long loops, after adjusting the data for measuring the third harmonic.

$$\begin{array}{c}\hfill B=C_{\mathrm{Kink}}\sqrt{\frac{{\mu }_0({\rho }_i+{\rho }_e)}{2}}=\frac{2L}{P_{\mathrm{Kink}}}\sqrt{\frac{{\mu }_0({\rho }_i+{\rho }_e)}{2}},\end{array}$$(21)

where μ₀ is the magnetic permeability. If we take for the observed P_Kink = P₃ instead of P₁, due to the described ambiguity, then for the average calculated magnetic field we have

$$\begin{array}{c}\hfill B_{\mathrm{calc}}=\frac{2L}{P_3}\sqrt{\frac{{\mu }_0({\rho }_i+{\rho }_e)}{2}}=\frac{P_1}{P_3}\left(\frac{2L}{P_1}\sqrt{\frac{{\mu }_0({\rho }_i+{\rho }_e)}{2}}\right),\end{array}$$(22)

and therefore B_calc = (P₁/P₃) B_real. For a mean background magnetic field of B_real = 30 G, as in our model, the field calculated from the observed quantities will be overestimated by a factor of P₁/P₃ (here 1.77), giving B_calc = 1.77 × 30 ∼ 53 G.

4. Discussion and conclusions

We explored the nature of the spectra in a gravitationally stratified loop with footpoints anchored in chromospheric plasma by performing transverse oscillations. Our base model consists of a 3D straight flux tube, which is perturbed in three different simulations by (a) a footpoint linearly polarised broadband velocity driver, (b) a symmetric sinusoidal velocity pulse, and (c) an asymmetric off-centre Gaussian pulse. We limited our 3D study mostly to the coronal part of our system and treated the chromosphere as a mass reservoir of near-constant temperature. We also performed a complimentary analysis, solving the generalised eigenvalue problem for the 1D wave equation with nodes at the footpoints, for a system with phase speed equal to the kink speed profile of the 3D loop. Finally, we performed forward modelling with our simulation data and then explored the seismological implications of our findings.

When tracking the centre of mass of our loop, the oscillation spectra reveal the existence of multiple harmonics when either the broadband driver or the off-centre velocity pulse are used to perturb the flux tube, as shown in Fig. 3. Our 1D analysis reveals that the observed signatures correspond to the frequencies of the system eigenmodes (Fig. 6), which have been deformed due to the large gradients in the values of density with height. The so-called half-harmonic that was reported in our past studies (Karampelas & Van Doorsselaere 2024; Karampelas et al. 2024) is now understood to consist of the deformed fundamental kink mode and the second harmonic, and in the case of the broadband driver it also has the signatures of driving frequencies. What was understood in past studies as the fundamental mode is in fact the third harmonic of the system. However, our analysis shows that this third harmonic of the full-length loop (and also the sixth, ninth, and twelfth harmonics) exhibits similar frequencies and spatial profiles in the corona as the fundamental mode (and also the second, third, and fourth harmonics) of the coronal part of that loop. This similarity is due to the nodes that these harmonics exhibit near the base of the corona for our system.

Creating synthetic observations of our 3D simulations and tracking the centre of the loop at different resolutions, we extracted the oscillation frequencies that would have been detected in real observations of a system matching our model. The broadband driver excites a standing wave with many frequencies present, consisting of multiple harmonics and frequencies associated with the driving. Most observations of decayless oscillations seem to have a single dominant period (e.g. Anfinogentov et al. 2015; Petrova et al. 2023), and it is rare to detect higher-order harmonics (Duckenfield et al. 2018). As such, our findings cast doubts on the idea that broadband drivers could be one of the possible mechanisms responsible for the excitation of decayless oscillations in coronal loops. For a conclusive answer on the nature of the drivers of these oscillations, realistic driving profiles that follow observational data need to be employed in future numerical studies.

Focusing on the oscillations excited by the initial velocity pulses, our spectra from the synthetic observations in the 171 Å line reveal only one dominant frequency band, which matches the third harmonic of the full-length loop. This mode, which also resembles the fundamental mode of the coronal part of the loop, can easily be misidentified as the fundamental mode of the full-length loop. The time-distance maps of the oscillating loops at the apex in Fig. 4 show signals of decaying amplitude over time with a single frequency. However, the oscillations get progressively harder to track after the first two to three periods due to the small initial amplitude. Considering the overlap of features in real observations, or the different levels of contrast of the intensity between the loop and the background, it is possible for these decaying oscillations to be categorised as decayless with a fundamental frequency at ∼4.75 mHz unless a longer part of the signal is clearly visible for a longer time. It is worth noting that some recent observations of decayless oscillations have shown them to have relatively short durations (up to four cycles; e.g. Shrivastav et al. 2024; Lim et al. 2024b), creating ambiguity as to their interpretation as decaying or decayless oscillations, similar to what was discussed for our data.

As mentioned, the third harmonic of our oscillating loop can be misidentified as the fundamental kink mode of the loop (see also Howson & Breu 2023). This misidentification is caused by the deformation of the harmonic modes due to the low values of the kink speed and the density scale height at the loop footpoints, below the transition region. However, our 3D numerical setup is characterised by a very wide chromosphere (Δz = 5 Mm) and an unphysically broadened transition region, which are not representative of the solar atmosphere. Therefore, our results are expected to differ from those for a more realistic profile of the solar atmosphere. Past studies have treated coronal loops as resonant cavities for low-frequency fluctuations (e.g. Hollweg 1984; Verdini et al. 2012), with resonant waves being confined in the coronal part and experiencing wave leakage through the sharp transition region. It is unclear, therefore, how a realistic chromospheric and transition region profile will affect the kink modes of our 3D model. Our 1D analysis predicts the existence of the previously mentioned deformed harmonics, as well as the respective values for the frequencies of these modes, as derived from the 3D simulations. However, as stated in Sect. 3.2, the thin tube approximation used to derive Eq. (12) is not valid below the transition region for our model. This means that the results of our 1D analysis can only be taken as an approximation of the actual transverse modes of the system. In addition, the results of the 1D analysis also depend on the profiles of the loop and the lower solar atmosphere considered in our model. In Appendix A we show the effects of different kink speed profiles on the solutions of Eq. (12), highlighting this further. From the above, it becomes clear that a proper parameter study based on realistic profiles of the lower solar atmosphere is required (a) to enable a detailed examination of the spatial profiles of the oscillation modes in the 3D simulations, (b) to examine the validity of using Eq. (12) for loops anchored in chromospheric plasma, and (c) to conclude on the possibility of misidentifying the order of the observed transverse oscillations in loops, and the applicability of our findings.

The presence of transverse oscillations has been confirmed in shorter loops (e.g. Gao et al. 2022, 2024; Li & Long 2023; Petrova et al. 2023; Shrivastav et al. 2024), with the recent study by Lim et al. (2024a) showing that they can be interpreted as standing waves with periods shorter than the temporal resolution of our current observatories. This suggests that the transverse oscillations observed in coronal loops might not necessarily be confined to the corona, similar to how the kink modes in our 3D numerical model also persist along the entire length of the loop. This seems to invalidate the treatment of loops as resonant cavities that confine waves in their coronal part while taking the effects of leakage from the footpoints into account (e.g. Hollweg 1984; Verdini et al. 2012). However, we again note that our models are characterised by an artificially broadened transition region. The effects of a sharp transition region on transverse standing waves in loops need to be properly explored in future numerical studies. Due to this uncertainty, observations at chromospheric and transition region lines at the footpoints of oscillating loops are necessary to identify the observed kink modes as either purely coronal modes or modes of the full-loop system before making any seismological estimations.

More advanced techniques of spatial coronal seismology, such as the χ² optimisation process shown in this study, produce results with a good agreement between full-loop harmonics and coronal loop harmonics. However, many data points along the entire length of the loop are required to properly constrain the estimated mean magnetic field values. Therefore, we conclude that the current approach in seismology, where only the standing modes in the coronal part of the loops are considered, can only be justified if there is definitive proof that the oscillatory motion is confined to the corona. Unveiling the dynamics of oscillating loops in the lower corona, the chromosphere and the transition region through simultaneous high-resolution and high-cadence observations from upcoming missions such as MUSE is essential for future seismological studies.

Acknowledgments

K.K. acknowledges support by an FWO (Fonds voor Wetenschappelijk Onderzoek – Vlaanderen) postdoctoral fellowship (1273221N). DL was supported by a Senior Research Project (G088021N) of the FWO Vlaanderen. TVD was supported by the C1 grant TRACEspace of Internal Funds KU Leuven, and a Senior Research Project (G088021N) of the FWO Vlaanderen. Furthermore, TVD received financial support from the Flemish Government under the long-term structural Methusalem funding program, project SOUL: Stellar evolution in full glory, grant METH/24/012 at KU Leuven. The research that led to these results was subsidised by the Belgian Federal Science Policy Office through the contract B2/223/P1/CLOSE-UP. It is also part of the DynaSun project and has thus received funding under the Horizon Europe programme of the European Union under grant agreement (no. 101131534). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union and therefore the European Union cannot be held responsible for them. Y.G. was supported by China Scholarship Council under file No. 202206010018. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation Flanders (FWO) and the Flemish Government – department EWI.

References

Appendix A: Solving the 1D wave equation eigenvalue problem for different phase speeds

In Sect. 3.2 we solve the Sturm-Liouville problem of Eq. 12 to find the modes and frequencies of non-axisymmetric, transverse oscillations of loops, considering the thin tube approximation. To further explore the effects that the profile of the lower solar atmosphere has on the solutions of Eq. 12, we performed a parameter study that considers different phase speed profiles. We used a set of five profiles (C₁ to C₅) that explores the effects of the width of the chromosphere and the transition region, as well as the values of the kink speed in the chromosphere with respect to the corona. The C₁ profile is the kink speed of our post-relaxation 2D loop profile, as seen in Fig. 1. C₂ is the kink speed derived from the initial conditions of our pre-relaxation 2D profile, and is characterised by a sharper transition region. Both C₁ and C₂ are derived for a loop with a wide chromosphere (Δz = 5 Mm). C₃ is a kink speed profile derived for a loop with a narrower chromosphere (Δz = 2 Mm). C₄ and C₅ are modified versions of C₃, for uniform chromospheric values of the kink speed.

The results of our parameter study can be seen in Fig. A.1. Comparing the harmonics for the C₁ and C₂, we see that the width of the transition region, for the range considered here, does not generate qualitatively different results regarding the shape of the kink modes or the values of the frequencies for each harmonic. Reducing the width of the chromosphere (C₃) can influence the shape of the kink modes and their frequencies substantially, but the results are highly affected by the properties of the chromosphere and the shape of the resulting kink speed profile there (C₄, C₅). For example, having chromospheric kink speeds values that are less than 10% of the coronal ones (C₄) can still lead to results similar to those presented for the C₁ profile, even for a narrow chromosphere with a sharp transition region.

Fig. A.1. Different kink speed profiles used as phase speeds in Eq. 12 and its respective solutions. Left: Kink profiles near one footpoint. Right: Corresponding first (solid blue lines) and third (dashed orange lines) harmonics along the entire length of the loop.

We conclude that the width of the solar chromosphere and the local properties that define the kink speed of a loop are of great importance when deriving the kink modes of the reduced 1D Sturm-Liouville problem of Eq. 12.

All Figures

Fig. 1. Post 2.5D MHD relaxation profiles for the temperature, density, Bz magnetic field, and kink speed of our stratified flux tube. Left: 2D temperature profile for a section of our model up to r = 3 Mm. Top right: 1D profiles along the z direction of the density (black lines) and Bz magnetic field (red lines) at r = 0 (solid lines) and at r = 8 Mm (dashed lines). Bottom right: Profile of the kink speed along the z-axis, with the coronal part (z ∈ [7, 193] Mm) highlighted (dashed orange line).

In the text

	Fig. 2. Profiles of the velocity driver and pulses used in the three 3D simulations. Left panels: Temporal profile of the vx driver at z = 0 for the first simulation (top) and its corresponding power spectrum (bottom). Right panel: Two initial velocity conditions used in the two simulations without the footpoint driver.
In the text

Fig. 3. Centre of mass displacement (left panels) and corresponding power spectra density profiles (right panels) for the three simulations of our 3D loop. The time-distance maps of the displacement along the x-direction show the magnified (×20) signal projected along the z-axis for visualisation purposes. Only the coronal part of the loop is depicted. Shown are the simulations with the driver (top), with the symmetric initial pulse (middle), and with the asymmetric pulse (bottom).

In the text

Fig. 4. Time-distance maps of the synthetic emission in the 171 Å line, for the three cases of oscillating loops, at the apex. Left column: SDO/AIA resolution. Middle column: Solar Orbiter/HRIEUV resolution. Right column: MUSE/SG resolution. Top panels: Loop with the broadband driver. Middle panels: Loop with the sinusoidal, symmetric pulse. Bottom panels: Asymmetric, off-centre Gaussian pulse. Overplotted in each panel are the oscillating signal of the loop centre (black asterisks) and a fitted sinusoidal function (solid red line) with the four prominent frequencies for each signal, as detected by Auto-NUWT.

In the text

Fig. 5. Distribution of oscillation frequencies, detected with Auto-NUWT at each height, z, from the synthesised 171 Å line emissions. The colour bars correspond to the normalised amplitude of the oscillation for each detected frequency. As in Fig. 4, the left, middle, and right columns correspond to SDO/AIA, Solar Orbiter/HRIEUV, and MUSE/SG resolution, while the top, middle, and bottom panels correspond to the simulations with the broadband driver, symmetric pulse, and asymmetric pulse, respectively.

In the text

	Fig. 6. Normalised eigenfunctions of the 1D oscillating string system and its corresponding frequencies. Left: Solution for the entire system. Right: Solution when only the coronal part of the system is used. The cyan shaded region highlights the coronal part in the left panel.
In the text

Fig. 7. Successful tests of χ2 goodness of fit between the eigenfunctions of the 1D oscillating string system (blue crosses, observations) and the eigenfunctions of the 1D system using a coronal model for the density and magnetic field distribution (solid orange line, prediction). Our models have a χ2 critical value of χcrit2 = 14 at α = 0.05 (95% significance level) with seven degrees of freedom. Bottom panels: Purely coronal system. Top and middle: Full-length 1D system. Shown are observations fit with the first harmonic of our predicted 1D model (left column) and the second harmonic (right). The first, second, third, and sixth harmonic are used as observations. The predicted scale height, magnetic field, scaling parameter, and χ2 test result are included in each panel.

In the text

	Fig. 8. Profiles of the kink speed along the z-axis (left) and the corresponding normalised eigenfunctions of the 1D system of Eq. (12) (right) for the short loop model used in Gao et al. (2023). The coronal part of the kink speed profile is highlighted in the left panel (z ∈ [4.2, 25.8] Mm, dashed orange line) and the eigenfunctions for the entire loop length in the right.
In the text

Fig. 9. Scatter plots for various parameters for observations of loops undergoing decayless oscillations, as reported in past works (Wang et al. 2012; Nisticò et al. 2013; Anfinogentov et al. 2013, 2015; Anfinogentov & Nakariakov 2019; Zhong et al. 2022a,b; Gao et al. 2022, 2024; Li & Long 2023; Petrova et al. 2023; Shrivastav et al. 2024). Left: Relation between the period and the loop length, as well as the best linear fit for the long loops. Right: Relation between the kink speed and the loop length. The dashed lines at 50 Mm divide the short- and long-loop data. The blue circles correspond to the short loops, the black circles to the long loops (original data), and the red circles to the long loops, after adjusting the data for measuring the third harmonic.

In the text

	Fig. A.1. Different kink speed profiles used as phase speeds in Eq. 12 and its respective solutions. Left: Kink profiles near one footpoint. Right: Corresponding first (solid blue lines) and third (dashed orange lines) harmonics along the entire length of the loop.
In the text