Modeling of Transverse Oscillations Driven by p-modes in Short Coronal Loops

Recent observations have revealed two types of decayless transverse oscillations in short coronal loops: one with short periods scaling with loop lengths, and the other with longer periods that exhibit a peak at around 5 minutes in the period distribution. To understand such a difference in period, we work in the framework of ideal MHD and model a short coronal loop embedded in an atmosphere with density stratification from the chromosphere to the corona. An inclined p-mode-like driver with a period of 5 minutes is launched at one loop footpoint. It is discovered that two types of decayless transverse oscillations can be excited in the loop. We interpret the 5 minutes periodicity as being directly driven by the footpoint driver, while the others, with periods of several tens of seconds, are regarded as kink eigenmodes of different harmonics. Therefore, our simulation shows that both types of decayless oscillations found in observations can be excited by p-modes in one short coronal loop. This study extends our understanding of ubiquitous decayless transverse oscillations in the corona. Furthermore, it suggests that p-modes could be an important energy source for coronal heating by driving decayless transverse oscillations.


Introduction
Transverse oscillations in coronal loops have been a subject of intensive investigations in solar physics since their first imaging observation in 1999 (Aschwanden et al. 1999;Nakariakov et al. 1999).They can be used as a powerful diagnostic tool of the coronal magnetic field (e.g., Nakariakov & Ofman 2001;Verwichte et al. 2013;Chen & Peter 2015;Su et al. 2018;Yang et al. 2020aYang et al. , 2020b)), and can also contribute to our understanding of coronal heating (see the review by Van Doorsselaere et al. 2020).External energy release processes can excite decaying kink oscillations in coronal loops (Zimovets & Nakariakov 2015;Zhang et al. 2020;Li et al. 2023).These oscillations begin with a large amplitude, and then rapidly decay over several periods (e.g., Nisticò et al. 2013;Goddard et al. 2016).Coronal loops can also support transverse oscillations without discernible excitation processes.Such oscillations were first discovered by Wang et al. (2012) and Tian et al. (2012) through imaging and spectroscopic observations respectively at nearly the same time.They have a smaller amplitude without any obvious damping, and are therefore called decayless oscillations (Nisticò et al. 2013).Both oscillations are believed to be standing kink waves, and their period scales linearly with the loop length in observations (e.g., Anfinogentov et al. 2015;Goddard et al. 2016).
Decayless oscillations are found to be ubiquitous in active region loops (Anfinogentov et al. 2015) and can be continually observed for a long time (e.g., Zhong et al. 2022a).Previous observational studies mainly focused on oscillations of long coronal loops with a length of several hundred megameters, and used data from the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) on board the Solar Dynamics Observatory (SDO; Pesnell et al. 2012).Their typical periods range from 1.5 to 10 minutes with an average displacement amplitude of ∼0.17 Mm (Anfinogentov et al. 2015).Recently, decayless oscillations have also been found in short coronal loops (50 Mm) with the high-resolution 174 A  images taken by the Extreme Ultraviolet Imager (EUI; Rochus et al. 2020) on board the Solar Orbiter (SolO; Müller et al. 2020).Petrova et al. (2023) reported high-frequency oscillations in two coronal loops with lengths of 4.5 Mm and 11 Mm located in a quiet-Sun region.The measured periods are 14 and 30 s, much shorter than in longer coronal loops.Zhong et al. (2022b) also identified such oscillations with an average period of 1.6 minutes in a 51 Mm long coronal loop using simultaneous observations from EUI and AIA.Furthermore, Li & Long (2023) conducted a statistical study of decayless short-period oscillations in 111 loops observed by EUI, and obtained a median period of 40 s and loop lengths of 10-30 Mm.The scaling law was found to be similar to that for longer loops, suggesting that all the oscillation events are of the same type, namely the (fundamental) kink eigenmode.Based on that, the authors estimated the magnetic field with the measured oscillation parameters following the seismology method used in previous studies (e.g., Nakariakov & Ofman 2001;Tian et al. 2012;Nisticò et al. 2013;Anfinogentov & Nakariakov 2019).Decayless transverse oscillations have also been detected in coronal bright points (CBPs; Tian et al. 2012;Gao et al. 2022a), which are generally composed of many small magnetic loops (see the review of Madjarska 2019).Gao et al. (2022a) investigated 31 decayless oscillation events in CBPs using EUV imaging observations obtained with AIA.The oscillation period is 1-8 minutes while the loop length ranges from 14 to 42 Mm.Unlike previous observations (Anfinogentov et al. 2015;Li & Long 2023;Petrova et al. 2023), no obvious linear correlation was found between the oscillation period and loop length.Meanwhile, the distribution of the period shows a peak at about 5 minutes (see Figure 3 in Gao et al. 2022a).The most straightforward interpretation is that the transverse oscillations are driven by a periodic external driver such as photospheric p-modes since their period is also around 5 minutes.Note that Figure 6 in Zhong et al. (2022b) also reveals a similar peak in the histogram of periods.
In fact, photospheric p-modes have been suggested to excite transverse waves and oscillations in the corona for a long time.Many observations show that an enhanced power occurs at ∼3 mHz in the coronal velocity power spectrum, implying that the observed transverse waves may be excited by p-modes (e.g., Van Doorsselaere et al. 2008;Tomczyk & McIntosh 2009;Morton et al. 2015Morton et al. , 2016Morton et al. , 2019)).Such a power spectrum has also been used to drive kink oscillations in simulated coronal loops (e.g., Pagano & De Moortel 2019;Howson & De Moortel 2023).Meanwhile, wave excitation by p-modes is also widely studied in a number of analytical and modeling works.Some previous works show that transverse waves can be generated from p-modes through mode conversion (e.g., Khomenko & Cally 2012;Cally 2017), while Riedl et al. (2019) and Skirvin et al. (2023) have found another scenario in which an inclined p-mode driver can directly excite kink oscillations in a flux tube.
In this study, we perform a 3D magnetohydrodynamic (MHD) simulation to study the relationship between p-modes and the excitation of two types of decayless transverse oscillations observed recently in short coronal loops.This paper is organized as follows: Section 2 describes our model and numerical setup.The simulation results are presented in Section 3. In Section 4, we further discuss the results and relate them to observations.Finally, Section 5 summarizes our findings.

Model
In this study, we aim to model a short coronal loop with a length of 30 Mm, which is close to the average loop length found in Gao et al. (2022a).Gravitational stratification and the lower atmosphere (i.e., chromosphere and transition region) are included. 4Therefore, we employ a magnetic flux tube similar to that used in Pelouze et al. (2023) and Guo et al. (2023).The simulation is divided into two steps: a 2D run first working in cylindrical coordinates (r, z) and then a 3D run in Cartesian coordinates.In both runs, the gravity is specified as along the loop, where L = 30 Mm is the loop length and g e = 274 m s −2 is the gravity at the solar surface.Such a distribution has been widely used in previous studies (e.g., Kohutova & Verwichte 2017;Karampelas et al. 2019a;Karampelas & Van Doorsselaere 2020;Riedl et al. 2021;Soler et al. 2021;Pelouze et al. 2023).The temperature profile from the chromosphere to the corona is given by which is derived from Aschwanden & Schrijver (2002).The temperature and height of the chromosphere are T ch = 2 × 10 4 K and z ch = 4 Mm, respectively.The chromosphere is set to be thicker than the real case, following some previous works (e.g., Howson & De Moortel 2022;Guo et al. 2023;Pelouze et al. 2023).Nevertheless, the physical properties of the chromosphere are appropriately described and the thick chromosphere can be regarded as an extended boundary.
The coronal temperature T co has a transverse distribution: Here the temperatures inside and outside the loop are set to be T i = 1.2 MK and T e = 3.6 MK, respectively.R represents the loop radius and b is a dimensionless number determining the thickness of the inhomogeneous layer.Initially, they are set to be 1 Mm and 10 respectively, corresponding to an inhomogeneous layer of 0.6 Mm.Similarly, the density profile at the bottom boundary is set to be where the internal (external) density is set to be ρ i,ch = 3.51 × 10 −8 kg m −3 (ρ e,ch = 1.17 × 10 −8 kg m −3 ).The density at greater heights is set to satisfy the field-aligned hydrostatic equilibrium, as described in Pelouze et al. (2023).The vertical profiles for temperature and density in the initial state are shown in Figure 1(A).
A uniform magnetic field of 42 G along the z-axis is adopted in the whole simulation domain.Given that short coronal loops have a lower height and are mainly distributed in active regions and CBPs, such a magnetic field strength can be reasonable (e.g., see Wang et al. 2007;Tian et al. 2008;Nisticò et al. 2013;Madjarska 2019).
As indicated by Pelouze et al. (2023) and Guo et al. (2023), such an initial setup is not in magnetohydrostatic (MHS) equilibrium.As a result of the nonuniform transverse structuring, it is not possible to obtain an analytical pressure balance equilibrium.Therefore, we relax the system numerically until the background velocities are small enough to be neglected.To save computational time, we first conduct the relaxation in the 2D simulation.To suppress the initial vertical flows in the loop, we multiply the three velocity components by a factor α(z, t) at each time step during the 2D simulation.Similar to Pelouze et al. (2023) The velocity is multiplied by the factor α < 1 between times t 0 and t 1 .In practice, we choose t 0 = 2τ and t 1 = 4500τ, where τ = 7.78 s is the unit of time.After t 1 , we continue to conduct the relaxation for another 1500τ, and succeed in obtaining an equilibrium state.Figure 1(B) illustrates temperature and density profiles along the loop axis after the relaxation process.The horizontal flows have been suppressed to a value below 0.01 km s −1 in the whole simulation domain, which indicates that a new MHS equilibrium has been achieved.The magnetic field inside the loop is now slightly different from that outside (see also Pelouze et al. 2023).Note that the transition region has been broadened to several megameters.It is a compromise hard to avoid since in reality the transition region is very thin and contains very large density and temperature gradients.Considering this, artificially broadening the transition region has become a widely used approach that allows a coarser resolution in the vertical direction (e.g., Lionello et al. 2009;Mikić et al. 2013).
We then rotate the 2D relaxed results to 3D Cartesian coordinates as the initial state for the next 3D MHD simulations, as shown in Figure 1(C).To ensure that there are no unphysical effects arising from the coordinate transition, we then let the 3D system evolve for another 300τ.After that, an inclined p-mode-like driver is employed at one footpoint (z = −15 Mm) to study waves excited in the loop.It is described by The driver has a velocity amplitude (V 0 ) of 500 m s −1 , a period (P) of 300 s, and an inclination angle (θ) of 15°relative to the loop axis.It is expected that the inclination can excite transverse Alfvénic motions/oscillations in the loop by breaking the azimuthal symmetry of the system (see Riedl et al. 2019;Skirvin et al. 2023).To ensure that the driver is located inside the loop, we multiply it with a Gaussian function of the form and a standard deviation σ = 2 Mm.Note that the loop radius has increased to about 2 Mm after relaxation.For both the 2D and 3D runs, we employ the PLUTO code (Mignone et al. 2007) to solve the ideal MHD equations.For spatial reconstruction, a piecewise total variation diminishing linear scheme is utilized, while numerical fluxes are computed using the Roe Riemann solver.A second-order characteristic tracing method is applied for time stepping.Anisotropic thermal conduction is also incorporated into our simulation while explicit resistivity and viscosity are not included.To maintain the divergence-free nature of the magnetic field, we employ the hyperbolic divergence cleaning method (e.g., Dedner et al. 2002).In the 2D case, the computational domain is [0, 6] Mm × [-15, 15] Mm, with a uniform grid of 128 × 512 cells, resulting in a resolution of 47 km in the rdirection and 59 km in the z-direction.In the 3D simulation, we use a uniform grid of 256 cells from −6.0 to 6.0 Mm in both the x-and y-directions, resulting in a resolution of 47 km.In the z-direction, 128 grid cells are adopted from [−7.5, 7.5] Mm and 512 grid points are uniformly distributed in the two left domains near the loop ends, namely [−15.0,−7.5] Mm and [7.5, 15.0] Mm.This ensures a higher resolution near the two footpoints to reveal the dynamics of the lower atmosphere.
The boundary conditions are described as follows.During the 2D relaxation process, we adopted an axisymmetric boundary condition at r = 0 (the loop axis) and an outflow boundary condition at r = 6 Mm.At two footpoints of the loop (z = ±15 Mm), we extrapolated the density and pressure from the hydrostatic equilibrium.Meanwhile, we followed the method used by Karampelas et al. (2019a) to set the magnetic field, which can ensure that the normal gradient of magnetic field equals zero.Additionally, the radial velocity v r and vertical velocity v z are set to be reflective.For the 3D simulation, we used outflow conditions for all the lateral boundaries, and similar conditions to the 2D case for the vertical boundaries.We only changed the velocity conditions according to Equation (5) at one footpoint (z = −15 Mm) to mimic the p-mode driver.

Results
In this section, we will present our numerical results.Note that a convergence study has been conducted to ensure that the current results are not influenced by the numerical resolution.As shown in previous studies (Riedl et al. 2019;Skirvin et al. 2023), an inclined p-mode driver can excite longitudinal waves, especially slow sausage modes, in the loop.Figure 2 presents the evolution of vertical/longitudinal velocity v z .The z-t distribution of v z is shown in Figure 2(A).Seen from the slope, the vertical velocity perturbation propagates along the loop with a speed of about 90 km s −1 in the corona, which is close to the local sound speed c s .The propagation is slower in the lower atmosphere, resulting from a smaller c s there (see also Figure 7(B) and the relevant discussion in Section 4.4).After about 500 s, the velocity distribution is affected by the reflection at another footpoint (z = 15 Mm), forming a pattern like the third-harmonic oscillation at the end of the simulation.Figures 2(B)-(D) show the time series for v z at the loop apex (z = 0) and the associated wavelet analysis results.It takes about 240 s for the wave signal to reach the apex, then the plasma oscillates there with a period of ∼280 s, which is close to the period of the driver.
In this study, we are more interested in the transverse waves/ oscillations excited in the loop.As Figure 3(A) illustrates, transverse oscillations are quickly excited and can reach the loop apex with a much higher propagation speed than slow waves (see also Section 4.4).This can also be seen in the time series of v x at the loop apex (Figure 3  driver, which means that the long-period oscillation is directly driven by the p-modes similar to the longitudinal motions.Specifically, it can be regarded as a propagating kink wave with a long wavelength that is longer than the loop length. 5As for the higher-frequency oscillations, we suggest that they correspond to the eigenmode kink oscillations of the loop, where the two different periods detected here should be related to different harmonics.Given that the wavelet power is derived from the transverse velocity at the apex, we expect to find only odd harmonic signals here.As a result, the period of 76.3 s should be related to the fundamental mode, while 41.6 s likely corresponds to the third harmonic.Interestingly, the global wavelet power at 41.6 s is much higher than that at 76.3 s, which is different from the expected case.This point will be further discussed in Section 4.3. Besides the oscillation period, the transverse velocity amplitude at the loop apex is also investigated.The apparent amplitude of the long-period oscillation at the loop apex is about 0.06 km s −1 at the beginning (before t ∼ 600 s) and then decreases to about 0.04 km s −1 .From the wavelet power spectrum (Figure 3(C)), the 5 minutes oscillation has a nearly constant power, indicating a non-decaying nature.In fact, it is as expected since the footpoint driver persistently offers energy to this periodicity.Therefore, we suppose that, from the perspective of observation, such an oscillation would appear like the decayless oscillation, although it is actually a driven motion.On the other hand, the short-period oscillations seem to experience a damping with time, which could be a result of resonant absorption (e.g., Goossens et al. 2011).However, as seen in Figure 3(B), a high-frequency oscillation without apparent damping can exist for more than ten cycles.As a comparison, decaying oscillations always damp quickly in less than six cycles (e.g., Nakariakov et al. 1999;Goddard et al. 2016).Hence, it is reasonable to regard both high-and low-frequency oscillations as decayless transverse oscillations in an observational sense.
However, the velocity amplitude in our simulation is smaller than in previous observations (e.g., Gao et al. 2022a).This could be because the simulation driver is different from the real case.Considering a linear scaling between the amplitude of the excited oscillation and that of the driver (see, e.g., Poedts et al. 1990;Karampelas et al. 2019b), we can expect that a larger V 0 or θ in Equation (5) can increase the transverse amplitude in the loop.Also, in such a closed system, the pileup of energy will gradually increase the amplitude (see, e.g., Karampelas et al. 2019a;Guo et al. 2019b).In addition, another possible reason for the small amplitude could be related to the cutoff in the transition region, which will be discussed in detail in Section 4.4.
In Figure 3, we choose the x-component of velocity (v x ) to illustrate the transverse oscillations.However, we need to be careful because sausage mode or fluting mode is also accompanied by the fluctuation of horizontal (radial) velocity.Therefore, in order to prove the existence of transverse kink mode with an azimuthal wavenumber m = 1, we need to further confirm whether the axisymmetry of the loop is broken (see Skirvin et al. 2023).In Figure 4(A), we plot v x and v y along the loop axis at several different time steps.We notice that v y , the horizontal velocity component in the direction perpendicular to the driver, is much smaller than v x , thus showing a non-axisymmetric property.We also determine the displacement of the loop by calculating the position of the center of mass.The time series of displacement at the loop apex (z = 0) is shown in Figure 4(B).Roughly speaking, the oscillation of displacement has a similar pattern to that of v x in Figure 3(B), but contains less detail.Nevertheless, the coexistence of two different period regimes can also be observed.In Figures 4(C) and (D), we plot two snapshots of the cross-section density profile (normalized) at the loop apex, with the horizontal velocity field overplotted.The distributions of the velocity vectors are consistent with the typical characteristics of the m = 1 kink mode, i.e., the "dipole-like" pattern (e.g., see Goossens et al. 2014;Guo et al. 2020;Skirvin et al. 2022).Based on these results, we can conclude that there are dominant transverse kink oscillations excited in the loop.

Two Types of Decayless Transverse Oscillations in Short
Coronal Loops In our simulation, we found that two types of transverse oscillations are excited by a velocity driver that is reminiscent of an inclined p-mode.The first one has periods of several tens of seconds, suggesting that it is of the same type as the highfrequency oscillations found in recent EUI observations (Li & Long 2023;Petrova et al. 2023).In other words, these oscillations are kink eigenmodes (this fact will be further verified in Section 4.3), and the period is determined by the physical parameters of the loop.Specifically, the period is proportional to L/c ph , where c ph is the phase speed.In the longwavelength limit, c ph is approximately equal to the kink speed c k (Edwin & Roberts 1983), which is further determined by the density and magnetic field inside and outside the loop.In the corona, the kink speed does not vary greatly with time and location, leading to a roughly linear scaling between L and P as observed in Anfinogentov et al. (2015), Goddard et al. (2016), and Li & Long (2023).To the best of our knowledge, our simulation shows for the first time that such oscillations in short coronal loops could be excited by the p-modes.
The second type of oscillation excited in our simulation has a period of about 5 minutes.It can be seen as a propagating wave whose period is modulated by an external driver, namely, the p-mode-like continuous velocity perturbation we add at the footpoint.We suggest that such a scenario may explain the observational results of decayless oscillations of short coronal loops in CBPs (Gao et al. 2022a).Moreover, it is easy to understand why no linear correlation between the loop length (wavelength) and period was found in this observation, because the frequency is not determined by the loop length L.
In fact, the decayless transverse oscillations found in both longer and shorter coronal loops could be a mixture of these two types of oscillations.Zhong et al. (2022b) summarized a number of previous observational results on oscillation periods and loop lengths (Wang et al. 2012;Anfinogentov et al. 2013Anfinogentov et al. , 2015;;Nisticò et al. 2013;Duckenfield et al. 2018;Anfinogentov & Nakariakov 2019;Mandal et al. 2021;Zhong et al. 2022a;Petrova et al. 2023), and made a scatter plot between these two parameters, as well as a histogram of the period.There is not only a linear scaling in the scatter plot, but also an obvious peak at about 300 s in the period distribution.Additionally, the scatter relative to the best linear fit between P and L is the largest at around 4-5 minutes.This may imply that the 5 minutes transverse oscillations related to the p-mode leakage are likely to occur in those observational studies, although these oscillations were all considered as the kink eigenmode.Furthermore, we note that, in reality, p-modes are known to be stochastically driven and have different periods, leading to a period distribution with a peak at ∼5 minutes, as illustrated in Tomczyk & McIntosh (2009, see their Figure 2).In our model, we just adopt a monoperiodic driver with a 5 minutes period for simplicity.However, if we consider the fact that p-modes can have varying periods, we could have a broadband period distribution of the driven transverse oscillations, which can be even closer to the observational results (e.g., Gao et al. 2022a;Zhong et al. 2022b).
In observations, the coexistence of two types of transverse oscillations in one oscillating loop as shown in Figure 3 has not been found yet.This could be attributed to the limitation of instruments.Most of the previous observational studies analyzed data from SDO/AIA (Wang et al. 2012;Anfinogentov et al. 2013Anfinogentov et al. , 2015;;Nisticò et al. 2013;Duckenfield et al. 2018;Anfinogentov & Nakariakov 2019;Mandal et al. 2021;Gao et al. 2022a;Zhong et al. 2022a).However, AIA's temporal resolution (12 s) of EUV bands limits our ability to observe high-frequency oscillations.This could explain why Gao et al. (2022a) failed to detect such oscillations in CBPs.On the other hand, the data from the High Resolution Imager (HRI) of EUI have both higher spatial and higher temporal resolution than AIA, allowing us to study short-period oscillations in shorter coronal loops (Zhong et al. 2022b;Li & Long 2023;Petrova et al. 2023).But unfortunately, EUI/HRI has only a limited field of view, and the continuous observation time is only several hours at most.Meanwhile, the shorter loops observed so far tend to be more dynamic, and most oscillation events last for just 2-4 cycles.Therefore, it is difficult to detect signals of the 5 minutes oscillation.Nevertheless, future EUI observations could still hopefully capture some short coronal loops that can exist stably for a long duration, and thus the simultaneous existence of these two types of oscillations might be observed.

Excitation Mechanism of the Decayless Oscillation
So far, the excitation mechanism of ubiquitous decayless oscillations is still under debate.Unlike decaying oscillations, which are excited by external pulses, decayless oscillations can persist for multiple cycles while retaining the oscillation amplitude.Previous studies have confirmed mechanisms that are responsible for the damping of kink oscillations in the coronal loop, such as resonant absorption (e.g., Goossens et al. 2011;Guo et al. 2020), formation of Kelvin-Helmholtz instability (e.g., Heyvaerts & Priest 1983;Terradas et al. 2008), andturbulence (Van Doorsselaere et al. 2021).These mechanisms may also be present in decayless oscillations, despite their lower transverse amplitudes (see, e.g., Antolin et al. 2016;Guo et al. 2019aGuo et al. , 2019b;;Karampelas et al. 2019aKarampelas et al. , 2019b;;Shi et al. 2021).In this case, however, a persistent external energy input is required to keep the amplitude from decaying.According to different types of drivers, current excitation models can be divided into three categories (e.g., see Nakariakov et al. 2021;Gao et al. 2022b): (1) a simple harmonic driver with a single frequency at the footpoint of the loop (e.g., Nisticò et al. 2013;Yuan et al. 2023); (2) an external quasi-steady flow such as supergranulation motions (also named the self-oscillatory model; see Nakariakov et al. 2016;Karampelas & Van Doorsselaere 2020, 2021); (3) a random driver with a broadband frequency (e.g., Afanasyev et al. 2020).
In this study, we use a model inspired by Pelouze et al. (2023) and Guo et al. (2023) but to study decayless oscillations driven by p-modes that are inclined to the vertical axis of the loop structure.Previous 3D MHD simulations often added a transverse driver with a period matching the eigenfrequency (e.g., Guo et al. 2019aGuo et al. , 2019bGuo et al. , 2023;;Karampelas et al. 2019aKarampelas et al. , 2019b)).However, our results show that a 5 minutes driver can also excite kink eigenmodes.To some extent, this mechanism may be similar to the self-oscillatory model, since the driver has a much longer period than that of the kink eigenmode in the short loop, and thus can be considered as a quasi-steady flow.Meanwhile, the driver can also excite another oscillation regime with a similar period, which distinguishes it from the typical self-oscillatory model.In any case, the oscillation-generating scenario in this study can be regarded as a new excitation mechanism for decayless oscillations, especially in short coronal loops, which is also consistent with recent observations.The link of widespread decayless oscillations to p-modes implies a potential energy source for wave heating of the corona.

Harmonics of the Eigenmode Kink Oscillation
Researchers are always interested in the harmonic characteristics of kink oscillations in coronal loops, since the period ratio can be used to diagnose the plasma structures in the corona (e.g., Andries et al. 2005;Van Doorsselaere et al. 2007;Duckenfield et al. 2018).Previous modeling works on decayless oscillations have revealed the excitation of high harmonics in loops (Afanasyev et al. 2020;Karampelas & Van Doorsselaere 2020, 2021).In this study, wavelet analysis of v x at the loop apex clearly revealed the existence of different harmonics of the eigenmode kink oscillation (see Figure 3(C)), and they are characterized as the fundamental mode with a period of 76.3 s and the third harmonic with a period of 41.6 s.
Here, Figure 5 demonstrates the presence of the second harmonic.We calculate the global wavelet power spectrum at all positions along the loop axis, and the results are shown in Figure 5(A).To highlight the harmonic signals, we take the logarithm of the power.For the vertical axis relating to the period, we only present a range of 30-90 s, corresponding to the range where the kink eigenmode periods lie.From the figure, we can see the second-harmonic oscillation with a period of 54.0 s has also been excited besides the fundamental mode and the third harmonic.The global wavelet power at 54.0 s for different z is also shown in panel (B).The location of the peak power can be determined as z = ±11 Mm.The time series for v x here are also analyzed, and the results are shown in Figures 5(C)-(E).The wavelet analysis revealed a new peak at about 54.0 s corresponding to the second harmonic, which is not present in Figure 3(D).In summary, we can conclude that there are at least three harmonics excited in our simulation.
Considering the density stratification and lower parts of the atmosphere in our model, it is difficult to analytically calculate the period of kink eigenmodes in the loop.Therefore, we choose to numerically determine this period by running our simulation with an initial velocity impulse transverse to the loop rather than a continuous footpoint driver, to excite kink eigenmodes (see, e.g., Guo et al. 2023).Specifically, we consider three different initial velocity perturbations.All of them are in the x-direction and restricted inside the loop, with an initial distribution given by v n z L sin 1 2 0 [ ( )] p + where v 0 = 5 km s −1 .The distribution has the form of standing waves with n antinodes, and thus the velocity perturbations can excite transverse kink oscillations of the nth harmonic.We conduct simulations for n = 1, 2, and 3, and successfully obtain the eigenmode periods for these three harmonics, i.e., P 1 = 71.5 s, P 2 = 53.7 s, and P 3 = 41.7 s.They are all close to the periods seen in the case of the p-mode-like driver (P 1 = 76.3 s, P 2 = 54.0 s, and P 3 = 41.6 s).Thus, we are convinced that our interpretation of the harmonic periods illustrated in Figures 3 and 5 is correct.
We note that these periods of different harmonics deviate significantly from the results under ideal conditions (P 1 /P 2 = 2 and P 1 /P 3 = 3).In fact, we have P 1 /P 2 = 1.33 and P 1 /P 3 = 1.72.Such a deviation can be interpreted as a consequence of a strong density stratification (e.g., Andries et al. 2005;Erdélyi & Verth 2007), and has also been discovered in previous observations (P 1 /P 2 = 1.4 obtained in Duckenfield et al. 2018) and simulations (P 1 /P 2 = 1.7-1.8obtained in Afanasyev et al. 2020;Karampelas & Van Doorsselaere 2020).Our model has a much shorter loop length and includes the lower atmosphere, so a larger deviation is expected.
From Figures 3(C) and (D), we can see one thing that is unusual: the power of the third harmonic is considerably higher than that of the fundamental mode.This is not only inconsistent with the traditional harmonic wave theory where more power in the fundamental mode is expected, but also different from the results in previous studies (Duckenfield et al. 2018;Afanasyev et al. 2020;Karampelas & Van Doorsselaere 2020).The exact cause of this result is unknown, but it might be associated with the longitudinal waves also produced by the driver.To demonstrate this, we run another simulation with only the driver's inclination angle θ changing from 15°to 90°.In other words, now the driver is fully horizontal, and longitudinal perturbations are not excited in the loop.In this case, we make a plot similar to Figure 3(D) and compare them in Figure 6.It is shown that for the case where θ = 90°, the fundamental mode has a larger power than the third harmonic, which is expected.Therefore, perhaps the longitudinal motions introduced by p-modes can influence the power of higher harmonics.However, the specific reason remains to be investigated in the future and is beyond the scope of this study.

Propagation and Cutoff of the Transverse Motions
In Figure 7(A), we plot the Alfvén speed c A and sound speed c s along the loop axis.The Alfvén speed reaches its peak value of about 1000 km s −1 at the height of the loop apex, much higher than the sound speed (90 km s −1 ).In Figures 7(B) and (C), some curves corresponding to c A (blue dotted lines) and c s (green dashed lines) are overplotted in space-time plots of v z and v x .The slope of these curves corresponds to the local Alfvén/sound speed.The propagation property for the excited longitudinal wave can be clearly seen in Figure 7(B), where the curves  corresponding to the sound speed have a good consistency with the longitudinal perturbations.This is as expected because the longitudinal wave is interpreted as the slow mode, which has a phase speed close to c s .When it comes to the propagation of transverse motions (v x ), we have a more complicated pattern due to the coexistence of short-period eigenmodes and long-period driven oscillations.Nevertheless, we can notice that the establishment of transverse oscillations in the whole loop is much faster than the sound speed.Some large-amplitude v x signals also have a good correspondence with the Alfvén speed.This is easy to understand since the transverse motions driven directly by the driver will travel along the loop with a speed close to c A (e.g., Magyar & Van Doorsselaere 2018).
The stratification in the transition region could influence the propagation of transverse motions along the loop.In Pelouze et al. (2023), the existence of the cutoff of kink waves has been confirmed.They claimed that kink waves with periods longer than the cutoff can reach the corona but with significant damping in amplitude.Although the physical parameters of our shorter loop are different, a similar scenario with a cutoff of about 100 s is still expected.In Figure 3(A), we can see the decrease in transverse velocity v x along the loop.The attenuation is especially significant in the lower atmosphere, as shown in Figure 4(A).Nevertheless, the kink waves can still be analyzed as mentioned above.To further confirm the existence of the cutoff, we conduct a new run with the same setup but a different driver that has a frequency higher than the cutoff.In this case, the driven kink waves are not subject to the cutoff effect.The result shows that the amplitude increases with height and reaches its peak at z = 0, rather than experiencing any attenuation.Thus, the long-period transverse motion in our simulation is indeed damped by the cutoff, while the standing kink waves with short periods seem to be less affected.
One limitation of this model is that the broadened transition region would reduce the reflection of kink waves compared to the real atmosphere.As revealed in Magyar et al. (2019), the vertical gradient in Alfvén speed determines the rate of reflection.A steeper transition region may thus prevent waves from propagating upward toward the corona.The observed presence of a 5 minutes peak in the coronal wave power spectrum (e.g., Tomczyk & McIntosh 2009;Morton et al. 2019;Gao et al. 2022a (3) waveguiding effects in the surroundings of structures (S.J. Skirvin et al. 2023, in preparation).The influence of these mechanisms on the wave driving requires further investigation.

Summary
In this study, we modeled a short coronal loop with two footpoints located in the lower atmosphere.An inclined p-modelike driver is employed at one footpoint to excite transverse oscillations in the loop.It is found that two types of decayless oscillations can be excited: the first has a period of about 5 minutes, close to the period of the driver, and the second has shorter periods of several tens of seconds.We interpreted the former as a propagating wave directly driven by the p-modes, and the latter as kink eigenmodes of different harmonics.Our results are consistent with recent observations of decayless oscillations in short coronal loops, including the high-frequency ones found with SolO/EUI data (Zhong et al. 2022b;Li & Long 2023;Petrova et al. 2023) and relatively low-frequency ones found inside CBPs with SDO/AIA (Gao et al. 2022a).
The photospheric p-modes have previously been believed to play a key role in the generation of propagating kink waves in the corona (Van Doorsselaere et al. 2008;Tomczyk & McIntosh 2009;Morton et al. 2015Morton et al. , 2016Morton et al. , 2019)), and our findings here suggest a possible relationship between these p-modes and decayless oscillations widely distributed in closed coronal magnetic structures.If a considerable amount of the energy from p-modes can be transported to the corona, it could contribute to balancing the radiative loss and could accelerate the solar wind.Given this, the energy analysis in this scenario will be the focus of our future work.Additionally, the decayless transverse oscillation has been seen as a useful seismological tool to diagnose physical properties in the corona (e.g., Tian et al. 2012;Nisticò et al. 2013;Anfinogentov & Nakariakov 2019;Li & Long 2023).In this regard, our results suggest that one needs to be careful when performing such coronal seismology, because some oscillations may be not kink eigenmodes but driven motions.In that case, the standard seismological scheme is not valid.
Finally, it is of great interest to observationally detect the coexistence of these two types of decayless oscillations in the future.Despite some limitations (as discussed in Section 4.1), with the unprecedented high spatial and temporal resolution of SolO/EUI, there is still the potential to observe such a phenomenon and thus further verify our conclusions.

Figure 1 .
Figure 1.The top two panels present the profiles of temperature (black) and number density (blue) along the loop axis (z-direction) before (A) and after (B) the 2D relaxation.The horizontal axis "Height" is actually the increase in z with respect to the footpoint.The solid and dashed lines represent the values inside and outside the loop.Panel (C) shows the 3D loop configuration.
(B)), as the wave signals arrive there much earlier than the longitudinal motions as shown in Figure 2(B).Meanwhile, the oscillations display both long-and short-period regimes.From the wavelet spectrum shown in Figures 3(C) and (D), we can find that there are actually three different periods, 305.2 s, 76.3 s, and 41.6 s.The longest period is close to 5 minutes, namely, the period of the

Figure 2 .
Figure 2. (A) Space-time plot for longitudinal velocity perturbation v z at the loop axis.(B) Time series of v z at the loop apex (z = 0), corresponding to the dashed black line in (A).(C) and (D) Wavelet analysis results of the time series shown in (B).A darker color represents a larger wavelet power.The dashed lines denote a 95% significance level.

Figure 3 .
Figure 3. Similar to Figure 2 but for the transverse velocity perturbation v x .

Figure 4 .
Figure 4. (A) Horizontal velocity along the loop axis at several time steps, with solid lines representing v x and dashed lines representing v y (multiplied by 10 because the value of v y is very low).(B) Transverse displacement of the loop obtained from considering the center of mass.(C) and (D) The cross-section profiles of the normalized density at z = 0 (loop apex) at t = 980.9s and t = 1121.0s, with the horizontal velocity marked by white arrows.

Figure 5 .
Figure 5. (A) The global wavelet power derived from time series of v x at different positions (z) along the loop axis.(B) The global wavelet power at the period of 54 s along the loop axis.(C)-(E) Similar to Figures 3(B)-(D) but for v x at the loop leg (z = 11 Mm) rather than the loop apex.
) could be explained by a number of alternative mechanisms: (1) mode conversion (e.g., Cally & Goossens 2008; Khomenko & Cally 2012), in which the leakage rate of the 5 minutes p-mode is not determined by the gradient in Alfvén speed; (2) reduced kink mode cutoffs in regions of inclined magnetic field (similar to slow waves, De Pontieu et al. 2005);

Figure
Figure (A) The Alfvén speed c A and sound speed c s along the loop axis.(B) and (C) Similar to panel (A) of Figures 2 and 3, but we use blue dotted lines and green dashed lines to show the Alfvén speed and sound speed, respectively.