Damping of nonlinear standing kink oscillations: a numerical study
Centre for mathematical Plasma Astrophysics,
Department of Mathematics, KU Leuven,
Celestijnenlaan 200B, bus 2400, 3001 Leuven, Belgium
email: norbert.magyar@wis.kuleuven.be
Received:
27
May
2016
Accepted:
11
September
2016
Aims. We aim to study the standing fundamental kink mode of coronal loops in the nonlinear regime, investigating the changes in energy evolution in the crosssection and oscillation amplitude of the loop which are related to nonlinear effects, in particular to the development of the KelvinHelmholtz instability (KHI).
Methods. We run ideal, highresolution threedimensional (3D) magnetohydrodynamic (MHD) simulations, studying the influence of the initial velocity amplitude and the inhomogeneous layer thickness. We model the coronal loop as a straight, homogeneous magnetic flux tube with an outer inhomogeneous layer, embedded in a straight, homogeneous magnetic field.
Results. We find that, for low amplitudes which do not allow for the KHI to develop during the simulated time, the damping time agrees with the theory of resonant absorption. However, for higher amplitudes, the presence of KHI around the oscillating loop can alter the loop’s evolution, resulting in a significantly faster damping than predicted by the linear theory in some cases. This questions the accuracy of seismological methods applied to observed damping profiles, based on linear theory.
Key words: Sun: corona / Sun: oscillations / magnetohydrodynamics (MHD)
© ESO 2016
1. Introduction
With the first observations of transverse oscillations in coronal loops by TRACE (Aschwanden et al. 1999; Nakariakov et al. 1999) a new era has begun in the exploration and understanding of the solar corona. These observations of waves allow the inference of different physical parameters which were previously largely unknown, a method called coronal seismology. It was quickly noted that the observed coronal loop oscillations were strongly damped, and if such a high dissipation was due to viscous or resistive damping, then the associated transport coefficients must be orders of magnitude higher than predicted by the classical theory. Nearly all of the observed, highamplitude, lower coronal eruptionrelated coronal loop oscillations (Zimovets & Nakariakov 2015) are strongly damped, with rare but nevertheless very intriguing exceptions (see, e.g. Aschwanden et al. 2002; Aschwanden & Schrijver 2011; Wang et al. 2012). Recently, small amplitude decayless kink oscillations were observed in coronal loops (Nisticò et al. 2013; Anfinogentov et al. 2013), and are thought to be ubiquitous in active regions (Anfinogentov et al. 2015). These oscillations appear to be constantly driven, and do not originate from an impulsive or eruptive event like the high amplitude oscillations. It is generally accepted that the mechanism responsible for the fast damping of transverse oscillations is the extensively studied resonant absorption or mode coupling (e.g. Ionson 1978; Hollweg 1984; Ruderman & Roberts 2002; Goossens et al. 2002). For resonant absorption to damp the oscillations, the existence of a surface within the limits of the loop, with Alfvén speed matching the global kink phase speed, is required. This is usually portrayed as a smooth layer around a homogeneous flux tube, but resonance has been proven to exist regardless of the geometrical shape of such a layer (Terradas et al. 2008b; Pascoe et al. 2011). The linear theory of kink oscillations is well studied (for a review, see Ruderman & Erdélyi 2009). The damping profile of oscillations, in the presence of a driver (or initial perturbation in the case of standing modes) has been shown to deviate from a purely exponential decay, being described rather by an initially Gaussian damping profile followed by an exponential damping profile (Pascoe et al. 2012; Hood et al. 2013; Pascoe et al. 2013, 2016b,a). Observed average periods, amplitudes and exponential damping times of kink oscillations in coronal loops were reported by, Aschwanden et al. (2002), for example. In their analysis of 26 loop oscillation events, they find average oscillation periods of 321 ± 140 s, oscillation amplitudes of 2200 ± 2800 km and damping times of 580 ± 385 s. The measured amplitudes correspond to relative amplitudes (to the loop length L) of approximately 1–5%. Arregui et al. (2007) and Goossens et al. (2008) used the analytical formula for the damping time (Ruderman & Roberts 2002) to construct a seismological method for inferring the loop parameters from the observed damping time and period.
Waves in the solar atmosphere often have high enough amplitudes to be considered nonlinear. The study of nonlinear waves in the solar atmosphere has been carried out especially in the context of chromospheric and coronal heating. For a review of earlier theoretical work on the subject, see Ruderman (2006). In particular, the analytical theory of nonlinear kink oscillations has been studied, both for propagating (Ruderman et al. 2010), and standing waves (Ruderman & Goossens 2014). In these studies it was shown that damping of propagating kink waves can be enhanced by the nonlinearity of an mmode resonance, where energy from the m = 1 kink mode is transferred to m ≥ 2 fluting modes, which, by resonant absorption, can damp faster due to shorter wavelengths. It was noted that the mmode resonance can damp the kink wave even in the absence of a resonant layer. However, in an axially inhomogeneous flux tube, due to density stratification, for example, the mmodes no longer have the same phase speed, thus the mmode resonance and enhanced damping disappears. These calculations were valid for weakly nonlinear oscillations, and the authors anticipated that the amplitude would be affected in the fully nonlinear case. In general, the study of fully nonlinear problems is only possible numerically. Numerical studies of nonlinear kink oscillations of coronal loops have been carried out by, Terradas & Ofman (2004), Ofman (2005), Terradas et al. (2008a), Ofman (2009b), Antolin et al. (2014, 2015), Magyar et al. (2015), Magyar & Van Doorsselaere (2016), for example. See also Ofman (2009a) for a review. Of particular and renewed interest in nonlinear evolution of transverse waves, even for low amplitudes, is the susceptibility of the resonant layer, due to its high velocity shear, to the KelvinHelmoltz instability (KHI; Heyvaerts & Priest 1983; Browning & Priest 1984; Karpen et al. 1994; Ofman et al. 1994; Poedts et al. 1997). This can enhance the wave energy dissipation via local turbulence, and thus is also relevant in the coronal heating problem. Direct observational evidence is lacking of the KHI in coronal loops. However, it has been observed in coronal mass ejections and quiescent prominences (Berger et al. 2010; Ryutova et al. 2010; Foullon et al. 2011; Ofman & Thompson 2011), and it was suggested that, in fact, the KHI may appear as strands of coronal loops as seen in EUV (Antolin et al. 2014). The previous statement is strengthened by the first observational evidence of resonant absorption in prominences (Okamoto et al. 2015; Antolin et al. 2015). Recently, the idea of observing the damping profile of kink oscillations of coronal loops in order to infer various local parameters such as density ratios and inhomogeneous layer widths in the context of coronal seismology was put forward (Pascoe et al. 2016a). Therefore, it is essential to know the effects that nonlinearity might have on damping profiles. New observational evidence suggests that the damping of the transverse oscillations of coronal loops is dependent on the amplitude (Goddard & Nakariakov 2016).
In this paper we investigate the nonlinear standing kink oscillation of coronal loops, for which the main effect of nonlinearity is the development of the KHI at the loop edges, altering the energy distribution in the loop cross section and ultimately leading to a change in the damping profile of the loop displacement.
2. Numerical model
2.1. Initial conditions
We model the coronal loop as a straight, densityenhanced magnetic flux tube. The loop consists of a homogeneous inner core, and an inhomogeneous annulus (transitional layer), in which the density varies sinusoidally from the core (ρ_{i} = 2.5 × 10^{12} kg m^{3}) to the background density value (ρ_{e} = 0.5 × 10^{12} kg m^{3}), where the subscripts i, e stand for internal and external, respectively. The numerical domain is permeated by a straight, homogeneous magnetic field, of 12.5 G strength. The distance between the central axis of the flux tube and the midpoint of the inhomogeneous layer defines its radius, R = 1.5 Mm. The thickness of the inhomogeneous layer is denoted by l. We neglect the effects of curvature, gravity (i.e. no stratification), energy sources and sinks (heating, thermal conduction, radiative processes), and we lack a realistic lower solar atmosphere (i.e. photosphere, chromosphere). The plasmaβ is constant throughout the domain, 0.06, corresponding to T_{e} = 4.5 MK and T_{i} = 0.9 MK. The resulting temperature profile is most probably not realistic. This simplistic profile was chosen in order to start with an initial equilibrium, that is, constant total pressure throughout the domain. Tests with isothermal models (T_{e} = T_{i} = 0.9 MK) show no significant difference in comparison with the present model, as expected (period, growth rates and damping rates not depending on the temperature profile).
Initially, we impose a perturbation in a component of the velocity transverse to the loop axis, of the form (1)V_{A,i} is the Alfvén speed inside the loop (V_{A,i} ≈ 0.7 Mm s^{1}), L = 120 Mm is the total loop length, and A is the perturbation amplitude. Thus the perturbation acts only inside the loop, including the inhomogeneous layer. This excites the fundamental kink mode of the loop, and we stop the simulation at t_{f} = 1500 s.
2.2. Boundary conditions
In order to excite a standing transverse oscillation, we fix the footpoint of the loop (z = 60 Mm), by setting the transverse velocities to antisymmetric, while the other variables are set to continuous (Neumanntype, zerogradient boundary condition) at this boundary. Exploiting the symmetric properties of the fundamental kink mode, we model only half of the loop in both the z and x directions. This reduces the computational time fourfold. The boundary conditions for these mirroring boundaries are described in Magyar et al. (2015). At the other, lateral boundaries we let any waves leave the domain freely by imposing an outflow (zerogradient) condition on all variables.
2.3. Numerical method and mesh
To solve the ideal 3D magnetohydrodynamic (MHD) problem, we use the MPIAMRVAC code (Keppens et al. 2012; Porth et al. 2014). We use the implemented secondorder “onestep” TVD method with the Roe solver and “Woodward” slope limiter. The constraint on the magnetic field divergence is maintained using Powell’s scheme. The numerical domain has dimensions of (0,6) × (−10,10) × (0,60) Mm. The base resolution is 24 × 80 × 32, and we use 4 levels of refinement, fully refining the region around the loop, resulting in an x − y plane resolution of 31.2 km per cell, or 0.02R, where R is the radius of the loop. The fullyrefined resolution in the z direction is 234 km per cell. Comparison with simulations with one more level of refinement shows no important differences in the dynamics, even though the KHI instability (and later turbulence) changes quantitatively.
3. Results and discussion
We ran a series of simulations which explored the parameter space, varying the initial velocity perturbation and thickness of the inhomogeneous layer of the loop. The chosen values are the following: A = { 0.005,0.01,0.02,0.035,0.05 } and l = { ≈ 0.0R,0.1R,0.33R,0.5R }, for a total of 20 different runs (Fig. 1).
Fig. 1 Density in the antinode crosssection of the loop (z = 0) at t = 540 s for the different parameters used (A = 0.005 not shown for brevity, looking almost identical to A = 0.01). Different columns have different initial inhomogeneous layer widths, shown at the top, and the different rows are for different amplitudes, shown in the left margin. Axis units are in Mm and density in units of 10^{12} kg m^{3}. 

Open with DEXTER 
The first value for l denotes an initially step profile which, after the start of the simulation, evolves into an inhomogeneous layer due to numerical diffusion.
Fig. 2 Density in the antinode crosssection of the loop (z = 0) at different times (in steps of approximately one halfperiod, shown in the upper part of each plot) for A = 0.035 and l = 0.33R. Axis units are in Mm and density in units of 10^{12} kg m^{3}. 

Open with DEXTER 
This layer is thus the thinnest possible layer in our simulation, corresponding to approximately 0.08R. The initial displacements of the loop are, in increasing order with the amplitude, { 0.079R,0.16R,0.32R,0.55R,0.78R }. Obtained oscillation periods are discussed later in the text. The nonlinearity parameter in the case of kink oscillations of a flux tube is given approximately by Ruderman & Goossens (2014): (2)The oscillations are weakly nonlinear for ν ≪ 1 and strongly nonlinear for ν ≥ 1. For the amplitudes chosen, the nonlinearity parameter varies between 0.4 and 4. Thus, all the simulations represent nonlinear and strongly nonlinear regimes of oscillations. As the simulations with A = 0.005 and A = 0.01 show almost similar evolution, lower amplitudes were not essential for this study. In Fig. 2, the evolution of the oscillation in the antinode crosssection can be followed, for a specific set of parameters. In this sequence, we can identify the principal signatures of nonlinear kink oscillations previously described in detail in numerical simulations by Terradas et al. (2008a), Antolin et al. (2014, 2015), Magyar et al. (2015). These works however, lack a description of one aspect which plays a role in the nonlinearly enhanced damping of the oscillations, namely the resonance between the kink and m ≥ 2 or fluting modes, which extracts energy from the kink mode, leading to damping of the transverse oscillation. This effect was studied analytically by Ruderman et al. (2010), and is present only in the absence of longitudinal inhomogeneities (e.g. gravitational stratification), which destroys the resonance between the modes (Ruderman & Goossens 2014). This effect can be seen clearly in the second snapshot of Fig. 2 (t = 135 s, approximately halfperiod), before the development of the KelvinHelmholtz instability (KHI), as a deformation of the initially circular crosssection in a manner resembling the m = 2 fluting mode. There are three main mechanisms acting to damp the transverse oscillation (aside from a very small numerical dissipation): resonant absorption or mode conversion, the previously mentioned m ≥ 2 resonance, and the development of the KHI (in the case shown in Fig. 2, in less than a period). These mechanisms act concomitantly and are coupled, resulting in the significance and magnitude of each mechanism varying in time, which is difficult to estimate. The kink mode presents a velocity shear near the lateral boundary of the loop, which is prone to the KHI. Resonant absorption, acting immediately after t = 0, further enhances the velocity shear, thus shortening the time after which the instability sets in. Note that, even in the absence of an inhomogeneous layer (which, due to numerical diffusivity, is not possible in our simulations), the KHI can still develop. In this sense, in the nonlinear regime, transverse oscillations will always be damped.
As previously stated, estimating the relative importance of the interacting damping mechanisms is not an easy task. By looking at the transverse velocity (v_{y}), in Fig. 3, we can appreciate the interaction of two mechanisms, resonant absorption and KHI.
Fig. 3 ycomponent of the velocity (v_{y}) at different times (selected times corresponding to plots in Fig. 2, shown in the upper part of each plot), for A = 0.035 and l = 0.33R. Axis units are in Mm and velocity in units of km s^{1}. 

Open with DEXTER 
Initially, the analytical kink eigenfunction for velocity develops after the perturbation at t = 0, representing a uniform speed within the boundaries of the loop and a dipolar field around it. As noted earlier, a velocity shear is present in this configuration, which, as shown in the t = 135 s plot, is later further enhanced by resonant absorption. Later on, the resonant layer is disrupted, its KE contributing to the development of the KHI. Due to the robust nature of resonant absorption (Terradas et al. 2008b), it is still present in patches, that is, on the surface of the rollups, but ceases to exist in the form which was present at t = 135 s. It is difficult to estimate the relative effectiveness and the contribution to the damping of the patchy resonant absorption compared to the predisruption phase. Ultimately, the rollups break up, resulting in a turbulent inhomogeneous layer, where most of the remaining wave energy is situated. The impact of the mmode resonance on the damping time is not clear. Magyar et al. (2015) studied nonlinear fundamental standing kink oscillations in a stratified atmosphere, in which the resonance is still noticeable, despite the longitudinal inhomogeneity. In Fig. 2 of Ruderman et al. (2010) the enhanced damping as a function of a nonlinearity parameter N is shown, which is roughly equal to the driving amplitude divided by the width of the inhomogeneous layer in units of R. For our parameters, this value is at most 0.5. This value corresponds to a very small effect on the oscillation amplitude. In this sense, we can be somewhat confident that this damping mechanism does not have a substantial impact on the resulting damping times, and in what follows we focus on the other damping mechanisms.
We studied the energy evolution in the antinode crosssection of the loop for the different parameters. For the case in Fig. 2, the average values of relevant variables are plotted as a function of simulation time in Fig. 4.
Fig. 4 Time evolution of different quantities in the antinode crosssection of the loop (z = 0) with A = 0.035 and l = 0.33R. Topleft: average kinetic energy (KE) density in the core and layer region of the loop. Topright: average internal energy density in the core and layer regions of the loop. Bottomleft: average temperature of the loop (average over both the core and layer). Bottomright: change in the area relative to the initial value, for the core region of the loop and the numerical domain outside the loop (outside both core and layer). 

Open with DEXTER 
Note that, in the antinode of the fundamental kink, all of the kink mode’s energy is in the form of KE. The oscillations present in this plot, with the double period of the kink mode, represent the exchange between KE and magnetic energy (ME), analogous to the exchange between kinetic and elastic potential energy in a classical harmonic oscillator (if we account only for the perturbation to the ME). To proceed further, we make the distinction between the core of the loop, defined as the region of the crosssection where ρ ≥ 0.98ρ_{i} and the ME is negligible (see Goossens et al. 2014, for example), and the inhomogeneous layer, defined as the region of the crosssection where 0.98ρ_{i} ≥ ρ ≥ 1.1ρ_{e}. Inspecting the figures, we observe that the average KE density in the core of the loop has a decreasing peak energy. At the same time the energy contained in the inhomogeneous layer is undergoing a different evolution: the peak values are decreasing, but its minimal values show a strong initial increase, which can be attributed to the increase in the energy of localized Alfvén waves, through resonant absorption (see, e.g. Poedts & Kerner 1991). Simultaneously, the KE in the layer is dissipated, resulting in an increase of internal energy (IE). The average IE density shows a longperiod oscillation, triggered by the perturbation at t = 0, corresponding to a longitudinal slow mode. The slow mode is present both in the core and the layer, however the background trend of average IE density shows a steady increase in the layer, in accordance with the KE loss. The IE density plot shows the average perturbation to its equilibrium value over time, and the rise due to the dissipation of the KE represents a ≈0.36% increase in IE density. This cannot account for the apparent heating of the loop (core with layer), from an average temperature of 1.17 MK to 1.47 MK: this would require around 50 times more energy converted into IE. Inspecting the change in the area of the core and the surrounding corona (outside the loop, area of the numerical domain), we notice that both shrink, resulting in an enlargement of the inhomogeneous layer. However, we also notice that the core area diminishes to a lesser extent than the corona. This implies that the higher average temperature of the loop is due to the enlarging inhomogeneous layer, mixing with the hot and rarefied plasma surrounding the loop, driven by the KHI. This effect is present only because of the chosen initial conditions, that is, with a corona five times hotter surrounding the loop. Note that the total IE in the layer (not shown here) presents a considerable increase (≈150%), but this is almost entirely (i.e. except the minute increase converted from KE) a result of the enlargement of the layer, caused by the KHI.
In the plots of Fig. 5 the evolution of the average KE density in the inhomogeneous layer for different initial amplitudes and layer thickness can be compared.
Fig. 5 Average kinetic energy density over time in the inhomogeneous layer of the antinode crosssection of the loop (z = 0). Left: l = 0.1R. Right: l = 0.5R. The different lines are for different initial amplitudes, from 0.01 to 0.05. 

Open with DEXTER 
These plots are normalized to the initial energy density contained in the layer, (i.e. of the initial perturbation), so that a comparison for different amplitudes is straightforward. We note that the relative increase of KE density in the layer is higher for smaller amplitudes, and it requires more time to dissipate. This demonstrates that the presence of KHI in the higher amplitude oscillations acts to enhance the conversion of KE to IE. This effect is more pronounced for wider initial inhomogeneous layers.
Fig. 6 Displacement amplitude over time in the loop antinode (z = 0). The values are normalized to the initial displacement. For the top plots (A = 0.01, left: l = 0.1R, right: l = 0.5R), the Gaussian and exponential profiles from Eq. (3) obtained by fitting Eq. (5) are plotted. For the bottom plots (A = 0.05, left: l = 0.1R, right: l = 0.5R), the Gaussian of Eq. (3) is fitted for the whole data and plotted, as well as the exponential resulting from fitting Eq. (6). 

Open with DEXTER 
In the following, we analyze the evolution of the kink oscillations triggered by the perturbation and their associated damping. For this, we track the centre of mass for the loop in the fundamental kink antinode crosssection where the displacement is maximum. The damping profile of impulsively triggered standing kink modes in a flux tube is not purely exponential (Hood et al. 2013; Pascoe et al. 2013, 2016a). Instead, initially the damping is approximated by a Gaussian profile, and after a time t_{s} by an exponential profile. This profile has the form: (3)and the switch time t_{s} is defined as: (4)Thus, by fitting observed standing kink oscillations of coronal loops to a function (5)one could seismologically estimate the density ratio of the loop and inhomogeneous layer thickness (Pascoe et al. 2016a). We fit Eq. (5), using mpfitfun.pro (Markwardt 2009), to the simulated damping profiles to obtain four parameters, A_{0},τ_{g},τ_{d}, and . Note that t_{s} is a constrained parameter (Eq. (4)) using the “.TIED” entry of the PARINFO structure, and φ is zero as it results from the initial condition. We find bestfittings for the period between 253.5 and 265.1 s; thicker inhomogeneous layers resulting in the shorter periods. The period for the simulation closest to the linear case (A = 0.005,l ≈ 0.0R) is 265.1 s. The theoretical fundamental kink period for a stepdensity flux tube with the same parameters is 263.6 s (Edwin & Roberts 1983). Some results from fitting to the simulated damping profiles are shown in Fig. 6. The amplitudes were normalized to the value of the maximal initial displacements. Note that for the lowamplitude case (A = 0.01, upper plots), the fit is much more accurate (lower χ^{2}) than for the highamplitude counterpart (A = 0.05, lower plots). Furthermore, we find that, in the runs which develop KHI (A ≥ 0.02), the last periods of oscillation are practically undamped. This can also be seen in Fig. 4 (top left), in the core average KE, which remains practically constant in maximal value after ≈1000 s. We interpret this as a consequence of the KHI: some of the kinetic energy percolates into the core region. Although a small effect, this unfortunately makes the fitting of Eq. (5) meaningless (i.e. a Gaussian profile will always fit better to the oscillations, for any value of t_{s}). Overcoming this problem by truncating the dataset where the constant KE period begins (e.g. 800 s) is unreliable as t_{s} varies greatly even for small changes in the truncation. Thus, great care should be taken when fitting Eq. (5) to observed oscillations for seismology. The oscillations resulting from A ≥ 0.02 are best described by a Gaussian profile. Therefore, we encourage the fitting of Gaussian damping profiles alongside the traditional exponential damping profile in future studies of transverse coronal loop oscillations. For our results to be readily comparable to available studies, we quantify the exponential damping time, and to obtain t_{s} for the high amplitude oscillations, we fit an exponentially damped sine; (6)for the whole dataset. The exponential damping times obtained in this way are plotted in Fig. 7.
Fig. 7 Exponential damping times for the different initial amplitudes, (A), and inhomogeneous layer widths, represented by diamond (⋄), triangle (△), square (□), and pentagon (), in increasing order. The stars (⋆) of the same colour as the symbols, represent the theoretically predicted damping time for the specific layer width and sinusoidal layer density profile. Analogously, the asterisks (E) represent the predicted damping for a linear layer density profile (Goossens et al. 2002). 

Open with DEXTER 
Note that for “wellbehaving” amplitude profiles (for A ≤ 0.01) we show τ_{d} obtained from Eq. (5). We also plotted the theoretical damping times due to resonant absorption for different inhomogeneous layer widths, for both linear and sinusoidal density profiles in the layer (Ruderman & Roberts 2002; Goossens et al. 2002). We find that, for the low amplitude cases (A ≤ 0.01) the simulated damping times are in between the predicted ones: the density profile in the inhomogeneous layer can greatly influence the resulting theoretical damping time, especially for thin layers (Soler et al. 2013). We would like to highlight three important characteristics of the nonlinear damping of kink oscillations from Fig. 7:

The damping time is no longer generally independent of amplitude, as in the linear case. There is a strong damping, even with initially thin inhomogeneous layers of the kink oscillations, for high enough amplitudes. In these cases, the damping time can be less than a third of the theoretically predicted one. Furthermore, the damping time appears to be saturated for amplitudes larger than some threshold value.

For thick enough inhomogeneous layers (in our case, for l = 0.5R), the damping time is independent of the amplitude, and thus coincidentally well approximated by the theory.

In the saturated high amplitude regime (here for A ≥ 0.035), the damping is only weakly dependent on the initial thickness of the inhomogeneous layer, and can thus not be used for seismological purposes.
For our initial conditions and measured oscillation periods, the time of transition between Gaussian and exponential damping profiles (Eq. (4)) is at t_{s} ≈ 400 s. In Fig. 8 we plotted the values for t_{s} obtained through the fitted Gaussian and exponential damping times.
Fig. 8 Switch time t_{s} (Eq. (4)) between the Gaussian and exponential damping, obtained through fitting as described in the article body. The symbols used for different layer widths are the same as in Fig. 7. The dashed line at 400 s is the theoretical value of t_{s}, for our parameters. 

Open with DEXTER 
As described above, we obtain τ_{d} and τ_{g} through fitting Eq. (5) for simulations with A ≤ 0.01. For the higher amplitude simulations, we fit the Gaussian of Eq. (5) to obtain τ_{g} (τ_{s} ≥ t_{f}) and Eq. (6) for τ_{d}. We can see that we recover the theoretical value for t_{s} well (within 10%) for A ≤ 0.01, but for higher amplitudes there are deviations, especially for the transitional amplitude of A = 0.02. Once again, this shows that we have to be careful when applying seismology on high amplitude, nonlinear oscillations.
Very recently, observational evidence has been produced showing that the damping time of transverse oscillations of coronal loops is a function of the initial displacement amplitude (Goddard & Nakariakov 2016). Unfortunately, a quantitative comparison of our results to observations is not possible given the large parameter space which would need to be covered. On the other hand, we can state that there are encouraging qualitative agreements between our Fig. 7 and Fig. 2 of Goddard & Nakariakov (2016): firstly, the damping time is decreasing as the initial amplitude of the oscillations increases, as is the case for our simulations with thin inhomogeneous layers. Secondly, this dependence weakens as the amplitude increases. The observational data also indicates that there is a considerable population between the values 1 and 2 of the damping time to the oscillation period ratio, which seem to be weakly dependent on the amplitude, similarly to our case of the loops with thick inhomogeneous layers. It is interesting to note that the transitional regime (defined earlier) appears to occur for much higher amplitudes for the observed oscillations than in our simulations. This could possibly imply a lower growth rate of the KHI for the observed loops.
4. Conclusions
In this study, we investigated the standing kink oscillations of a straight flux tube, aiming to model transverse oscillations of coronal loops. We ran 3D ideal MHD simulations, exploring the parameter space of initial amplitudes and inhomogeneous layer thicknesses. The chosen amplitudes cover the weakly linear and fully nonlinear regimes. The resulting oscillation periods are well described by the linear theory. Furthermore, for low amplitudes, we find that the resulting damping of the oscillation is close to the damping times computed by linear theory. However, for higher amplitudes, nonlinear effects have a definitive impact on the resulting oscillation characteristics, especially the development of the KHI around the loop edges where the velocity shear is highest, coinciding with the layer where resonant absorption is taking place. The development of KHI (and ultimately the threshold amplitude at which its effects become important) is dictated by its growth rate, which in turn depends on the ratio of loop radius to length, oscillation amplitude, Alfvén speeds, inhomogeneous layer thickness and numerical dissipation. We show that even for initially thin inhomogeneous layers, the oscillations undergo rapid damping due to the presence of KHI. For high amplitudes (with KHI developing in under one oscillation period), the damping time is almost independent of the initial thickness of the inhomogeneous layer. On the other hand, for thick inhomogeneous layers, the damping time does not seem to depend on the initial amplitude of the perturbation. To put things in perspective, the highest amplitude perturbation used in our simulations initially displaces the loop by less than its radius, and is thus at the lower boundary of observed flarerelated coronal loop oscillation displacements (Aschwanden et al. 2002). Studying the energy distribution in the antinode crosssection of the loop, we arrive at the conclusion that the kinetic energy in the inhomogeneous layer is converted to plasma internal energy more quickly in the presence of KHI. The increase in average internal energy is less than one percent, thus the energy budget of the wave is not enough to cause any significant temperature change in the 0.9 MK plasma. However, this heating might be significant for prominences (Antolin et al. 2015). Even if the dissipated wave energy is not enough to cause significant heating, we show that if the loop is surrounded by hotter plasma, mixing induced by the KHI can increase the average temperature of the loop. In the presence of KHI, the peak value of average kinetic energy deposited in the inhomogeneous layer is lower than in simulations without KHI. This is a consequence of the accelerated conversion of kinetic to internal energy in the presence of KHI, which cascades energy to smaller scales where it can be dissipated (by numerical dissipation) more efficiently. The disruption of the resonant layer may also contribute to the reduction in the peak value of average kinetic energy, though it is unclear how effective the resulting “patchy” resonant absorption is.
According to the present study, it becomes uncertain whether seismology schemes based on the linear theory for the damping rates of coronal loops are valid for highamplitude, nonlinear transverse oscillations. Using the observed switch between Gaussian and exponential damping profiles of transverse coronal loop oscillations for coronal seismology has recently been suggested. However, because of the nature of nonlinear kink oscillations it is questionable how accurately one can infer parameters such as inhomogeneous layer thickness and density ratio from the observed damping profiles, given that the nonlinear damping times are nearly insensitive to them. However, it is important that future studies also include Gaussian damping profiles in their analyses, as this profile seems to better describe the damping of nonlinear transverse oscillations of coronal loops. The new observations of amplitudedependent damping times qualitatively support our conclusions.
Acknowledgments
The authors would like to thank the referee for valuable comments which helped to improve the manuscript. N.M. acknowledges the Fund for Scientific ResearchFlanders (FWOVLaanderen). T.V.D. was supported by an Odysseus grant, the Belspo IAP P7/08 CHARM network and the GOA2015014 (KU Leuven). Inspiration for this research was found during ISSI and ISSIBJ workshops. Visualization was done with the help of VisIt software (Childs et al. 2012).
References
 Anfinogentov, S., Nisticò, G., & Nakariakov, V. M. 2013, A&A, 560, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Anfinogentov, S. A., Nakariakov, V. M., & Nisticò, G. 2015, A&A, 583, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Antolin, P., Yokoyama, T., & Van Doorsselaere, T. 2014, ApJ, 787, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Antolin, P., Okamoto, T. J., De Pontieu, B., et al. 2015, ApJ, 809, 72 [NASA ADS] [CrossRef] [Google Scholar]
 Arregui, I., Andries, J., Van Doorsselaere, T., Goossens, M., & Poedts, S. 2007, A&A, 463, 333 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aschwanden, M. J., & Schrijver, C. J. 2011, ApJ, 736, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Aschwanden, M. J., Fletcher, L., Schrijver, C. J., & Alexander, D. 1999, ApJ, 520, 880 [NASA ADS] [CrossRef] [Google Scholar]
 Aschwanden, M. J., de Pontieu, B., Schrijver, C. J., & Title, A. M. 2002, Sol. Phys., 206, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Berger, T. E., Slater, G., Hurlburt, N., et al. 2010, ApJ, 716, 1288 [NASA ADS] [CrossRef] [Google Scholar]
 Browning, P. K., & Priest, E. R. 1984, A&A, 131, 283 [NASA ADS] [Google Scholar]
 Childs, H., Brugger, E., Whitlock, B., et al. 2012, in High Performance VisualizationEnabling ExtremeScale Scientific Insight (Chapman & Hall/CRC Computational Science), 357 [Google Scholar]
 Edwin, P. M., & Roberts, B. 1983, Sol. Phys., 88, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Foullon, C., Verwichte, E., Nakariakov, V. M., Nykyri, K., & Farrugia, C. J. 2011, ApJ, 729, L8 [NASA ADS] [CrossRef] [Google Scholar]
 Goddard, C. R., & Nakariakov, V. M. 2016, A&A, 590, L5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goossens, M., Andries, J., & Aschwanden, M. J. 2002, A&A, 394, L39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goossens, M., Arregui, I., Ballester, J. L., & Wang, T. J. 2008, A&A, 484, 851 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Goossens, M., Soler, R., Terradas, J., Van Doorsselaere, T., & Verth, G. 2014, ApJ, 788, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Heyvaerts, J., & Priest, E. R. 1983, A&A, 117, 220 [NASA ADS] [Google Scholar]
 Hollweg, J. V. 1984, ApJ, 277, 392 [NASA ADS] [CrossRef] [Google Scholar]
 Hood, A. W., Ruderman, M., Pascoe, D. J., et al. 2013, A&A, 551, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ionson, J. A. 1978, ApJ, 226, 650 [NASA ADS] [CrossRef] [Google Scholar]
 Karpen, J. T., Dahlburg, R. B., & Davila, J. M. 1994, ApJ, 421, 372 [NASA ADS] [CrossRef] [Google Scholar]
 Keppens, R., Meliani, Z., van Marle, A. J., et al. 2012, J. Comput. Phys., 231, 718 [NASA ADS] [CrossRef] [Google Scholar]
 Magyar, N., & Van Doorsselaere, T. 2016, ApJ, 823, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Magyar, N., Van Doorsselaere, T., & Marcu, A. 2015, A&A, 582, A117 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Markwardt, C. B. 2009, in Astronomical Data Analysis Software and Systems XVIII, eds. D. A. Bohlender, D. Durand, & P. Dowler, ASP Conf. Ser., 411, 251 [Google Scholar]
 Nakariakov, V. M., Ofman, L., Deluca, E. E., Roberts, B., & Davila, J. M. 1999, Science, 285, 862 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Nisticò, G., Nakariakov, V. M., & Verwichte, E. 2013, A&A, 552, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ofman, L. 2005, Space Sci. Rev., 120, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Ofman, L. 2009a, Space Sci. Rev., 149, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Ofman, L. 2009b, ApJ, 694, 502 [NASA ADS] [CrossRef] [Google Scholar]
 Ofman, L., & Thompson, B. J. 2011, ApJ, 734, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Ofman, L., Davila, J. M., & Steinolfson, R. S. 1994, Geophys. Res. Lett., 21, 2259 [NASA ADS] [CrossRef] [Google Scholar]
 Okamoto, T. J., Antolin, P., De Pontieu, B., et al. 2015, ApJ, 809, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Pascoe, D. J., Wright, A. N., & De Moortel, I. 2011, ApJ, 731, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Pascoe, D. J., Hood, A. W., de Moortel, I., & Wright, A. N. 2012, A&A, 539, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pascoe, D. J., Hood, A. W., De Moortel, I., & Wright, A. N. 2013, A&A, 551, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pascoe, D. J., Goddard, C. R., Nisticò, G., Anfinogentov, S., & Nakariakov, V. M. 2016a, A&A, 589, A136 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pascoe, D. J., Goddard, C. R., Nisticò, G., Anfinogentov, S., & Nakariakov, V. M. 2016b, A&A, 585, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poedts, S., & Kerner, W. 1991, Phys. Rev. Lett., 66, 2871 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Poedts, S., Toth, G., Belien, A. J. C., & Goedbloed, J. P. 1997, Sol. Phys., 172, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Porth, O., Xia, C., Hendrix, T., Moschou, S. P., & Keppens, R. 2014, ApJS, 214, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Ruderman, M. S. 2006, Philos. Trans. Roy. Soc. London Ser. A, 364, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Ruderman, M. S., & Erdélyi, R. 2009, Space Sci. Rev., 149, 199 [NASA ADS] [CrossRef] [Google Scholar]
 Ruderman, M. S., & Goossens, M. 2014, Sol. Phys., 289, 1999 [NASA ADS] [CrossRef] [Google Scholar]
 Ruderman, M. S., & Roberts, B. 2002, ApJ, 577, 475 [NASA ADS] [CrossRef] [Google Scholar]
 Ruderman, M. S., Goossens, M., & Andries, J. 2010, Phys. Plasmas, 17, 082108 [NASA ADS] [CrossRef] [Google Scholar]
 Ryutova, M., Berger, T., Frank, Z., Tarbell, T., & Title, A. 2010, Sol. Phys., 267, 75 [NASA ADS] [CrossRef] [Google Scholar]
 Soler, R., Goossens, M., Terradas, J., & Oliver, R. 2013, ApJ, 777, 158 [NASA ADS] [CrossRef] [Google Scholar]
 Terradas, J., & Ofman, L. 2004, ApJ, 610, 523 [NASA ADS] [CrossRef] [Google Scholar]
 Terradas, J., Andries, J., Goossens, M., et al. 2008a, ApJ, 687, L115 [NASA ADS] [CrossRef] [Google Scholar]
 Terradas, J., Arregui, I., Oliver, R., et al. 2008b, ApJ, 679, 1611 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, T., Ofman, L., Davila, J. M., & Su, Y. 2012, ApJ, 751, L27 [NASA ADS] [CrossRef] [Google Scholar]
 Zaqarashvili, T. V., Zhelyazkov, I., & Ofman, L. 2015, ApJ, 813, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Zimovets, I. V., & Nakariakov, V. M. 2015, A&A, 577, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Figures
Fig. 1 Density in the antinode crosssection of the loop (z = 0) at t = 540 s for the different parameters used (A = 0.005 not shown for brevity, looking almost identical to A = 0.01). Different columns have different initial inhomogeneous layer widths, shown at the top, and the different rows are for different amplitudes, shown in the left margin. Axis units are in Mm and density in units of 10^{12} kg m^{3}. 

Open with DEXTER  
In the text 
Fig. 2 Density in the antinode crosssection of the loop (z = 0) at different times (in steps of approximately one halfperiod, shown in the upper part of each plot) for A = 0.035 and l = 0.33R. Axis units are in Mm and density in units of 10^{12} kg m^{3}. 

Open with DEXTER  
In the text 
Fig. 3 ycomponent of the velocity (v_{y}) at different times (selected times corresponding to plots in Fig. 2, shown in the upper part of each plot), for A = 0.035 and l = 0.33R. Axis units are in Mm and velocity in units of km s^{1}. 

Open with DEXTER  
In the text 
Fig. 4 Time evolution of different quantities in the antinode crosssection of the loop (z = 0) with A = 0.035 and l = 0.33R. Topleft: average kinetic energy (KE) density in the core and layer region of the loop. Topright: average internal energy density in the core and layer regions of the loop. Bottomleft: average temperature of the loop (average over both the core and layer). Bottomright: change in the area relative to the initial value, for the core region of the loop and the numerical domain outside the loop (outside both core and layer). 

Open with DEXTER  
In the text 
Fig. 5 Average kinetic energy density over time in the inhomogeneous layer of the antinode crosssection of the loop (z = 0). Left: l = 0.1R. Right: l = 0.5R. The different lines are for different initial amplitudes, from 0.01 to 0.05. 

Open with DEXTER  
In the text 
Fig. 6 Displacement amplitude over time in the loop antinode (z = 0). The values are normalized to the initial displacement. For the top plots (A = 0.01, left: l = 0.1R, right: l = 0.5R), the Gaussian and exponential profiles from Eq. (3) obtained by fitting Eq. (5) are plotted. For the bottom plots (A = 0.05, left: l = 0.1R, right: l = 0.5R), the Gaussian of Eq. (3) is fitted for the whole data and plotted, as well as the exponential resulting from fitting Eq. (6). 

Open with DEXTER  
In the text 
Fig. 7 Exponential damping times for the different initial amplitudes, (A), and inhomogeneous layer widths, represented by diamond (⋄), triangle (△), square (□), and pentagon (), in increasing order. The stars (⋆) of the same colour as the symbols, represent the theoretically predicted damping time for the specific layer width and sinusoidal layer density profile. Analogously, the asterisks (E) represent the predicted damping for a linear layer density profile (Goossens et al. 2002). 

Open with DEXTER  
In the text 
Fig. 8 Switch time t_{s} (Eq. (4)) between the Gaussian and exponential damping, obtained through fitting as described in the article body. The symbols used for different layer widths are the same as in Fig. 7. The dashed line at 400 s is the theoretical value of t_{s}, for our parameters. 

Open with DEXTER  
In the text 