Amplitudes and energy fluxes of simulated decayless kink oscillations

Recent observations with the Atmospheric Imaging Assembly (AIA) instrument on the SDO spacecraft have revealed the existence of decayless coronal kink oscillations. These transverse oscillations are not connected to any external phenomena like flares or coronal mass ejections, and show significantly lower amplitudes than the externally excited decaying oscillations. Numerical studies have managed to reproduce such decayless oscillations in the form of footpoint driven standing waves in coronal loops, and to treat them as a possible mechanism for wave heating of the solar corona. Our aim is to investigate the correlation between the observed amplitudes of the oscillations and input the energy flux from different drivers. We perform 3D MHD simulations in single, straight, density-enhanced coronal flux tubes for different drivers, in the presence of gravity. Synthetic images at different spectral lines are constructed with the use of the FoMo code. The development of the Kelvin-Helmholtz instability leads to mixing of plasma between the flux tube and the hot corona. Once the KHI is fully developed, the amplitudes of the decayless oscillations show only a weak correlation with the driver strength. We find that low amplitude decayless kink oscillations may correspond to significant energy fluxes of the order of the radiative losses for the Quiet Sun. A clear correlation between the input energy flux and the observed amplitudes from our synthetic imaging data cannot be established. Stronger drivers lead to higher vales of the line width estimated energy fluxes. Finally, estimations of the energy fluxes by spectroscopic data are affected by the LOS angle, favoring combined analysis of imaging and spectroscopic data for single oscillating loops.


INTRODUCTION
Over the past twenty years, observations of the Sun have shown the existence of waves and oscillations throughout the solar corona (Aschwanden, 2006;De Moortel and Nakariakov, 2012). The discovery of transverse magnetohydrodynamic (MHD) standing (Aschwanden et al., 1999;Nakariakov et al., 1999) waves in coronal loops, and propagating waves in open magnetic field structures (Verwichte et al., 2005) has lead to many observational and numerical studies. The ubiquity of such waves has also been established in prominence threads (Okamoto et al., 2007), coronal loops (McIntosh et al., 2011), as well as greater areas of the corona (Tomczyk et al., 2007;Tomczyk and McIntosh, 2009;Thurgood et al., 2014;Morton et al., 2016), renewing the interest on the effects of these waves in the solar atmosphere.
Analytical studies on the nature of the transverse oscillations in inhomogeneous plasmas (Zajtsev and Stepanov, 1975;Ryutov and Ryutova, 1976;Edwin and Roberts, 1983;Allcock and Erdélyi, 2017) have described the different surface waves expected in a non-uniform plasma. In order to explain the observed damping of such oscillations (Tomczyk and McIntosh, 2009;Verth et al., 2010;Terradas et al., 2010;Pascoe et al., 2016aPascoe et al., , 2017Pascoe et al., 2019), extensive theoretical and numerical work has also been performed. The mechanisms of resonant absorption and mode coupling (Sakurai et al., 1991;Goossens et al., 1992;Ruderman and Roberts, 2002;Arregui et al., 2005;Goossens et al., 2011;Pascoe et al., 2012Pascoe et al., , 2016bDe Moortel et al., 2016;Yu et al., 2017;Pascoe et al., 2018) are considered the reason behind this spatial and temporal attenuation of the oscillations, by transferring the energy of the global kink mode to local azimuthal Alfvén modes. Through phase mixing (Heyvaerts and Priest, 1983;Soler and Terradas, 2015), the energy is then transferred to ever decreasing smaller scales until it gets dissipated by resistivity and viscosity (Ofman et al., 1994a,b;Poedts and Boynton, 1996;Ofman et al., 1998). Observational studies of waves in the solar chromosphere and corona (De Pontieu et al., 2007;Tomczyk and McIntosh, 2009;Morton et al., 2012) suggest the existence of enough energy flux to sustain the radiative losses of ∼ 100 W m −2 for the non-active region corona (Withbroe and Noyes, 1977). In De Moortel and  it was shown that LOS integration of footpoint driven multistrand coronal loop oscillations leads to an underestimation of the wave energy. In Antolin et al. (2017) it was found for a non-driven oscillating loop that the wave energy ends up underestimated through the localisation of the energy by resonant absorption. Recent simulations of coronal loop waves (Pagano and De Moortel, 2017;Pagano et al., 2018;Pagano and De Moortel, 2019) have not reported sufficient heating to balance the radiative losses.
Alongside the high amplitude, externally initiated, decaying transverse oscillations in coronal loops (Verwichte et al., 2009(Verwichte et al., , 2010 a new group of smallamplitude decayless transverse oscillations have been identified in coronal loops (Nisticò et al., 2013;Anfinogentov et al., 2015;Duckenfield et al., 2018), with amplitudes ∼ 0.1-0.4 Mm. These decayless oscillations have been interpreted in different ways over the years. They have been treated as continuously driven kink waves with a footpoint driver Karampelas et al., 2019;Afanasyev et al., 2019), as a self-oscillatory process due to the interaction of the loops with quasi-steady flows (Nakariakov et al., 2016), or as a line of sight (LOS) effect from the development of the Kelvin-Helmholtz instability (KHI)  in impulsive standing loop oscillations.
The development of the KH instability has been theorized in plasma structures were standing surface waves are observed (Heyvaerts and Priest, 1983;Browning and Priest, 1984;Zaqarashvili et al., 2015;Barbulescu et al., 2019;Hillier et al., 2019), caused by the strong shear velocities generated by the azimuthal Alfvén waves. Recent numerical studies (Terradas et al., 2008;Antolin et al., 2014;Magyar and Van Doorsselaere, 2016;Howson et al., 2017;Terradas et al., 2018;Karampelas et al., 2017;Antolin et al., 2018) have confirmed the development of transverse wave induced Kelvin-Helmholtz (TWIKH) rolls for standing kink waves. Additional work has been performed in order to develop methods of identifying the effects of KHI in oscillating loops Van Doorsselaere et al., 2018). Spatially extended TWIKH rolls have been found in the case of continuously driven standing waves, which fully deform the initial monolithic loop cross-section into a fully turbulent one Karampelas et al., 2019). These spatially extended TWIKH rolls in simulations of continuously driven loops have also been reported as sites of mixing of plasma and heating in the solar corona, by effectively spreading the effects of phase mixing across the cross-section of loops (Karampelas et al., 2017;Guo et al., 2019;Karampelas et al., 2019;Afanasyev et al., 2019).
One of the main challenges in creating an efficient wave heating model is to find a way of providing a high enough energy flux over large periods of time, while still remaining within the constraints imposed by the observational data. The decayless oscillations seem to be ideal candidates for prolonged energy input in such a model, and also pose the constraint of having relatively small oscillation amplitudes. In Guo et al. (2019), the use of a mixed Alfvén and transverse footpoint driver suggested as a way of increasing the total energy input while retaining small amplitudes in the synthetic images. The authors in Afanasyev et al. (2019) suggest the use of broadband drivers as a way to include additional energy for a decreased observational impact, while still mantaining sufficient energy input. In the current work, we study the effects of transverse footpoint drivers of different strengths on the observational signatures of the induced oscillations in synthetic images. The aim is to detect the existence of a correlation between the input energy flux and the oscillation amplitudes and to determine the observed fluxes from synthetic spectroscopic and imaging data, and set possible constraints for future work.

Basic setup
For our 3D simulations, we use straight, density-enhanced magnetic flux tubes in a low-β coronal environment, following Karampelas et al. (2019). This setup models gravitationally stratified, active region coronal loops in ideal MHD, in the presence of numerical resistivity. Each loop has a full length (L) of 200 Mm and an initial minor radius (R) of 1 Mm, which is constant with height. The radial density profile of our cylindrical flux tubes at the footpoint is We denote the basic values of our physical parameters with the index i (e) for internal (external) values, with respect to our tube. The external and internal density at the footpoint are equal to ρ e = 10 9 µ m p cm −3 = 0.836 × 10 −12 kg m −3 (µ = 0.5 and m p is the proton mass) and ρ i = 3 × ρ e . We denote the coordinates in the plane perpendicular to the loop axis by x and y , and by z the coordinate along its axis. b = 20 sets the width of the boundary layer to ∼ 0.3R. For the models studied, we set the temperature to be constant with height and equal to T e = 3 × T i = 2.7 MK. Finally, we consider an initial uniform magnetic field parallel to the flux tube axis (along the z-axis) equal to B z = 22.8 G. This gives us a plasma β = 0.018.
Gravity varies sinusoidally along the flux tube, corresponding to a semi-circular loop with major radius Lπ −1 , and takes a zero value at the loop apex (z = 0) and a maximum absolute value at the footpoints (z = ±100 Mm). Due to the hydrostatic equilibrium we have stratification of pressure and temperature along the loop. To counteract the initial pressure imbalance at the loop boundary caused by the gravitational stratification, the magnetic field is restructuring inside the flux tube, once we let the system reach a semiequillibrium state after one period. By the end of the relaxation period, neither of temperature, pressure, nor density deviate significantly from their initial state, as in Karampelas et al. (2019). The radial profiles of density and temperature at several positions along the z-axis are shown in Fig. 1.

Boundary conditions and driver
Our tubes are driven from the footpoint (z = 100 Mm), using a continuous, monoperiodic 'dipole-like' driver (Karampelas et al., 2017), inspired by that used by Pascoe et al. (2010). The period of the driver is P 2L/c k , coinciding with the corresponding fundamental eigenfrequency for our gravitationally stratified flux tube (Edwin and Roberts, 1983;Andries et al., 2005). For a kink speed of c k = 1167 km s −1 we obtain a period of P = 171 s. The driver at the bottom boundary has a uniform and time varying velocity inside the loop, where υ 0 (in km s −1 ) is the peak velocity amplitude. Outside the loop, the velocity follows the relation where α(t) = υ 0 (0.5 P/π) sin(2πt/P ) is a function that recentres the driver, tracking the footpoint. A transition region following the density profile exists between the two areas, in order to avoid any numerical instabilities due to jumps in the velocity.
At the bottom boundary, we also keep the velocity component parallel to the z-axis (v z ) antisymmetric (z = 100 Mm) and we extrapolate the values for density and pressure, using the equations for hydrostatic equilibrium. A zero normal gradient condition is used to extrapolate the values of each magnetic field component through the bottom boundary. Studying the fundamental standing kink mode for an oscillating flux tube allows us to take advantage of the inherent symmetries of this mode, as well as the symmetric nature of our driver. In the top boundary at z = 0, the location of the loop apex, we kept v z , B x , and B y antisymmetric in the x − y plane at the apex. The rest of the quantities on that boundary are defined as symmetric. Thus, only half the loop is simulated along the loop axis. All side boundaries are set to outflow (Neumann-type, zero-gradient) conditions for all variables.

Numerical method and grid
The 3D ideal MHD problem is solved using the PLUTO code (Mignone et al., 2018), with the extended GLM method from Dedner et al. (2002) keeping the solenoidal constraint on the magnetic field. We use the finite volume piecewise parabolic method (PPM) and the characteristic tracing method for the timestep. The domain dimensions are (x, y, z) = (16, 6, 100) Mm for models D1, D2 and D4, and (x, y, z) = (16, 10, 100) Mm for models D6 and D8. The latter was chosen in order to properly resolve the expanded turbulent loop cross-section from the stronger drivers. We have cell dimensions of 40×40×1562.5 km for all models. The resolution is higher in the x − y plane, to better resolve the small-scale structures that appear in the loop cross section. The density scale height of our setup is 28 Mm, and is close to the value of the loop major radius Lπ −1 ∼ 64 Mm. Alongside the lack of radiation or thermal conduction, this allows for a coarser resolution on the z-axis, that can still sufficiently model the density stratification along the flux tube. In all of our models, we also have the inevitable numerical dissipation effects, which lead to values of effective resistivity and viscosity many orders of magnitude larger than the expected ones in the solar corona.

Forward modelling
We use the FoMo code  to create synthetic images and compare to real observations. We focus on the use of the Fe XII 193Å emission line, with maximum formation temperature of log T = 6.19, in order to create spectroscopic data of our flux tubes. The 193Å line tracks the warmer flux tube boundary layer  and is better suited for our studies of driven oscillations with turbulent cross-sections . The same analysis has been performed for the Fe IX 171Å line. This line has a maximum formation temperature of log T = 5.93, and is more sensitive to the colder plasma from the loop interior. We also use the 'Imaging'functionality of 'fomo-c'from FoMo, in order to create emission images for the corresponding AIA channels.
We create time-distance maps of the loop apex for each model, by placing a slit at LOS angles equal to 0 • (perpendicular to the axis of the oscillation, or x axis), 45 • and 90 • (parallel to the axis of the oscillation).To better compare to the observations of decayless coronal loop oscillations, we degrade the original spatial resolution of the AIA synthetic images to that of SDO/AIA (0 .6). Then we resample the data to match the pixel size of the target instrument. A similar process is used for the spectroscopic data, where we try to mimic Hinode/EIS by taking a plate-scale of 1 and a Gaussian PSF of 3 FWHM Guo et al., 2019), while also adding poisson distributed photon noise. For the spectroscopic data we also reduce the initial spectral resolution from 3 km s −1 to 36 km s −1 for the degraded spectrum. The temporal resolution is kept at ∼ 11 s, which is very close to both that of AIA and EIS instruments.

RESULTS
We drive our loops for fifteen cycles. As in Karampelas et al. (2017), the first waves to reach the apex (z = 0) are the azimuthal Alfvén waves at the boundary layer of our tube, thanks to their higher propagation speed, followed by the propagating kink waves. Because of the symmetry at the apex, the driver induced waves lead to the initiation of a standing oscillation resembling the fundamental kink mode for gravitationally stratified loops Magyar et al., 2015;Karampelas et al., 2019).
Once the standing mode is initiated, we have the gradual development of the KH instability and the expansion of the flux tube cross-section (Karampelas et al., 2017). The values of the effective numerical dissipation coefficients are low enough not to hinder the development of the instability (see Karampelas Frontiers et al., 2019). By performing forward modeling for our setups, we construct time distance maps of our oscillations at the apex, for the different spectral lines. Two examples of these forward modeling timedistance maps are presented in Fig. 2. There we see the apex displacement over time for the D2 and D4 model in the 171Å and 193Å AIA channels. As expected for loops with turbulent cross-sections , the development of the KH instability and the spatially extended TWIKH rolls lead to extensive mixing of the loop with the surrounding plasma. For our setup of a cold loop embedded in a warmer corona, this leads to a big emission drop in the Fe IX 171Å line.
By saturating the image of the 171Å emission line for model D4, we can see that this drop occurs faster for stronger drivers. A smaller drop is also observed in the 193Å line, for stronger drivers, and is caused by the drop of the average density of the emitting plasma, and its mixing with the surrounding plasma (Karampelas et al., 2017), due to the spatially extended TWIKH rolls. In coronal observations we expect to mostly see the cooling stage of loops (Viall and Klimchuk, 2012). The loops observed performing these decayless observations are expected to be driven for long periods of time, while cooling from higher initial temperatures. Hence, the observed decayless regime should correspond to the later stages of our simulations, where we observe the mixing induced emission drop.
In Fig. 3 we present additional synthetic images for the D2 model, both at numerical and instrument resolution. In that figure, we see the 193Å AIA channel for the emission intensity at 0 • LOS angle, as well as the Doppler velocities for the Fe XII 193Å line at 0 • and 90 • angles. We observe the development of emission peaks prior to the development of the TWIKH rolls, due to the excitation of higher order harmonics and the deformation of the flux tube from the combined effect of the inertia and fluting modes Ruderman et al., 2010;Yuan and Van Doorsselaere, 2016;Antolin et al., 2017;Terradas et al., 2018). Once the KHI develops, we see the emergence of stronger emission peaks in the 193 A channel, stronger Doppler velocities perpendicular to the direction of the oscillation, and the gradual expansion of the resonant flows across the loop cross-section (Antolin et al., 2015;Karampelas and Van Doorsselaere, 2018). The latter is evident in the bottom panels of Fig. 3, where we look at the loop from the side. These out of phase movements from the TWIKH rolls are present in all of our models and get more intense and spatially extended the stronger each driver is.
We track the oscillations in the constructed time-distance maps for our models by setting a threshold value off I thr = I min + 0.1(I max − I min ), where I is the intensity across the loop on each frame, I max is its maximum value and I min is the background intensity. The point with this intensity is defined as the 'edge'of the loop. A spline interpolation for the position is used to inhibit tracking errors, with the number of interpolation points being half the number of frames in the TD map. We present the displacement of the apex for each loop in Fig. 4, where panels for all the models are shown at a 0 • LOS angle. An initial increase of the amplitude is followed by saturation and an eventual drop in the second half of the simulations. This is the stage of the simulations where the decayless regime should correspond to. In the last stages, we see that the D2 and D4 models show similar values of displacement, while the D6 and D8 models show lower amplitudes than the former ones.
We see that similar profiles are obtained for both emission lines used in the present work. The synthetic images in the 171Å line show higher initial amplitudes than those in the 193Å line. This is because the loop core has the highest transverse velocities before the development of the KH instability. The loop boundary layer has radially decreasing υ x velocities, as result of the driver shape that we use at the footpoint, leading to smaller displacements than the loop core, as seen in the two emission lines, each being more sensitive to the corresponding loop region. A similar thing was observed in Antolin et al. (2016). In that work, the loop core showed a higher oscillation amplitude than the loop boundary layer, before the initiation of the KHI.
In our models, the constant driving leads to overall stronger displacements, while the deformation of the loop cross-section leads to higher amplitudes in the 193Å line for the given temperature distribution.
In Fig. 4 we get a small phase difference between the D1 model and all the rest. This small change in the phase between the observed displacement and the driving frequency is a LOS effect and caused by the integration of the emission across TWIKH rolls with varying phases. A similar phase shift was observed for the turbulent layer of oscillating loops performing in Antolin et al. (2016). For the D1 model, the slower development of the KHI leads to a delay in the manifestation of that phase difference. Similar trends have been observed for synthetic images at 45 • , with the difference that the obtained displacements are smaller due to the projection of the oscillation plane on the plane of sky (POS).
In order to determine whether the obtained amplitude drop in Fig. 4 is caused purely by the development of strong azimuthal resonant flows and the KH eddies, we plot in Fig We can see that the efficiency of the drivers change over time, with the strongest changes happening for the stronger drivers. These variations with time are the result of the boundary conditions employed here. The zero normal gradient condition for the magnetic field allowed it to freely evolve with time, instead of fixing its value at the footpoint. This is affecting the Poynting flux from the driver in the bottom boundary, causing the evolution of the energy input, similarly to Guo et al. (2019) and Karampelas et al. (2019). In the scope of this study, we do not focus on creating a realistic model for a driver, but we focus on the effects of the driver strength on the synthetic data. Therefore, tracking the evolution of the driver efficiency is sufficient within the context of the current study. The initial increase in the energy flux is followed by a drop, which is again more pronounced in the models with the stronger drivers. By comparing the curves for the models D2 and D4 with the synthetic images of Fig. 2, we see place that this drop takes place when the loops become turbulent due to the KHI.
Focusing at the last 6 cycles of the the simulation in Fig. 5, we see that D4 shows a stronger input than D6, while still remaining weaker than D8. D2 on the other hand remains weaker from all the aforementioned ones, providing almost half or less as much input flux. This is in contrast to the results of Fig. 4, where the oscillation amplitude of the D2 model was comparable to that of the D4 model, and higher than the respective amplitudes of models D6 and D8. This shows that the previously obtained amplitudes are affected by the development of the out of phase movements of the spatially extended TWIKH rolls. All the drivers provide fluxes within the range of the radiative losses of the non-active region corona (Withbroe and Noyes, 1977), with the exception of D1, which is significantly weaker than the rest.
After calculating the input energy flux and the observed loop displacements, we plot amplitude-flux diagrams in Fig. 6. For these diagrams, we took the average of the (absolute values) of the maximum displacements and fluxes, for each model, at time intervals 19.95 − 28.5 minutes, 28.5 − 34.2 minutes and 34.2 − 42.75 minutes. The diagrams contain the amplitudes traced from AIA-resolution synthetic images of 171 and 193 AIA channels at both 0 • and 45 • LOS angle. Taking this and the results from Fig. 4 into account, we explain the differences between the 171 and 193Å AIA channels as LOS effects of our models. The highest values of the amplitudes are obtained for the highest driver input fluxes. However, a clear correlation is absent for fluxes ≤ 110 W m −2 . Instead, we see a relatively uniform distribution of amplitudes among the different fluxes in that region. The corresponding amplitudes are below 1 Mm as we look perpendicular to the oscillation direction (0 • .) At a LOS angle of 45 • , the value of these amplitudes is between 0.15 and 0.55 Mm, placing them within the range of the observed decayless oscillations (Nisticò et al., 2013;Anfinogentov et al., 2015).
Finally, we want to calculate the observed kinetic energy flux (W m −2 ) across the line of sight, for a specific spectral line. Based on the work of , we calculate the Doppler energy flux as: and the non-thermal energy flux as: With c k = 1167 km s −1 we denote the kink speed. L ⊥LOS is the length of the domain across the LOS, υ Dop,λ the Doppler velocity (in km s −1 ) for the spectral line, ξ λ the corresponding Doppler width in km s −1 , and the thermal width (in km s −1 ): With c we denote the speed of light, λ 0 the wavelength at rest for the spectral line, k B the Boltzmann constant, m p the proton mass and with µ the atomic weight in proton masses of the emitting element. The emitting temperature T λ is taken as an approximation from the ξ λ,th in the beginning of the simulations, where the velocities along the LOS are practically zero. The quantity < ρ λ (l, t) > is based on the emissivity-weighted average density from Antolin et al. (2017): where we instead divide the integral of density times emissivity with the maximum value of the integrated emissivity along the LOS at each frame. The fluxes from all setups are then normalized to the domain dimensions of models D6 and D8 for all setups, to make direct comparisons possible between the different models.
We plot our results for the Fe XII 193Å spectral line in Fig. 7, following Hinode/EIS spatial, spectral and temporal resolution. For a LOS angle of 90 • , we observe an initial increase in the energy flux estimated by the Doppler velocities, which persisted for the duration of the simulations for models D1 and D2. For the models D4, D6 and D8 we observe a drop in the estimated Doppler energy flux, in the last cycles of the simulations. This drop appears once the turbulent loop cross-section shows very strong out of phase flows, due to the TWIKH rolls. These azimuthal flows do not register in the Doppler velocities at this angle, leading to the aforementioned drop. A similar behavior was observed in Antolin et al. (2017) for a non-driven transverse oscillation. We also see a correlation with the drop of the input energy flux, as observed for the the models D4, D6 and D8. However, a similar correlation is not observed for the D1 and D2 models. This shows that the observed drop in the Doppler energy flux is the combined result of the input drop and the loop deformation by the TWIKH rolls.
The diagram of the non-thermal energy flux in Fig. 7 shows a saturation once the KHI fully deforms the loop cross-section. This agrees with the results of Antolin et al. (2017), where a similar saturated domain was obtained after the manifestation of the KHI. For the models D4, D6 and D8, this saturation takes place before the drop in the Doppler energy flux. The small increase of the non-thermal energy flux observed in the later parts of the simulation, does not fully compensate for the drop in the Doppler energy flux. This shows that the aforementioned Doppler flux follows closely the input energy flux, in the case of continuously driven oscillating loops. We also observe that the non-thermal energy fluxes are increasing for stronger drivers. This suggests that the observed correlation between line widths and Doppler shifts reported by McIntosh and De Pontieu (2012) can correspond to a continuous input of transverse MHD waves, supporting their hypothesis. However, the use of a single loop for each driver excludes the LOS effects of integrating over many structures, and cannot give a safe statistical result to be compared with the observations (McIntosh and De Pontieu, 2012).
We see a disagreement between the total observed (Doppler + non-thermal) energy flux, for a LOS angle of 90 • , and the input in Fig.5 for models D1 and D2. The former, as seen in Fig. 8 is overestimated with respect to the input energy flux during the later part of the respective simulations. This comes in contrast with the results of Antolin et al. (2017), where the estimated total energy flux was lower than the available energy flux calculated directly from the simulation data. This disagreement is due to the inaccuracies introduced through the method of estimating the density and temperature of the emitting plasma from the synthetic data. The latter appears in the correction for the thermal terms in the non-thermal energy fluxes. Additional uncertainty can be caused by the angle of observation. As seen in Fig. 8, the total observed energy for a LOS angle of 0 • follows the development of the TWIKH rolls and leads to an underestimation of the energy flux.
The uncertainty revolving around the LOS angle can be bypassed by a combined approach of multiple observations. Once a small angle of observation is identified by combining the energy diagrams and the spectra, like those in Fig. 4, we can then approximate the total energy fluxes. We achieve this by adding an additional term to the sum of the Doppler and non-thermal energy fluxes. This additional term is derived from eq. 4, by replacing the Doppler velocities with the the velocity derived from the displacement of the oscillating loop, like in the time-distance maps of Fig. 4. The results of these calculation are presented in Fig. 9, and show a better agreement with the input energy than the underestimated fluxes in panel (B) of Fig. 8.

DISCUSSION
In the present work we aimed to study the effects of different transverse drivers in the dynamics and energetics of decayless oscillations in coronal loops, through the use of synthetic images. We modeled the decayless oscillations as the standing waves from continuously footpoint-driven straight flux tubes, in the presence of gravity. We performed MHD simulations for footpoint drivers of different strengths and performed forward modeling to our data, constructing time-distance maps at the loop apex. We then degraded the spatial and spectral resolution of our signal, targeting different instruments, namely the SDO/AIA and Hinode/EIS. From the created acquired synthetic images for the 193Å and the 171Å lines we studied the oscillation amplitudes in the plane of sky (POS) at the apex, generated by the five different drivers. To do so, we tracked the 'edge'of the loop from the emission images, which is defined by threshold value for the intensity, set slightly higher than the background intensity.
In Fig. 2 we observe the drop of the emission in synthetic images for both lines employed here. This drop is the result of the spatially extended TWIKH rolls, expanding the loop cross-section and causing mixing with the surrounding plasma. This drop in the emission is stronger and occurs faster for stronger drivers, as we can see in the example of Fig. 2 for the D2 and D4 models. In observations we mostly expect to observe loops during their cooling stage (Viall and Klimchuk, 2012), starting from unknown higher temperatures of formation. For decayless oscillations, we therefore expect them to be already in a turbulent state, because of the continuous driving. Any change in the strength of the driver could then lead to sudden changes in the observed emission in various lines, depending on the temperature gradient between the loop and the surrounding plasma, as well as the temperature of said coronal plasma.
In Fig. 4 we see that the final values of the amplitudes are comparable between the different drivers. We are mostly interested in the second half of the simulations, during which the KHI has deformed the cross-section of each loop. This is because we expect the loops performing decayless oscillations to be already in a turbulent state. During that phase, we can observe higher oscillation amplitudes from weaker drivers. This becomes obvious by comparing these time distance maps and the input energy flux diagrams of Fig. 5 for the D4 and D8 models. The D4 model showed bigger displacements towards the end of the simulations than D8 model, despite having a less efficient driver at that stage of the simulations. In addition to that, the D2 model showed very similar amplitudes to the D4, D6 and D8 models, despite providing less than half the energy flux of these drivers. This result is a clear indication of the importance of the out of phase flows, developed by the TWIKH rolls, in the characteristics of the observed oscillations.
Following the results of Fig. 7, we observe a drop in the Doppler energy flux for the models D4, D6 and D8. This drop is primarily related to the respective drop of the input energy flux, and secondarily to the development of the TWIKH rolls, since the drop in the input energy flux of models D1 and D2 does is not followed by a drop in their Doppler flux. This is related to the less deformed loop cross-sections expected for weaker drivers in continuously driven oscillations of coronal loops . The connection between the Doppler energy flux and the resonant flows has also been established in Antolin et al. (2017), for a small-amplitude non-driven transverse oscillation of coronal loop.
From the diagram of the non-thermal energy flux over time in Fig. 7, we observe the saturation of the energy flux, attributed to the turbulent motions present across the loop cross-section. We also observe a positive correlation between the average values of the saturated non-thermal energy fluxes and the driver efficiency from Fig. 5. This suggests that a continuous input of transverse MHD waves could lead to the observed correlation between line widths and Doppler shifts reported by McIntosh and De Pontieu (2012). However, the small number of measurements (one observation per loop per driver) does not allow us to acquire safe statistical results. Furthermore, our study has excluded the LOS effects of integrating over many structures, which prevents us from making a direct comparison with the observations (McIntosh and De Pontieu, 2012).
By taking into account the results of Fig. 9, we see that the sum of the non-thermal and Doppler energy flux can lead to different values, depending on the angle of observation and the methods used in the calculations. By adding an additional term, derived by estimating the oscillation velocity from the timedistance maps, we can compensate for the effects of the observation angle when calculating the energy fluxes.
We need to stress here that the observed saturation in the amplitudes is not caused by numerical dissipation, but is instead the combined result of the strong out of phase flows generated by the TWIKH rolls and the effects of plasma mixing in the emission. Scaling tests with setups of higher resolution, as well as past results Karampelas et al., 2019) have revealed a similar saturation in the oscillation amplitudes. In Karampelas et al. (2019), we have seen that the input energy from our drivers turns into kinetic energy (strong out of phase flows), magnetic energy and internal energy of the plasma (wave heating). The derived spectra in that study reveal a turbulent profile, represented by the inertial range. This inertial range is connected to the existence of smaller scales and is present as long as the KHI is not suppressed by the very large dissipation parameters. The development of these smaller scales could be observed in the Doppler velocities and non-thermal line-widths from spectroscopic observations. These observations, however, would be dependent upon the angle of observation, the instrumental resolution, as well as the assumptions used for the initial temperature and density profile.
Finally, from Fig. 6 we conclude that we cannot obtain a clear correlation between the driver input and the observed oscillation amplitudes just from the synthetic imaging data. In our setups, the amplitudes between 0.15 and 0.55 Mm were evenly distributed over a wide range of input fluxes, up to 110 W m −2 . Small amplitude oscillations can potentially be hiding enough energy to sustain the non-active region corona (∼ 100 W m −2 ) (Withbroe and Noyes, 1977). The development of spatially extended TWIKH rolls can mask its observational signatures, by causing out of phase movements of the loop plasma and affecting the emission at various lines. This is relevant when explaining the observed decayless oscillations (Nisticò et al., 2013;Anfinogentov et al., 2015) as the result of continuous footpoint driving of a coronal loop, and stresses the need for combined instrument observations.
To summarize our results, the mixing of plasma caused by the developed KHI affects the emission in synthetic images, leading to a drop once the loop cross-section becomes turbulent. The developed TWIKH rolls show strong out of phase flows, which lead to the saturation of the oscillation amplitudes for each model. For different angles, we obtain oscillation amplitudes near the observed ones, which do not show a clear correlation with the input energy from the driver. These observed amplitudes can potentially carry enough energy to sustain the radiative losses from the Quiet Sun. A better correlation can be obtained between the input energy and the spectroscopic results from the Doppler energy flux and the non-thermal energy fluxes. All of our results, however, are dependent upon the angle of observations and the approximations used, hinting towards the need for combined spectroscopic and imaging observations. Apart from our suggestion that small-amplitude decayless oscillations contain enough energy flux to support the QS, there are other ways to introduce additional energy within the context of such observed waves. In De Moortel and  it was shown that the wave energy is significantly underestimated by integration of multiple loops along the LOS. In that work, the kinetic energy estimated by the LOS Doppler velocities was ∼ 5% − 20% of the energy in the domain. A combination of a kink and Alfvén driver  can also provide more energy into the system, while retaining a similar profile in AIA synthetic images. Alternatively, the use of a broadband driver could also increase the driver efficiency, for only small deviations from the expected oscillation amplitudes. The detection of such higher harmonics in decayless oscillations by Duckenfield et al. (2018) suggests that the latter method should be studied further in future works.

CONFLICT OF INTEREST STATEMENT
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

AUTHOR CONTRIBUTIONS
K.K. has performed numerical simulations and forward modeling for the models in this project, and set the basic ideas for the analysis of the results. D.J.P. has performed specialized data extraction and analysis of the synthetic images. T.V.D. has set the basic guidelines for this project, consulted during the analysis and supervised the completion of the project. M.G. and P.A. have contributed on the methods of the analysis and have given feedback during the discussion of the results. K.K. has written the first draft of the current paper, received input from all co-authors, and written the final version of the manuscript.