Horizontally Polarized Kink Oscillations Supported by Solar Coronal Loops in an Asymmetric Environment

Kink oscillations are ubiquitously observed in solar coronal loops, their understanding being crucial in the contexts of coronal seismology and atmospheric heating. We study kink modes supported by a straight coronal loop embeded in an asymmetric environment using three-dimensional magnetohydrodynamic (MHD) simulations. We implement the asymmetric effect by setting different exterior densities below and above the loop interior, and initiate the simulation using a kink-like velocity perturbation perpendicular to the loop plane, mimicking the frequently measured horizontally polarized kink modes. We find that the external velocity fields show fan blade structures propagating in the azimuthal direction as a result of the successive excitation of higher azimuthal Fourier modes. Resonant absorption and phase mixing can still occur despite an asymmetric environment, leading to the development of small scales at loop boundaries. These small scales nonetheless develop asymmetrically at the upper and lower boundaries due to the different gradients of the Alfven speed. These findings enrich our understanding of kink modes in coronal loops embedded within an asymmetric environment, providing insights helpful for future high-resolution observations.


Introduction
Magnetohydrodynamic (MHD) waves abound in the solar atmosphere (see reviews by, e.g., Li et al. 2020;Van Doorsselaere et al. 2020;Banerjee et al. 2021;Nakariakov et al. 2021;Wang et al. 2021).In particular, kink waves supported by the structured corona have attracted considerable attention since their first imaging detection (Aschwanden et al. 1999;Nakariakov et al. 1999).Further observations show that kink oscillations in coronal loops may be decaying (e.g., Zimovets & Nakariakov 2015;Goddard et al. 2016) or decayless (Tian et al. 2012;Wang et al. 2012;Nisticò et al. 2013;Anfinogentov et al. 2015).Ubiquitous propagating kink modes were also detected in the outer corona (Tomczyk et al. 2007;Van Doorsselaere et al. 2008).Studies of kink modes are often in the context of either coronal heating or coronal magnetic field diagnosis (e.g., Nakariakov & Ofman 2001;Yang et al. 2020).The two contexts are not mutually exclusive.As an example, Goddard & Nisticò (2020) connected the evolution of kink waves with the thermal evolution of their host loops.
A thorough understanding of kink modes is crucial for their related applications.In classical theories in which a straight coronal loop and a piecewise constant loop interior and loop exterior are assumed, kink modes are a collective eigenmode with azimuthal wavenumber q = 1 (e.g., Edwin & Roberts 1983).When a continuous transition layer connects the loop interior and exterior, kink modes are subject to resonant absorption (Goossens et al. 2011).During Movie associated to Fig. 3 is available at https://www.aanda.orgthis process, the coordinated wave energy is transferred to localized Alfvénic motions, which subsequently phase-mix to generate small-scale structures.Along this line of thinking, a variety of modeling works have been performed to explore either the properties of kink modes or their heating effects (e.g., Antolin et al. 2015;Magyar & Van Doorsselaere 2016;Karampelas et al. 2017;Guo et al. 2019Guo et al. , 2020;;Shi et al. 2021a,b;De Moortel & Howson 2022).
In reality, however, the environment of coronal loops is unlikely to be uniform, as is the case with the much-studied horizontally polarized kink modes (Wang et al. 2008).The upper and lower ambient coronae of the loop apex may be asymmetric due to density stratification.While the properties of MHD waves in an asymmetric slab configuration have been examined (e.g., Allcock & Erdélyi 2017;Zsámberger et al. 2018;Chen et al. 2022), there has been no study dedicated to collective waves supported by an asymmetric tube.This work examines the property of kink modes supported by a straight coronal loop in an asymmetric environment using 3D MHD simulations.We see new features in comparison with the case in which a symmetric ambient corona is implemented.Section 2 provides the numerical setup and is followed by the simulation results in Sect.3. Section 4 provides a summary of this work.

Numerical setup
Our equilibrium configuration, illustrated in Fig. 1b, is a straightened version of a curved tube anchored in the photosphere (Fig. 1a).We see coronal loops as field-aligned magnetic tubes.The equilibrium magnetic field is z-directed, and all equilibrium  quantities are z-independent.The equilibrium density is prescribed by ρ(x, y) = ρ e (θ) where (r, θ) denotes a standard polar coordinate system in the x−y plane, and the azimuthal coordinate, θ, is measured counterclockwise from the positive direction of the x-axis.A continuous profile, f (r) = exp[−(r/R) α ], is used to connect the internal density, ρ i = 5 × 10 8 m p cm −3 (m p is the proton mass), to some external density, ρ e (θ), with α = 6 and the nominal loop radius R = 1 Mm.This setup normalizes the length of the transition layer, l, to the loop radius as l/R ≈ 0.6, which is within the range of some seismological estimations based on resonant absorption (e.g., Arregui et al. 2007).We constructed an asymmetric background (hereafter the asymmetric case) by setting where [ρ e1 , ρ e2 ] = [1, 3.5] × 10 8 m p cm −3 .Here ρ e1 represents the density at (x, y) = (0, +∞) and ρ e2 at (x, y) = (0, −∞).
Figures 2a2 and b2 show the equilibrium density profiles along the x-axis (i.e., y = 0) and y-axis (i.e., x = 0) for this case.One sees that the loop environment is asymmetric in the y direction but symmetric in the x direction.The x (y) direction is shown as horizontal (vertical), in accordance with the convention in naming the polarization properties of kink motions (see Fig. 1a).The magnetic field was set to be uniform (15 G).The interior temperature at (x, y) = (0, 0) is T i = 1 MK.A uniform thermal pressure was implemented to enforce transverse force balance, and the resulting plasma 1463,3272,1749] km s −1 .The loop length is L = 50 Mm.A further simulation with the symmetric background (ρ e1 = ρ e2 = 1 × 10 8 m p cm −3 , hereafter the symmetric case) was conducted for comparison (Figs.2a1 and b1).
For both cases, the equilibrium is perturbed by a horizontal kick, where v 0 = 10 km s −1 is the amplitude, and σ = R characterizes the spatial extent.The z dependence in Eq. ( 3) dictates that only axial fundamentals are of interest here.The set of ideal MHD equations was evolved using the PLUTO code (Mignone et al. 2007).We used the piecewise parabolic method for spatial reconstruction, and the second-order Runge-Kutta method for time-stepping.The Harten-Lax-van Leer discontinuities approximate Riemann solver (Miyoshi & Kusano 2005) was used to compute inter-cell fluxes, and a hyperbolic divergence cleaning method was used to keep the magnetic field divergence-free.We simulated only half of the loop, taking advantage of the symmetry with respect to the apex plane (z = L/2).The simulation domain is [−20, 20] Mm × [−20, 20] Mm × [0, 25] Mm in x, y, and z.In the x and y directions, we used 600 uniform grids for [−4, 4] Mm and 200 stretched grids for the other domains.A uniform grid with 50 cells was adopted in the z direction.We used the outflow boundary condition for all primitive variables at the lateral boundaries (x and y).The symmetric boundary condition was applied at the top boundary (z = L/2).For the bottom boundary (z = 0), the transverse velocities (v x and v y ) were fixed at zero.The density, pressure, and B z were fixed at their initial values.The other variables were set as outflows.

Simulation results
Figure 3 shows the density distributions together with the velocity fields in the apex plane at t = 130 s for the symmetric (top panels) and asymmetric (bottom) cases.The left column displays these quantities for a large domain, while the right column zooms in on the inner portion.From Fig. 3 and the associated animation, one sees that the velocity field for the symmetric case is typical of azimuthally standing kink motions.The development of vortical motions, a result of resonant absorption and phase-mixing, is clear at the loop boundary.For the asymmetric case, one sees the development of fan-blade-like motions propagating azimuthally.Furthermore, the velocity fields at the upper and lower loop boundaries are quite different, with small-scale structures much easier to identify at the upper boundary.
We can further see the difference between the velocity fields of the two cases by examining their radial velocities, v r .Figure 4 shows the temporal evolution of v r along a circle of radius r = 2 Mm in the apex plane (the dashed lines in Fig. 3).For the symmetric case (a), azimuthally standing motions are clearly shown.For the asymmetric case (b), however, v r is seen to propagate in the azimuthal direction as the velocity stripes become increasingly inclined with time.Put another way, Fourier modes with increasingly high azimuthal wave numbers appear as time proceeds.
A154, page 2 of 5  We then decomposed the radial velocity, v r , shown in Fig. 4 into different Fourier modes.Following Terradas et al. (2018), this decomposition was performed using the discrete Fourier transform, where G v r (q) represents the coefficient of Fourier mode q, N is the length of the discrete v r sequence, and v r (n) is the nth element.The real (imaginary) part of G v r nearly vanishes for even (odd) q values given that v r is almost symmetric about the yaxis.The Fourier coefficients are displayed in Fig. 5 as functions of q and time.For the symmetric case, only a single Fourier mode of q = 1 is recognizable, with G v r (q = 1) being almost purely real (Fig. 5a).For the asymmetric case, however, Fourier modes with higher q values are excited.Figure 5b displays the real or imaginary parts of G v r (q) for four selected Fourier modes.One finds that the Fourier coefficients of higher q modes reach their maximum at later stages, indicating a coupling among different Fourier modes.The modulus |G v r (q)| as a function of q and time is displayed in the bottom panels for the symmetric (c) and asymmetric (d) cases.Similarly, one sees a single Fourier mode of q = 1 in the symmetric case.For the asymmetric case,  multiple Fourier modes coexist and higher q modes are successively excited as time proceeds.These results demonstrate that the fan-blade-like velocity fields propagating in the azimuthal direction are attributable to the continuous excitation of higher Fourier modes and the coupling among these modes.We verified that wave energy is well contained in the system after a brief transient phase.The kink mode is then subject to damping due to resonant absorption and the subsequent phasemixing (Goossens et al. 2011).This process manifests itself as the vortical motions and the development of small-scale structures at loop boundaries.The top panels of Fig. 6 show the temporal evolution of v x along the y-axis in the apex plane (z = L/2).For the symmetric case (a1), small-scale spikes are symmetrically developing at both sides of the loop boundary.For the asymmetric case, however, spikes develop more rapidly at the upper boundary than at the lower one.This difference is further demonstrated in the bottom panels of Fig. 6, where two snapshots of v x (t = 200 s and 400 s) are displayed.We compared our results with the theoretical prediction given by A154, page 3 of 5 Fig. 7. Temporal profiles of v x (averaged over the loop center with ρ > 0.9ρ i ) for the symmetric (a) and asymmetric (b) cases.The dashed red lines are the best-fit exponential damped sinusoids, with the fitted period, P, and damping time, τ, displayed in each panel.Mann et al. (1995), and our phase-mixing length, L ph , is estimated as where ω A = (π/L)v A (x, y) is the local Alfvén frequency.That small-scale structures develop more rapidly at the upper boundary is therefore due to the stronger gradient in the Alfvén frequency therein.We proceeded more quantitatively by evaluating the Alfvén frequencies [ω Ai , ω Ae1 , ω Ae2 ] ≈ [0.09, 0.21, 0.11] rad s −1 with the corresponding densities.We first examined the Alfvén continuum pertinent to the upper loop boundary.Equation ( 20) from Mann et al. (1995) indicates that the instantaneous phase variation across the continuum is ∆ϕ 1 = (ω Ae1 − ω Ai )t, meaning that the number of spikes is approximately N spikes,1 ≈ ∆ϕ 1 /π.Our setup then yields that ∼7.6 spikes are accumulated every 200 s, which is in line with Fig. 6b2.Likewise, one expects that ∼1.3 extrema will be accumulated every 200 s close to the lower loop boundary, which is also quantitatively reproduced.
The development of small-scale structures is accompanied by the damping of coordinated kink motions.The damping envelope of kink modes has been extensively examined given their seismological potential (e.g., Ruderman & Roberts 2002;Pascoe et al. 2012;Hood et al. 2013).Here we only focused on the difference in terms of the period and damping times between asymmetric and symmetric cases.Figure 7 presents the temporal evolution of v x averaged over the region of ρ > 0.9ρ i for the symmetric (a) and asymmetric cases (b; labeled "sym" and "asym").We fit the profiles using an exponentially damped sinusoid, printing the best-fit periods (P) and damping times (τ).One sees that P asym > P sym , which is unsurprising given the overall heavier ambient corona and hence a larger inertia experienced by the kink motions in the asymmetric case.The fact that τ asym > τ sym can be understood by drawing an analogy with the damping of kink motions in an axisymmetric setup, for which the ideal quasi-mode (QM) theory predicts that τ/P in the thinboundary-thin-tube limit decreases with the ratio of the internal density, ρ i , to the ambient value, ρ amb (e.g., Goossens et al. 1992;Ruderman & Roberts 2002).While it is not straightforward to quantify, ρ amb in the asymmetric case is expected to exceed ρ e1 , which can be unambiguously identified as ρ amb in the symmetric case.That said, an exponential damping envelope does not apply to the symmetric case at late times (Fig. 7a), reinforcing the notion that ideal QMs offer only a partial picture of the timedependent coordinated motions (see Soler & Terradas 2015, for details).The profile of v x in the symmetric case (Fig. 7a) shows a beat at t ≈ 160 s.One possible cause of this beat is the interplay between kink motions and Alfvénic motions (e.g., Fig. 9 of Soler & Terradas 2015).It is known that the presence of similar beats can lead to an underestimation of the actual damping time (Nisticò et al. 2013).We therefore also computed the period (P QM ) and damping time (τ QM ) expected for kink QMs resonantly absorbed in the Alfvén continuum (see the review by Goossens et al. 2011 for conceptual clarifications).We specifically adapted our previous resistive eigenmode code (Chen et al. 2021) to address the present equilibrium configuration, finding P QM ≈ 51.6 s and τ QM ≈ 67.1 s.These QM values are rather close to the best-fit values in Fig. 7a (P sym and τ sym ), suggesting that the presence of a beat-like behavior may not significantly affect our fitting approach for evaluating the kink periods and damping times.
In this work we paid attention to the effect of an asymmetric loop environment on kink oscillations, and thus neglected both gravity and the gravitationally stratified atmosphere to avoid further complexity.The density profile we adopted (Eq.( 2)) may not exactly represent a realistic stratified density profile.In terms of modeling, however, our model is different from most previous simulations, in which a uniform background was assumed.The results shown here are relevant for further exploring more realistic models.
It is necessary to mention that for the asymmetric case, the density variation in the environment (ρ e2 /ρ e1 = 3.5) across the loop width is more significant than the density variations over a scale height (Λ ≈ 50 Mm) of a 1 MK corona.While we did not attempt to model a stratified density profile, we nonetheless performed an additional simulation with ρ e2 = 1.2 × 10 8 m p cm −3 .This density ratio of ρ e2 /ρ e1 = 1.2 is in line with some loops observed by the Transition Region and Coronal Explorer (TRACE) or Solar and Heliospheric Observatory (SOHO), whose width, w, can be as large as 8 Mm (i.e., exp(w/Λ) ≈ 1.2 1 ; e.g., Schrijver 2007).We performed similar analyses, and their results are shown in Figs. 8 and 9.We find that the key features we previously obtained, the azimuthally propagating velocity field and the excitation of higher azimuthal Fourier modes, still exist, though not as strong as in the case with ρ e2 /ρ e1 = 3.5.Our model thus provides a proof-of-concept test of the effect of an asymmetric environment on kink oscillations.More realistic models that take both the stratified atmosphere and the gravity effect into account are called for.

Summary
We studied kink modes supported by a straight coronal loop in an asymmetric environment using 3D MHD simulations.The new features compared with the case of symmetric ambients are summarized as follows: (1) The velocity fields in the exterior show the development of fan-blade features that propagate azimuthally.This is attributed to the successive excitation of higher azimuthal Fourier modes.(2) Resonant absorption and phase-mixing can still occur in the asymmetric case.Velocity spikes are generated nonsymmetrically, with small-scale structures developing more rapidly at the upper boundary than at the lower.This work provides a simplified model that studies horizontally polarized kink modes supported by coronal loops in an asym-metric ambient corona.These findings enrich our understanding of kink modes in more realistic configurations and are helpful for future high-resolution observations. 1

Fig. 1 .
Fig. 1.Illustrations of a curved tube anchored in the photosphere (a) and its straightened version (b), which is used in our configuration.

Fig. 3 .
Fig. 3. Density distributions and velocity fields at the loop apex for the symmetric (top) and asymmetric (bottom) cases at t = 130 s.The dashed green lines mark a circle of radius r = 2 Mm.The right panels show zoomed-in views of the left ones.An animated version of this figure is available online; it has the same layout and runs from 0 to 400 s.

Fig. 4 .
Fig. 4. Temporal evolution of the radial velocity, v r , along the circle delineated in Fig. 3 for the symmetric (a) and asymmetric (b) cases.

Fig. 5 .
Fig. 5. Fourier coefficients, G vr , for the symmetric (left) and asymmetric (right) cases.Top panels show the real or imaginary part of G vr for some selected q as a function of time.Bottom panels show the modulus |G vr | as a function of time and q.