© ESO, 2014

1. Introduction

Since the first observations of transverse oscillations of coronal loops by the Transition Region And Coronal Explorer (TRACE) in 1998 (Aschwanden et al. 1999; Nakariakov et al. 1999), theorists have been striving to improve the simple models of coronal loops as one-dimensional single magnetic flux tubes, developed in the late seventies (e.g. Edwin & Roberts 1983), to more realistic models for coronal loops. Building better models of coronal loop oscillations leads to better diagnostics of the local plasma parameters by coronal seismology (see e.g. the review papers of Nakariakov & Verwichte 2005; Goossens 2008; Verwichte et al. 2009; De Moortel & Nakariakov 2012).

The observed rapid damping of the oscillations (Schrijver et al. 2002; Aschwanden et al. 2002) is thought to be caused by resonant absorption. The theory of plasma heating by resonant absorption has been developed by Chen & Hasegawa (1974) and was applied to coronal loop oscillations for instance by Hollweg & Yang (1988), Sakurai et al. (1991), and Goossens et al. (1995). Comparing the observed wave periods and damping times with theoretical models of magnetohydrodynamic (MHD) waves, it has been possible to estimate the magnetic field strength (Nakariakov & Ofman 2001), the thickness of the inhomogeneous layer (Goossens et al. 2002), the density contrast between the loop and background plasma (Aschwanden et al. 2003), and the density scale height (Andries et al. 2005), amongst others. Arregui et al. (2007) and Goossens et al. (2008) showed that due to the underdetermination of the system of governing equations, it is not possible to define a unique equilibrium state given the observed values of wave periods and damping times. Instead, the possible equilibrium states define a one-dimensional curve in a three-dimensional space. If more information about the density contrast between the loop and the surrounding plasma were present, Arregui & Asensio Ramos (2011) have shown how to constrain the estimates for the magnetic field strength, thickness of inhomogeneous layer, and density contrast using Bayesian inference on the observed damping times and frequencies. Bayesian inference also provides confidence intervals for these quantities.

Because of the multitude and detail of recent observations of transverse loop oscillations using the Coronal Multichannel Polarimeter (CoMP; Tomczyk et al. 2007) and especially the Atmospheric Imaging Assembly on board the Solar Dynamics Observatory (AIA/SDO; Aschwanden & Schrijver 2011; Wang et al. 2012; White & Verwichte 2012; Nisticò et al. 2013), it is now possible to find the distribution of equilibrium parameters using a statistical approach, either using the observed frequencies and damping times (Verwichte et al. 2013), or the displacement vectors themselves (see Asensio Ramos & Arregui 2013). Alternatively, Arregui et al. (2013) used the ratio of the period of the fundamental mode to the first overtone to estimate the parameters for density stratification along the magnetic field lines and magnetic field divergence, and formulated how a Bayesian model comparison can determine which of these two models is most plausible for a given measured period ratio and the associated uncertainty.

Observations show that often a coronal loop does not oscillate in isolation, but as part of a system of coronal loops that exhibits collective behaviour (Verwichte et al. 2004). Therefore, it is desirable to model the interaction within such a system. The oscillations of a system of two homogeneous parallel magnetic cylinders have first been studied numerically by Luna et al. (2008) and later analytically by Van Doorsselaere et al. (2008) in the long-wavelength limit. More recently, Luna et al. (2009, 2010) extended the model to systems of more than two coronal loops using the T-matrix method (Waterman 1969). This method has first been used in the context of solar physics to model the interaction of waves with sunspot structures (Bogdan & Zweibel 1985; Bogdan & Cattaneo 1989; Keppens et al. 1994) and describes the action of different magnetic cylinders on one another.

The effects of damping in systems of coronal loops have also been studied. Arregui et al. (2008) were the first to study the damped oscillations of a system of two coronal loops using Cartesian geometry. Later, Terradas et al. (2008) studied the damped oscillations of a multistranded loop using a 2D numerical code, and Ofman (2009) modelled a system of four parallel coronal loops using a 3D numerical code. The effects of line-of-sight integration on the estimate of the energy content and the wave mode identification in configurations of multiple loops have been discussed by De Moortel & Pascoe (2012).

The results of simulations of initial value problems should be interpreted as a superposition of eigenmodes. Therefore it is useful to study analytically the damped eigenmodes of a system of coronal loops. Keppens et al. (1994) already developed analytic expressions for resonantly damped oscillations in the context of the interaction of sunspots with acoustic waves using the T-matrix method. Their results were adapted to coronal conditions by Soler (2010), who solved the resulting system of equations numerically. Recently, Robertson & Ruderman (2011) developed a fully analytic model for the damped oscillations of a system of two identical parallel coronal loops.

This paper aims to generalise the results of Robertson & Ruderman (2011) to a system of two not necessarily identical parallel loops. The density contrast between the two loops is an important parameter that characterises the interaction between the loops; hence this parameter cannot be ignored. We also investigate the difference between the damping rates of so-called standard and anomalous systems (Van Doorsselaere et al. 2008). Finally, we perform an extensive parametric study to examine the dependence of the damping rate on separation distance, thickness of inhomogeneous layer, and relative loop density in detail and provide explanations. Our analysis gives a quantitative point of view to complement the results obtained previously by Soler (2010).

The paper is organised as follows: in the next section, the equations for the damping rate for a system of two arbitrary parallel coronal loops are derived. Sections 2.1 and 2.2 briefly review the common foundations of the studies of Robertson & Ruderman (2011) and ours for completeness, while Sect. 2.3. is concerned with the generalisation to arbitrary two-loop systems. In Sect. 3 we study the dependence of the damping rate on the equilibrium parameters. Our conclusions are summarised in Sect. 4.

Fig. 1 Sketch of the equilibrium configuration (from Robertson & Ruderman 2011).

2. Derivation of the damping rate

2.1. Equilibrium configuration and governing equations

The system of two coronal loops is modelled as a pair of parallel magnetic cylinders of length L and radii R_L and R_R, as seen in Fig. 1. The loops consist of a uniform plasma with densities ρ_L and ρ_R with transition layers of width l_L and l_R, respectively, in which the density adjusts itself continuously and monotonically to the value ρ_e of the exterior plasma. We adopt the cold plasma approximation. A constant magnetic field B lies parallel to the loop axes; we align a Cartesian coordinate system such that the z-axis points in this direction. We study the problem in bicylindrical (σ,τ,z) coordinates. A sketch of the bicylindrical coordinate system is given in Fig. 2. This coordinate system has two points x =  ± a such that all coordinate surfaces σ = const. pass through these two points. The relation between Cartesian and bicylindrical coordinates is given by $\mathit{x}\mathrm{=}\frac{\mathit{a}\sinh \mathit{\tau }}{\cosh \mathit{\tau }\mathrm{-}\cos \mathit{\sigma }}\mathit{,} \mathit{y}\mathrm{=}\frac{\mathit{a}\sin \mathit{\sigma }}{\cosh \mathit{\tau }\mathrm{-}\cos \mathit{\sigma }}\mathrm{·}$(1)The tube boundaries can be expressed by τ =  − τ_L and τ = τ_R, where τ_L,R > 0. The tube radii and the distance d between the two tube centres can be written in terms of the parameters of the coordinate system as ${\mathit{R}}_{\mathrm{L}}\mathrm{=}\frac{\mathit{a}}{\sinh {\mathit{\tau }}_{\mathrm{L}}}\mathit{,} {\mathit{R}}_{\mathrm{R}}\mathrm{=}\frac{\mathit{a}}{\sinh {\mathit{\tau }}_{\mathrm{R}}}\mathit{,} \mathit{d}\mathrm{=}\mathit{a}\mathrm{(}\coth {\mathit{\tau }}_{\mathrm{L}}\mathrm{+}\coth {\mathit{\tau }}_{\mathrm{R}}\mathrm{)}\mathit{.}$(2)

Fig. 2 Sketch of a plane z = const. in the bicylindrical (σ,τ) coordinate system. The lines σ = const. are circles passing through the points x =  ± a, the lines τ = const. are nested circles around these points, whose radius tends to zero and whose centre tends to ±a as |τ| → ∞ (from Van Doorsselaere et al. 2008).

Two geometric limits are relevant throughout the paper. When |τ| → ∞, the circles τ = const., which are nested around the points x =  ± a, will become smaller and their centres will converge to the points x =  ± a. When |τ| ≪ 1, the two circles become very large (filling the entire left and right half-plane in the limit τ → 0). A layer defined by τ₀ < τ < τ₀ + l will become very asymmetric for low values of τ₀. The parameters l_L and l_R are essentially parameters of the coordinate system. A physical quantity related to l_L and l_R is the mean thickness of the inhomogeneous layer, defined in bicylindrical coordinates by ${\mathit{l}}_{\mathrm{avg,i}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{R}}_{\mathrm{i}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\hspace{0.17em}\frac{{\mathit{l}}_{\mathrm{i}}\mathrm{\left(}\mathit{\sigma }\mathrm{\right)}\hspace{0.17em}{\mathit{a}}^{\mathrm{2}}}{\mathrm{(}\cosh {\mathit{\tau }}_{\mathrm{i}}\mathrm{-}\cos \mathit{\sigma }{\mathrm{)}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{d}\mathit{\sigma }\mathit{.}$(3)The integration in (3) can be carried out and yields ${\mathit{l}}_{\mathrm{avg,i}}\mathrm{=}{\mathit{l}}_{\mathrm{i}}{\mathit{R}}_{\mathrm{i}}\coth {\mathit{\tau }}_{\mathrm{i}} \stackrel{{\mathit{\tau }}_{\mathrm{i}}\mathrm{\to }\mathrm{\infty }}{\mathrm{-}\mathrm{\to }} {\mathit{l}}_{\mathrm{i}}\mathit{R}i$(4)(Robertson & Ruderman 2011), where the subscript i can be L or R. We impose that the ratios of the average thickness of the inhomogeneous layers to the tube radii are small parameters, ${\mathit{\epsilon }}_{\mathrm{L}}\mathrm{:}\mathrm{=}\frac{{\mathit{l}}_{\mathrm{avg,L}}}{{\mathit{R}}_{\mathrm{L}}}\mathrm{\ll }\mathrm{1}\mathit{,} {\mathit{\epsilon }}_{\mathrm{R}}\mathrm{:}\mathrm{=}\frac{{\mathit{l}}_{\mathrm{avg,R}}}{{\mathit{R}}_{\mathrm{R}}}\mathrm{\ll }\mathrm{1}\mathit{,}$(5)in order to conserve the ring-like shape of the inhomogeneous layers. This is the thin boundary (TB) assumption. The TB assumption breaks down when τ_L,R is not large enough.

Averaging over the different directions has consequences for the damping properties of the two-loop system, especially for small loop distances. We address them in Sect. 3.1.

Following Ruderman & Roberts (2002) and Robertson & Ruderman (2011) we have used the linearised viscous MHD equations to describe the plasma motions. They are given by $\begin{array}{ccc}\mathit{\rho }\frac{{\mathit{\partial }}^{\mathrm{2}}{\xi }}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{\mu }}_{\mathrm{0}}}\hspace{0.17em}\mathrm{(}\mathrm{\nabla }\mathrm{\times }{b}\mathrm{)}\mathrm{\times }{B}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial t}}\mathrm{\left[}\mathrm{\nabla }\mathrm{\right(}\overline{\mathit{\nu }}\hspace{0.17em}\mathrm{\nabla }\mathrm{·}{\xi }\mathrm{)}\mathrm{-}\mathrm{\nabla }\mathrm{\times }\mathrm{(}\overline{\mathit{\nu }}\hspace{0.17em}\mathrm{\nabla }\mathrm{\times }{\xi }\mathrm{\left)}\mathrm{\right]}\mathit{,}\\ {b}& \mathrm{=}& \mathrm{\nabla }\mathrm{\times }\mathrm{(}{\xi }\mathrm{\times }{B}\mathrm{)}\mathit{.}\end{array}$Here $\overline{\mathit{\nu }}$ denotes the coefficient of shear viscosity, ξ is the Lagrangian displacement vector, μ₀ the permeability of free space, and b denotes the perturbation of the magnetic field. Viscosity is important only in a narrow dissipative layer about the Alfvén resonant surface. This is the magnetic surface where the equation Re ω = ω_A is satisfied, where ω_A is the Alfvén frequency defined by ${\mathit{\omega }}_{\mathrm{A}}\mathrm{=}{\mathit{k}}_{\mathit{z}}{\mathit{v}}_{\mathrm{A}}\mathrm{=}\mathit{B}{\mathit{k}}_{\mathit{z}}\mathit{/}\sqrt{\mathit{\mu \rho }}$, v_A is the Alfvén speed, and k_z the parallel wavenumber. In previous studies that used cylindrical coordinates, the damping rate has been found to be independent of $\overline{\mathit{\nu }}$ when this parameter is small (e.g. Poedts & Kerner 1991; Goossens et al. 1992; Ruderman & Roberts 2002) because of the ideal nature of the damping mechanism. We can assume a priori that this is also the case in the two-loop system. Hence it does not matter whether we use the viscous or resistive MHD equations to describe resonant damping. The eigenfunctions and location of heating are dependent on the dissipative mechanism (Van Doorsselaere et al. 2007), but we are not concerned with these problems in this paper.

Finally, to obtain analytical progress, we assume that the loop length L is much larger than the size of the system in the transverse direction d (the long-wavelength approximation): $\frac{\mathit{d}}{\mathit{L}}\mathrm{\ll }\mathrm{1.}$(8)Here we will study the eigenvalue problem. The boundary condition at the ends of the loops is the frozen-in condition at the photosphere at z =  ± L/2, that is ${\xi }\mathrm{=}\mathrm{0}\mathit{,} \mathit{z}\mathrm{=}\mathrm{\pm }\frac{\mathit{L}}{\mathrm{2}}\mathrm{·}$(9)This quantises the axial wavenumber k_z as ${\mathit{k}}_{\mathit{z}}\mathrm{=}\frac{\mathit{n\pi }}{\mathit{L}}$. The boundary conditions across the layers is given in the next subsection.

2.2. Boundary conditions across the dissipative layer

We are not interested in the behaviour of the solutions in the dissipative layer itself. By solving the linearised MHD equations separately in the core and outer region and in the dissipative layer, we can match the local solutions to obtain boundary conditions across the dissipative layer in the form of jump conditions. This approach has often been used in the study of resonant absorption in cylindrical geometry, for example by Sakurai et al. (1991), and Goossens et al. (1995) for the driven problem and Tirry & Goossens (1996) for the eigenvalue problem, and has recently been adapted to a bicylindrical coordinate system by Robertson & Ruderman (2011). Hence we limit ourselves in this subsection to citing the major results needed in the remainder of the paper.

Van Doorsselaere et al. (2008) have computed the eigenmodes of the two-loop system analytically in the long-wavelength limit if we replace the inhomogeneous layers by discontinuous jumps in density from ρ_L,R to ρ_e. We introduce the dispersion function $\begin{array}{ccc}{\mathit{D}}_{\mathrm{0}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{:}\mathrm{=}{\mathit{F}}^{\mathrm{2}}{\mathit{\omega }}^{\mathrm{4}}\mathrm{(}{\mathit{\rho }}_{\mathrm{L}}\mathrm{-}{\mathit{\rho }}_{\mathrm{e}}\mathrm{)}\mathrm{(}{\mathit{\rho }}_{\mathrm{R}}\mathrm{-}{\mathit{\rho }}_{\mathrm{e}}\mathrm{)}& & \mathrm{-}\mathrm{[}\mathrm{(}{\mathit{\rho }}_{\mathrm{L}}\mathrm{+}{\mathit{\rho }}_{\mathrm{e}}\mathrm{)}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{\rho }}_{\mathrm{e}}{\mathit{v}}_{\mathrm{Ae}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}\mathrm{]}\\ & & \mathrm{\times }\mathrm{[}\mathrm{(}{\mathit{\rho }}_{\mathrm{R}}\mathrm{+}{\mathit{\rho }}_{\mathrm{e}}\mathrm{)}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}\mathrm{2}{\mathit{\rho }}_{\mathrm{e}}{\mathit{v}}_{\mathrm{Ae}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}\mathrm{]}\end{array}$(10)to write the dispersion relation of Van Doorsselaere et al. (2008) as ${\mathit{D}}_{\mathrm{0}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\mathrm{0}$(11)with solutions ${\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}\mathrm{\{}{\mathit{\zeta }}_{\mathrm{L}}\mathrm{+}{\mathit{\zeta }}_{\mathrm{R}}\mathrm{+}\mathrm{2}\mathrm{\pm }\sqrt{\mathrm{(}{\mathit{\zeta }}_{\mathrm{L}}\mathrm{-}{\mathit{\zeta }}_{\mathrm{R}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{4}\mathrm{(}{\mathit{\zeta }}_{\mathrm{L}}\mathrm{-}\mathrm{1}\mathrm{\left)}\mathrm{\right(}{\mathit{\zeta }}_{\mathrm{R}}\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{F}}^{\mathrm{2}}}\mathrm{\}}}{\mathrm{(}{\mathit{\zeta }}_{\mathrm{L}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}{\mathit{\zeta }}_{\mathrm{R}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{(}{\mathit{\zeta }}_{\mathrm{L}}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{(}{\mathit{\zeta }}_{\mathrm{R}}\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{F}}^{\mathrm{2}}}\mathrm{·}$(12)In these equations ζ_L,R = ρ_L,R/ρ_e are the normalised densities and F = exp [−(τ_L + τ_R)] is a geometric factor.

The normal modes of the homogeneous two-loop system have previously been studied numerically by Luna et al. (2008). They found four eigenfrequencies that can be split into two pairs of frequencies that lie very close to one another, instead of the two eigenfrequencies described by Eq. (12). The four different eigenmodes can be classified as symmetric (S) or antisymmetric (A), and the displacement can be either in the x- or the y-direction. The lowest frequencies are found for the S_x and A_y eigenmodes. The long-wavelength assumption ignores the dispersion between the two eigenmode pairs and reduces the number of normal modes of the system from four to two.

Note that ω₊ is a monotonically decreasing function of τ_L + τ_R and ω₋ is a monotonically increasing function of τ_L + τ_R, with limits $\begin{array}{ccc}& & \underset{{\mathit{\tau }}_{\mathrm{L}}\mathrm{+}{\mathit{\tau }}_{\mathrm{R}}\mathrm{\to }\mathrm{\infty }}{\lim }\hspace{0.17em}{\mathit{\omega }}_{\mathrm{+}}\mathrm{(}{\mathit{\tau }}_{\mathrm{L}}\mathrm{+}{\mathit{\tau }}_{\mathrm{R}}\mathrm{)}\mathrm{=}\begin{array}{c}\max \end{array}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{kL}}\mathit{,}{\mathit{\omega }}_{\mathrm{kR}}\mathrm{\right)}\mathit{,}\\ & & \underset{{\mathit{\tau }}_{\mathrm{L}}\mathrm{+}{\mathit{\tau }}_{\mathrm{R}}\mathrm{\to }\mathrm{\infty }}{\lim }\hspace{0.17em}{\mathit{\omega }}_{\mathrm{-}}\mathrm{(}{\mathit{\tau }}_{\mathrm{L}}\mathrm{+}{\mathit{\tau }}_{\mathrm{R}}\mathrm{)}\mathrm{=}\begin{array}{c}\min \end{array}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{kL}}\mathit{,}{\mathit{\omega }}_{\mathrm{kR}}\mathrm{\right)}\mathit{.}\end{array}$(13)Here ω_ki is the kink frequency defined by ${\mathit{\omega }}_{\mathrm{ki}}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{\rho }}_{\mathrm{i}}{\mathit{\omega }}_{\mathrm{Ai}}^{\mathrm{2}}\mathrm{+}{\mathit{\rho }}_{\mathrm{e}}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}}{{\mathit{\rho }}_{\mathrm{i}}\mathrm{+}{\mathit{\rho }}_{\mathrm{e}}}\mathit{,}$(14)where again the subscript i can be either L or R. Equation (13) shows that for large loop separations, the system of coronal loops becomes decoupled. The collective eigenfrequencies tend to the kink eigenfequencies of the individual loops. A consequence of the limits (13) is that for two nonidentical loops (by convention, we set ζ_L < ζ_R in this case) a region in the plane determined by ζ_L/ζ_R and τ_L + τ_R (i.e. density contrast and distance between the loops) always exists such that ω₋ < ω_AL (Van Doorsselaere et al. 2008). Systems of two coronal loops in which this situation occurs are called anomalous.

The introduction of a thin inhomogeneous layer leads to a continuous Alfvén spectrum. At the Alfvén resonant magnetic surface, defined by Re ω = ω_A, the ideal MHD equations become “quasi-singular” (Soler et al. 2013). For the two-loop system, we calculated that just as in cylindrical geometry, the eigenfrequencies are complex and hence the resonant position does not lie on the real axis in the complex plane. This means there are no actual singularities in the ideal MHD equations. In resistive MHD, the quasi-singularities are lifted and the quasi-modes found in ideal MHD become proper resistive eigenmodes.

When a smooth initial perturbation is applied to the system, each surface resonates at its own Alfvén frequency. Energy pumped into the inhomogeneous layer is converted into local Alvén wave energy propagating along the magnetic field lines. After a while, the neighbouring magnetic flux surfaces drift out of phase with one another, creating steep gradients on small length scales. These gradients would grow unboundedly in time in ideal MHD, but in resistive MHD the Ohmic dissipation term will dominate locally due to the decrease of the magnetic Reynolds number (phase mixing, see e.g. Heyvaerts & Priest 1983). This converts the wave energy into heat. The global phenomenon of transfer of incident wave energy into energy of local Alfvén waves and its subsequent dissipation on small scales due to finite resistivity is what is meant by resonant absorption (Andries 2003).

The qualitative description of resonant absorption in coronal cylinders (loops, jets, etc.) is independent of the coordinate system and has amassed a substantial literature (starting with Sakurai et al. 1991; Goossens et al. 1995). Improvements on the theory of resonant absorption consist mainly of extensions of the theory to more realistic equilibria and the effects on measurable quantities relevant to coronal seismology such as wave periods and damping times.

Recently, the applicability of the “classical theory of resonant damping” (M. Ruderman, priv. comm.) to the initial value problem as outlined above has been challenged by Pascoe et al. (2013) and Hood et al. (2013) for propagating kink waves and by Ruderman & Terradas (2013) for standing kink modes. The time (resp. length) scales required to create the large gradients near the resonant layer are of the order of the wave period (wavelength) times Re^1/3, where Re is the magnetic Reynolds number. Because Re ≫ 1 in the solar corona, these scales are much larger than the typical damping times (lengths) of the oscillations.

It was found numerically by Pascoe et al. (2013) and analytically by Hood et al. (2013) that the initial spatial damping profile of propagating kink oscillations is better described by a Gaussian damping profile. After a few wavelengths, the exponential damping profile predicted by the classical theory dominates. The same results were found by Ruderman & Terradas (2013) for temporal damping profiles of standing kink waves. These results indicate that the classical theory of resonant absorption underestimates the damping time (length) of the oscillations. The error increases with increasing thickness of the inhomogeneous layer. Ruderman & Terradas (2013) found that for thin boundary layers such as considered here, the error is still lower than 20% for l/R = 0.3. This justifies our use of the classical theory for the two-loop system.

Since ω₋ lies outside of the Alfvén spectrum of the left loop in anomalous systems, there is only a resonant surface for ω₊ in this loop. In all other cases resonant surfaces exist for both eigenfrequencies in both loops.

In the outer and core regions we can use an ideal MHD description. Van Doorsselaere et al. (2008) calculated that the displacement vector and total pressure perturbation in bicylindrical coordinates there is given by $\begin{array}{ccc}& & {\mathit{P}}^{\mathrm{\prime }}\mathrm{=}\mathrm{\Theta }\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\cos \hspace{0.17em}\mathrm{(}\mathit{\sigma }\mathrm{-}{\mathit{\sigma }}_{\mathrm{0}}\mathrm{)}\mathrm{-}\mathrm{\Theta }\mathrm{\left(}\mathrm{0}\mathrm{\right)}\cos {\mathit{\sigma }}_{\mathrm{0}}\mathit{,}\end{array}$(15)$\begin{array}{ccc}& & \mathrm{\Theta }\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \\ \\ \end{array}\end{array}$(16)$\begin{array}{ccc}& & {\mathit{\xi }}_{\mathit{\tau }}\mathrm{=}\mathit{\xi ̂}\mathit{\tau }\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{(}\cosh \mathit{\tau }\mathrm{-}\cos \mathit{\sigma }\mathrm{)}\cos \mathrm{(}\mathit{\sigma }\mathrm{-}{\mathit{\sigma }}_{\mathrm{0}}\mathrm{)}\mathit{,}\end{array}$(17)$\begin{array}{ccc}& & \mathit{\xi ̂}\mathit{\tau }\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \\ \\ \end{array}\end{array}$(18)in which C_L,C₁,C₂ and C_R are constants to be determined by the boundary conditions across the inhomogeneous layers.

In the dissipative layer we use the viscous MHD equations. Robertson & Ruderman (2011) found that, just as in cylindrical coordinates, the τ-dependence of the total pressure Θ can be taken constant, $\mathrm{\Theta }\mathrm{=}{\mathrm{\Theta }}_{\mathrm{A}}\mathit{,}$(19)without further loss of accuracy. For the τ-component of the displacement vector we have $\mathit{\xi ̂}\mathit{\tau }\mathrm{=}\mathrm{-}\frac{{\mathrm{\Theta }}_{\mathrm{A}}}{\mathit{a}{\mathit{\rho }}_{\mathrm{A}}\mathrm{\Delta }}\hspace{0.17em}{\mathit{G}}_{\mathrm{\Lambda }}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{+}\mathit{C,}$(20)where C is an integration constant and the function ${\mathit{G}}_{\mathrm{\Lambda }}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\hspace{0.17em}\frac{{\mathit{e}}^{\mathrm{-}{\mathit{u}}^{\mathrm{3}}\mathit{/}\mathrm{3}}}{\mathit{u}}\mathrm{\left\{}\exp \mathrm{\right(}\mathit{iu\tau }\hspace{0.17em}\mathrm{sgn}\mathrm{\left(}\mathrm{\Delta }\mathrm{\right)}\mathrm{)}\mathrm{+}\mathrm{\Lambda }\mathit{u}\mathrm{-}\mathrm{1}\mathrm{\}}\hspace{0.17em}\mathrm{d}\mathit{u}$(21)also appears in studies of eigenmodes of dissipative MHD in Cartesian (Ruderman et al. 1995) and cylindrical (Tirry & Goossens 1996) geometry. These functions depend on the density ρ_A and the gradient of Alfvén frequency squared | Δ | at the resonant surface. In the eigenvalue problem the former only depends on the real part of the eigenfrequency of the two-loop system by definition of the resonant surface. Here ρ_A and | Δ | are given by $\begin{array}{ccc}\mathrm{\Delta }& \mathrm{=}& \\ {\mathit{\rho }}_{\mathrm{A}}& \mathrm{=}& \end{array}$The jump of any quantity f across the dissipative layer is defined as $\mathrm{\left[}\mathit{f}\mathrm{\right]}\mathrm{=}\underset{\mathit{h}\mathrm{\to }\mathrm{0}}{\lim }\hspace{0.17em}\mathrm{\left\{}\mathit{f}\mathrm{\right(}{\mathit{\tau }}_{\mathrm{i}}\mathrm{+}\mathit{h}\mathrm{)}\mathrm{-}\mathit{f}\mathrm{(}{\mathit{\tau }}_{\mathrm{i}}\mathrm{-}\mathit{h}\mathrm{\left)}\mathrm{\right\}}\mathit{.}$(24)Asymptotic matching of the local ideal and viscous equations leads to the following jump conditions across the dissipative layer: for standard systems $\begin{array}{ccc}\mathrm{\left[}\mathrm{\Theta }\mathrm{\right]}& \mathrm{=}& \mathrm{0}\mathit{,}\\ \mathrm{\left[}\mathit{\xi ̂}\mathit{\tau }\mathrm{\right]}& \mathrm{=}& \mathrm{-}\frac{\mathit{i\pi }\mathrm{\Theta }\mathrm{\left(}{\mathit{\tau }}_{\mathrm{A}}\mathrm{\right)}}{\mathit{a}{\mathit{\rho }}_{\mathrm{A}}\mathrm{\left|}\mathrm{\Delta }\mathrm{\right|}}\end{array}$holds for both loops. In anomalous systems the absence of a resonant surface for the lower eigenfrequency in the left loop means that only one resonant surface in the left loop remains, defined by ω₊ = ω_AL. We can use the same formulae for standard and anomalous systems if we redefine the jump condition (26) for anomalous systems as follows: $\mathrm{\left[}\mathit{\xi ̂}\mathit{\tau ,}\mathrm{L,}\mathrm{-}\mathrm{\right]}\mathrm{=}\mathrm{0}\mathit{,} \mathrm{\left[}\mathit{\xi ̂}\mathit{\tau ,}\mathrm{L,}\mathrm{+}\mathrm{\right]}\mathrm{=}\mathrm{-}\frac{\mathit{i\pi }\mathrm{\Theta }\mathrm{\left(}{\mathit{\tau }}_{\mathrm{AL}}\mathrm{\right)}}{\mathit{a}{\mathit{\rho }}_{\mathrm{A}}\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{L}}\mathrm{\right|}}\mathit{,}\mathrm{\left[}\mathit{\xi ̂}\mathit{\tau ,}\mathrm{R}\mathrm{\right]}\mathrm{=}\mathrm{-}\frac{\mathit{i\pi }\mathrm{\Theta }\mathrm{\left(}{\mathit{\tau }}_{\mathrm{AR}}\mathrm{\right)}}{\mathit{a}{\mathit{\rho }}_{\mathrm{A}}\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{R}}\mathrm{\right|}}\mathrm{·}$(27)

2.3. Dispersion relation and damping rate

We use the jump conditions (25), (26), and (27) to derive an expression for the damping rate γ for two loops with arbitrary densities. One approach is to integrate these jump conditions across a nonuniform equilibrium to obtain boundary conditions across the entire inhomogeneous layer. This leads to a description in the form of principal value integrals (Robertson & Ruderman 2011).

Another, simpler, approach can be motivated as follows: the dissipative solutions (19) and (20) are local in the sense that they were obtained using a Taylor expansion of the Alfvén frequency at the resonant layer. If we assume that the density profile is almost linear in the inhomogeneous layer, Eqs. (19) and (20) can be considered valid across the entire inhomogeneous layer. Then the jump conditions (25), (26), and (27) can be used across this layer as well.

We do not lose any generality compared to almost all previous studies using the TB assumption, which assumed linear or sinusoidal density profiles in the nonuniform layer. Furthermore, in cylindrical coordinates, the only effect of the density profile is a constant factor appearing in solutions for the damping rate (see Eq. (77) of Goossens et al. 2002, which needs to be corrected for a typo: the factor $\mathrm{(}{\mathit{\omega }}_{\mathrm{Ai}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}\mathrm{)}$ needs to be squared). This expression only depends in fact on the gradient of the Alfvén speed at the resonance through the factor 1/Δ. One needs determine a posteriori whether this dependence still holds in bicylindrical coordinates. It is important to note that the density profile does matter when we consider thick inhomogeneous layers (Soler et al. 2013).

A second argument in favour of our approach is that ultimately, retention of the principal value integrals in the calculations only leads to a phase shift of the eigenfunctions and does not influence the imaginary part γ of the eigenfrequency (Robertson & Ruderman 2011). By dropping the principal value integral from the start, we can calculate more general expressions for the damping rate than were found by Robertson & Ruderman (2011).

We focus on the case of standard systems. The derivation for standard and anomalous systems is similar; the results for anomalous systems are given at the end of the derivation. The jump conditions across the inhomogeneous layer are $\begin{array}{ccc}\mathit{\xi ̂}\mathit{\tau }\mathrm{(}\mathrm{-}{\mathit{\tau }}_{\mathrm{L}}\mathrm{)}\mathrm{-}\mathit{\xi ̂}\mathit{\tau }\mathrm{(}\mathrm{-}{\mathit{\tau }}_{\mathrm{L}}\mathrm{-}{\mathit{l}}_{\mathrm{L}}\mathrm{)}\mathrm{=}\mathrm{\left[}\mathit{\xi ̂}\mathit{\tau ,}\mathrm{L}\mathrm{\right]}& \mathrm{=}& \mathrm{-}\frac{\mathit{i\pi }\mathrm{\Theta }\mathrm{\left(}{\mathit{\tau }}_{\mathrm{AL}}\mathrm{\right)}}{\mathit{a}{\mathit{\rho }}_{\mathrm{A}}\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{L}}\mathrm{\right|}}\mathit{,}\\ \mathrm{\Theta }\mathrm{(}\mathrm{-}{\mathit{\tau }}_{\mathrm{L}}\mathrm{-}{\mathit{l}}_{\mathrm{L}}\mathrm{)}\mathrm{-}\mathrm{\Theta }\mathrm{(}\mathrm{-}{\mathit{\tau }}_{\mathrm{L}}\mathrm{)}& \mathrm{=}& \mathrm{0}\mathit{,}\\ \mathit{\xi ̂}\mathit{\tau }\mathrm{(}{\mathit{\tau }}_{\mathrm{R}}\mathrm{+}{\mathit{l}}_{\mathrm{R}}\mathrm{)}\mathrm{-}\mathit{\xi ̂}\mathit{\tau }\mathrm{\left(}{\mathit{\tau }}_{\mathrm{R}}\mathrm{\right)}\mathrm{=}\mathrm{\left[}\mathit{\xi ̂}\mathit{\tau ,}\mathrm{R}\mathrm{\right]}& \mathrm{=}& \mathrm{-}\frac{\mathit{i\pi }\mathrm{\Theta }\mathrm{\left(}{\mathit{\tau }}_{\mathrm{AR}}\mathrm{\right)}}{\mathit{a}{\mathit{\rho }}_{\mathrm{A}}\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{R}}\mathrm{\right|}}\mathit{,}\\ \mathrm{\Theta }\mathrm{(}{\mathit{\tau }}_{\mathrm{R}}\mathrm{+}{\mathit{l}}_{\mathrm{R}}\mathrm{)}\mathrm{-}\mathrm{\Theta }\mathrm{\left(}{\mathit{\tau }}_{\mathrm{R}}\mathrm{\right)}& \mathrm{=}& \mathrm{0.}\end{array}$Equations (28)–(31) are the appropriate boundary conditions to be used in Eqs. (16) and (18). Since Θ is constant across the annuli, we can take Θ_AL = Θ(−τ_L − l_L) = C₁e^− τ_L + C₂e^τ_L in the left annulus and Θ_AR = Θ(τ_R + l_R) = C₁e^τ_R + C₂e^− τ_R in the right annulus. This yields the following system of equations: $\begin{array}{ccc}{\mathit{C}}_{\mathrm{L}}{\mathit{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{L}}\mathrm{-}{\mathit{l}}_{\mathrm{L}}}& \mathrm{=}& {\mathit{C}}_{\mathrm{1}}{\mathit{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{L}}}\mathrm{+}{\mathit{C}}_{\mathrm{2}}{\mathit{e}}^{{\mathit{\tau }}_{\mathrm{L}}}\mathit{,} \\ {\mathit{C}}_{\mathrm{R}}{\mathit{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{R}}\mathrm{-}{\mathit{l}}_{\mathrm{R}}}& \mathrm{=}& {\mathit{C}}_{\mathrm{1}}{\mathit{e}}^{{\mathit{\tau }}_{\mathrm{R}}}\mathrm{+}{\mathit{C}}_{\mathrm{2}}{\mathit{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{R}}}\mathit{,} \\ \frac{{\mathit{C}}_{\mathrm{1}}{\mathit{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{L}}}\mathrm{-}{\mathit{C}}_{\mathrm{2}}{\mathit{e}}^{{\mathit{\tau }}_{\mathrm{L}}}}{{\mathit{\rho }}_{\mathrm{e}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}\mathrm{)}}\mathrm{-}\frac{{\mathit{C}}_{\mathrm{L}}{\mathit{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{L}}\mathrm{-}{\mathit{l}}_{\mathrm{L}}}}{{\mathit{\rho }}_{\mathrm{L}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{AL}}^{\mathrm{2}}\mathrm{)}}& \mathrm{=}& \frac{\mathrm{-}\mathit{i\pi }\mathrm{(}{\mathit{C}}_{\mathrm{1}}{\mathit{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{L}}}\mathrm{+}{\mathit{C}}_{\mathrm{2}}{\mathit{e}}^{{\mathit{\tau }}_{\mathrm{L}}}\mathrm{)}}{{\mathit{\rho }}_{\mathrm{A}}\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{L}}\mathrm{\right|}}\mathit{,} \\ \frac{\mathrm{-}{\mathit{C}}_{\mathrm{R}}{\mathit{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{R}}\mathrm{-}{\mathit{l}}_{\mathrm{R}}}}{{\mathit{\rho }}_{\mathrm{R}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{AR}}^{\mathrm{2}}\mathrm{)}}\mathrm{-}\frac{{\mathit{C}}_{\mathrm{1}}{\mathit{e}}^{{\mathit{\tau }}_{\mathrm{R}}}\mathrm{-}{\mathit{C}}_{\mathrm{2}}{\mathit{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{R}}}}{{\mathit{\rho }}_{\mathrm{e}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}\mathrm{)}}& \mathrm{=}& \frac{\mathrm{-}\mathit{i\pi }\mathrm{(}{\mathit{C}}_{\mathrm{1}}{\mathit{e}}^{{\mathit{\tau }}_{\mathrm{R}}}\mathrm{+}{\mathit{C}}_{\mathrm{2}}{\mathit{e}}^{\mathrm{-}{\mathit{\tau }}_{\mathrm{R}}}\mathrm{)}}{{\mathit{\rho }}_{\mathrm{A}}\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{R}}\mathrm{\right|}}\mathrm{·}\end{array}$This system has a nontrivial solution for C₁,...,C₄ only if its coefficient matrix has a determinant equal to zero. We now introduce the notation $\begin{array}{ccc}\mathit{L}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{L}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{{\mathit{\omega }}_{\mathrm{AL}}^{\mathrm{2}}}^{\mathrm{)}} \mathit{F}\mathrm{=}{\mathit{e}}^{\mathrm{-}\mathrm{(}{\mathit{\tau }}_{\mathrm{L}}\mathrm{+}{\mathit{\tau }}_{\mathrm{R}}\mathrm{)}}\\ \mathit{E}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{e}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}}^{\mathrm{)}} {\mathit{I}}_{\mathrm{L}}\mathrm{=}\frac{\mathit{i\pi }}{{\mathit{\rho }}_{\mathrm{A}}\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{L}}\mathrm{\right|}}\\ \mathit{R}& \mathrm{=}& {\mathit{\rho }}_{\mathrm{R}}\mathrm{(}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{{\mathit{\omega }}_{\mathrm{AR}}^{\mathrm{2}}}^{\mathrm{)}} {\mathit{I}}_{\mathrm{R}}\mathrm{=}\frac{\mathit{i\pi }}{{\mathit{\rho }}_{\mathrm{A}}\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{R}}\mathrm{\right|}}\end{array}$(36)to write the dispersion relation for ω² as $\begin{array}{ccc}{\mathit{D}}_{\mathrm{0}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}& \mathrm{=}& \mathrm{-}\mathit{ER}{\mathit{I}}_{\mathrm{R}}\mathrm{\left[}{\mathit{F}}^{\mathrm{2}}\mathrm{\right(}\mathit{L}\mathrm{-}\mathit{E}\mathrm{)}\mathrm{+}\mathrm{(}\mathit{L}\mathrm{+}\mathit{E}\mathrm{\left)}\mathrm{\right]}\\ & & \mathrm{-}\mathit{EL}{\mathit{I}}_{\mathrm{L}}\mathrm{\left[}{\mathit{F}}^{\mathrm{2}}\mathrm{\right(}\mathit{R}\mathrm{-}\mathit{E}\mathrm{)}\mathrm{+}\mathrm{(}\mathit{R}\mathrm{+}\mathit{E}\mathrm{\left)}\mathrm{\right]}\mathrm{+}{\mathit{E}}^{\mathrm{2}}\mathit{LR}{\mathit{I}}_{\mathrm{L}}{\mathit{I}}_{\mathrm{R}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{F}}^{\mathrm{2}}\mathrm{)}\mathit{,}\end{array}$(37)where D₀(ω) was introduced in (10). By Eqs. (36), the final term on the right-hand side of (37) is proportional to 1/|Δ_L|·1/|Δ_R |. For quasi-linear density profiles, we see in Eq. (52) that this term is proportional to ε_L·ε_R and hence can be ignored with respect to the other terms, which are of the order of ε_L,R. We can then isolate the terms ε_L,R in the right-hand side of Eq. (37) and write the right-hand side compactly as εD_ε(ω). This yields the dispersion relation ${\mathit{D}}_{\mathrm{1}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}$(38)where ${\mathit{D}}_{\mathrm{1}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{=}{\mathit{D}}_{\mathrm{0}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{-}\mathit{\epsilon }{\mathit{D}}_{\mathit{\epsilon }}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathit{.}$(39)Since the introduction of a thin inhomogeneous layer only represents a weak perturbation to the system, it is logical to search for solutions in the neighbourhood of ω_±. Hence we develop D₁(ω) in a Taylor series about ω = ω_±: $\begin{array}{ccc}{\mathit{D}}_{\mathrm{1}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathrm{\approx }{\mathit{D}}_{\mathrm{1}}\mathrm{(}{\mathit{\omega }}_{\mathrm{\pm }}\mathrm{+}\mathit{\epsilon \delta \omega }\mathrm{)}& \mathrm{=}& {\mathit{D}}_{\mathrm{0}}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{\pm }}\mathrm{\right)}\mathrm{+}\mathit{\epsilon }{\mathit{D}}_{\mathit{\epsilon }}\mathrm{\left(}{\mathit{\omega }}_{\mathrm{\pm }}\mathrm{\right)}\\ & & \mathrm{+}\mathit{\epsilon }\hspace{0.17em}\mathit{\delta \omega }{\left[\frac{\mathit{\partial }{\mathit{D}}_{\mathrm{0}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}}{\mathit{\partial \omega }}\right]}_{\mathit{\omega }\mathrm{=}{\mathit{\omega }}_{\mathrm{\pm }}}\mathrm{+}\mathit{O}\mathrm{(}{\mathit{\epsilon }}^{\mathrm{2}}\mathrm{)}\mathit{.}\end{array}$(40)Using Eqs. (11) and (12) for the frequencies of the dissipationless system in Eq. (40) yields $\mathit{\delta \omega }\mathrm{=}\frac{{\mathit{D}}_{\mathit{\epsilon }}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}}{\mathrm{\left[}\mathit{\partial }{\mathit{D}}_{\mathrm{0}}\mathrm{\right(}\mathit{\omega }\mathrm{)}\mathit{/}\mathit{\partial \omega }{\mathrm{]}}_{\mathit{\omega }\mathrm{=}{\mathit{\omega }}_{\mathrm{\pm }}}}\mathrm{·}$(41)Equations (36) and (37) imply that D_ε(ω) is purely imaginary. This means that δω is also imaginary, hence we write δω = iγ, where γ is the damping rate we aim to determine. In conclusion, introducing a thin inhomogeneous layer only introduces damping of the eigenmodes and does not result in a frequency shift of the oscillations.

After some algebra, we obtain the following expression of the damping rate for the transverse oscillations in a system of two parallel coronal loops: $\begin{array}{ccc}\frac{{\mathit{\gamma }}_{\mathrm{\pm }}}{{\mathit{\omega }}_{\mathrm{\pm }}}& \mathrm{=}& \frac{\mathrm{-}\mathit{\pi }\mathrm{(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}\mathrm{)}}{\mathrm{\mp }\mathrm{4}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{e}}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{A}}\sqrt{\mathrm{(}{\mathit{\zeta }}_{\mathrm{L}}\mathrm{-}{\mathit{\zeta }}_{\mathrm{R}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{4}\mathrm{(}{\mathit{\zeta }}_{\mathrm{L}}\mathrm{-}\mathrm{1}\mathrm{\left)}\mathrm{\right(}{\mathit{\zeta }}_{\mathrm{R}}\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{F}}^{\mathrm{2}}}}\\ & & \hspace{0.17em}\\ & & \mathrm{\times }(\frac{{\mathit{\rho }}_{\mathrm{R}}\mathrm{(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{AR}}^{\mathrm{2}}\mathrm{)}\mathrm{\left[}\mathrm{\right(}\mathrm{1}\mathrm{+}{\mathit{F}}^{\mathrm{2}}\mathrm{\left)}{\mathit{\rho }}_{\mathrm{L}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{AL}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{F}}^{\mathrm{2}}\mathrm{\left)}{\mathit{\rho }}_{\mathrm{e}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}}{\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{R}}\mathrm{\right|}}\\ & & \hspace{0.17em}\\ & & \mathrm{+}\frac{{\mathit{\rho }}_{\mathrm{L}}\mathrm{(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{AL}}^{\mathrm{2}}\mathrm{)}\mathrm{\left[}\mathrm{\right(}\mathrm{1}\mathrm{+}{\mathit{F}}^{\mathrm{2}}\mathrm{\left)}{\mathit{\rho }}_{\mathrm{R}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{AR}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{F}}^{\mathrm{2}}\mathrm{\left)}{\mathit{\rho }}_{\mathrm{e}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}}{\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{L}}\mathrm{\right|}})\mathrm{·}\end{array}$(42)The two last terms describe the contribution of the resonances in both loops (through the factors 1/|Δ_L,R |) to the total damping rate. The meaning of these terms is discussed in detail in Sect. 3. Formally, the equations for the damping rate in the left loop in an anomalous system can be obtained from the equations of a standard system by letting | Δ_L| → ∞. For future reference, we will write down the equation for the damping rate in anomalous systems, which is $\begin{array}{ccc}\frac{{\mathit{\gamma }}_{\mathrm{-}\mathit{,}\mathrm{a}}}{{\mathit{\omega }}_{\mathrm{-}}}& \mathrm{=}& \frac{\mathrm{-}\mathit{\pi }\mathrm{(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}\mathrm{)}}{\mathrm{\mp }\mathrm{4}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{e}}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{A}}\sqrt{\mathrm{(}{\mathit{\zeta }}_{\mathrm{L}}\mathrm{-}{\mathit{\zeta }}_{\mathrm{R}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{4}\mathrm{(}{\mathit{\zeta }}_{\mathrm{L}}\mathrm{-}\mathrm{1}\mathrm{\left)}\mathrm{\right(}{\mathit{\zeta }}_{\mathrm{R}}\mathrm{-}\mathrm{1}\mathrm{)}{\mathit{F}}^{\mathrm{2}}}}\\ & & \hspace{0.17em}\\ & & \mathrm{\times }\frac{{\mathit{\rho }}_{\mathrm{R}}\mathrm{(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{AR}}^{\mathrm{2}}\mathrm{)}\mathrm{\left[}\mathrm{\right(}\mathrm{1}\mathrm{+}{\mathit{F}}^{\mathrm{2}}\mathrm{\left)}{\mathit{\rho }}_{\mathrm{L}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{AL}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{F}}^{\mathrm{2}}\mathrm{\left)}{\mathit{\rho }}_{\mathrm{e}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}}{\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{R}}\mathrm{\right|}}\mathrm{·}\end{array}$(43)Expression (42) generalises several results found previously in studies of two-loop oscillations. Firstly, if we remove dissipation, this means we set I_L = 0,I_R = 0 in Eq. (37), then substituting the relations (36) in this reduced dispersion relation leads back to Eq. (11) found by Van Doorsselaere et al. (2008).

Secondly, we consider what happens when the tubes are identical. Such systems exhibit the standard behavior independent of the separation between the loops. The two terms in the second factor of (42) are then also identical. After some algebra, we find that the expression for the damping rate becomes ${\mathit{\gamma }}_{\mathrm{\pm }}\mathrm{=}\frac{\mathrm{-}\mathit{\pi }{\mathit{\omega }}_{\mathit{Ae}}^{\mathrm{4}}{\mathit{\rho }}_{\mathrm{e}}\mathrm{(}\mathit{\zeta }\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{F}}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right(}\mathrm{1}\mathrm{\pm }\mathit{F}\mathrm{)}}{\mathrm{2}{\mathit{\omega }}_{\mathrm{\pm }}{\mathit{\rho }}_{\mathrm{A}}\mathrm{\left|}\mathrm{\Delta }\mathrm{\right|}\mathrm{(}\mathit{\zeta }\mathrm{+}\mathrm{1}\mathrm{\mp }\mathrm{(}\mathit{\zeta }\mathrm{-}\mathrm{1}\mathrm{)}\mathit{F}{\mathrm{)}}^{\mathrm{3}}}\mathrm{·}$(44)This equation is identical to Eq. (82) of Robertson & Ruderman (2011), up to a sign, since these authors write ω ≈ ω_± − iγ, while we used ω ≈ ω_± + iγ.

Finally, we consider the case d/R → ∞ in which the tubes are placed far away from each other. This situation corresponds mathematically to τ_L + τ_R → ∞ or F → 0. Even though the loop centres converge to x =  ± a and hence remain at a finite distance from one another, compared to their radius (which is the typical length scale of the loops), the loops will diverge infinitely far apart from each other.

We note that in the long-wavelength approximation and using bicylindrical coordinates, it is in fact impossible to realise a configuration of loops that are far apart without shrinking to a zero thickness. This is because the Eqs. (2) imply $\frac{{\mathit{R}}_{\mathrm{L,R}}}{{\mathit{L}}_{\mathrm{L,R}}}\mathrm{=}\frac{\mathit{ak}}{\mathit{n\pi }\sinh \mathrm{\left(}{\mathit{\tau }}_{\mathrm{L,R}}\mathrm{\right)}}\mathit{,}$(45)while the long-wavelength assumption together with (2) implies that $\frac{\mathit{d}}{\mathit{L}}\mathrm{=}\frac{\mathit{a}\mathrm{(}\coth {\mathit{\tau }}_{\mathrm{L}}\mathrm{+}\coth {\mathit{\tau }}_{\mathrm{R}}\mathrm{)}}{\mathit{L}} \stackrel{\mathit{\tau }\mathrm{\to }\mathrm{\infty }}{\mathrm{-}\mathrm{\to }} \frac{\mathrm{2}\mathit{a}}{\mathit{L}}\mathrm{=}\frac{\mathrm{2}\mathit{ak}}{\mathit{\pi }}\mathrm{\ll }\mathrm{1.}$(46)Equation (46) demands that ak ≪ 1 such that from Eq. (45) it follows that R_L,R/L_L,R → 0 as τ_L + τ_R → ∞.

When the loops are placed far away from one another, we expect that the loops oscillate independently from one another. Indeed, using the limits (13) on Eq. (42) (in this limiting case, the system always exhibits the standard behaviour) and substituting ω_ki in the resulting expression by (14), we recover after some algebra the formula $\frac{\mathit{\gamma }}{{\mathit{\omega }}_{\mathrm{ki}}}\mathrm{=}\frac{\mathrm{-}\mathit{\pi }}{\mathrm{2}{\mathit{\omega }}_{\mathrm{ki}}^{\mathrm{2}}}\hspace{0.17em}\frac{{\mathit{\rho }}_{\mathrm{i}}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{e}}^{\mathrm{2}}\mathrm{(}{\mathit{\omega }}_{\mathrm{Ai}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{(}{\mathit{\rho }}_{\mathrm{i}}\mathrm{+}{\mathit{\rho }}_{\mathrm{e}}{\mathrm{)}}^{\mathrm{3}}\mathit{\rho }\mathrm{(}{\mathit{\tau }}_{\mathrm{Ai}}\mathrm{\left)}\mathrm{\right|}\mathrm{\Delta }\mathrm{\left(}{\mathit{\tau }}_{\mathrm{Ai}}\mathrm{\right)}\mathrm{|}}$(47)for the damping rate. This result is valid for both loops. We compared this equation to Eq. (77) of Goossens et al. (1992). Since the thicknesses of the inhomogeneous layers can vary with position (see the remarks concerning Fig. 2), we compared Eq. (47) with the damping rates of equivalent tubes described in cylindrical geometry whose inhomogeneous layers are rings whose thickness equals the average thickness of the original layers, described in bicylindrical coordinates. If we assume that the density in the inhomogeneous layer varies almost linearly from its value in the core to the value of the outside plasma, then Eq. (4) implies for the loops described in bicylindrical coordinates that $\begin{array}{ccc}\mathit{\rho }\mathrm{\left(}{\mathit{\tau }}_{\mathrm{Ai}}\mathrm{\right)}\hspace{0.17em}\mathrm{\left|}\mathrm{\Delta }\mathrm{\right|}\mathrm{=}\mathit{\rho }\mathrm{\left(}{\mathit{\tau }}_{\mathrm{Ai}}\mathrm{\right)}{\left|\frac{\mathrm{d}{\mathit{\omega }}_{\mathrm{Ai}}^{\mathrm{2}}}{\mathrm{d}\mathit{\tau }}\right|}_{\mathit{\tau }\mathrm{=}{\mathit{\tau }}_{\mathrm{Ai}}}& \mathrm{=}& {\mathit{\omega }}_{\mathrm{ki}}^{\mathrm{2}}\hspace{0.17em}\frac{{\mathit{\rho }}_{\mathrm{i}}\mathrm{-}{\mathit{\rho }}_{\mathrm{e}}}{{\mathit{l}}_{\mathrm{i}}}\\ & \mathrm{=}& {\mathit{\omega }}_{\mathrm{ki}}^{\mathrm{2}}\hspace{0.17em}\frac{\mathrm{(}{\mathit{\rho }}_{\mathrm{i}}\mathrm{-}{\mathit{\rho }}_{\mathrm{e}}\mathrm{)}{\mathit{R}}_{\mathrm{i}}}{{\mathit{l}}_{\mathrm{avg,i}}}\mathrm{·}\end{array}$(48)This is equal to R_i times the corresponding value of ρ_Ai |Δ | for the equivalent loops in cylindrical coordinates. Substituting Eq. (48) into (47) exactly gives Eq. (77) of Goossens et al. (1992). This perfect agreement between the two coordinate systems is possible only because the geometric term cothτ_i vanishes in Eq. (4) in the limit τ_i → ∞. For loops closer to one another, geometric effects due to the variation of annulus thickness will play a role.

3. Parametric study of the damping rate

In this section we study the dependence of the damping rate γ, given in Eqs. (42) and (43), on the separation distance d/R between the loops and the density contrast between the two loops. We assume that both loops have the same radius, since the numerical calculations of Luna et al. (2008) of eigenoscillations of a two-loop system – without restriction to long wavelengths – showed that the eigenfrequencies only weakly depend on the tube radius. For analytical simplicity we suppose that the density profile is linear throughout the inhomogeneous layers τ ∈  { −τ_i − l, − τ_i } and τ ∈  { τ_i,τ_i + l }: $\overline{\mathit{\rho }}\mathrm{\left(}\mathit{\tau }\mathrm{\right)}\mathrm{=}{\mathit{\rho }}_{\mathrm{e}}\left(\frac{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\zeta }}_{\mathrm{L,R}}\mathrm{)}\mathrm{(}{\mathit{\tau }}_{\mathrm{i}}\mathrm{-}\mathrm{|}\mathit{\tau }\mathrm{\left|}\mathrm{\right)}}{\mathit{l}}\mathrm{+}\mathrm{1}\right)\mathit{.}$(49)Here τ_i = τ_L = τ_R and the subscripts L and R are dropped for clarity of notation where possible. In Eq. (49) the relative density ζ_L was used for τ < 0 and ζ_R was used for τ > 0. With this density profile we have $\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{L}}^{\mathrm{+}\mathit{,}\mathrm{-}}\mathrm{\right(}{\mathit{\tau }}_{\mathrm{A}}\mathrm{\left)}\mathrm{\right|}\mathrm{=}\frac{{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{4}}\mathrm{(}{\mathit{\zeta }}_{\mathrm{L}}\mathrm{-}\mathrm{1}\mathrm{)}}{\mathit{l}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}}\mathrm{·}$(50)We need the thin boundary assumption to be valid throughout the entire parameter space. We know, however, that the tube radius and the shape of the inhomogeneous layer strongly depend on the value of τ_i. The most straightforward solution is to fix the average relative thickness of the nonuniform layer l_avg/R (defined in (3)) to a small number. This also ensures that the expressions for the damping rate become independent of the radius. Finally, all length scales were normalised with respect to the loop length L and the density was normalised with respect to the density ρ_e of the exterior plasma, that is, ζ = ρ/ρ_e. As a consequence, all frequencies are normalised with respect to the Alfvén speed of the exterior medium, ω → ω/ω_Ae.

3.1. Dependence on the separation between loops

Fig. 3 Plot of the lower eigenfrequency ω− (full line) as a function of tube separation d/R, together with the Alfvén frequency of the left tube (dotted line), marking the transition point between standard and anomalous systems, and the kink frequency of the right tube (dashed line).

First we investigated the dependence of the damping rate on the separation distance between the two loops. From (2) we know that d/R = 2coshτ_i. We investigated two cases. In the first case, the loop densities are equal: ζ_L = ζ_R = 3. This means we can use the standard expression (42) for the damping rate everywhere. In the second case, we set ζ_L = 2.5, ζ_R = 2.9. As Fig. 3 shows, only for very small distances (d/R ≲ 2.05) do we need to use the anomalous expression (43) to determine γ₋. It can be read off from Eqs. (42) and (43) that in the limit ω₋ → ω_AL the damping rates of the transverse oscillations in standard and anomalous systems converge to the same value, such that γ₋ is a continuous function of both density and distance at this frequency.

Fig. 4 Dependence of the damping rate γ/ω0 (where ω0 can be either ω+ or ω−) on the separation between the loops. a) Equal loop densities ζL = ζR = 3. b) Unequal loop densities ζL = 2.5, ζR = 2.9. The damping rate corresponding to ω− is shown as full lines, while the damping rate corresponding to ω+ is shown as thick dashed lines. The damping rates for the kink oscillations of the individual loops (thin dashed and dotted lines) are also shown.

Figure 4 shows how the damping properties of the two-loop system vary with increasing distance between the loops. The damping of the kink frequencies of the individual flux tubes are given by $\frac{{\mathit{\gamma }}_{\mathrm{kL,R}}}{{\mathit{\omega }}_{\mathrm{kL,R}}}\mathrm{=}\frac{\mathrm{-}\mathit{\pi }}{\mathrm{8}}\hspace{0.17em}\frac{{\mathit{l}}_{\mathrm{avg}}}{\mathit{R}}\hspace{0.17em}\frac{{\mathit{\zeta }}_{\mathrm{L,R}}\mathrm{-}\mathrm{1}}{{\mathit{\zeta }}_{\mathrm{L,R}}\mathrm{+}\mathrm{1}}$(51)(Goossens et al. 1992). For both eigenfrequencies the sign of γ is negative so that we are dealing with an actual damping mechanism. As is the case in cylindrical coordinates, the results using the TB assumption predict long damping times in contradiction to observed damping times (Aschwanden et al. 2002).

In the previous subsection, we analytically showed that the damping rates of the two eigenfrequencies converge to the damping rates of the individual loop oscillations because they are placed farther apart.

When the separation between the loops becomes smaller, they interact more strongly because the collective homogeneous eigenfrequencies differ more from the individual kink frequencies, as can be seen for the lower eigenfrequency in Fig. 3. Just as in Robertson & Ruderman (2011), we found that the interaction between the loops reduces the efficiency of resonant damping. The damping disappears completely in the limit τ_i → 0. However, this last fact is due to geometric effects corresponding to the bicylindrical coordinate system and is not caused by physical effects. The main effect lies in the factor 1/Δ, which appears in the expressions (42) and (43) of the damping rate. Equations (4) and (50) can be combined to yield $\frac{\mathrm{1}}{\mathrm{\left|}\mathrm{\Delta }\mathrm{\right|}}\mathrm{~}\mathit{l}\mathrm{=}\frac{{\mathit{l}}_{\mathrm{avg}}}{\mathit{R}}\tanh {\mathit{\tau }}_{\mathrm{i}}\mathit{.}$(52)The quantity l is linked to the coordinate system, while we fixed the physical quantity l_avg/R in the calculations. The limit τ_i → 0 hence forces 1/|Δ| → 0, even though the inhomogeneous layer does not disappear here (which is the usual interpretation of this limit). Instead, the inhomogeneous layer becomes very asymmetric because both loops expand when taking τ_i → 0. Hence the limit 1/|Δ| → 0 is due to geometrical effects, and not because of the disappearance of the resonant damping mechanism.

Aside from that, one can calculate that the real part of the eigenfrequencies tends to ${\mathit{\omega }}_{\mathrm{+}}^{\mathrm{2}}\mathrm{\to }{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}\mathit{,} {\mathit{\omega }}_{\mathrm{-}}^{\mathrm{2}}\mathrm{\to }\frac{\mathrm{2}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}}{{\mathit{\zeta }}_{\mathrm{L}}\mathrm{+}{\mathit{\zeta }}_{\mathrm{R}}}$(53)as τ_i → 0. This causes the term ${\mathit{\omega }}_{\mathrm{+}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}$ in Eqs. (42) and (43) to vanish for the damping rate γ₊ corresponding to ω₊. The factor (1 ± F) in Eq. (44) should be identified with ${\mathit{\omega }}_{\mathrm{+}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}$ in Eqs. (42) and (43).

When the two loops have equal densities, Eq. (53) yields that ${\mathit{\omega }}_{\mathrm{-}}^{\mathrm{2}}\mathrm{\to }{\mathit{\omega }}_{\mathrm{Ai}}^{\mathrm{2}}$ with ${\mathit{\omega }}_{\mathrm{Ai}}^{\mathrm{2}}$ the interior Alfvén frequency of both loops. One can see from the second factor in Eq. (42) that because (1 − F²) = 0 a factor $\mathrm{(}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ai}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}$ can be factored out, which vanishes because of (53). These factors correspond with the term (1 − F²) that appears in (44). The limits (53) are themselves influenced by the description in bicylindrical coordinates because of the term F² appearing in (12). However, in anomalous systems these terms do not play a role for the lower frequency as Eq. (43) shows, so the essence of the problem lies in the terms 1/Δ_L,R, as explained before.

In short, the analysis of the expressions for the damping rate shows that the limit 1/|Δ| → 0 as τ_i → 0, caused by the geometry of our description, causes the damping rate to go to zero. It follows that the result that the damping rate decreases with decreasing loop distance is also influenced by the geometrical description in bicylindrical coordinates. Indeed, Fig. 9.4 (a) in Soler (2010), who used the T-matrix method, shows that in cylindrical geometry the damping rate of the high-frequency eigenoscillation increases as the loops move closer together. It is not clear how one can separate the geometry from the physics for small loop distances in a description in bicylindrical coordinates.

Because we assumed a linear or quasi-linear density profile, the average thickness of the inhomogeneous layers l_avg/R only appears as a constant of proportionality in Eq. (42) (with 1/|Δ | substituted from (52)). The damping rate is proportional to the average thickness of the inhomogeneous layers, as is the case for cylindrical structures in the TB limit. (e.g. Goossens et al. 2009). It is reasonable to assume that the results obtained with the TB assumption can be extrapolated to larger thicknesses of the inhomogeneous layers up to l_avg/R ≈ 0.3, as in Van Doorsselaere et al. (2004).

We note from Fig. 4 that the high-frequency oscillations (corresponding to ω₊) are more strongly damped than the low-frequency oscillations (corresponding to ω₋) for equal loop densities, while for unequal loop densities the faster oscillations are more strongly damped only up to a certain separation distance between the loops. In general, the high-frequency oscillations are more strongly damped when the interaction between the loops is strong. A physical explanation for this can be found when we consider the motions of the intermediate fluid. As Luna et al. (2008) explained, the plasma between the loops moves in such a way that it acts as an additional restoring force for the oscillatory motions of the two-loop system for the eigenmodes corresponding to ω₊, while they support the collective motions of the two-loop system for the eigenmodes corresponding to ω₋. (In their formulation, there were four eigenfrequencies since they dropped the long-wavelength assumption, but because these frequencies were split into two pairs lying very close to ω_±, their argument can be used here as well.) The same holds for the damping rate, which is enhanced by the oppositely directed motions at both sides of the resonant layers for the forced eigenmodes corresponding to ω₊. These eigenmodes can be more efficiently damped than the unforced eigenmodes corresponding to ω₋. This effect persists until the system decouples at larger separation distances, when the individual kink properties of the loops dominate.

Fig. 5 Relative difference between the damping rate of transverse oscillations in standard systems with the equivalent expression for anomalous systems for loop densities ζL = 2.5, ζR = 2.9. (For a given separation distance, only one expression is valid.)

Finally, Fig. 5 shows the difference between the standard and anomalous expressions (42) and (43) as a function of separation distance. Of the two expressions, only one can be valid for a specified separation distance given the nature of the system as shown in Fig. 3. We see, however, that the relative error one makes by ignoring the standard/anomalous dichotomy is not larger than 15%. The physical reason for this is the following: in anomalous regimes, the less dense loop cannot support the global oscillations of the system anymore (Van Doorsselaere et al. 2008). Hence, the contribution of the damping in this loop to the total damping rate becomes negligible in anomalous systems.

3.2. Dependence on loop densities

Fig. 6 Dependence of the real part ω of the eigenfrequencies and the corresponding damping rate γ/ω in the two-loop system as a function of the density in the right tube ζR for fixed ζL = 3 and d/R = 3. a) Eigenfrequencies ω+ (thick dashed line) and ω− (full line) as a function of ζR, together with the Alfvén frequencies of the left and right tubes (dashed-dotted lines), marking two transition points between standard and anomalous systems (indicated by saR and saL on the horizontal axis, referring to the loop that cannot support the global oscillation), and the individual kink frequencies of the two tubes (dotted and dashed lines). b) Ratio of damping rate to real part of the eigenfrequency of the two-loop system γ/ω0, where ω0 can be either ω+ (thick dashed line) or ω− (full line), and of the left (dotted line) and right (dashed line) loop as a function of ζR.

In this subsection we investigate the influence of the different tube densities on the damping rate. We set d/R = 3, which is a realistic distance between the loops according to the observations by Aschwanden et al. (2003). We fixed the density of the left loop ζ_L = 3 and let ζ_R vary between 1 and 6. This means that either loop can be the less dense one, and that an anomalous system exists at both ends of the density range ζ_R. This is illustrated in Fig. 6a; the two intersections where such a transition happens are marked on the horizontal axis. A similar dependency of the eigenmodes on the loop densities has been obtained by Luna et al. (2009) using the T-matrix method (see their Fig. 1). For low values of ζ_R, the right loop cannot support the global oscillations while for high values of ζ_R, the left loop will be unable to support them. It can be shown that here as well, the contribution of the anomalous regime is negligible because only one loop is oscillating. When the loop densities are equal, the interaction between the loops is maximal since the frequencies ω₊ and ω₋ differ most from the individual kink frequencies.

Figure 6b illustrates the dependence of the observational quantity γ/ω₀ (where ω₀ can be ω₊ or ω₋) on the density contrast between the two loops when the distance between them is kept fixed. When the loop densities are equal, the higher eigenfrequency is damped more strongly than for unequal loop densities. The lower eigenfrequency is least damped in this case. Moreover, the damping of the high-frequency oscillation is higher than that of the low-frequency oscillation, |γ₊/ω₊| > |γ₋/ω₋|, when the loops have a similar density. The situation is reversed when the relative density contrast between the loops is large.

The reason for this can be found when we analyse Eq. (42) as a function of the two loop densities ζ_L and ζ_R. We isolated the contribution of each loop to the total damping rate. The only terms that are not invariant under the permutation ζ_L ↔ ζ_R are the two terms in the second factor of (42). The contribution of these terms is plotted in Fig. 7. We see that the contribution of the term corresponding to the denser loop dominates the expression for γ₋, while the contribution of the less dense loop dominates in the expression for γ₊. Because the inequality $\frac{{\mathit{\zeta }}_{\mathrm{L}}\mathrm{(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{AL}}^{\mathrm{2}}\mathrm{)}}{\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{L}}^{\mathrm{+}\mathit{,}\mathrm{-}}\mathrm{\right|}}\mathit{<}\frac{{\mathit{\zeta }}_{\mathrm{R}}\mathrm{(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{AR}}^{\mathrm{2}}\mathrm{)}}{\mathrm{\left|}{\mathrm{\Delta }}_{\mathrm{R}}^{\mathrm{+}\mathit{,}\mathrm{-}}\mathrm{\right|}}$ holds when ζ_L < ζ_R for linear density profiles, this means that the factor $\mathrm{\left[}\mathrm{\right(}\mathrm{1}\mathrm{+}{\mathit{F}}^{\mathrm{2}}\mathrm{\left)}{\mathit{\rho }}_{\mathrm{L,R}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{AL,R}}^{\mathrm{2}}\mathrm{)}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{F}}^{\mathrm{2}}\mathrm{\left)}{\mathit{\rho }}_{\mathrm{e}}\mathrm{\right(}{\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right]}$ dominates in the terms describing the contribution from each loop to the total damping rate.

Fig. 7 Contribution of the two terms in the last factor of Eq. (42) to the total damping rate of the higher (+) and lower (–) eigenfrequency. The symbols L and R refer to the terms containing 1/|ΔL,R |, respectively.

The sharper the density contrast, the more substantial the contribution of the denser loop becomes for γ₋, and conversely, the more substantial the contribution of the less dense loop becomes for γ₊. In these regimes the damping rate γ₋ can be associated with the damping of the kink oscillation of the denser loop. From this and Eq. (51) for the damping of the individual kink modes by resonant absorption it follows that the low-frequency eigenoscillation is more strongly damped when ζ_L ≠ ζ_R.

Soler (2010) noted that even if the real part ω₋ of the collective oscillation eigenfrequency falls outside of the Alfvén spectrum for one of the loops, the collective eigenmodes can still be efficiently damped due to the resonance in the other loop. This can now be explained for the zero-beta model since the global damping rate γ₋ depends more on the damping properties of the denser loop. It is reasonable to assume that similar results continue to hold for a finite plasma beta. An interesting question is whether the denser loops also contribute most to the total damping rate of kink-like collective oscillations for systems of more than two coronal loops. This we cannot model using bicylindrical coordinates. We conjecture that the answer in general is affirmative, since Luna et al. (2010) have shown that in a collective oscillation of multiple nonidentical strands, the largest amplitudes occur in the most dense strands for the lowest eigenfrequencies corresponding to kink-like modes and vice versa.

The quantitative results of the previous paragraphs depend on the specific form of the density profile through the gradients of the Alfvén speeds at the resonant surface in the two tubes. About these we do not have much physical information. In particular, if the density profile were different in the two loops, this would affect the relative contribution of both resonances to the global damping of the oscillation. Hence our results are only as accurate as the assumption of quasi-linearity.

When the loop densities are similar, ζ_L ≈ ζ_R, the two terms of the second factor in (42) become about equally important because this expression reduces to (44), found by Robertson & Ruderman (2011), for identical loops. From this equation it is evident that | γ₊/ω₊| > |γ₋/ω₋ |. This is due to two factors: on the one hand (1 ± F), whose spatial dependence we identified in subsection 3.1. with that of the factor ${\mathit{\omega }}_{\mathrm{+}}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{Ae}}^{\mathrm{2}}$ in Eq. (42). The factor ζ + 1 ∓ (ζ − 1)F that appears in the denominator of (44) comes from the reduction of ${\mathit{\omega }}_{\mathrm{\pm }}^{\mathrm{2}}$ when ζ_L = ζ_R. Hence, when the loops interact more strongly with one another, the difference between γ₊ and γ₋ is explained mathematically by the collective oscillation properties of the loops through the factor ω_±. Therefore, we believe that the explanation in terms of the flow of the intermediate plasma of Luna et al. (2008), cited when considering the effect of distance between the loops in the previous subsection, can also physically explain the predicted damping times for the high- and low-frequency eigenoscillations.

Fig. 8 Complex eigenspectrum of the two-loop system and of the left and right loops as a function of ζR. a) Complex spectrum for the higher (crosses) and lower (diagonal crosses) eigenfrequencies of the two-loop system, and the eigenfrequencies of the equivalent left (solid diamond) and right (boxes) loop in cylindrical coordinates as a function of ζR. The damping rates of the left and right loop have been corrected for geometric effects (see text). The parameter ζR increases from 1 to 9 in steps of 0.2. The direction in which the eigenfrequencies move when ζR increases are indicated with arrows. b) Damping rates of the two-loop system γ+ (thick dashed line) and γ− (full line) and of the left (dotted line) and the right (dashed line) loop as a function of ζR. The damping rates of the left and right loop have been corrected for geometric effects (see text). c) The density ζR of the right loop as a function of the eigenfrequencies ω+ (thick dashed line) and ω− (full line) of the two-loop system, and ωkL (dotted line) and ωkR (dashed line) of the left and right loop, respectively.

When the densities of the two loops approach each other, the eigenoscillation corresponding to ω₊ gradually changes its character from being related to the less strongly damped oscillations of the low-density loop to the more heavily damped forced collective oscillations of the two-loop system. The global emerging phenomenon which arises from this is an avoided crossing of the kink modes of the individual loops, as was previously noted by Luna et al. (2009). The characteristics of this phenomenon are now analysed from the same point of view as in the analysis performed by Van Doorsselaere & Poedts (2007) for the complex spectrum of resistive MHD in cylindrical coordinates.

In the complex spectrum shown in Fig. 8a, the complex eigenfrequencies of the two-loop system and of the right loop travel to the left in the complex plane (indicating a decrease in frequency) as ζ_R increases. The eigenfrequency of the viscous eigenmode of the left loop remains constant since ζ_L is fixed. The interaction between the loops is the strongest for ζ_L ≈ ζ_R with the turning points of the two-loop eigenfrequencies at ζ_R = 3. A qualitatively similar picture emerges when we fix another value for ζ_L; the only difference is that the point (ω_kL,γ_kL) moves to the left along the curve defined by (ω_kR,γ_kR) as ζ_L increases, because expressions (12) and (47) with (48) substituted are the same for both loops.

The effects of avoided crossing on the imaginary part of the eigenfrequencies are illustrated in Fig. 8b. The standard and anomalous expressions (42) and (43) for γ were used where appropriate. In Figs. 8a and 8b the geometrical effects associated with the limits (ζ_R − 1)/(ζ_L − 1) → 0 and (ζ_R − 1)/(ζ_L − 1) → ∞ calculated when discussing Fig. (6b) were eliminated artificially by multiplying the damping rates of the kink oscillations by tanhτ_i. The transfer of damping properties of the individual loops from one eigenmode to the other is visible in Fig. 8b, but the convergence of γ_+ , − to γ_kL,R when ζ_R → ∞ is slower than the convergence of ω_+ , − to ω_kL,R.

From Fig. 8c we infer that as (ζ_R − 1)/(ζ_L − 1) increases, the eigenmode corresponding to ω₊ and γ₊ transfers some of its properties of the left loop to the eigenmode corresponding to ω₋ and γ₋ and inherits the properties of the right loop from it, and vice versa. Near the region where ζ_L = ζ_R the collective properties of the two-loop system are dominant.

When the density of the right tube decreases to the density of the exterior fluid, the eigenfrequency ω₊ tends to the Alfvén frequency of the exterior fluid (Fig. 6a), since it is the highest frequency present in the system. In this limiting case there remains only one loop, hence γ₊ converges to zero in Fig. 6b. The damping rate γ₋ does not converge to γ_kL but rather to γ_kLtanhτ_i = γ_kLtanh(arccosh(d/2R)) because of Eqs. (2) and (52). Similarly, we have γ₊ → γ_kLtanhτ_i and γ₋ → γ_kRtanhτ_i as (ζ_R − 1)/(ζ_L − 1) → ∞.

Finally, we compared our results with those of Fig. 9.10 (b) of Soler (2010). The behaviour of damping times as a function of relative loop density for the same equilibrium parameters as considered by Soler (2010) is shown in Fig. 9. The behaviour of the damping times near equal loop densities is the same as we found, and has been explained earlier in this subsection. In Fig. 9.10 (b) of Soler (2010) we have τ₋/P₋ < τ₊/P₊ for ζ_R > ζ_L in the range of values of ζ_R/ζ_L shown. This is probably because this range is too narrow, since the limits lim_{(ζR−1)/(ζL−1)→∞} τ₊/P₊ = τ_kL/P_kL > τ_kR/P_kR = lim_{(ζR−1)/(ζL−1)→∞} τ₋/P₋ hold in cylindrical geometry, as is also visible from Fig. 9. Soler (2010) noted that the damping times of the low-frequency eigenoscillations are more weakly dependent on the relative loop densities than the damping times of high-frequency oscillations. This is not the case in our analytic description. The difference can be probably explained by considering that the geometric factors describing the dependence of the damping rate on the separation distance for similar loop densities are slightly different in the description with the T-matrix method than in the description in bicylindrical coordinates.

The results in this section can be compared with those of an equivalent cylindrical loop with a symmetric inhomogeneous layer of width l_avg/R (in the sense that geometrical effects do not hide physical effects such as interaction between the loops) when d/R ≳ 4.58 because then tanh(arccoshd/2R) > 0.90. This is a rather large separation distance compared with the observations of Aschwanden et al. (2003). As Fig. 4 of Van Doorsselaere et al. (2008) shows, for such separation distances the loops are already mostly decoupled from one another.

4. Conclusions

We have studied the damped transverse oscillations of a system of two parallel not necessarily identical coronal loops. We restricted ourselves to long wavelengths and thin inhomogeneous layers. For linear or quasi-linear density profiles, we derived analytic expressions for the damping rate γ of the two-loop system that generalise the results found for single cylinders and systems of two identical loops, as we have shown. A subdivision between standard and anomalous systems must be made here.

Fig. 9 Ratio of damping time to oscillation period of the two-loop system τ/P0, where P0 can be either P+ (thick dashed line) or P− (full line), and of the left (dotted line) and right (dashed line) loop as a function of ζR. The equilibrium parameters have been chosen identical to those of Soler (2010).

A parametric study to investigate the dependence of γ on the dimensionless parameters led to two important conclusions. Firstly, we saw that the interaction between the loops increases the damping rate of the high-frequency oscillation, while the damping rate of the low-frequency oscillation is reduced by it. The analysis of Eq. (42) showed that for a strong relative density contrast the denser loop contributes most to the damping rate γ₋ and the less dense loop contributes most to the damping rate γ₊, resulting in the inequality |γ₊/ω₊| < |γ₋/ω₋| due to the dominating kink-like properties of the loops. On the other hand, for similar loop densities the terms related to ω_± dominate the expression for the damping rate. This means that the explanation of Luna et al. (2008) for the different eigenfrequencies of the two-loop system, in terms of the motions of the intermediate fluid, is very likely valid and applicable to the damping rate as well. The intermediate plasma motions act as an additional restoring force for the high-frequency eigenoscillation, which results in shearing motions and an increased damping rate for this collective oscillation, while they support the motions of the low-frequency eigenoscillation, which leads to a decreased damping rate for this collective oscillation. Mathematically, the behaviour of the eigenmodes of the two-loop system as a function of ζ_R was explained by an avoided crossing of the kink eigenmodes of the individual loops in the complex plane. Our analysis complements the findings of Luna et al. (2008, 2009), Soler (2010) and quantitatively justifies the physical explanation in terms of the motions of the intermediate fluid provided by these authors.

Secondly, we found that in our description, the interaction between the loops reduces the efficiency of resonant damping. However, as conjectured in the final paragraph of Robertson & Ruderman (2011) and confirmed by Eq. (52), this effect is attributed at least partly to the change of shape of the inhomogeneous layer linked to the bicylindrical coordinate system and not to physical effects. Equation (53) shows which terms in the equations for ω_± and γ_± force the damping rate to vanish in the limit of touching loops.

Even though we also used the bicylindrical coordinate system in this paper, we made clear that there are several problems with this formalism. Summarising, the problems are that

• 1.

  it is not possible to preventR → 0 if d/R → ∞.

• 2.

  The change of shape of inhomogeneous layer cannot be dealt with in a straightforward way within the bicylindrical coordinate system to investigate the limit of small separation distances between the two loops.

• 3.

  The results cannot be generalised to systems of more than two coronal loops.

Therefore we believe additional analytical studies should make use of other mathematical techniques. For a system of two parallel coronal loops, one could look for analytic expressions of the eigenfrequencies and eigenfunctions using the T-matrix method in a local coordinate system, both for undamped and damped oscillations of coronal loops. This would facilitate the comparison with known results from cylindrical geometry. Alternatively, it would be interesting to extend the analysis of Soler (2010) to systems of more than two coronal loops. This approach could lead to a unified understanding of the damped oscillations of systems of coronal loops.

Acknowledgments

The authors thank the referee and Michael Ruderman for valuable comments, which helped improve the content of this paper considerably. S.E.G. acknowledges Arno Kuijlaars and Marcel Goossens for useful discussions. T.V.D. and S.E.G. have received funding from an Odysseus grant of the FWO Vlaanderen and FP7 under grant agreement 276808. We acknowledge the support from Belspo under IAP P7/08 CHARM.

References

All Figures

	Fig. 1 Sketch of the equilibrium configuration (from Robertson & Ruderman 2011).
In the text

	Fig. 2 Sketch of a plane z = const. in the bicylindrical (σ,τ) coordinate system. The lines σ = const. are circles passing through the points x =  ± a, the lines τ = const. are nested circles around these points, whose radius tends to zero and whose centre tends to ±a as |τ| → ∞ (from Van Doorsselaere et al. 2008).
In the text

	Fig. 3 Plot of the lower eigenfrequency ω− (full line) as a function of tube separation d/R, together with the Alfvén frequency of the left tube (dotted line), marking the transition point between standard and anomalous systems, and the kink frequency of the right tube (dashed line).
In the text

Fig. 4 Dependence of the damping rate γ/ω0 (where ω0 can be either ω+ or ω−) on the separation between the loops. a) Equal loop densities ζL = ζR = 3. b) Unequal loop densities ζL = 2.5, ζR = 2.9. The damping rate corresponding to ω− is shown as full lines, while the damping rate corresponding to ω+ is shown as thick dashed lines. The damping rates for the kink oscillations of the individual loops (thin dashed and dotted lines) are also shown.

In the text

	Fig. 5 Relative difference between the damping rate of transverse oscillations in standard systems with the equivalent expression for anomalous systems for loop densities ζL = 2.5, ζR = 2.9. (For a given separation distance, only one expression is valid.)
In the text

Fig. 6 Dependence of the real part ω of the eigenfrequencies and the corresponding damping rate γ/ω in the two-loop system as a function of the density in the right tube ζR for fixed ζL = 3 and d/R = 3. a) Eigenfrequencies ω+ (thick dashed line) and ω− (full line) as a function of ζR, together with the Alfvén frequencies of the left and right tubes (dashed-dotted lines), marking two transition points between standard and anomalous systems (indicated by saR and saL on the horizontal axis, referring to the loop that cannot support the global oscillation), and the individual kink frequencies of the two tubes (dotted and dashed lines). b) Ratio of damping rate to real part of the eigenfrequency of the two-loop system γ/ω0, where ω0 can be either ω+ (thick dashed line) or ω− (full line), and of the left (dotted line) and right (dashed line) loop as a function of ζR.

In the text

	Fig. 7 Contribution of the two terms in the last factor of Eq. (42) to the total damping rate of the higher (+) and lower (–) eigenfrequency. The symbols L and R refer to the terms containing 1/|ΔL,R |, respectively.
In the text

Fig. 8 Complex eigenspectrum of the two-loop system and of the left and right loops as a function of ζR. a) Complex spectrum for the higher (crosses) and lower (diagonal crosses) eigenfrequencies of the two-loop system, and the eigenfrequencies of the equivalent left (solid diamond) and right (boxes) loop in cylindrical coordinates as a function of ζR. The damping rates of the left and right loop have been corrected for geometric effects (see text). The parameter ζR increases from 1 to 9 in steps of 0.2. The direction in which the eigenfrequencies move when ζR increases are indicated with arrows. b) Damping rates of the two-loop system γ+ (thick dashed line) and γ− (full line) and of the left (dotted line) and the right (dashed line) loop as a function of ζR. The damping rates of the left and right loop have been corrected for geometric effects (see text). c) The density ζR of the right loop as a function of the eigenfrequencies ω+ (thick dashed line) and ω− (full line) of the two-loop system, and ωkL (dotted line) and ωkR (dashed line) of the left and right loop, respectively.

In the text

	Fig. 9 Ratio of damping time to oscillation period of the two-loop system τ/P0, where P0 can be either P+ (thick dashed line) or P− (full line), and of the left (dotted line) and right (dashed line) loop as a function of ζR. The equilibrium parameters have been chosen identical to those of Soler (2010).
In the text