© The Authors 2024

1. Introduction

In recent years, increasing numbers of observations have revealed that magnetic structures in the solar corona are dynamic and support a rich variety of low-frequency waves in the magnetohydrodynamic (MHD) regime (see, e.g., Li et al. 2020; Wang et al. 2021; Nakariakov et al. 2021; Banerjee et al. 2021, for recent reviews). As energy carriers, these waves may be of potential importance for tackling the coronal heating problem (see the reviews by, e.g., De Moortel & Browning 2015; Arregui 2015; Van Doorsselaere et al. 2020). Likewise, these waves themselves carry information about their hosts, and thus have been considered useful for probing the physical conditions in the solar corona (e.g., Nakariakov & Verwichte 2005; De Moortel & Nakariakov 2012; Nakariakov & Kolotkov 2020). Furthermore, coronal waves are believed to be related to the energy release processes in solar flares as well (e.g., McLaughlin et al. 2018; Zimovets et al. 2021; Inglis et al. 2023).

Canonically accepted as the only motions that displace structural axes, kink modes have frequently been imaged in the solar corona (see, e.g., Nakariakov et al. 2021 for a recent review; see also Nakariakov et al. 1999; Aschwanden et al. 1999 for early dectections). They have been routinely employed in coronal seismology, which combine wave theories and observations to deduce those coronal parameters that are difficult to directly measure. For instance, the measured periods of kink motions have been in routine use for inferring the coronal magnetic field strength (e.g., Nakariakov & Ofman 2001; Yang et al. 2020). Likewise, transverse inhomogeneity lengthscales can be seismologically deduced with the measured damping times of decaying kink motions, provided that this damping is attributable to such mechanisms as resonant absorption (e.g., Goossens et al. 2002; Aschwanden et al. 2003; Arregui & Asensio Ramos 2014; Arregui et al. 2019). Regarding this aspect, we note that resonant absorption is an ideal process that transfers kink energy into localized Alfvénic motions (see the review by Goossens et al. 2011). This process has been shown to be robust for a rich set of configurations (e.g., Van Doorsselaere et al. 2004a; Terradas et al. 2008; Pascoe et al. 2011; Howson et al. 2017; Guo et al. 2020; Shi et al. 2024a).

Solar coronal loops have been customarily modeled as straight magnetic cylinders in theoretical studies on kink motions (e.g., Edwin & Roberts 1983; see also Zajtsev & Stepanov 1975; Wentzel 1979; Spruit 1982; Cally 1986). For instance, resonantly damped kink motions have been extensively studied from the ideal quasi-mode perspective with the Frobenius approach (e.g., Soler et al. 2013; Soler 2019; Geeraerts et al. 2022) or more frequently dissipative eigenmode computations (e.g., Poedts & Kerner 1991; Van Doorsselaere et al. 2004b; Terradas et al. 2006b; Guo et al. 2016). In particular, concise expressions are available for the oscillation frequencies and damping rates of kink quasi-modes when the thin-tube-thin-boundary (TTTB) approximation applies (e.g., Goossens et al. 1992, 2009; Soler et al. 2013). From the initial value problem perspective, it has been shown that the theoretical expectations for kink quasi-modes tend to well describe the temporal evolution of kink motions (e.g., Ruderman & Roberts 2002), particularly when loop boundaries are thin (Soler & Terradas 2015).

Closed loops, believed to be typical magnetic structures in the corona (Reale 2014), serve as the main waveguides for kink motions. A question then arises of how the loop curvature influences the eigenfunctions of kink modes. Previous wave-related investigations focused on two-dimensional (2D) loop models (e.g., Smith et al. 1997; Brady & Arber 2005; Selwa et al. 2005, 2006; Verwichte et al. 2006a,b,c). In particular, Verwichte et al. (2006a,b,c) analyzed fast magnetoacoustic modes and their seismological potential in detail in a curved slab model under the plasma β = 0 limit. Lateral leakage of kink motions via wave tunneling was examined. However, these two-dimensional studies are constrained to vertically polarized kink modes only. While vertical polarizations have been reported (e.g., Wang & Solanki 2004), a more comprehensive examination, including the more frequently observed horizontal polarization (e.g., Zhong et al. 2023), should also be taken into account. This can only be achieved in three-dimensional (3D) examinations. Initial analytical progress was made by Van Doorsselaere et al. (2004a, 2009) in a toroidal model embedded in a radially dependent force-free field, showing that the curvature of a loop has little influence on the kink eigenmodes in the absence of wave leakage. Beyond the consideration of a constant loop cross section, Ruderman (2009) found that the difference in eigenfrequencies of vertical and horizontal oscillations is proportional to the tube expansion parameter under the thin tube (TB) approximation. Numerically, the semi-torus model with a constant loop cross section considered by Terradas et al. (2006a) distinguished horizontally and vertically polarized oscillations. The frequency difference between these two polarizations is not observable. Meanwhile, lateral leakage that is naturally induced in curved loops has also been included, and it is found to be less efficient than resonant absorption. The response of such a semi-circular loop model to external perturbations was also numerically examined by Pascoe & De Moortel (2014). In addition, kink oscillations in a force-free magnetic field with a dipole configuration have been considered by, e.g., McLaughlin & Ofman (2008), Magyar & Nakariakov (2020). In particular, strong damping of kink oscillations is illustrated by McLaughlin & Ofman (2008), and it was attributed to the curvature of their loop model. It seems that the efficiency of lateral leakage and resonant damping in curved loops needs to be clarified further. Regarding the frequency of kink modes in curved loops, Magyar & Nakariakov (2020) found that the WKB approximation can reasonably describe the kink period in their planar loop model without considering a sigmoid geometry. However, the influence of a potential field configuration on the different polarizations of kink modes still remains unclear.

In the current work we consider a curved loop model embedded in a potential magnetic field and examine its response to initial velocity perturbations from different directions. The influence of loop curvature on the eigenfrequency and damping rate of the excited kink polarizations is examined. Our study differs from the available studies in the following aspects. First, both horizontal and vertical polarizations are considered and their frequencies are compared with the WKB approximation. Second, lateral leakage is investigated in 3D simulations, and the efficiency of this damping mechanism is examined and compared with resonant absorption. This paper is organized as follows. Section 2 describes the equilibrium configuration and the numerical setup. The simulation results are presented in Sect. 3, followed by a summary and discussion of the present study in Sect. 4.

2. Numerical model

We considered a potential magnetic field to mimic the magnetic structures in the solar corona. A similar configuration was considered by Oliver et al. (1993), Selwa et al. (2005), and Rial et al. (2013), among others. We let (x, y, z) be a Cartesian coordinate system. The equilibrium magnetic field is given by

$$\begin{array}{cc}& B_x(x,z)=B_{0}cos\left(\frac{x}{{\mathrm{\Lambda }}_B}\right)exp(-\frac{z}{{\mathrm{\Lambda }}_B}),\hfill \end{array}$$(1)

$$\begin{array}{cc}& B_z(x,z)=-B_{0}sin\left(\frac{x}{{\mathrm{\Lambda }}_B}\right)exp(-\frac{z}{{\mathrm{\Lambda }}_B}),\hfill \end{array}$$(2)

where B₀ = 20 G, Λ_B = 2L_B/π with L_B = 45 Mm. This magnetic field can be alternatively expressed via the vector potential

$$\begin{array}{c}\hfill \mathit{A}=\psi \widehat{y}=B_0{\mathrm{\Lambda }}_{B}cos\left(\frac{x}{{\mathrm{\Lambda }}_B}\right)exp(-\frac{z}{{\mathrm{\Lambda }}_B})\widehat{y},\end{array}$$(3)

where ψ is the magnetic flux function whose contours delineate magnetic lines of force. A set of right-handed orthonormal basis vectors $({\widehat{\mathit{e}}}_{\mathrm{t}},{\widehat{\mathit{e}}}_{\mathrm{h}},{\widehat{\mathit{e}}}_{\mathrm{v}})$ can be defined, where

$$\begin{array}{cc}& {\widehat{\mathit{e}}}_{\mathrm{t}}=cos(x/{\mathrm{\Lambda }}_B)\widehat{x}-sin(x/{\mathrm{\Lambda }}_B)\widehat{z},\hfill \\ \hfill & {\widehat{\mathit{e}}}_{\mathrm{h}}=\widehat{y},\hfill \\ \hfill & {\widehat{\mathit{e}}}_{\mathrm{v}}=sin(x/{\mathrm{\Lambda }}_B)\widehat{x}+cos(x/{\mathrm{\Lambda }}_B)\widehat{z}.\hfill \end{array}$$(4)

Clearly, ${\widehat{\mathit{e}}}_{\mathrm{t}}$ (${\widehat{\mathit{e}}}_{\mathrm{v}}$) is locally parallel (perpendicular) to the magnetic field. Thus, we further identify ${\widehat{\mathit{e}}}_{\mathrm{h}}$ and ${\widehat{\mathit{e}}}_{\mathrm{v}}$, which characterize the horizontal and vertical polarizations, respectively.

Density-enhanced loops are constructed as follows. Consider, for now, a uniform-density loop, by which we mean the density is uniform along the magnetic field. The density distribution is given by

$$\begin{array}{c}\hfill \rho (\overline{r})={\rho }_{\mathrm{e}}+({\rho }_{\mathrm{i}}-{\rho }_{\mathrm{e}})f(\overline{r}),\end{array}$$(5)

with

$$\begin{array}{c}\hfill f(\overline{r})=\frac{1}{2}\{1-tanh\left[(\frac{\overline{r}}{R}-1)b\right]\},\end{array}$$(6)

where

$$\begin{array}{c}\hfill \overline{r}=\sqrt{{\left(\frac{\psi -{\psi }_0}{\mathrm{\Delta }\psi }\right)}^2+{\left(\frac{y}{\mathrm{\Delta }y}\right)}^2},\end{array}$$(7)

with ψ₀ = ψ(x = 30 Mm, z = 0), Δψ = B₀Λ_B, Δy = 2Λ_B. Here ρ_i (ρ_e) represents the internal (external) density, and R = cos(30/Λ_B)−cos(30.5/Λ_B) prescribes the x-extent of the loop at its footpoints. The parameter b determines the width of the boundary layer. The 3D view of the loop can be seen in Fig. 1a.

Fig. 1. Three-dimensional rendering of the density-enhanced loop. Several arbitrarily selected magnetic field lines are also shown and color-coded by the local magnetic field strength. The left and right panels correspond to our “uniform-density” and “stratified” computations (see text for more details).

We proceed to specify the rest of the physical parameters. The internal (external) density is ρ_i = 1.17 × 10⁻¹⁴ g cm⁻³ (ρ_e = 2.34 × 10⁻¹⁵ g cm⁻³). The thickness of the boundary layer is determined by b = 10, which gives the layer width l ≈ 0.3a, where a is the half-width of the cross section at the loop apex. The temperature of the background corona is T_e = 1 MK. The pressure is uniform in the entire computational domain, which is set to be p = 2.3k_Bρ_eT_e, with k_B the Boltzmann constant.

To distinguish between different polarizations of kink modes, we employ two initial velocity perturbations. To excite the horizontal polarization, we introduce an initial velocity perturbation of the form

$$\begin{array}{c}\hfill \mathit{v}(x,y,z;t=0)=v_{0}sin\theta f(\overline{r}){\widehat{\mathit{e}}}_{\mathrm{h}},\end{array}$$(8)

where v₀ represents the amplitude of the velocity perturbation, $sin\theta =\frac{z}{\sqrt{x^2+z^2}}$ (with θ the angle along the loop from the footprint), and $f(\overline{r})$ is given by Eq. (6). It ensures maximum velocity perturbation at loop apex and zero velocity at loop footpoints in the y-direction. Similarly, we introduce the other velocity perturbation to excite the vertical polarization, which reads

$$\begin{array}{c}\hfill \mathit{v}(x,y,z;t=0)=v_{0}sin\theta f(\overline{r}){\widehat{\mathit{e}}}_{\mathrm{v}}.\end{array}$$(9)

To limit the simulations to the linear regime, we set the velocity amplitude to v₀ = 5 km s⁻¹. This choice ensures that the maximum displacement of the loop is much less than the loop radius. The velocity fields of the initial perturbations are shown in Figs. 2a and 3a. Although these initial perturbations cannot exactly match the form of the eigenmodes, which are not straightforward to obtain beforehand in the current model, kink eigenmodes are subsequently excited in our loop models as a response to these initial velocity perturbations.

Fig. 2. Velocity fields (arrows) for the horizontal polarization at (a) t = 0, (b) t = 34 s, and (c) t = 661 s. The left and right columns correspond to the y = 0 and x = 0 cuts, respectively. Overplotted are the filled contours of vh = vy. The green curves in the left column represent the magnetic lines of force; the dashed lines outline the density-enhanced loop. These snapshots are from the online animations, which run from t = 0 to 1270 s.

Fig. 3. Similar to Fig. 2, except that the filled contours are for $v_{\mathrm{v}}:=\mathit{v}·{\widehat{\mathit{e}}}_{\mathrm{v}}$. An animation is available online and runs from t = 0 to 1270 s.

The boundary conditions are specified as follows. Reflective boundary conditions are applied for the three components of velocity at z = 0 to mimic the footpoints of the magnetic structures anchored in the lower solar atmosphere. All other variables at z = 0 are set to have zero gradients. Asymmetric boundary conditions are used for v_x, B_y, and B_z to follow the oscillating properties of kink modes at x = 0. All other variables there are set to be continuous. Zero gradient boundary conditions are employed at the other side boundaries and the top boundary.

We evolved the set of 3D ideal MHD equations with the finite-volume code MPI-AMRVAC (Xia et al. 2018; Keppens et al. 2021). Gravity was neglected throughout. As described in Xia et al. (2017), we split the magnetic field into B = B₀ + B₁, taking B₀ as the equilibrium one given by Eq. (2) and actually solving for the time-dependent component (B₁). An approximate HLL Riemann solver and a third-order “cada3” slope limiter (Čada & Torrilhon 2009) were employed for evaluating inter-cell fluxes. We used the third-order TVD Runge-Kutta method for time marching with a Courant number of 0.5. The computational domain was chosen to be [0, 40]×[−10, 10]×[0, 40] Mm, which incorporates only the x ≥ 0 portion given the symmetric properties of our numerical implementation. Adaptive mesh refinement (AMR) was employed to further reduce the computational cost. We adopted a base grid of 40 × 60 × 120. Four levels of AMR were implemented, leading to the highest resolution of 41.7 km across the loop region.

3. Results

3.1. General snapshots

We first examine the response of the uniform-density loop to the initial perturbations by analyzing the distribution of velocity in the loop. Figures 2 and 3 show the velocity field at y = 0 and x = 0 (the loop apex) for horizontal and vertical polarizations at different times. The density-enhanced loop is outlined by dashed lines, and the magnetic lines of force are shown as green lines. We start by noting that wave reflection can be seen in Fig. 3b due to the asymmetric boundary at the bottom surface. However, the reflected fast waves quickly leave the computational domain and do not influence the following analysis. In both Figs. 2b and 3b we can clearly see the well-known dipole-like velocity fields that are characteristic of kink motions (e.g., Goossens et al. 2014; Guo et al. 2020) around the loop boundary. From the related animations, we can observe damping kink oscillations in the loop region for both polarizations. A signature of resonant absorption is identified as an increase in localized Alfvénic motions around the loop boundary. For the horizontal polarization, Alfvénic motions can be found around z = 20.5 Mm and z = 19.5 Mm at the loop apex, namely the upper and lower boundaries of the loop cross section. Likewise, such Alfvénic motions can be observed around y = −0.5 Mm and y = 0.5 Mm at the loop apex for the vertically polarized mode. Resonant absorption is thus one of the damping mechanisms in the current model. This mechanism was investigated in many previous works (e.g., Guo et al. 2020). Here we demonstrate that resonant absorption is a robust damping mechanism that manifests itself in curved loops as well. In addition to resonant absorption, lateral leakage of the wave modes can also be observed in the animations available online, and in Figs. 2c and 3c as selected snapshots. Oscillations due to the wave leakage can be seen from the upper region outside the loop for both polarizations. A detailed discussion of damping mechanisms is presented below.

3.2. Oscillating frequency

We examined the oscillation frequency in both polarizations. Figure 4 shows the evolution of velocity sampled at the loop apex for both polarizations. We located the local extrema in both oscillation profiles, thereby evaluating the oscillation period as twice the average temporal spacing between two adjacent extrema. Consequently, we obtained two different periods for the horizontal polarization (P_h = 326.3 s) and the vertical case (P_v = 401.8 s).

Fig. 4. Oscillation profiles of the horizontally (black) and vertically (blue) polarized modes. The velocity is sampled at x = 0, y = 0, z = 19.85 Mm. The dots represent the extrema of each curve.

We can compare these derived periods with analytical predictions. For straight cylinders that are structured only in the radial direction, a simple expression is available for the oscillation frequency of kink modes in the TTTB limit (e.g., Goossens et al. 1992, 2002), writing

$$\begin{array}{c}\hfill {\omega }_{\mathrm{k}}=\sqrt{\frac{{\rho }_{\mathrm{i}}{\omega }_{\mathrm{Ai}}^2+{\rho }_{\mathrm{e}}{\omega }_{\mathrm{Ae}}^2}{{\rho }_{\mathrm{i}}+{\rho }_{\mathrm{e}}}},\end{array}$$(10)

where ω_Ai (ω_Ae) represents the internal (external) Alfvén frequency. Equation (10) becomes even simpler when thermal pressure is negligible, reading

$$\begin{array}{c}\hfill {\omega }_{\mathrm{k}}={\omega }_{\mathrm{Ai}}\sqrt{\frac{2{\rho }_{\mathrm{ie}}}{{\rho }_{\mathrm{ie}}+1}},\end{array}$$(11)

where ρ_ie is the ratio of the internal to the external density. Equation (11) serves as a reasonable starting point for our further analysis, given the small value of the plasma β in our equilibrium setup (β ≲ 0.07 in the loop region). The internal Alfvén frequency changes with height, and we thus consider a WKB approximation to estimate the eigenfrequency

$$\begin{array}{c}\hfill P_{\mathrm{k}}=2\sqrt{\frac{{\rho }_{\mathrm{ie}}+1}{{\rho }_{\mathrm{ie}}}}{\int }_0^{\mathrm{L}}\frac{\mathrm{d}s}{v_{\mathrm{Ai}}(s)},\end{array}$$(12)

with

$$\begin{array}{c}\hfill v_{\mathrm{Ai}}(s)=\frac{B(s)}{\sqrt{{\mu }_0{\rho }_{\mathrm{i}}}},\end{array}$$(13)

where L represents the loop length, s is the distance from a footpoint along the loop axis, and B(s) is the magnetic field strength along s. Regardless of the damping first, we can roughly estimate the eigenperiod of the current loop by considering Eq. (12). The period can thus be readily obtained as P_k = 369 s. One can find that the horizontal period P_h has a deviation of about 11.6%, compared with P_k, while the vertically polarized mode has a deviation of about 8.9%. It should be noted that these deviations are acceptable, given the assumptions considered in deriving Eq. (11). Detailed discussions regarding the period discrepancy can be found in Sect. 3.4. Here we show that the frequency given by WKB theory is a good approximation to estimate the eigenfrequency of different kink polarizations in a curved loop. Such a small frequency deviation from the WKB approximation in the curved loop model can be regarded as a positive result for coronal seismology since it gives an acceptable error of around 10%, for example in probing the coronal magnetic field based on Eq. (12).

3.3. Resonant damping and lateral leakage

Resonant absorption can be observed in the current oscillating loop. It is well known that the oscillation amplitude decreases due to the resonant energy transfer from the collective kink modes to the local Alfvén waves (e.g., Guo et al. 2020). In Fig. 4, we can indeed see the damping of oscillations in the loop region. Meanwhile, wave leakage can also be captured, as shown in the bottom panels of Figs. 2 and 3. Thus a comparison between the damping efficiency of resonant absorption and lateral leakage is necessary.

Before proceeding, we first revisit the physics of lateral leakage in a curved loop. As discussed in Verwichte et al. (2006a), under a slab model, the oscillation property depends on the slope of the Alfvén frequency profile. The local Alfvén frequency can change with height, due to the variation in magnetic field strength or vertical stratification of density. We now consider the upward leakage. Once the local Alfvén frequency decreases with height and becomes smaller than the wave frequency at a certain region above an oscillating loop, the wave energy leaks out from the loop region via tunneling the evanescent barrier, leading to an oscillatory solution at a larger distance above the oscillating loop. In our current model, the local Alfvén frequency along the z-direction is of the form

$$\begin{array}{c}\hfill {\omega }_{\mathrm{A}}(z)=k(z)v_{\mathrm{Ai}}(z)=\frac{\pi }{L(z)}\frac{B(z)}{\sqrt{{\mu }_0\rho (z)}},\end{array}$$(14)

where the local wavenumber k(z) changes due to the variation of field line length L(z). Here, B(z) = B₀exp( − z/Λ_B) represents the magnetic field strength at a given height z. Given that the magnetic field strength also drops with z, so dω_A/dz < 0. Therefore, we should expect an evanescent barrier and the wave tunneling effect. Figure 5 shows the variation in Alfvén frequency ω_A along the z-direction at the loop apex. The green dashed line labels the predicted oscillating frequency of the loop according to Eq. (12). It reveals an evanescent barrier around 21 Mm ≲ z ≲ 26 Mm and lower Alfvén frequencies than the eigenfrequency when z ≳ 26 Mm. The location of the evanescent barrier in the numerical models may slightly differ from this prediction, due to the deviation of the eigenfrequencies of kink polarizations from the analytical value.

Fig. 5. Alfvén frequency variation along the z-direction at x = 0, y = 0. The cross represents the predicted eigenfrequency of the kink mode in the curved loop, which is calculated from Eq. (12).

The oscillatory patterns shown in the upper loop region are now understandable with the aid of Fig. 5. For a clear illustration, we also considered a time-height map of velocity at y = 0 of the loop apex, as shown in Fig. 6. Oblique stripes above z ∼ 23 Mm can be seen in both polarizations. In Fig. 5 the Alfvén frequency decreases and becomes smaller than the predicted eigenfrequency when z ≳ 26 Mm (it should be z ≳ 23 Mm in the horizontal case due to the eigenfrequency deviation), leading to a longer period with height in the upper loop region. Therefore, we could see the oblique stripes become less vertical with time. Similar cross-field wave propagation can also be observed in, e.g., Kaneko et al. (2015), Raes et al. (2017) due to phase mixing. For the horizontal polarization, the evanescent barrier where the oscillations should be evanescent is mixed with the velocity signals of resonant Alfvén waves. Nonetheless, we can still observe a non-velocity region around z ∼ 23 Mm after about t = 500 s. In the case of vertical polarization, however, the evanescent patterns are not visible until t ≳ 900 s. The evanescent barrier is probably covered by the extension of the external velocity field, namely the wing of kink oscillations in the z-direction.

Fig. 6. Time-distance map of velocity for (a) horizontal and (b) vertical polarization. The evolution of vy and vz over time is sampled at x = 0, y = 0 along the z-direction.

To quantitatively examine the effect of lateral leakage in the current model, a straightforward idea is to isolate the resonant absorption by considering a non-leakage model. To this end, a vertically stratified density distribution is considered. The Eq. (5) is then replaced by

$$\begin{array}{c}\hfill \rho \prime (\overline{r})=\rho {\prime }_{\mathrm{e}}+(\rho {\prime }_{\mathrm{i}}-\rho {\prime }_{\mathrm{e}})f(\overline{r}),\end{array}$$(15)

where

$$\begin{array}{c}\hfill \rho {\prime }_{\mathrm{i}}=2{\rho }_{\mathrm{i}}exp\left[-\frac{{(z/L_0)}^{\mu }}{{\mathrm{\Lambda }}_B/L_0}\right],\rho {\prime }_{\mathrm{e}}=2{\rho }_{\mathrm{e}}exp\left[-\frac{{(z/L_0)}^{\mu }}{{\mathrm{\Lambda }}_B/L_0}\right]\end{array}$$(16)

with L₀ = 10 Mm being the unit of length, and the index μ = 1.8. The density structure of this stratified loop is shown in Fig. 1b. For a more reasonable comparison, this stratified model ensures the same eigenfrequency as the uniform density model. The eigenperiod can be estimated using the WKB approximation given by Eq. (12), but with the internal Alfvén speed defined by

$$\begin{array}{c}\hfill v_{\mathrm{Ai}}(s)=\frac{B(s)}{\sqrt{{\mu }_0\rho (s)}},\end{array}$$(17)

where ρ(s) is the density variation along s. In this model, the local Alfvén frequency along the z-direction is illustrated in Fig. 7. We can see the Alfvén frequency outside the loop is always higher than the eigenfrequency (dashed line), resulting in the absence of the evanescent barrier in this model. In an intuitive manner, we can compare the time-distance map of velocity for both polarizations with Fig. 6, as shown in Fig. 8. Not surprisingly, the patterns in the upper external loop region in Fig. 6 disappear in Fig. 8 for both polarizations, indicating that no wave leakage occurs in this stratified loop. We can thus attribute any damping in this model to resonant absorption.

Fig. 7. Similar to Fig. 5, but for the density stratified loop. The cross marker represents the predicted eigenfrequency of the kink mode according to Eq. (12).

Fig. 8. Similar to Fig. 6, but for the density stratified loop.

It should be noted that we only examine the local Alfvén frequency profile at the loop apex (x = 0). In Eq. (17), we find that the term $B(s)/\sqrt{\rho (s)}$ monotonically decreases with height lower than the loop apex. This implies that the Alfvén frequency at other heights in the loop region is always higher than the eigenfrequency shown in Fig. 7. Thus, the evanescent barrier is indeed absent in the entire loop region.

Figure 9 displays the fitting results for oscillating profiles in all models. Using the similar procedure described in Sect. 3.2, we obtain a period of $P_{\text{h}}^{\prime }=274.8$ s for the horizontal polarization in the stratified loop. This period shows a deviation of about 25.8%, compared with the predicted eigenfrequency P_k. This discrepancy is discussed in detail in Sect. 3.4. An exponential function of the form Aexp( − t/τ) is employed for the fit, represented by the black (leakage model) and green (non-leakage model) dashed lines in Fig. 9. The fitting procedure yields a damping-time-to-period ratio τ_h/P_h = 1.5 (τ_v/P_v = 1.5) for the oscillating loop with lateral leakage and ${\tau }_{\text{h}}^{\prime }/P_{\text{h}}^{\prime }=1.5$ (${\tau }_{\text{v}}^{\prime }/P_{\text{v}}^{\prime }=1.5$) for the loop without wave leakage. For comparison, we use the expression for the damping time due to resonant absorption under the thin-tube-thin-boundary (TTTB) approximation. For a transversely linear density distribution (e.g., Goossens et al. 2009; Soler et al. 2013), it gives,

$$\begin{array}{c}\hfill \tau =\frac{2d}{{\pi }^{2}l}\frac{{\rho }_{\mathrm{ie}}+1}{{\rho }_{\mathrm{ie}}-1}{P}_{\mathrm{k}},\end{array}$$(18)

Fig. 9. Oscillation profiles of the (a) horizontally and (b) vertically polarized modes in all models. The velocity is sampled at loop apex (x = 0, y = 0, z = 19.85 Mm). Oscillations with and without wave leakage are shown in black and green, respectively. The dashed lines represent fitting results with Ah = 3 km s−1 and Av = 4.5 km s−1 in the exponential functions.

where d represents the width of the loop. If we consider the loop width varying from d = a at the footpoint to d = 2a at the loop apex, we can readily estimate the damping-time-to-period ratio to be τ/P_k = 1.01 − 2.03. This implies that our fitted damping time is very close to the analytically derived value for the damping time due to resonant absorption. It should be noted that this is a rough estimation, as the damping rate may vary with different transverse density profiles, as discussed in many previous studies (e.g., Ruderman & Roberts 2002; Goossens et al. 2009). Nonetheless, the same damping-time-to-period ratio in the uniform-density model and the stratified model indicates that the damping effect due to lateral leakage is not significant. This probably means that the damping time of lateral leakage is much longer than that of resonant absorption. This quantitatively confirms that lateral leakage is less efficient as a damping mechanism in curved coronal loops, compared with resonant absorption.

In addition, the same damping-time-to-period ratio for both polarizations indicates that the damping caused by resonant absorption does not have a preferred direction. Thus, the significant difference in the number of observational events of horizontal and vertical kink polarizations may just be due to the preference of possible external exciters (coronal eruptions or ejections as suggested by Zimovets & Nakariakov 2015). In other words, the exciters of kink modes primarily interact with loops in horizontal directions rather than vertical directions. Additionally, the same damping-time-to-period ratio in the stratified model indicates that the rare occurrence of vertical polarizations should not be attributed to the density stratification.

It should be noted that our current analysis only considers the exponential damping stage. Previous studies (e.g., Pascoe et al. 2012, 2016; Hood et al. 2013; Magyar & Van Doorsselaere 2016; Guo et al. 2020) have confirmed that the damping profile of kink modes consists of two stages: the Gaussian stage and the subsequent exponential stage. We can see that the first period of the oscillation profiles in Fig. 9 cannot be perfectly described by exponential functions, probably suggesting a Gaussian damping stage, which is known to robustly manifest in non-axisymmetric loops (Guo et al. 2020). In the current analysis, however, we only focus on the exponential damping profile for a more direct comparison with the thin-tube-thin-boundary (TTTB) expression, as discussed in Guo et al. (2020).

We now compare the efficiency of resonant absorption and lateral leakage from the energetics perspective. As shown above, we would expect a scenario in which the damping time of lateral leakage is much longer than that of resonant absorption. In practice, these timescales can be obtained by comparing the energy and energy flux densities in the computational domain. As discussed in Guo et al. (2020), among others, we can define wave-related energy and energy flux densities for linear perturbations

$$\begin{array}{cc}& E={\int }_V(\frac{1}{2}{\rho }_0{\mathit{v}}_1^2+\frac{{\mathit{B}}_1^2}{2{\mu }_0}+\frac{p_1^2}{2\gamma p_0})\mathrm{d}V,\hfill \end{array}$$(19)

$$\begin{array}{cc}& F={\oint }_S[p_1{\mathit{v}}_1+\frac{{\mathit{B}}_1\times ({\mathit{v}}_1\times {\mathit{B}}_0)}{{\mu }_0}]·\mathrm{d}\mathit{S}.\hfill \end{array}$$(20)

Here the subscript “1” represents the linear perturbation to the equilibrium. V represents the whole computational domain and S is the surface of it. We note that the source term defined by Eq. (B.8) in Guo et al. (2020) is not relevant in the current potential field given an initially constant pressure throughout. In Fig. 10, we illustrate the evolution of total energy and energy flux densities in terms of the second-order small quantities given by Eqs. (19) and (20). For a straightforward comparison we consider the amplitude of energy flux density variations, as shown by the red and blue curves in Fig. 10. In particular, the energy fluxes through open boundaries (except the left symmetric and the bottom reflective boundaries) should induce energy losses associated with lateral leakage. In Fig. 10 the acceptable small difference between the red and blue curves reflects the efficiency of the symmetric and reflective boundaries in practice. We find, for both polarizations, that the total energy E is much larger than the amplitude of the energy flux |F|, and the timescale E/|F| is longer than 10³ s. This means that the damping time associated with lateral leakage is about one order of magnitude longer than that of resonant absorption. Therefore, lateral leakage is less efficient as a damping mechanism, compared with resonant absorption.

Fig. 10. Energy and energy flux densities given by Eqs. (19) and (20) for (a) horizontal and (b) vertical polarizations. Here the subscript “h” (“v”) represents the horizontal (vertical) case. The energy densities for both polarizations are depicted with black lines. Energy flux densities are computed in two ways: by considering all boundaries and by excluding the left symmetric boundary and bottom reflective boundary. These two cases are shown as red and blue curves, respectively.

The decrease in total energy in Fig. 10 is probably numerical. We recomputed the energy E shown in Fig. 10 by considering a lower AMR level (lower numerical resolution) and found a faster decrease with time. This hints that a larger grid spacing leads to an underestimation of the energy in the computational domain, particularly in the non-uniform boundary layer as discussed in Shi et al. (2024b). Nonetheless, the statement that the damping time of lateral leakage is much longer than that of resonant absorption remains unchanged.

3.4. Frequency discrepancies

We now attempt to understand period discrepancies shown up in different models. Seen from Fig. 4, we can find a frequency discrepancy between the two polarizations. Likewise, the oscillation periods for the uniform density model and the stratified model are different as well, as shown in Fig. 9a. These differences can be understood intuitively. For instance, the internal density decreases in the z-direction in the stratified loop, leading to different dynamics across the loop in the ${\widehat{\mathit{e}}}_{\mathrm{v}}$ direction. While the uniform density loop has the same density distribution across the loop in the same direction. Therefore, the period difference between the uniform density model and the stratified model is reasonable since the density distributions perpendicular to the oscillating direction are different. Following this line of thought, we would expect no difference between the two periods of vertical polarization, given that the density has no variation in the y-direction across the loop. This means that the vertically oscillating loop cannot “feel” any difference in the environment from the left to the right in the y-direction. Indeed, Fig. 9b shows that the period of the vertical polarization is $P_{\text{v}}^{\prime }=401.8$ s, which is the same as the period of vertical polarization in the uniform density loop model. This hints that the WKB approximation can only reasonably describe the eigenfrequency of kink modes in a loop without an asymmetric cross-loop density profile perpendicular to the oscillation direction. The asymmetric cross-loop density profile can result in different dynamics from one side of the cross-section to the other, leading to a deviation of the eigenfrequency from analytical predictions.

To further quantitatively understand the period discrepancies, we consider the force analysis by employing the linearized momentum equation,

$$\begin{array}{c}\hfill {\rho }_0\frac{\partial \mathit{v}}{\partial t}=F_{\mathrm{p}}+F_B,\end{array}$$(21)

where

$$\begin{array}{c}\hfill F_{\mathrm{p}}=-\mathrm{\nabla }{p}_1,F_B=\frac{(\mathrm{\nabla }\times {\mathit{B}}_1)\times {\mathit{B}}_0}{{\mu }_0}+\frac{(\mathrm{\nabla }\times {\mathit{B}}_0)\times {\mathit{B}}_1}{{\mu }_0}·\end{array}$$(22)

Here p₁ represents the linear perturbation to the equilibrium pressure, and ${\rho }_0=\rho (\overline{r})$. We calculate the average force density in a large circular region ($\sqrt{{\mathit{y}}^2+{(z-19.85)}^2}<5$ Mm) at the apex (x = 0). Figures 11 and 12 illustrate the profiles of the pressure gradient force density F_p (blue), the magnetic force density F_B (green), and the sum of these two. We observe that for all cases, especially for the leakage model (uniform density loop), the resultant force is very close to the magnetic force until the amplitude of the forces becomes too small to recognize. This means that the magnetic force¹ is still the dominant restoring force. This agrees with what we already know about kink modes. We now focus on the resultant forces that are practically relevant in the current case. By fitting all the resultant force profiles (black curves) in Figs. 11 and 12, we can obtain the corresponding amplitudes. In the uniform density loop, the amplitude of the resultant force for the horizontal polarization is F_h = 2.0 × 10⁻¹¹ N m⁻³, while the corresponding amplitude for the vertical polarization is F_v = 2.3 × 10⁻¹¹ N m⁻³. Considering that the amplitudes of the velocity profiles are A_h = 3 km s⁻¹ and A_v = 4.5 km s⁻¹, respectively, we can readily deduce the frequency difference ω_h/ω_v = 1.29, according to a Fourier decomposition in time of Eq. (21). Therefore, the period difference between the horizontal and vertical polarizations seems understandable. Likewise, we obtain the corresponding amplitude of the fitting for the horizontal polarization is $F_{\text{h}}^{\prime }=2.0\times {10}^{-11}$ N m⁻³ in the density stratified loop. However, we note that the average density at the loop apex surface decreases by about 25%. Therefore, the period deviation of about 25% can be observed for the horizontal polarization in the stratified loops according to Eq. (21). In addition, the amplitude of the resultant force for the vertical polarization is $F_{\text{v}}^{\prime }=1.6\times {10}^{-11}$ N m⁻³ for the stratified loop, and thus the periods show almost no difference between the uniform density model and the stratified model.

Fig. 11. Evolution profiles of the magnetic force density (blue), the pressure gradient force density (green), and the resultant force density of these two (black). The corresponding forces are calculated by Eqs. (21) and (22). The dashed red lines represent the fitting results with Fh = 2.0 × 10−11 N m−3 and Fv = 2.3 × 10−11 N m−3 in the exponential functions.

Fig. 12. Similar to Fig. 11, but for the density-stratified model. The dashed red lines represent the fitting results with $F_{\text{h}}^{\prime }=2.0\times {10}^{-11}$ N m−3 and $F_{\text{v}}^{\prime }=1.6\times {10}^{-11}$ N m−3 in the exponential functions.

4. Discussion and conclusions

We investigated the horizontally and vertically polarized kink eigenmodes excited by initial velocity perturbations in curved coronal loops embedded in a potential magnetic field. In the case of a uniform-density loop, we found that the eigenfrequencies of both polarizations deviate from the WKB approximation by about 10%. This implies that the WKB approximation effectively describes both horizontal and vertical polarizations of the kink eigenmodes in this uniform-density loop. Lateral leakage of both polarizations was clearly observed in the upper loop region as expected. For comparison, we also considered a density-stratified loop without an evanescent barrier, resulting in the absence of wave leakage. The damping-time-to-period ratios of both polarizations in the density-stratified loop were similar to those in the uniform-density loop, indicating that wave leakage does not significantly influence the current damping of both polarizations. All damping rates closely aligned with the predictions of the thin-tube-thin-boundary (TTTB) approximation, indicating that resonant absorption is the dominant damping mechanism in the current curved loops. Examining the oscillating frequency of both polarizations in the stratified loop, we found that the eigenfrequency of the vertical polarization can be described by the WKB approximation, while the frequency of the horizontal polarization deviates from the WKB approximation by around 25%. This implies that the WKB approximation is effective only in describing kink modes in a loop without an asymmetric cross-loop density variation in the direction perpendicular to the polarization, the vertical polarization for instance.

The deviations from the WKB approximation in the frequencies are reasonable. Magyar & Nakariakov (2020) point out that the WKB approximation can well describe the eigenfrequency of horizontal kink oscillations when taking into account the expansion of the loop due to the chosen potential magnetic field and the density stratification. In our current model, however, a deviation of about 11% was found in the horizontal polarization, even if the changes in the loop cross section are included, as seen in Eqs. (11) and (12). This deviation may have several causes. One is that our current loop is not long enough, and thus leads to a deviation from the long wavelength limit (kR ≪ 1 with k being the axial wave number and R being the loop radius) considered in deriving Eq. (10). In addition, we assume a uniform magnetic field across the loop (i.e., pressureless condition) when performing the WKB approximation. Even so, the deviation is only around 10%. This suggests that estimating parameters via seismology may not require measuring the external magnetic field. However, it is worth noting that the measurement of the internal Alfvén speed is still necessary for seismology practices, as indicated in Eq. (12). Directly using Eq. (11) based on the thin-tube (TB) theory by assuming an average Alfvén frequency profile is still possible to induce significant errors.

The frequency discrepancy between the two polarizations in the uniform loop model seems to contradict previous analytical and numerical findings (Ruderman 2003; Guo et al. 2020). In previous studies, a scenario was discussed where a larger (smaller) axis in the loop cross section results in a lower (higher) frequency of kink modes, given that more contribution of the external (internal) Alfvén frequency is involved in the minor (major) axis. In the present model, the radius of the loop cross section in the horizontal direction is fixed, while the radius of the loop cross section in the vertical direction increases from the bottom to the loop apex. However, the frequency of the horizontal polarization is higher than that of the vertical polarization. This contradiction is probably induced by the changing magnetic field strength with height in the present model. In this case, the Alfvén frequency (or Alfvén speed if we neglect the axial wave number variation since the loop is thin) decreases from the lower loop region to the upper loop region, which is the main difference from previous studies (e.g., Ruderman 2003; Guo et al. 2020). In the current model, we expect the variation in Alfvén speed in the cross section to influence the dynamics of the loop. At a given height, the Alfvén speed changes along the minor axis of the cross section (i.e., the ${\widehat{e}}_{\mathrm{v}}$ direction). For the vertical polarization, the Alfvén speed variation in the oscillation direction induces the contraction or extension of the cross section in the minor direction. In the horizontally polarized case, the variation in Alfvén speed is perpendicular to the oscillation direction, thus inducing an oblique minor axis. In addition, the variation in the Alfvén speed along the loop also affects the eigenfunctions of kink modes. Therefore, the dynamics of the present loop are different from canonical elliptical loops. It is not straightforward to make a direct comparison with the scenario discussed in Ruderman (2003) and Guo et al. (2020).

The density-stratified loop is artificially designed for comparison with the uniform-density loop. To isolate the effect of resonant absorption, the density-stratified loop is modified by adjusting the scale height determined by the index μ, as defined in Eq. (16). This modification ensures the absence of an evanescent barrier above the loop due to a higher Alfvén frequency, compared with the uniform-density loop. Meanwhile, the density-stratified loop maintains the same predicted eigenfrequency as given by the WKB approximation. This makes the comparison between the two types of oscillating loops more reasonable. Even though lateral leakage is absent in the density-stratified model, the damping-time-to-period ratio remains the same as in the uniform-density loop. This clearly shows that lateral leakage is less effective than resonant absorption.

A nonlinear regime has not been examined in the current study. Given that both resonant absorption and lateral leakage are linear processes, we thus consider a linear regime with small amplitude perturbations to excite different polarizations, ensuring amplitudes significantly smaller than the loop width. In this manner, the damping of kink polarizations is not influenced by the nonlinearity discussed for example in Goddard & Nakariakov (2016), Magyar & Van Doorsselaere (2016), and Van Doorsselaere et al. (2021), thereby making the comparisons with linear analytical theories more straightforward.

Movies

Movie 1 associated with Fig. 2 (fig2_vh_uni) Access here

Movie 2 associated with Fig. 3 (fig3_vh_uni) Access here

Acknowledgments

The authors acknowledge the funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No. 724326). T.V.D. was also supported by the C1 grant TRACEspace of Internal Funds KU Leuven and a Senior Research Project (G088021N) of the FWO Vlaanderen. Furthermore, TVD received financial support from the Flemish Government under the long-term structural Methusalem funding program, project SOUL: Stellar evolution in full glory, grant METH/24/012 at KU Leuven. The research that led to these results was subsidised by the EU through the DynaSun project (number 101131534 of HORIZON-MSCA-2022-SE-01) and the Belgian Federal Science Policy Office through the contract B2/223/P1/CLOSE-UP. M.G. and B.L. acknowledge support from the National Natural Science Foundation of China (12203030, 12373055, 41974200). The authors thank Norbert Magyar and Konstantinos Karampelas for helpful discussions.

References

All Figures

	Fig. 1. Three-dimensional rendering of the density-enhanced loop. Several arbitrarily selected magnetic field lines are also shown and color-coded by the local magnetic field strength. The left and right panels correspond to our “uniform-density” and “stratified” computations (see text for more details).
In the text

Fig. 2. Velocity fields (arrows) for the horizontal polarization at (a) t = 0, (b) t = 34 s, and (c) t = 661 s. The left and right columns correspond to the y = 0 and x = 0 cuts, respectively. Overplotted are the filled contours of vh = vy. The green curves in the left column represent the magnetic lines of force; the dashed lines outline the density-enhanced loop. These snapshots are from the online animations, which run from t = 0 to 1270 s.

In the text

	Fig. 3. Similar to Fig. 2, except that the filled contours are for $v_{\mathrm{v}}:=\mathit{v}·{\widehat{\mathit{e}}}_{\mathrm{v}}$. An animation is available online and runs from t = 0 to 1270 s.
In the text

	Fig. 4. Oscillation profiles of the horizontally (black) and vertically (blue) polarized modes. The velocity is sampled at x = 0, y = 0, z = 19.85 Mm. The dots represent the extrema of each curve.
In the text

	Fig. 5. Alfvén frequency variation along the z-direction at x = 0, y = 0. The cross represents the predicted eigenfrequency of the kink mode in the curved loop, which is calculated from Eq. (12).
In the text

	Fig. 6. Time-distance map of velocity for (a) horizontal and (b) vertical polarization. The evolution of vy and vz over time is sampled at x = 0, y = 0 along the z-direction.
In the text

	Fig. 7. Similar to Fig. 5, but for the density stratified loop. The cross marker represents the predicted eigenfrequency of the kink mode according to Eq. (12).
In the text

	Fig. 8. Similar to Fig. 6, but for the density stratified loop.
In the text

	Fig. 9. Oscillation profiles of the (a) horizontally and (b) vertically polarized modes in all models. The velocity is sampled at loop apex (x = 0, y = 0, z = 19.85 Mm). Oscillations with and without wave leakage are shown in black and green, respectively. The dashed lines represent fitting results with Ah = 3 km s−1 and Av = 4.5 km s−1 in the exponential functions.
In the text

Fig. 10. Energy and energy flux densities given by Eqs. (19) and (20) for (a) horizontal and (b) vertical polarizations. Here the subscript “h” (“v”) represents the horizontal (vertical) case. The energy densities for both polarizations are depicted with black lines. Energy flux densities are computed in two ways: by considering all boundaries and by excluding the left symmetric boundary and bottom reflective boundary. These two cases are shown as red and blue curves, respectively.

In the text

	Fig. 11. Evolution profiles of the magnetic force density (blue), the pressure gradient force density (green), and the resultant force density of these two (black). The corresponding forces are calculated by Eqs. (21) and (22). The dashed red lines represent the fitting results with Fh = 2.0 × 10−11 N m−3 and Fv = 2.3 × 10−11 N m−3 in the exponential functions.
In the text

	Fig. 12. Similar to Fig. 11, but for the density-stratified model. The dashed red lines represent the fitting results with $F_{\text{h}}^{\prime }=2.0\times {10}^{-11}$ N m−3 and $F_{\text{v}}^{\prime }=1.6\times {10}^{-11}$ N m−3 in the exponential functions.
In the text