Decayless oscillations in 3D coronal loops excited by a power-law driver

Aims. We studied the manifestation of decayless oscillations in 3D simulations of coronal loops, driven by random motions. Methods. Using the PLUTO code, we ran magnetohydrodynamic (MHD) simulations of a straight gravitationally stratified flux tube, with its footpoints embedded in chromospheric plasma. We consider transverse waves drivers with a horizontally polarised red noise power-law spectrum. Results. Our broadband drivers lead to the excitation of standing waves with frequencies equal to the fundamental standing kink mode and its harmonics. These standing kink oscillations have non-decaying amplitudes, and spectra that depend on the characteristics of the loops, with the latter amplifying the resonant frequencies from the drivers. We thus report for the first time in 3D simulations the manifestation of decayless oscillations from broadband drivers. The spatial and temporal evolution of our oscillation spectra reveals the manifestation of a half harmonic, which exhibits half the frequency of the identified fundamental mode with a similar spatial profile. Our results suggest that this mode is related to the presence of the transition region in our model and could be interpreted as being the equivalent to the fundamental mode of standing sound waves driven on pipes closed at one end. Conclusions. The potential existence of this half harmonic has important implications for coronal seismology, since misinterpreting it for the fundamental mode of the system can lead to false estimations of the average kink speed profile along oscillating loops. Finally, its detection could potentially give us a tool for distinguishing between different excitation and driving mechanisms of decayless oscillations in observations.


Introduction
The term decayless transverse (or kink) oscillations has been used over the past decade to describe a category of generally small-amplitude transverse oscillations in coronal loop observations in the extreme ultraviolet (EUV).The name is inspired by the near constant value of their amplitude, persisting over a large number of oscillation periods (Nisticò et al. 2013).This comes in stark contrast to the other known regime of large amplitude fast decaying transverse loop oscillations (Aschwanden et al. 1999;Nakariakov et al. 1999).Both decaying and decayless transverse oscillations are treated as standing kink oscillations in cylindrical loops (e.g.Edwin & Roberts 1983;Van Doorsselaere et al. 2008;Anfinogentov et al. 2015), and alongside propagating kink waves they are some of the most studied and ubiquitous phenomena in the solar corona (see Nakariakov et al. 2021, for a review).Ever since their discovery in coronal loops (Wang et al. 2012;Tian et al. 2012), the ubiquitous nature of decayless oscillations has been confirmed by a large number of observational studies in coronal loops (e.g.Anfinogentov et al. 2013Anfinogentov et al. , 2015;;Zhong et al. 2022a,b), in short coronal loops with length of a tens of Mm (e.g.Petrova et al. 2023;Shrivastav et al. 2023), and in coronal bright points (CBPs) (e.g.Gao et al. 2022).
Unlike their decaying counterparts, which are excited by nearby transients (see Nechaeva et al. 2019), the constant ampli-tude of these decayless oscillations suggests a continuous supply of energy that will counteract the damping from effects like resonant absorption (e.g.Goossens et al. 2002), phase mixing (e.g.Heyvaerts & Priest 1983), and the Kelvin-Helmholtz instability (KHI; Van Doorsselaere et al. 2021a,b).This makes non-damping kink oscillations a potential heating mechanism for the solar atmosphere (e.g.Lim et al. 2023; see also Van Doorsselaere et al. 2020 for a review).In addition, decayless oscillations can be used in coronal seismology, as in Anfinogentov & Nakariakov (2019) where they were used to determine the distribution of kink and Alfvén speeds in active regions.
However, the excitation mechanism of decayless oscillations is still under debate.In Antolin et al. (2016), they are hypothesised to be the result of (line of sight) LOS effects due to KHI.Nisticò et al. (2013) proposed instead a harmonic driver in resonance with the loop, an idea used extensively in simulations (e.g.Karampelas et al. 2017;Guo et al. 2019).However, harmonic drivers lead to equal manifestation of horizontally and vertically polarised oscillations, while observations favour horizontal polarisations (Anfinogentov et al. 2015).Excitation via supergranulation flows as a self-oscillation was proposed (Nakariakov et al. 2016(Nakariakov et al. , 2022)), and was shown numerically (Karampelas & Van Doorsselaere 2020) to create horizontally linearly polarised oscillations.However, this can not explain decayless oscillations of loops rooted in sunspots, where supergranulation flows are absent (Mandal et al. 2022).Excitation by vortex shedding (Nakariakov et al. 2009) can also generate non-damping oscillations (Karampelas & Van Doorsselaere 2021) and could explain observations of flare-induced decayless oscillations in Mandal et al. (2022).Afanasyev et al. (2020) explored the spectrum of excited oscillations from a red noise driver in a 1D study, providing the first proof of concept for broadband drivers exciting multiple harmonics in decayless oscillations (see also Ruderman & Petrukhin 2021;Ruderman et al. 2021).The contribution of footpoint-driven transverse waves with a multi-frequency spectrum in the energy evolution and heating of 3D coronal loops has also been studied in recent years (e.g.Pagano & De Moortel 2019;Pagano et al. 2020;Howson & De Moortel 2023).However, no numerical study of 3D coronal loops perturbed by broadband drivers has reported on the excitation and evolution, spatial and temporal, of standing waves matching the observed decayless oscillations, in the presence or absence of a transition region and a lower chromosphere.
In this letter we study the spectra of decayless oscillations in a 3D gravitationally stratified loop, generated by a driver with a power-law spectrum.In Sect. 2 we describe our initial and boundary conditions and the numerical scheme used throughout our simulations.Our findings are presented in detail in Sect.3. Finally, a thorough discussion of our results and their implications for the topic focused on here takes place in Sect. 4.

Numerical set-up
Initial conditions.We modelled a coronal loop with a chromospheric part and a transition region, following Pelouze et al. (2023) and Guo et al. (2023).We calculated the hydrostatic equilibrium in a 2.5D slice with r ∈ [0, 8] Mm and z ∈ [0, 200] Mm, a uniform grid (200 × 2048 points), and sinusoidal gravity (g z = 274 cos(π z/200) m s −2 ) along the vertical magnetic field lines B = B z k = 30 G, using the following profiles: (5) The subscripts Ch, Cor, and C refer to the chromospheric (z = 0 and z = 200 Mm), coronal, and apex (z = 100 Mm) values, respectively, inside (i) and outside (e) of the loop.The function ζ(r) gives us the radial loop profile, with a crosssection radius of R = 1 Mm.Finally, we have T Ch = 0.02 MK, T C,i = 1 MK, T C,e = 1.5 MK, ρ Ch,i = 3.51 × 10 −8 kg m −3 , and ρ Ch,e = 1.17 × 10 −8 kg m −3 .We let the system evolve for 3890 s, reducing each velocity component per iteration (v i = v i /1.0001) for 2334 s, reaching a semi-equilibrium state with maximum residual velocity values of ∼4 km s −1 along the z direction and ∼0.05 km s −1 along the x and y directions.The initial conditions after the 2.5D MHD relaxation are shown in Boundary conditions.In the 2.5D set-up, we consider axisymmetry at r = 0 and open boundaries at r = 8 Mm.At the z = 0 and z = 200 Mm, we take zero-gradient conditions for the magnetic field, antisymmetry for the velocity, symmetry for the density, and constant temperature T Ch .For the 3D setup, we have open boundaries at x = −6 Mm, x = 6 Mm, and y = 4 Mm and a reflective boundary at y = 0, simulating half the loop.At z = 0 and z = 200 Mm, the boundaries are the same as in the 2.5D case, but with added symmetry for the pressure.At t 0 = 202 s, we apply a velocity driver ({v x , v y } = {V(t) ζ(r, t), 0}) at z = 0, where we take R = R d = 2.5 Mm in ζ(r, t), to ensure that the loop always remains inside the area of the driver.For the velocity signal V(t), we take a red noise power-law spectrum (S ∝ f −1.66 , with f being the frequency), using the colorednoise 2.1.0Python package.Our driver is shown in Fig. 2; the orange line depicts its background trend, calculated with a low pass filter.Our driver also tracks the location of the footpoint (r(t)) by numerically integrating the velocity signal over time.
Numerical scheme.We solve the full magnetohydrodynamic (MHD) equations for a hydrogen plasma with the PLUTO code (Mignone et al. 2007), using the parabolic method in 2.5D and the MP5 method in 3D, the Roe solver, the third-order Runge-Kutta method, and the extended GLM formulation.The effective numerical diffusivity (η) is estimated at ∼10 −5 −10 −4 in units of the inverse magnetic Reynolds number.In the 2.5D set-up, we also add explicit magnetic diffusivity η = 10 −4 .We include thermal conduction (κ = 9 × 10 −12 T 5/2 in J s −1 K −1 m −1 ) and saturation effects for large temperature gradients.For T ≤ T cut = 0.25 MK, we apply the correction by Linker et al. (2001), Lionello et al. (2009), and Mikić et al. (2013), to treat the transition region with a coarser grid.

Results
We limit our analysis to the coronal part of the loop, while using the lower chromospheric part as a mass reservoir.As is shown below, our driver displaces the loop while also exciting kink oscillations.To study these oscillations, we track the horizontal position of the loop centre of mass for each z plane, by calculating the average displacement, using as weight the quantity (ρ(z) − ρ e (z)) 2 , where ρ e (z) is the density outside the loop.The averaging takes place across the entire domain at each z-plane.This gives us the transverse displacement of the loop along x over time, for each height.
The top left panel of Fig. 3 shows the displacement of the coronal part of the loop for waves excited by the red noise driver.
L6, page 2 of 6  We considered a number of 'slits' (21 for z ∈ [25, 175]) and projected the signal along z.For better visualisation we also magnified the signal of each slit by simply multiplying by a factor of 20.This multiplication is not employed anywhere else in our analysis, but is only used to visually show the loop motion.The right panel shows the normalised Fourier power spectral density (PSD), along the loop, with resolution along z corresponding to the number of slits used.The low frequency components of the red noise driver lead to large a displacement of the driven footpoint, diminishing as we travel along z towards the anchored footpoint.This aperiodic displacement also manifests in the spectra shown in the top right panel of Fig. 3, for frequencies close to 0. We also see the spatial profile of the corresponding frequencies, being stronger closer to the driven footpoint.
In order to focus on the periodicities of the observed oscillating patterns in the top left panel of Fig. 3, we removed the background trend of the aperiodic displacement from our signal by applying a high pass Gaussian filter, similarly to what could have been applied on a signal derived from EUV observations of decayless oscillations.This process is also the same as that used to calculate the background trend of our driver, shown in Fig. 2. The now detrended displacement and respective Fourier spectra are shown in the middle panels of Fig. 3.In the middle right panel, we now clearly see two bands of frequencies, one centred at ∼2.5−3 mHz and one at ∼4−5 mHz.The nature of these bands, which are also visible in the spectrum of the full signal shown in the top right panel of the same figure, can be understood by identifying the fundamental kink mode of our loop.
The top left panels of Fig. 4 show the non-magnified detrended displacement signal at the apex and its correspond-ing wavelet spectrum over time, for our loop.We see that the detrended signal depicts a transverse oscillation with amplitudes of the order of ∼0.2 Mm.This amplitude shows no significant decay over time, and is typical of the values found in decayless oscillations of coronal loops (e.g.Anfinogentov et al. 2015).The two frequency bands at ∼2.5−3 mHz and ∼4−5 mHz can also be seen here, persisting over time, with the former exhibiting an overall stronger signal.As a control, we also performed a simulation of the same loop oscillating with a non-driven decaying fundamental standing kink mode, excited by an initial velocity field.The signal and wavelet spectrum is shown in the top right panels of Fig. 4, from which we confirm that the ∼4−5 mHz frequency band corresponds to the fundamental standing kink mode for our oscillating loop.Looking again at the middle right panel of Fig. 3, we see that the spatial profile of this frequency band indeed matches the expected one for the fundamental standing kink mode.Having oscillations with the period of the fundamental standing kink mode and a non-decaying amplitude, we are led to describe these excited standing waves as decayless oscillations.
Alongside the fundamental mode, in the middle right panel of Fig. 3 we can also see patterns of additional higher harmonics: the second (∼7.5−8 mHz), third (∼12 mHz), fourth (∼16 mHz), and even fifth (∼20 mHz) harmonic.These were made clear after the removal of the contribution from the lower frequencies.However, they can also be detected in the PSD contours in the top right panel of Fig. 3, although their signal is much weaker than that of the low frequency loop displacement of the original time series.
The ∼2.5−3 mHz frequency band, which we descriptively refer to as the 'half harmonic', has frequencies around half that L6, page 4 of 6 of the identified fundamental mode.This frequency band is also present within the 99% significance level contours in the wavelet spectra of the decaying fundamental kink mode (see top right panels of Fig. 4).As can be seen by the wavelets, this mode is only weakly excited by the initial velocity field.Finally, both the half harmonic and the fundamental frequency bands have finite widths in the PSD graph, due to the limited sample size of the signal, as well as the development of instabilities and the restructuring of the loop density.
To ensure that the identified frequency bands are the result of the loop filtering and magnifying its natural frequencies, and not exclusively the result of the power distribution of our driver, we performed an additional set of simulations, one with the red noise driver and one with an impulsive oscillation, for a different loop, with initial coronal magnetic field of B z = 20 G.This loop has temperature and density profiles very similar to those shown in Fig. 1 for the loop with B z = 30 G. The detrended displacement and corresponding Fourier spectrum is shown in the bottom panels of Fig. 3.The bottom panels of Fig. 4 again show the signal and wavelet spectrum at the apex for the detrended displacement and for an impulsive decaying oscillation.By studying the Fourier and wavelet PSDs, we see that the new fundamental mode (∼3 mHz) and its corresponding half harmonic (∼1.5 mHz) and overtones (second at ∼6 mHz, third at ∼9 mHz) are excited.
Finally, we performed a simulation of our original loop with a coronal magnetic field of B z = 30 G, where we simultaneously employ two drivers: (a) the original one (see Fig. 2) at z = 0 and (b) a new driver at z = 200 Mm.The second driver, shown in the left panel of Fig. 5, was created from a new red noise spectrum, different from that of the original driver, with the additional step of removing the low frequency background trend.Driving the loop from both footpoints led again to the superposition of a decayless oscillation and a low frequency displacement.Taking the PSD of the detrended displacement signal in the right panel of Fig. 5, we see its similarity to that in the middle right panel of Fig. 3.The main difference comes from the second driver, which increases the power of the half harmonic near the second footpoint.

Discussion and conclusions
In this study we report for the first time in 3D MHD simulations the excitation of decayless kink oscillations driven by a footpoint power-law driver, for a gravitationally stratified coronal loop with footpoints embedded in chromospheric plasma.The amplitudes of our oscillations are of the order of 0.2 Mm and do not show any signs of decay for the duration of our simulations (∼5−7 oscillation periods, depending on the characteristics of each loop used).Figures 3 and 4 also show a modulation of the oscillation amplitudes, similar to what was reported in Nakariakov et al. (2022) for a model of self-oscillations with additional random-motion terms.However, further exploration of this effect with respect to past studies of self-oscillations is required, which falls outside the scope of this present work.The low frequency component in our driver does not suppress the decayless oscillations when the fully compressible 3D MHD equations are considered, but instead adds an overlaying displacement that appears aperiodic.Expanding upon this, driving a loop at any frequency band would manifest in the oscillation spectrum, for example in simulations of decayless oscillations of short coronal loops driven by p-modes (Gao et al. 2023).Using a high pass Gaussian filter to remove this aperiodic, very low frequency motion, we again end up with a signal of a decayless oscillation.
Using loops of different initial conditions, we see that the same driver leads to the excitation of standing modes for the resonant frequencies corresponding to each loop.This further supports the findings of the 1D study in Afanasyev et al. (2020), proving that coronal loops can filter the driving frequencies in broadband drivers by responding to the resonant ones.We also see our loops manifesting the fundamental standing kink mode and its overtones and we identify modes as high as the fifth harmonic.Such overtones, mainly the second and third harmonic, have also been observed in kink oscillations of loops in EUV, both for decayless (e.g.Duckenfield et al. 2018) and decaying (e.g.Duckenfield et al. 2019) oscillations, with our model also reproducing them for random motion footpoint driving in the case of decayless oscillations.
In addition to the fundamental mode and its overtones, our model also predicts the excitation of a half harmonic for each loop, with wavelength λ half = 4L, where L is the loop length.This frequency band scales with the fundamental, with faint traces of it detected in the spectra of decaying kink oscillations.Its scaling with the fundamental mode is difficult to explain, as it does not agree with the expected solutions of the wave equation in a stratified medium with a transition region present (see Howson & Breu 2023, for the harmonics in a 1D model using Alfvén waves).A possible interpretation can be given by considering the effects of the transition region, which is not a perfectly reflective boundary (Pelouze et al. 2023).The wave reflections there will set up standing waves, establishing the L6, page 5 of 6 fundamental standing kink modes and its overtones that we detect.At the same time, the movement of the footpoints in the transition region will create a velocity antinode there.As the loops oscillate, this will lead to a superposition of waves, each seeing one 'closed' and one 'open' footpoint, leading to the manifestation of the half harmonic with a frequency that is half that of the fundamental.The equivalent mechanical analogue is that of the fundamental mode of standing sound waves in cylindrical pipes closed at only one end (like pan flutes) having half the frequency of the fundamental in open-ended (or closed-ended) pipes.The lower velocities in the transition region for footpoints anchored at the chromosphere, could explain why our simulations of decaying oscillations also show faint traces of this half harmonic.
The exact nature of this half harmonic needs to be further explored in a dedicated study.However, this predicted half harmonic can have profound effects in coronal seismology.The ratio of the periods from multiple harmonics detected in standing loop oscillations can be used to calculate the profile of the average kink speed (C K ) along the loop (Jain & Hindman 2012).For example, the models by Andries et al. (2005) and Safari et al. (2007) attribute deviations from the expected values of these ratios to density stratification, calculating the periods of the modified first, second, and third harmonics as Here P kink = 2L/C K is the period of the fundamental standing kink mode and H the density scale height.A half harmonic with a similar spatial profile to the first harmonic, but with significantly more power, as our Fourier analysis suggests, could be misidentified as the P 1 mode, leading to false estimations of the density scale height for the model previously mentioned, the magnetic field or the C K profile along the loop in general.
Misidentifying this half harmonic could also lead to errors in the slopes of a factor of 2, when studying the correlation of the oscillation period and the loop length (e.g.Anfinogentov et al. 2015;Shrivastav et al. 2023).Therefore, this predicted half harmonic, if observed, can be of great importance when using decayless oscillations as tools for coronal seismology.Finally, we note that such a half harmonic would be more prominent if the loop is driven by displacing its footpoints with such broadband drivers.Models of decayless oscillations that require both footpoints to be anchored or that they move primarily with very small amplitudes at dissonant frequencies, will exhibit fainter traces of this mode.It remains to be tested whether the detection of the half harmonic, or lack thereof, could give us a tool to distinguish between footpoint driving and alternative excitation mechanisms for decayless oscillations in loops, such as driving by vortex shedding, or generating selfoscillations through the loop interacting with supergranulation flows.

Fig. 1 .
Fig. 1.Temperature, density, and B z magnetic field along z inside (solid line) and outside (dashed line) of the flux tube, at the end of the 2.5D MHD relaxation.

Fig. 2 .
Fig. 2. Velocity signal (left panel) and spectrum (right panel) for our red noise driver.The orange line shows its velocity trend.

Fig. 3 .Fig. 4 .
Fig. 3. Centre of mass x displacement and corresponding spectra for two different loops.Left panels: magnified (×20) x displacement of the loop centre of mass per height (projected on z) from slits placed along the coronal part of the loop.Right panels: fourier power spectra of the displacement taken from the slits on the left.The top and middle panels correspond to the original and detrended displacement signals for a loop with a coronal magnetic field B z ∼ 30 G. The bottom panels correspond to the detrended displacement signal for a loop with B z ∼ 20 G.

Fig. 5 .
Fig. 5. Signal of the second velocity driver at z = 200 Mm (left panel) and PSD of the respective centre of mass displacement (right panel) for a loop with a coronal magnetic field B z ∼ 30 G, driven from both footpoints.The original driver is also active at z = 0.