Coronal loops are modeled as density structures in the solar atmosphere which act as waveguides for several types of magnetohydrodynamic (MHD) modes. One of the most readily detectable are standing kink oscillations, which causes periodic transverse perturbations of the loop axis. These oscillations were first detected using the Transition Region And Coronal Explorer (TRACE; Aschwanden et al., 1999, 2002; Nakariakov et al., 1999). The number of observations has greatly increased since the launch of the Atmospheric Imaging Assembly (AIA; Lemen et al., 2012) onboard the Solar Dynamics Observatory (SDO). Catalogs of observations have been produced by Zimovets and Nakariakov (2015) and Goddard et al. (2016). Kink oscillations attract a lot of attention due to their potential for seismological studies of the coronal plasma (e.g., reviews by De Moortel and Nakariakov, 2012; Stepanov et al., 2012; Pascoe, 2014; De Moortel et al., 2016), in particular to estimate the strength of the coronal magnetic field (e.g., Nakariakov et al., 1999; Nakariakov and Ofman, 2001; Van Doorsselaere et al., 2007; White and Verwichte, 2012; Pascoe et al., 2016; Sarkar et al., 2016).

Kink oscillations are observed to be strongly damped, only having a detectable amplitude for a few cycles, which is attributed to resonant absorption (e.g., Chen and Hasegawa, 1974). Coronal loops are modeled as having a higher density than the surrounding plasma, and for resonant absorption to take place only requires that the transition between the higher and lower density plasma is smooth. Inside this inhomogeneous layer, wave energy is transferred from kink to Alfvén modes where the local Alfvén speed matches the kink speed C_k. The timescale of this process is comparable to the period of oscillation (e.g., Hollweg and Yang, 1988; Goossens et al., 1992, 2002; Ruderman and Roberts, 2002).

Resonant absorption is a robust mechanism which occurs even in loops which are not cylindrically symmetric (Terradas et al., 2008b; Pascoe et al., 2011). Furthermore, numerical studies demonstrate that the Kelvin-Helmholtz instability (e.g., Terradas et al., 2008a; Antolin et al., 2015; Okamoto et al., 2015; Magyar and Van Doorsselaere, 2016; Hillier et al., 2019) is most efficient in loops with thin inhomogeneous layers. This instability leads to a mixing of plasma and effective widening of the inhomogeneous layer (e.g., Goddard et al., 2018; Karampelas and Van Doorsselaere, 2018) in addition to increased heating due to phase mixing (e.g., Heyvaerts and Priest, 1983; Karampelas et al., 2017; Pagano et al., 2018; Guo et al., 2019) which can account for observations of broad differential emission measures in coronal loops (Van Doorsselaere et al., 2018).

The transverse density profile can be described by two dimensionless parameters; the density contrast ratio being the ratio of the density at the center of the loop to the density far from it ζ = ρ₀/ρ_e, and the width of the inhomogeneous layer normalized to the minor loop radius ϵ = l/R. The damping rate due to resonant absorption depends on both of these parameters. For this reason the problem is underdetermined when trying to infer ζ and ϵ from the (exponential) damping time alone and some additional information is required (e.g., Goossens et al., 2008; Arregui and Asensio Ramos, 2014; Arregui and Goossens, 2019). Pascoe et al. (2013) produced a more accurate description of the damping profile due to resonant absorption, which includes the initial Gaussian damping regime of kink oscillations in addition to the later exponential damping regime. This damping profile is characterized by two damping times and so allows both ζ and ϵ to be estimated. This method was first applied by Pascoe et al. (2016) and later extended by Pascoe et al. (2017a,c) to include additional effects such as a time-dependent period of oscillation, the presence of additional parallel harmonics, and the use of Bayesian analysis (e.g., Arregui et al., 2013; Arregui, 2018) to improve the estimation of uncertainties. This seismological method requires both the Gaussian and exponential damping regimes to be accurately detected in the data and so depends on the oscillation data having a sufficiently high quality.

Another method for estimating ϵ is by forward modeling the appearance of the density profile for direct comparison with the transverse EUV intensity profile of the loop (Goddard et al., 2017; Pascoe et al., 2017b). It is possible to apply both of these methods simultaneously to observational data. This was recently demonstrated by Pascoe et al. (2018) for a loop for which the oscillation data alone was too noisy to allow strong constraint of ζ and ϵ. However, the spatial information from the intensity profile produced a strong constraint on ϵ, such that the oscillation data was only required to infer ζ when both methods were applied simultaneously. On the other hand, the value of ϵ ≈ 0.9 for this observation is significantly outside the validity of the thin boundary (TB) approximation. To correct for this effect, Pascoe et al. (2018) performed a narrow parametric study using the TB estimate as a starting point. The result of this study was a change in the estimated value of ζ from the TB value of 2.3 to a value of 2.8 based on numerical simulations for ϵ = 0.9. This case demonstrates the need for a seismological method that can account for the behavior of kink oscillations in loops with wide inhomogeneous layers without the need for separate studies and corrections applied afterwards. The use of a self-consistent seismological method is particularly important for future development of techniques for data analysis where multiple observational signatures are forward modeled and a systematic error arising from the TB approximation would have a deleterious influence on other observables. For example, the EUV intensity is I∝ρ² and so a change in inferred density contrast from 2.3 to 2.8 in the example above corresponds to an intensity change by a factor of approximately 1.5.

In this paper we study the behavior of kink oscillations of coronal loops for various transverse density profiles. Our aim is to provide a simple method of estimating the damping profile for a chosen profile which may be used for seismological investigations. The damping profiles for resonant absorption used in previous studies and this one are described in section 2. The results of our parametric study and the generation of a lookup table (LUT) are presented in section 3. In section 4 we present examples of the application of our method to synthetic test cases and observational data. Conclusions are presented in section 5.

Initial applications of resonant absorption to explain the strong damping of kink oscillations (Goossens et al., 1992, 2002; Ruderman and Roberts, 2002) considered only the asymptotic state of the system. The damping profile was an exponential of the form

where τ_d is the exponential damping time, P is the period of oscillation, κ = (ρ₀−ρ_e)/(ρ₀+ρ_e) and the factor 4/π² corresponds to a linear transition between ρ₀ and ρ_e. In this paper we consider a linear profile for the density in the inhomogeneous layer since it is the simplest smooth profile, can describe the widest range of possible structures, and is the only profile for which the analytical solution for all times is presently known (see discussion section 6.2 of Pascoe et al., 2018).

Numerical simulations by Pascoe et al. (2012) demonstrated that the damping profile of strongly damped propagating kink oscillations is more accurately described by a Gaussian damping profile rather than an exponential one. The existence of these two regimes was reconciled by the analytical description derived by Hood et al. (2013) which described the damping profile for all times and demonstrates that Gaussian and exponential profile can be obtained in the limits of small and large time, respectively. The derivation by Hood et al. (2013) was performed for propagating kink waves with the damping rate expressed in terms of damping length scales but we may consider the case of standing kink waves with a damping time using the long wavelength limit λ = C_kP, giving

We note that the relationship between damping length scales (propagating waves) and timescales (standing waves) has been demonstrated explicitly for the exponential regime (e.g., comparing the derivations of Goossens et al., 2002; Terradas et al., 2010) but presently not the Gaussian regime. Nonetheless we demonstrate the applicability of this relationship (proposed by Pascoe et al., 2010) by comparison with our numerical simulations. The applicability of this relationship has also previously been demonstrated in numerical simulations by Ruderman and Terradas (2013) and Magyar and Van Doorsselaere (2016).

Pascoe et al. (2013) proposed a general damping profile (GDP) that combined both of these damping regimes into a single approximation. This is described in terms of a switch from the Gaussian damping profile that applies at the start of the oscillation to the exponential profile which applies later.

where the switch from the Gaussian to exponential damping regime occurs at t_s and As=D(t=ts).

The above damping profiles are also based on the thin tube approximation. In this limit the period of the fundamental standing kink mode is

where L is the loop length, and the kink speed for a low-β plasma (uniform magnetic field) is

where C_A0 is the internal Alfvén speed. The thin tube thin boundary (TTTB) approximation for the period of oscillation of kink modes therefore depends on ζ but not ϵ (Goossens et al., 2008). However, the parametric studies by Van Doorsselaere et al. (2004) and Soler et al. (2014) find that P_k does depend on ϵ outside of the applicability of the TB approximation.

To illustrate the different damping profiles, Figure 1 shows the results of our numerical simulations (described in section 3) for three values of ϵ with ζ = 2. For this value of density contrast the GDP suggests a switch from the Gaussian to exponential damping regime at t = 3P. For kink oscillations in low density contrast loops such as these, the Gaussian damping profile (blue curves) provides a much better description than the exponential damping profile (red curves), and the general damping profile (green curves) further improves the description for later times. As expected, all three analytical profiles become poorer as ϵ increases and the TB approximation becomes less appropriate.

In this paper we wish to characterize the damping behavior of kink oscillations as accurately as possible, and without prior assumption of the form of the damping profile. The plus symbols in Figure 1 represent the amplitudes A_i which we use to characterize the oscillation. These amplitudes are defined at every half cycle of the oscillation. The LUT damping profile (dashed lines) is constructed from these amplitudes by spline interpolation. The dashed lines indicate that this method allows us to accurately describe the damping of the kink oscillation, albeit at a cost of requiring more information. The exponential and Gaussian damping profiles can each be characterized by a single parameter, i.e., the damping time τ_d or τ_g. The GDP combines both the Gaussian and exponential damping regimes and so is characterized by both these damping times (with the switch time given in terms of these two parameters in Equation 6). For the six cycles indicated in Figure 1 our interpolation method uses 13 parameters, or more generally 2n+1 parameters for n cycles. On the other hand, this number is still sufficiently small that the results of hundreds of numerical simulations can be compiled into a lookup table of minimal size.

For each of the three simulations in Figure 1, and the additional simulations presented in section 3, the amplitudes A_i used to describe the LUT damping profile were found by a least-squares fit to the data using MPFIT (Markwardt, 2009). Simulation data covered at least six cycles of the oscillation, with a greater number for simulations with smaller values of ζ or ϵ for which the damping is weaker and so there are a greater number of cycles at an observable amplitude. This is also intended to ensure that the exponential damping regime is present in the data, which occurs later for smaller ζ. This allows the damping profile to be extrapolated beyond the number of simulated cycles using an exponential profile with the last measured damping time. We require that the LUT damping profile function is able to return results for arbitrary values of time in order for it to appropriately handle requests from the user or from a fitting routine transparently i.e., without needing to take the details of individual simulations into account.

In this section we study the behavior of kink oscillations in coronal loops for various combinations of ζ and ϵ. Soler et al. (2014) and Van Doorsselaere et al. (2004) performed similar numerical studies investigating the damping time for the exponential regime. These studies demonstrated that the thin boundary approximation produces poor estimates of the damping behavior in loops with wider inhomogeneous layers.

As in Soler et al. (2014), we consider L/R = 100 which is typical for observations of standing kink oscillations. Weak dependence of the damping on the longitudinal wavenumber k_z has been demonstrated by Van Doorsselaere et al. (2004). Numerical simulations were performed using a Lax-Wendroff code to solve the linear MHD equations in cylindrical coordinates (r, θ, z) for m = 1 symmetry corresponding to kink oscillations (and the Alfvén waves generated by resonant absorption). The magnetic field is constant and aligned with the z-direction. This code was previously used in Pascoe et al. (2012, 2013, 2015) and Hood et al. (2013) to study the spatial damping of propagating kink waves and here is applied to the case of the temporal damping of standing kink waves by appropriate choice of boundary and initial conditions. The boundary conditions are line-tied to simulate the loop footpoints being fixed in the photosphere, while the boundary in the r-direction is placed sufficiently far away to not affect the results. The fundamental longitudinal harmonic of the standing kink mode is excited by a perturbation to the radial and azimuthal velocities of the form

In the following simulations, the numerical domain covers r = [0, 10R] and z = [0, 100R], with a resolution of 1, 000 × 1, 000 grid points. Each of our 300+ simulations took approximately 1 h to run using 80 × 2.8 GHz processor cores. We study the standing mode by considering the variation in the amplitude of the transverse velocity measured at the loop apex z = L/2.

Figure 2 shows the simulations performed in our parametric study which were used to generate the first version of our lookup table. The solid, dashed, and dotted lines correspond to the damping rates shown in the right panel. These damping profiles are based on the TB approximation for the Gaussian damping profile and so will not accurately describe the behavior for large ϵ but serve as an indication of the range of parameters we are interested in considering with respect to observational studies. The curves demonstrate that we are not equally interested in all regions of the ζ-ϵ parameter space. Large amplitude standing kink oscillations are typically observed for fewer than 6 cycles (e.g., Goddard and Nakariakov, 2016) and so we are mainly interested in parameters which produce this level of damping. The dotted, solid, and dashed lines correspond to τ_g/P = 1, 2, and 5, respectively, and indicate the parameters we expect to be relevant to observations. For large values of both ζ and ϵ, kink oscillations would be very strongly damped and hence unlikely to be reliably detected. The shorter time series available would also generally limit the seismological information that could be obtained. Unlike ϵ which has the defined upper limit of 2, there is no upper limit for ζ. However, the damping rate is asymptotic in this limit and so we can consider a reasonable upper limit, which is taken to be 7 in this study but may be extended in the future (e.g., high contrast filament threads considered by Arregui et al., 2008). For small values of both ζ and ϵ the oscillations would be weakly damped. Such oscillations are not typically observed, although we are still interested in the behavior for small values of ϵ with regard to checking convergence of our results to the analytical profiles based on the TB approximation.

The distribution of the simulations we have performed reflects these areas of interest. When compiling the results of these simulations into a lookup table, it is therefore not convenient to use a 2D array to describe the ζ-ϵ plane. Instead, the results are considered as a list of scattered positions in the parameter space using the IDL routine GRIDDATA which can be used to interpolate our simulation results to return the damping profile for an arbitrary value of ζ and ϵ. Additional simulations were performed for testing purposes, including generating synthetic observational data used in section 4. These additional simulations are not included version 1.0 LUT but may be in future applications of the LUT to actual observational data. The use of a LUT and interpolation methods for scattered data allows the method to be improved over time by incorporating additional results as they are obtained. Other examples of a LUT strategy for solar applications include the CHIANTI emission database (Del Zanna et al., 2015) used as part of EUV forward modeling codes such as FoMo (Van Doorsselaere et al., 2016), and the CAISAR code for inversions of solar Ca II spectra (Beck et al., 2015).

The LUT and the corresponding IDL code are available at https://github.com/djpascoe/kinkLUT. The routine requires as input the values of ϵ, ζ, and the normalized times t_n = t/P at which the damping profile is desired. The value of each amplitude A_i is determined by 2D interpolation of the simulation results using the IDL routine GRIDDATA. (In this paper we use the linear method, chosen as the simplest method with fewest assumptions, for which requested grid points are linearly interpolated from triangles formed by Delaunay triangulation. These triangles were constructed with the TRIANGULATE routine and are included in the LUT save file). The damping profile is then returned by spline interpolation of these amplitudes A_i for the user-requested times t_n. This procedure can be used within a forward modeling function used for comparing a model with data. For example, a simple model for an oscillation with a single harmonic and no background trend, as considered in this paper, is

where A₀ is the initial amplitude of a sinusoidal oscillation with period P and phase shift ϕ, and the damping profile DLUT is based on our lookup table. We demonstrate the results of such a method in section 4.

In this section we demonstrate the application of our LUT as a seismological tool to use the observed damping of a kink oscillation to infer information about the transverse density profile of the oscillating loop.

Figure 8 shows the results of a test of our method for a kink oscillation simulated in a loop with ζ = 2.15 and ϵ = 0.75. This data point is not included in version 1.0 of our LUT used in the following analysis. The top panel shows the analyzed oscillation which includes uniformly distributed random noise with a maximum amplitude of 5% of the initial kink oscillation amplitude. The middle and bottom panels show 2D histograms approximating the marginalized posterior probability density function for ζ and ϵ based on 10⁵ Markov chain Monte Carlo (MCMC) samples of the GDP and LUT models, respectively (see also Pascoe et al., 2017a, 2018). This data comes from a simulation with an inhomogeneous layer width that is sufficiently large for the TB approximation to produce inaccurate results. The GDP approach overestimates the value of ϵ, and correspondingly underestimates the value of ζ.

The seismological estimates in Figure 8 demonstrate that it is important to consider the dependence of the damping on ζ and ϵ, aside from the quality of the observational data. Due to the asymptotic nature of the inversion curve, the extent to which ζ and ϵ are constrained depends on whether the density contrast is high or low; near the low-ζ asymptote ζ is well-constrained but ϵ is poorly constrained, and vice versa for the high-ζ asymptote. This is demonstrated in Figure 8 where the GDP inversion underestimates the actual value of ζ = 2.15 and produces very low estimates of ζ. The red bars indicate the 95% credible intervals and show the value of ζ is well-constrained (but excludes the actual value of ζ due to the systematic error from the TB approximation). In comparison, the inversion results using the LUT are less constrained for ζ but the maximum a posteriori probability (MAP) estimate of ζ ≈ 2.11 is consistent with the actual value for the synthetic data. The nature of these constraints makes 2D histograms such as those in Figure 9 the simplest way of representing the available information for the transverse density profile. Quoting confidence intervals alone can be misleading if the asymptotic behavior is not kept in mind. Furthermore, 1D histograms for ζ and ϵ would generally not be well described using a normal distribution (they have been approximated by the exponentially modified Gaussian function by Pascoe et al., 2017a) and so estimates of uncertainties based on assuming this distribution would not be accurate.

Figure 9 shows the results of a test of our method for a loop with a thin boundary layer ζ = 3.75 and ϵ = 0.15 (this data point is again not included in version 1.0 of our LUT). In this limit we see both the GDP and LUT inversions producing good estimates of the actual value of ϵ. The asymptotic behavior means ζ is far less constrained than ϵ, and less constrained for the LUT results since the inferred value of ζ is larger (it is underestimated using the GDP).

In this paper we have performed a parametric study for the dependence of the damping of standing kink oscillations by resonant absorption on the density contrast ratio ζ and the width of the inhomogeneous layer ϵ. Previous studies (Van Doorsselaere et al., 2004; Soler et al., 2014) have demonstrated the inaccuracy of the classical analytical description for the damping rate τ_d, which is based on the TB approximation, when applied to loops with large values of ϵ. We have expanded on these studies by developing a description of the damping profile which includes the initial non-exponential regime, and summarized our results into a lookup table which may be used to produce improved estimates of loop parameters.

The method presented in this paper makes analysis simpler in terms of not requiring an additional step to correct for the effect of a wide inhomogeneous layer, as in Pascoe et al. (2018), and automatically propagates the corresponding uncertainties. Being self-consistent also makes the method more suitable than the TB approximation as a starting point to include other effects which depend on the density profile parameters. On the other hand, this work only accounts for the influence of wide inhomogeneous layers on the kink oscillation damping profile, with other assumptions and approximations needing to be considered.

The simulations performed in this paper use a linear profile for the density in the inhomogeneous layer. In terms of the damping rate due to resonant absorption, the key physical parameter is the gradient of the Alfvén speed at the resonance location. Arregui and Goossens (2019) discuss how the choice of density profile (e.g., linear or sinusoidal in the inhomogeneous layer) affects the seismologically-inferred loop parameters, since they have different constants of proportionality (e.g., Goossens et al., 2002; Roberts, 2008; Soler et al., 2014) between the damping time and the loop parameters (ζ, ϵ). Arregui and Goossens (2019) note that seismological applications have so far used the linear profile and that choosing a different profile would return a different answer. However, presently no such choice exists since the linear profile is the only one for which a full analytical solution is known (Hood et al., 2013; Ruderman and Terradas, 2013). The linear profile is therefore the natural choice of profile for use in this study, to allow a comparison of numerical results with the TB approximations in the limit ϵ → 0. Pascoe et al. (2018) estimate that the difference in the seismological results using the linear and sinusoidal density profiles is ≲10%. However, these two profiles are not the only possible choices. If some profile (other than the linear one used) should be shown to be more appropriate for coronal loops then it would be simple (though computationally expensive) to reproduce our parametric study for that profile. If the spatial resolution of loops is sufficiently high then it is possible to test different profiles by forward modeling the EUV profile (e.g., Aschwanden et al., 2007; Goddard et al., 2017; Pascoe et al., 2017b, 2018). Such methods still typically involve further assumptions or approximations (e.g., loops having a cylindrical cross-section and azimuthal symmetry) however the aim of these investigations should not be considered the impossible task of inferring the precise density profile of coronal loops but rather estimating the most appropriate characteristic scales that influence the relevant physical processes. Resonant absorption, phase mixing, and the Kelvin-Helmholtz instability all depend on the transverse loop structure. In this context it is not required that the coronal loop density profile is exactly linear for the analysis based on this profile (or any other simplified model) to provide accurate and useful estimates of parameters such as the size of the transverse inhomogeneity l = ϵR.

Our LUT is also based on results of numerical simulations for the linear regime and so excludes the effect of the Kelvin-Helmholtz instability (KHI). The current technique is therefore more suitable for loops with wider inhomogeneous layers in which KHI develops at a slower rate. Simulations by Magyar and Van Doorsselaere (2016) suggest KHI is weak for ϵ≳0.5 and the weaker development for larger values of ϵ is also evident in simulations by Goddard et al. (2018) and Pagano et al. (2018). In principle the same LUT strategy could be extended to include effects such as KHI by performing a similar parametric study with a third parameter being the amplitude of the oscillation, A₀, with the study in this paper representing the ζ-ϵ plane in the limit A₀ → 0. However this wider parametric study (ζ, ϵ, A₀) where each data point is based a 3D nonlinear simulation would be orders of magnitude more computationally expensive than our current one. A wider parametric study which describes the evolution of the loop, e.g., the development of fine structure by KHI, would also require additional observable signatures to be modeled and compared with data. In the current study we focus on the shape of the damping profile due to resonant absorption. Additional observational signatures have recently been considered by Antolin et al. (2017) based on forward modeling techniques (Van Doorsselaere et al., 2016) to describe the appearance of loops observed simultaneously in multiple EUV channels.

Our parametric study was also performed for thin loops which satisfy the long wavelength approximation k_zR≪1. As a first approximation the effect of a large k_z may be considered in terms of reducing the kink mode phase speed from C_k, but otherwise not influencing the damping per cycle (Van Doorsselaere et al., 2004). Variations in the loop length in time can account for a time-dependent period of oscillation (e.g., Nisticò et al., 2013; Russell et al., 2015; Pascoe et al., 2017c) but should also remain well-approximated by considering the damping per cycle provided by our LUT as unaffected. The period of oscillation (and ratios of periods for different longitudinal harmonics) can also be affected by longitudinal structuring due to gravity (e.g., Andries et al., 2005). Cooling of coronal loops causes the scale height to vary in time, accompanied by decreases in the period and period ratio (Morton and Erdélyi, 2009) and an amplification of kink oscillation which acts against the damping due to resonant absorption (Ruderman, 2011). Numerical simulations suggest this effect is typically small but could be approximately 10% for larger amplitude oscillations (Magyar et al., 2015). Similarly, the period of oscillation and period ratio can be affected by expansion of the loop at the apex (Verth and Erdélyi, 2008) and a time-dependent expansion can also reduce the damping of kink oscillations (Shukhobodskiy et al., 2018). It may be possible to include effects such as these which modify the period of oscillation relative to P_k and introduce additional damping/amplification terms into a more general oscillation model alongside our LUT damping profile for resonant absorption. However, including loop evolution which affects the transverse structuring of the loop, such as KHI, would be more challenging.

DP and TVD contributed conception and design of the study. AH and DP developed the code for the numerical simulations. DP wrote the first draft of the manuscript, and received input from all co-authors.