Excitation and evolution of coronal oscillations in self-consistent 3D radiative MHD simulations of the solar atmosphere

Coronal loops, the basic building blocks of solar corona, are commonly subject to oscillatory behaviour. Due to their magnetic structure they act as wave guides and can support a variety of wave modes (see e.g. reviews by Nakariakov & Verwichte 2005; Nakariakov & Kolotkov 2020). Various types of wave modes have been reported at coronal heights; these include nearly incompressible transverse oscillations in the direction perpendicular to the loop axis (Aschwanden et al. 1999; Nakariakov et al. 1999; Verwichte et al. 2013), compressible longitudinal oscillations along the loop axis (Berghmans & Clette 1999; De Moortel et al. 2000) and incompressible torsional oscillations (Kohutova et al. 2020b). Transverse oscillations of coronal loops are by far the most commonly observed type. They were first discovered in the TRACE observations of coronal loops following a solar flare (Aschwanden et al. 1999; Nakariakov et al. 1999) and since then have been studied extensively through high-resolution solar observations, numerical simulations and analytical theory. The most commonly reported oscillation periods lie in the range of 1 to 10 min, and the oscillation amplitudes range from few hundred km to several Mm. Most of the observed transverse coronal loop oscillations have been identified as corresponding to either standing or propagating kink modes. This mode identification is based on modelling individual coronal loops as cylindrical flux tubes (Edwin & Roberts 1983). It should however be noted that the density structuring in the corona is more complex than this and assumptions about the individual mode properties based on idealised modelling of homogeneous plasma structures do not always apply in a non-homogeneous and dynamic environment that is real solar corona (e.g. Goossens et al. 2019). The fundamental harmonic of standing kink mode is the most commonly observed, where the footpoints of the coronal loop act as nodes and the point of maximum displacement amplitude lies at the loop apex. Despite fundamental harmonic being the most intuitive resonant response of a line-tied magnetic fieldline to an external perturbation, excitation of higher-order harmonics is also perfectly viable; higher order harmonics have also been identified in coronal loop oscillations (Verwichte et al. 2004; Duckenfield et al. 2018).

Transverse oscillations of coronal loops are commonly used for coronal seismology (see e.g. reviews by De Moortel 2005; De Moortel & Nakariakov 2012; Nakariakov & Kolotkov 2020). Coronal seismology is a powerful method that relies on using oscillation parameters that are readily observable in the coronal imaging and spectral data to deduce plasma properties such as density, temperature, magnetic field strength or density structuring, which are otherwise very challenging to measure directly. The accuracy of the seismologically deduced quantities is however limited by somewhat simplifying assumptions of modelling coronal structures as homogeneous plasma cylinders and is subject to uncertainties caused by projection and line-of-sight effects.

Transverse coronal loop oscillations are often classified as occurring in two regimes: large-amplitude rapidly damped oscillations and small amplitude seemingly ’decayless’ oscillations. Large amplitude oscillations are typically triggered by a external impulsive event, which is usually easy to identify, such as blast wave following a solar fare (e.g. White & Verwichte 2012). Amplitudes are of the order of 1 Mm and the oscillations are rapidly damped within 3-4 oscillation periods.

Small amplitude oscillations are observed without any clearly associated driver and can persist for several hours (Wang et al. 2012; Nisticò et al. 2013; Anfinogentov et al. 2015). Their amplitudes are of the order of 100 km and they show no observable damping. Even though it is clear they must be driven by some small-scale abundant process, the exact excitation mechanism of this oscillation regime is unclear. Several excitation mechanisms have been proposed to explain the persistent nature of these oscillations. These include mechanisms acting at coronal heights such as onset of thermal instability in the corona (Kohutova & Verwichte 2017; Verwichte & Kohutova 2017), presence of siphon flows in the coronal loop (Kohutova & Verwichte 2018), self-oscillation of coronal loops due to interaction with quasi-steady flows (Nakariakov et al. 2016; Karampelas & Van Doorsselaere 2020) as well as footpoint-concentrated drivers. Commonly assumed excitation mechanism is associated with turbulent vortex flows which are ubiquitous in the lower solar atmosphere (Carlsson et al. 2010; Shelyag et al. 2011; Liu et al. 2019) and the resulting random footpoint buffeting. This is typically modelled as a stochastic driver at the footpoints (e.g. Pagano & De Moortel 2019). Such driver however does not reproduce the observational characteristics of this type of oscillations, in particular the stability of the oscillation phase and lack of any observable damping (Nakariakov et al. 2016). It has been proposed that a footpoint driver of the form of a broadband noise with a power-law fall-off can in principle lead to excitation of eigenmodes along coronal loops (Afanasyev et al. 2020), this has however not been tested through MHD simulations. Finally, excitation of transverse oscillations through solar p-mode coupling has also been proposed in several studies (Tomczyk & McIntosh 2009; Morton et al. 2016, 2019; Riedl et al. 2019). Since p-modes correspond to acoustic waves generated by turbulent flows in the solar convection zone, they are subject to acoustic cut-off at the $\beta=1$ equipartition layer. High magnetic field inclinations can however lead to p-mode leakage into chromosphere even for frequencies below the cut-off frequency (De Pontieu et al. 2004), and to subsequent coupling to different magnetoacoustic wave modes (Santamaria et al. 2015). A 3D MHD simulations of coronal fluxtubes with gravity-acoustic wave driver at the footpoints however did not find any evidence of the acoustic driving leading to significant transverse displacement of the fluxtubes at coronal heights (Riedl et al. 2019).

It is in fact very likely that the oscillatory behaviour observed in the corona is caused by a combination of several different excitation mechanisms. 3D MHD simulations of solar atmosphere spanning from convection zone to corona in principle include all of the potential mechanisms discussed above. However, so far only impulsively generated rapidly damped coronal loop oscillations have been studied in such 3D MHD simulations (Chen & Peter 2015), which show good match of damping rates compared with observed impulsively generated oscillations. Such numerical works also highlight the implications of observational limitations for coronal seismology, as they enable direct comparison of seismologically deduced and real values of physical quantities. Detailed properties of transverse coronal loop oscillations have been studied numerically in great detail in recent years, including their evolution, damping and development of instabilities and resulting nonlinear dynamics (e.g. Antolin et al. 2014; Magyar & Van Doorsselaere 2016; Karampelas et al. 2017). Most of these studies however employ simplified geometries that model coronal structures as straight uniform fluxtubes. Additionally, such studies also rely on artificially imposed harmonic footpoint drivers (Karampelas et al. 2017; Pagano & De Moortel 2017), without discussing where do these drivers originate from. Such approach makes it possible to isolate individual effects of interests, and focus on evolution at small scales through techniques such as adaptive mesh refinement (e.g. Magyar et al. 2015) due to obvious computational limitations that come with studying physical mechanisms operating at vastly different spatial scales. There is, however, need for more complex approach to see how such effects manifest in a realistic solar atmosphere setup and to identify their observational signatures.

In this work we investigate the presence of oscillations of coronal loops in a simulation spanning from convection zone to corona, with a realistic lower solar atmosphere dynamics, magnetic field geometry and density structuring. We for the first time analyse excitation of transverse coronal oscillations in a self-consistent simulation of the solar atmosphere. We use combination of forward modelling of coronal EUV emission and evolution of physical quantities along magnetic field to determine oscillation characteristics.

We analyse coronal oscillations in the numerical simulation of magnetic enhanced network spanning from the upper convection zone to the corona (Carlsson et al. 2016) using 3D radiation MHD code Bifrost (Gudiksen et al. 2011). This simulation has been previously used for analysing oscillatory behaviour of chromospheric fibrils by Leenaarts et al. (2015) and for studying signatures of propagating shock waves by Eklund et al. (2020).

Bifrost solves resistive MHD equations and includes non-LTE radiative transfer in the photosphere and low chromosphere and parametrized radiative losses and heating in the upper chromosphere, transition region and corona. The simulation includes the effects of thermal conduction parallel to the magnetic field and the non-equilibrium ionization of hydrogen.

The physical size of the simulation box is 24 $\times$ 24 $\times$ 16.8 Mm. The vertical extent of the grid spans from 2.4 Mm below the photosphere to 14.4 Mm above the photosphere, which corresponds to $z=0$ surface and is defined as the (approximate) height where the optical depth $\tau_{500}$ is equal to unity. The simulation is carried out on a 504 $\times$ 504 $\times$ 496 grid. The grid is uniform in the $x$ and $y$ direction with grid spacing of 48 km. The grid spacing in $z$ direction is non-uniform in order to resolve steep gradients in density and temperature in the lower solar atmosphere and varies from 19 km to 98 km. The domain boundaries are periodic in the $x$ and $y$ direction and open in the $z$ direction. The top boundary uses characteristic boundary conditions such that any disturbances are transmitted through the boundary with minimal reflection. Tests of the characteristic boundaries implemented in Bifrost suggest that the reflection of the energy from the upper boundary is 5% or less (Gudiksen et al. 2011). At the bottom boundary the flows are let through and the magnetic field is passively advected while keeping the magnetic flux through the bottom boundary constant. The average horizontal pressure is driven towards a constant value with a characteristic timescale of 100 s. This pressure node at the bottom boundary gives rise to acoustic wave reflection which resembles the refraction of waves in the deeper solar atmosphere. This leads to global radial box oscillations with a period of 450 s, which are a simulation counterpart of Solar p-modes (Stein & Nordlund 2001; Carlsson et al. 2016) (radial modes in this context correspond to oscillations in horizontally averaged quantities polarised along z-axis in a cartesian simulation box).

The photospheric magnetic field has average unsigned value of about 50 G and is concentrated in two patches of opposite polarity about 8 Mm apart; this leads to development of several magnetic loops at coronal heights (Fig. 1). Upon initialisation of the simulation, the large-scale magnetic field configuration was determined by the use of potential extrapolation of the vertical magnetic field, specified at the bottom boundary. After inserting the magnetic field into the domain it is quickly swept around by convective motions. These lead to magnetic field braiding and associated Ohmic and viscous heating which together maintain high temperature in the chromosphere and corona (see e.g. Carlsson et al. (2016); Kohutova et al. (2020a) for further details of the numerical setup).

In order to compare oscillatory behaviour seen in the numerical simulation with the commonly observed transverse coronal loop oscillations seen by SDO/AIA, we forward model the EUV emission using the FoMo code (Van Doorsselaere et al. 2016). The FoMo code is capable of modelling optically thin coronal emission in individual spectral lines as well as in the broadband SDO/AIA channels. The instrumental response function of a given SDO/AIA bandpass $\kappa_{\alpha}$ is calculated for a grid of electron densities $n_{\mathrm{e}}$ and temperatures $T$ as

where $G(\lambda_{\alpha},n_{\mathrm{e}},T)$ is the contribution function calculated from the CHIANTI database, $R_{\alpha}(\lambda_{\alpha})$ is the wavelength-dependent response function of the given bandpass, and the wavelength integration is done over all spectral lines within the given bandpass. $\kappa_{\alpha}$ is then integrated along the line-of-sight parallel to z-axis, which corresponds to looking straight down on the solar surface.

Figure 2 shows forward modelled emission in the SDO/AIA 171 Å bandpass, which is dominated by the Fe IX line with formation temperature $\log T=5.9$. In order to highlight the oscillatory motion, a time-distance plot is created by taking the intensities along a slit across the loop structures and stacking them in time. Several oscillating structures can be seen in such plot, suggesting transverse coronal loop oscillations are abundant in the model. The forward-modelled EUV emission in the model is more diffuse (Fig. 3) and subject to less observable structuring across the magnetic field than coronal emission typically observed at similar wavelengths, where the individual coronal loops can appear very threadlike (e.g. Peter et al. 2013; Williams et al. 2020). We note that the simulation does not include magnetic flux emergence which means the loops are not mass-loaded and pushed into the corona from the lower solar atmosphere. Instead, the dense loops are filled via chromospheric evaporation caused by localised heating (Kohutova et al. 2020a). Some transverse structuring is however still present thanks to the temperature variation across different magnetic fieldlines. This makes it possible to see oscillatory behaviour of the individual coronal strands in the intensity time-distance plot.

Coronal loops in our study are represented by closed magnetic fieldlines extending into the corona. Several of those have enhanced densities compared to the surroundings. Density-enhanced coronal loops with well-defined boundaries in the simulation are however less common than what is usually observed in the corona. As there is no flux emergence mechanism included in the simulation, most of the dense loops are filled by chromospheric evaporation instead of being lifted up from the lower atmosphere.

The magnetic field in the simulation domain is constantly evolving and undergoing complex motions which include sideways displacement, oscillatory motion, torsional and upward/downward motion resulting in the change of the total coronal loop length. Therefore, in order to investigate the oscillatory behaviour of coronal loops in the simulation, it is necessary to trace the evolution of the corresponding magnetic field lines through both time and space. To do this, we use a field tracing method previously used by Leenaarts et al. (2015) and Kohutova et al. (2020a). A magnetic fieldline is defined as a curve in 3D space $\@vec{r}(s)$ parametrised by the arc length along the curve $s$, for which $\mathrm{d}\@vec{r}/\mathrm{d}s=\@vec{B}/|\@vec{B}|$.

The tracing of the evolution of magnetic fieldlines is done by inserting seed points (one seedpoint per tracked fieldline) into the locations in the simulation domain which show oscillatory behaviour in the velocity time-distance plots. This is equivalent to tracking the evolution of magnetic conectivity of an oscillating plasma element. The seed points are then passively advected forward and backward in time using the velocity at the seed point position. At every time step the magnetic field is then traced through the instantaneous seed point position in order to determine the spatial coordinates of the traced fieldline. This is done using the Runge–Kutta–Fehlberg integrator with an adaptive step size. Even though the accuracy of this method is limited by the size of the time step between the two successive snapshots (i.e. 10 seconds), it works reasonably well provided that the magnetic field evolution is smooth and there are no large amplitude velocity variations occurring on timescales shorter than the size of the time step. We note that this method leads to a jump in the fieldline evolution in the instances where magnetic reconnection occurs.

We investigate the evolution of 5 different fieldlines labelled F1…F5 (Fig. 4). We note that the footpoints of chosen magnetic loops lie in the regions of enhanced vorticity. We analyse the evolution of the temperature, density, and 3 velocity components, $v_{\mathrm{L}}$, $v_{\mathrm{V}}$ and $v_{\mathrm{T}}$ along each loop. The longitudinal velocity $v_{\mathrm{L}}=\@vec{v}\cdot\@vec{T}$ corresponds to the velocity component aligned with the tangent vector of the magnetic fieldline given by $\@vec{T}=\@vec{B}/|\@vec{B}|$. The vertical velocity $v_{\mathrm{V}}=\@vec{v}\cdot\@vec{N}$ is the velocity component along the fieldline normal vector given by $\@vec{N}=\frac{\mathrm{d}\@vec{T}}{\mathrm{d}s}/|\frac{\mathrm{d}\@vec{T}}{\mathrm{d}s}|$ and corresponds to the motion in the plane of the coronal loop. Finally the transverse velocity component along the binormal vector is given by $v_{\mathrm{T}}=\@vec{v}\cdot\@vec{R}$ where $\@vec{R}=\@vec{T}\times\@vec{N}$ and corresponds to transverse motion perpendicular to the plane of the coronal loop. Unit vectors $\@vec{T}$, $\@vec{N}$ and $\@vec{R}$ together form an orthogonal coordinate system equivalent to a Frenet frame of reference. Such coordinate frame is well suited for analysing oscillations in complex 3D magnetic field geometries and is commonly used in such studies (e.g. Felipe 2012; Leenaarts et al. 2015; González-Morales et al. 2019). We further calculate a wavelet power spectrum (Torrence & Compo 1998) for all three velocity components at 3 different locations along the fieldline; close to the loop footpoints 1 Mm above the transition region (height of which is tracked with time) at the beginning and end of the fieldline (labelled FP1 and FP2 for left and right footpoint respectively) and at the loop apex halfway between FP1 and FP2.

An important point for consideration is the applicability of the fieldline tracing method that is used to study the fieldline evolution in the instances of magnetic reconnection. Tracing of the magnetic fieldlines is equivalent to tracing the grid points in the simulation domain which are magnetically connected. Since in the corona the matter and energy transport happens along the magnetic field, tracing the magnetic fieldlines is therefore the closest we can get to studying the true evolution of coronal structures without being affected by line of sight effects, which can significantly influence determination of oscillation parameters (see e.g. De Moortel & Pascoe 2012). When a magnetic loop ‘reconnects’, the magnetic connectivity of a certain part of the loop changes. However, advecting the seed points initially placed into regions where oscillations were detected ensures we keep tracing the magnetic connectivity of the oscillating part even if reconnection occurs in other parts of the loop. We therefore argue that this approach remains the most feasible way of tracing the true evolution of coronal structures.

We note that evolution in line-of-sight integrated emission does not necessarily copy the evolution of traced magnetic fieldlines with the same $x,y$ coordinates. This is likely due to line-of-sight superposition of emission from several magnetic structures, which can vary in emissivity and in presence or lack of oscillatory behaviour. This effect of line-of-sight superposition on oscillation analysis has been discussed by previous studies (De Moortel & Pascoe 2012; Leenaarts et al. 2015). Our work further highlights the limitations posed by line-of-sight effects which should be taken into account when analysing observations of coronal oscillations.

The identification of observed transverse oscillations as either standing and propagating is complicated by the relatively short length of coronal loops in the simulation and by high typical values of kink speed in the corona which translates to high phase speed of the oscillation. Assuming a mean kink speed in the range 100 - 1000 km s^-1, the expected propagation time along a 20 Mm long loop varies from 200 s to 20 s, corresponding to 20 and 2 time steps of the simulation output, respectively. The position of the maximum oscillation amplitude higher up along the loop rather than at the loop footpoints suggests the oscillations are in fact standing with the loop footpoints acting as oscillation nodes. The term nodes is however used loosely here, since the footpoints are in fact not static, but continuously moving (Fig. 10). There are no velocity nodes clearly observable in the coronal sections of studied loops and the longitudinal profile of $V_{\mathrm{T}}$ component matches fundamental oscillation harmonic. We note that similar oscillations in chromospheric fibrils analysed in synthetised H-alpha emission have been identified as propagating by Leenaarts et al. (2015). The method commonly used to investigate presence of oscillation phase shift and the corresponding phase speed is based on detecting phase shift in time-distance plots taken at several different positions along the studied loop, using observed or forward-modelled emission. It should be noted however, that such method only works for static loops. It is likely to produce spurious results in the presence of the bulk motion of the fieldline across of the slits in addition to the oscillatory motion, which seems to be the case for all 5 fieldines studied in this work. Instead we focused of the evolution of the oscillation phase at different locations along each loop that have been traced with time, while accounting for the loop 3D structure and evolution. No persistent phase shift has been identified between the loop apex and the footpoints in any of the analysed cases, that would suggest observed transverse oscillations are propagating.

On the other hand, such distinction is easily made for longitudinal oscillations as the typical values of the sound speed in the corona are order of magnitude smaller compared to Alfvén or kink speed (see Table 1). Longitudinal oscillations propagating from one footpoint to another with local sound speed were detected in F2. The mean time delay between the three observable successive wave trains is around 250 s. Such propagating waves can be a result of an acoustic driver (potentially linked to the global p-mode equivalent oscillations) which drives compressive longitudinal waves that steepen into shocks in the chromosphere and propagate through the corona. Such propagating oscillations were however not universally present for all of the studied fieldlines, hence we refrain from drawing any conclusions about the periodicity (or lack thereof) of the driver. We note that the period of the global p-mode oscillation seen in the simulation is greater, around 450 s, as discussed in section 2. Standing longitudinal oscillations sustained for few oscillation periods were detected in F1, F3, F4 and F5, with relatively short periods between 100-200 s, depending on the loop. The longitudinal velocity profiles showing an oscillation node at the apex of the loops suggest these correspond to a second harmonic and likely represent a response of the loop to external perturbations occurring at coronal heights (Nakariakov et al. 2004). In the loops that show the presence of second harmonic of a longitudinal standing mode the onset of the oscillation follows events associated with increase of temperature of the plasma, and hence likely linked to small scale reconnection events/other impulsive heating events. Oscillation damping is observable for most cases, however, detailed analysis of the damping profiles and associated dissipation is beyond the scope of this work.

We also highlight the evolution of $v_{\mathrm{T}}$ in F3, which shows clear oscillation with 180 s period that is sustained over time and does not show any observable damping (Fig. 11). There is no consistent phase shift between the oscillatory behaviour at the loop apex and in the loop legs, suggesting the oscillation is standing. This is further supported by the fact that the location of the maximum velocity amplitude lies close to the loop apex. This oscillation pattern is very similar to the commonly observed regime of decayless oscillations (Wang et al. 2012; Nisticò et al. 2013; Anfinogentov et al. 2015). Conversely, an example of the damped impulsively generated (or ”decaying”) oscillation regime is observable in F2 following an event associated with impulsive temperature increase in the loop (Fig. 11). The oscillation velocity amplitude is largest at the apex and the evolution of the longitudinal $V_{\mathrm{T}}$ profile matches fundamental harmonic of a standing transverse oscillation. This suggests it is a natural response of the loop to a perturbation accompanied by an impulsive energy release. We also note that for several cases of oscillations in loops F1, F4 and F5, the classification into ”decaying” and ”decayless” is not as clear (due to lack of clear steady damping pattern such as shown in Fig. 11). F3 is the only loop that shows persistent undamped oscillation present for the whole duration of simulation.

The oscillatory motions seen in the simulation in magnetic structures, both ”decaying” and ”decayless” regimes correspond to standing transverse oscillations with the oscillation antinodes lying at coronal heights and are therefore consistent with the observations of transverse oscillations of coronal loops as seen in TRACE and SDO/AIA data (e.g. Aschwanden et al. 1999; White & Verwichte 2012), where the standing transverse oscillation are the most commonly observed regime.

We further note that the detected oscillation periods vary between different fieldlines and can also change with time, particularly in loops that are subject to change in physical properties, or with increasing/decreasing loop length, as would be expected for standing oscillations (this is true for all oscillation modes visible in the simulation). This spread and variability of the oscillation periods as well as lack of oscillation coherence in different parts of the simulation domain suggests that they are not a result of a coherent global oscillation in the simulation. Our analysis therefore does not agree with the premise that the transverse oscillations in the corona are driven by the global p-mode oscillation which is based on comparison of peaks in the power spectra from velocity measurements in the corona and the p-mode power spectra (Morton et al. 2016, 2019) (it should however be noted that analysis in such studies focuses on very long coronal loops, i.e. the spatial scales are very different from the short low-lying loops studied in this work). Such mechanism requires coronal loops to oscillate with the same period regardless of their length, as they are driven at their footpoints by a harmonic driver. The variability of the oscillation period in the simulation therefore excludes the presence of a harmonic driver, or at the very least suggests that such driver is not a dominant mechanism responsible for the excitation of transverse oscillations in this work. We stress that care should be taken when drawing conclusions from the global power spectra due to the temporal variability of the oscillation periods of individual magnetic structures which has been shown both observationally (Verwichte & Kohutova 2017) and by numerical simulations (Kohutova & Verwichte 2017).

Excitation mechanism proposed by Nakariakov et al. (2016) is based on identifying the decayless oscillations as self-oscillations of the coronal loop, that are excited by interaction with quasi-steady flows, similar to a bow moving across a violin string. The flows in the simulation domain are too dynamic and not stable enough to drive persistent self-oscillation of the fieldlines; this is evidenced by the motion of loop footpoints that are dragged around by turbulent flows in the lower solar atmosphere shown in Fig. 10 and show variability on the order of minutes. 3D numerical models however seem to suggest that excitation of self-oscillations requires persistent flows that are steady on timescales of the order of an hour (Karampelas & Van Doorsselaere 2020).

The absence of increased oscillatory power in the chromosphere at the positions of the footpoints of the oscillating fieldlines also suggests that a harmonic footpoint driver is not a prerequisite for the excitation of coronal oscillations. Open questions addressing the excitation of the decayless oscillations therefore still remain. If their excitation mechanism is global, why are they not observed in all coronal loops? The decayless oscillations are sufficiently abundant to not be driven by an isolated event (Anfinogentov et al. 2015), but not enough to classify global process as the driver.

Finally, the differences between the magnetic and density structure of the simulated and real solar corona, and hence the limitations of the model used for our analysis need to be addressed. The forward modelled EUV emission from the simulated corona is comparatively more diffuse and lacks a lot of fine-scale structure seen in real coronal observations, especially fine-scale coronal strands revealed by observations from the second Hi-C flight (Williams et al. 2020), which are not resolved in SDO/AIA observations. The simulated corona also suffers from a lack of coronal loops understood in the traditional sense as distinct structures denser than the surroundings with distinct boundaries. As the simulation does not include any flux emergence, the coronal loops in the simulation do not correspond to overdense flux tubes lifted into the corona from the lower solar atmosphere but are instead filled with evaporated plasma due to impulsive heating events. Due to enhanced heating regions having an irregular shape (Kohutova et al. 2020a) the structures filled with evaporated plasma have greater spatial extent and lack well-defined boundaries. Simulation resolution might also be a limiting factor, eventhough the characteristic transverse scales seen in time-distance plots are well above the horizontal grid size of 48 km. Distinct oscillating strands are however still observable in Fig. 2. Furthermore, recent numerical studies of evolution of initially homogeneous coronal loops as a response to transverse motions suggest that our highly idealised picture of coronal loops as monolithic plasma cylinders is not very realistic in the first place (Magyar & Van Doorsselaere 2018; Karampelas et al. 2019; Antolin & Van Doorsselaere 2019).

We have studied the excitation and evolution of coronal oscillations in 3D self-consistent simulations of solar atmosphere spanning from convection zone to solar corona using radiation-MHD code Bifrost. We have combined forward-modelled EUV emission with 3D-tracing of magnetic field though both space and time in order to analyse oscillatory behaviour of individual magnetic loops, while accounting for their dynamic evolution. We have analysed evolution of different velocity components, using wavelet analysis to capture changes in the oscillatory periods due to evolution of properties of the oscillating magnetic loops. Various types of oscillations commonly observed in the corona are reproduced in such simulations. We detected standing oscillations in both transverse and longitudinal velocity components, including higher order oscillation harmonics. We have further shown that self-consistent simulations reproduce existence of two distinct regimes of transverse coronal oscillations, that is rapidly decaying oscillations triggered by impulsive events and sustained small-scale oscillations showing no observable damping. Damped transverse oscillations were found to be associated with instances of impulsive energy release, such as small scale reconnection events, or external perturbations leading to sudden changes in the loop properties, in agreement with coronal observations. Persistent transverse oscillations on the other hand were not linked to any such impulsive events. We did not find any evidence for this oscillation regime being driven by a global (simulation) p-mode. Lack of enhanced oscillatory power near the footpoint regions of the studied loops together with variability of oscillation periods between different coronal loops and lack of oscillation coherence across the simulation domain exclude any type of harmonic driver as being responsible for excitation of the oscillations.

Our work therefore highlights the complexity of coronal oscillations in simulations with realistic magnetic field configurations which include the complex dynamics of the lower solar atmosphere. Care needs to be taken when translating findings from highly idealised models into real solar observations, as there are several limitations to treating coronal loops as static structures. We have shown that individual fieldlines are very dynamic, their footpoints migrate and their overall length changes significantly over realistic timescales. This might have non-negligible consequences for accuracy of estimates of coronal plasma parameters deduced using coronal seismology. The oscillating coronal structures we analysed are far from idealised plasma cylinders. The focus of modelling work should therefore be shifting from rigid cylindrical models towards more realistic descriptions that account for rapid variability, complex morphology and presence of nonlinearities. There are obviously natural limitations to this approach resulting from the computational expense of building such models and also from the associated limits on physical size of the simulation domain. Work is currently underway to investigate the evolution of coronal loops in 3D self-consistent simulations that span high into the corona. We would hence like to highlight the potential of using self-consistent simulations of the solar atmosphere as a laboratory for testing assumptions made by coronal seismology and models of various damping and dissipation mechanisms in an environment with realistic density structure and magnetic field geometry.