11institutetext:  Shandong Provincial Key Laboratory of Optical Astronomy and
Solar-Terrestrial Environment, Institute of Space Sciences, Shandong
University, Weihai 264209, China
11email<EMAIL_ADDRESS>22institutetext: Centre for Mathematical Plasma
Astrophysics, Department of Mathematics, KU Leuven, Celestijnenlaan 200B,
B-3001 Leuven, Belgium

# Damped kink motions in a system of two solar coronal tubes with elliptic
cross-sections

Mijie Shi 11 Bo Li 11 Shaoxia Chen 11 Hui Yu 11 Mingzhe Guo 1122

(Received ……… / Accepted ………)

###### Abstract

Aims. This study is motivated by observations of coordinated transverse
displacements in neighboring solar active region loops, addressing
specifically how the behavior of kink motions in straight two-tube equilibria
is impacted by tube interactions and tube cross-sectional shapes.

Methods. We work with linear, ideal, pressureless magnetohydrodynamics.
Axially standing kink motions are examined as an initial value problem for
transversely structured equilibria involving two identical, field-aligned,
density-enhanced tubes with elliptic cross-sections (elliptic tubes).
Continuously nonuniform layers are implemented around both tube boundaries. We
numerically follow the system response to external velocity drivers, largely
focusing on the quasi-mode stage of internal flows to derive the pertinent
periods and damping times.

Results. The periods and damping times we derive for two-circular-tube setups
justify available modal results found with the T-matrix approach. Regardless
of cross-sectional shapes, our nonuniform layers feature the development of
small-scale shears and energy accumulation around Alfv$\acute{\rm e}$n
resonances, indicative of resonant absorption and phase-mixing. As with two-
circular-tube systems, our configurational symmetries make it still possible
to classify lower-order kink motions by the polarization and symmetric
properties of the internal flows; hence such mode labels as $S_{x}$ and
$A_{x}$. However, the periods and damping times for two-elliptic-tube setups
further depend on cross-sectional aspect ratios, with $A_{x}$ motions
occasionally damped less rapidly than $S_{x}$ motions. We find uncertainties
up to $\sim 20\%$ ($\sim 50\%$) for the axial Alfv$\acute{\rm e}$n time (the
inhomogeneity lengthscale) if the periods (damping times) computed for two-
elliptic-tube setups are seismologically inverted with canonical theories for
isolated circular tubes.

Conclusions. The effects of loop interactions and cross-sectional shapes need
to be considered when the periods and particularly the damping times are
seismologically exploited for coordinated transverse displacements in adjacent
coronal loops.

###### Key Words.:

magnetohydrodynamics (MHD) — Sun: corona — Sun: magnetic fields — waves

##  1 Introduction

Cyclic transverse displacements of solar coronal loops are arguably the most
extensively observed collective motion in modern solar coronal seismology (see
e.g., Nakariakov & Kolotkov, 2020; Nakariakov et al., 2021, for recent
reviews). Two regimes have been established. The decayless regime is such that
the displacements show little damping and their magnitudes are usually
substantially smaller than visible loop widths, with both features already
clear when this regime was first identified in measurements with Hinode (Tian
et al., 2012) and the Solar Dynamics Observatory/Atmospheric Imaging Assembly
(SDO/AIA; Wang et al., 2012; Nisticò et al., 2013; Anfinogentov et al., 2013).
This regime is known not to be connected with eruptive events but ubiquitous
in active region (AR) loops, as evidenced by a statistical survey of the
SDO/AIA data (Anfinogentov et al., 2015) and by the recent series of analyses
of the measurements with the Extreme Ultraviolet Imager (EUI) on board the
Solar Orbiter (e.g., Zhong et al., 2022; Petrova et al., 2023; Zhong et al.,
2023). Decaying loop displacements, on the other hand, typically damp over
several cycles and are of larger amplitudes as revealed by their first imaging
observations by the Transition Region and Coronal Explorer (TRACE; Schrijver
et al., 1999; Aschwanden et al., 1999; Nakariakov et al., 1999). These
decaying motions are usually associated with lower coronal eruptions (Zimovets
& Nakariakov, 2015). Regardless, there exist ample detections of decaying
displacements, an inexhaustive list being those by Hinode (Ofman & Wang, 2008;
Van Doorsselaere et al., 2008a; Erdélyi & Taroyan, 2008), the Solar
TErrestrial RElations Observatories (STEREO; Verwichte et al., 2009), and by
SDO/AIA (Aschwanden & Schrijver, 2011; White & Verwichte, 2012). Statistical
studies therefore prove possible either through compiling published results
(e.g., Verwichte et al., 2013) or via directly cataloging the SDO/AIA events
(Goddard et al., 2016; Nechaeva et al., 2019). Decaying cyclic displacements
were specifically established to be nearly exclusively axial fundamentals
(e.g., Fig.9 in Goddard et al. 2016 and Fig.5 in Nechaeva et al. 2019).

Seismological applications of cyclic transverse displacements typically start
with their identification as trapped fast kink motions. Practically, this
identification largely relies on the scheme for classifying collective motions
in straight, static, field-aligned configurations where isolated density-
enhanced tubes with circular cross-sections (circular tubes hereafter) are
embedded in an otherwise uniform corona (Edwin & Roberts 1983, ER83; also
Zajtsev & Stepanov 1975; Cally 1986). That the configuration is structured
only transversely and in a one-dimensional (1D) manner means that the
azimuthal wavenumber $m$ makes physical sense, with kink motions corresponding
to $m=1$ (see the textbooks by Roberts 2019 and Goedbloed et al. 2019). Let
“ER83 equilibria” label specifically those where the transverse structuring is
piecewise constant. The relevant theories then enable the measured periods of
axially standing kink motions to be employed to infer the axial
Alfv$\acute{\rm e}$n time and hence the coronal magnetic field strength for
decayless (e.g., Anfinogentov & Nakariakov, 2019; Li & Long, 2023) and
decaying regimes alike (Nakariakov & Ofman 2001; also the reviews by e.g.,
Nakariakov & Verwichte 2005; De Moortel & Nakariakov 2012; Nakariakov et al.
2021). While undamped for ER83 equilibria, trapped fast kink motions are in
general resonantly absorbed in the Alfv$\acute{\rm e}$n continuum when the 1D
structuring is allowed to be continuous (Ruderman & Roberts, 2002; Goossens et
al., 2002, and references therein). Theoretically, the concept of kink quasi-
modes then arises and the internal kink motions are damped in conjunction with
the accumulation and phase-mixing of localized Alfv$\acute{\rm e}$nic motions
(e.g., Poedts & Kerner 1991; Tirry & Goossens 1996; Soler et al. 2013; also
the review by Goossens et al. 2011). Seismologically, resonant absorption has
been customarily invoked to interpret the decay of large-amplitude loop
displacements, thereby allowing the measured damping times to be employed to
deduce such key parameters as the transverse inhomogeneity lengthscales (e.g.,
Aschwanden et al., 2003; Goossens et al., 2008; Arregui & Asensio Ramos, 2011;
Arregui et al., 2015; Arregui, 2022).

Deviations from the canonical ER83 equilibria are known to impact the behavior
of collective motions, and we choose to focus on two geometrical properties
that render the meaning of the azimuthal wavenumber not as clear (e.g., the
review by Li et al., 2020). One concerns loop cross-sections, which may
actually be tied to coronal heating via their key role in determining the
morphology of coronal loops in, say, soft X-ray (Klimchuk et al., 1992) and
EUV (Watko & Klimchuk, 2000) in the first place. Recent imaging observations
with Hi-C (Klimchuk & DeForest, 2020) and Hi-C 2.1 (Williams et al., 2021)
suggest that AR loops may maintain a circular cross-section throughout their
visible segments. However, there also exist suggestions that favor elliptic
cross-sections to better account for the morphology of AR loops, as deduced
with the aid of coronal magnetic field modeling (e.g., Wang & Sakurai, 1998;
Malanushenko & Schrijver, 2013) and/or multi-vantage-point measurements with
STEREO (e.g., McCarthy et al., 2021). In particular, the derived aspect ratios
may readily attain $\sim 1.5-5$ (Malanushenko & Schrijver, 2013), a range also
compatible with the spectroscopic measurements with Hinode/EUV Imaging
Spectrometer (EIS; Kucera et al., 2019). Restrict ourselves to straight
configurations where coronal loops preserve a constant elliptic cross-section
along their axes (elliptic tubes hereafter). Substantial attention has been
paid to collective perturbations in such equilibria from both the modal (e.g.,
Ruderman, 2003; Erdélyi & Morton, 2009; Morton & Ruderman, 2011; Aldhafeeri et
al., 2021) and initial-value-problem (IVP) perspectives (Guo et al., 2020).
Despite the lack of axisymmetry, kink motions remain identifiable as those
that transversely displace the tube axes, and in general they remain subjected
to resonant absorption for continuous transverse structuring. However, as was
first shown by Ruderman (2003), one now needs to discriminate two distinct
polarizations, where the internal flows are primarily directed along the major
(“major-polarized” for brevity) and minor axes (“minor-polarized”),
respectively. In addition, the periods of major-(minor-) polarized kink
motions increase (decrease) with the major-to-minor-axis ratio, while the
damping times for both polarizations tend to increase with this ratio. It was
further demonstrated that the density contrast between elliptic tubes and
their surrounding fluids play a subtle role in mediating the differences
between the dispersive properties of the differently polarized motions (Guo et
al., 2020).

Another geometrical property pertains to systems involving multiple tubular
structures. Observationally, it has long been known that perturbations in such
systems may evolve collectively, an incomplete list of examples being the
coordinated transverse displacements detected either in groups of prominence
threads (e.g., Yi et al. 1991; Lin et al. 2003; Okamoto et al. 2007; also the
review by Arregui et al. 2018) or in neighboring AR loops (e.g., Schrijver &
Brown, 2000; Schrijver et al., 2002; Verwichte et al., 2004; Wang et al.,
2012; White et al., 2013). Kink motions in straight multi-circular-tube
systems have also been extensively examined theoretically. We choose to
concentrate on the modal analyses where a harmonic time-dependence is assumed
a priori, noting that multi-dimensional IVP studies prove equally informative
(e.g., Terradas et al., 2008b; Ofman, 2009; Magyar & Van Doorsselaere, 2016;
Guo et al., 2019a). The first modal examination was due to Luna et al. (2008,
hereafter L08), who numerically solved the pertinent eigenvalue problem (EVP)
for undamped kink modes in a two-tube system. Broadly speaking, two different
approaches have been employed in further modal studies. The ensuing EVPs turn
out to be analytically tractable in the thin-tube (TT) limit for two-tube
equilibria when formulated in bi-cylindrical coordinates, for both undamped
(Van Doorsselaere et al., 2008b; Robertson et al., 2010; Ruderman & Petrukhin,
2023) and damped kink motions (Robertson & Ruderman, 2011; Gijsen & Van
Doorsselaere, 2014). Similar EVPs have also been formulated with the T-matrix
formalism of scattering theory, originally introduced to solar contexts by
Bogdan & Zweibel (1985, 1987). This formalism is sufficiently general for
examining kink motions in composite equilibria comprising an arbitrary number
of circular tubes, with specific investigations available for undamped motions
in systems with two or three (Luna et al., 2009) up to tens of tubes (Luna et
al., 2010, 2019). Addressing continuous transverse structuring, the T-matrix
approach proves capable of handling the resonant absorption of collective
motions in general (Keppens et al., 1994) and particularly that of kink
motions in two-tube equilibria (Soler & Luna, 2015, SL15 hereafter). The
damping is nonetheless assumed to be weak by construction, given that it was
incorporated in the T-matrix framework exclusively via the thin-boundary (TB)
connection formulae (see SL15 for more details; also see Sakurai et al. 1991
for the first derivation of the TB formulae). Regardless, kink motions are
known to be much more complicated in multi-tube equilibria than for isolated
circular tubes, and we restrict ourselves to two-identical-tube systems. It
was shown by L08 (see Fig.2 therein) that kink motions, namely those that
displace both tube axes, need to be classified according to both the
polarization and symmetric properties of the two internal flows. The notions
of $S_{x}$, $A_{x}$, $S_{y}$, and $A_{y}$ then ensue, one subtlety being that
the $S_{x}-A_{y}$ and $A_{x}-S_{y}$ pairs become indistinguishable in terms of
frequencies (Van Doorsselaere et al., 2008b) and damping rates (SL15) when the
tubes are sufficiently thin. This subtlety notwithstanding, $A_{x}$ motions
always turn out to oscillate and damp more rapidly than $S_{x}$ motions: the
differences between the frequencies and particularly the damping rates may be
substantial enough to impact seismological deductions (see Fig.4 in SL15).

This study is intended to examine damped lower-order kink motions in a two-
identical-elliptic-tube system, thereby addressing the situation where neither
the equilibrium configuration as a whole nor an individual tube allows the
azimuthal wavenumber to make exact sense. We will adopt linear, pressureless,
ideal magnetohydrodynamics (MHD) throughout, given that quiescent AR loops are
of primary interest. We will adopt an IVP perspective, paying particular
attention to the frequencies and damping rates of axially standing motions by
largely focusing on the duration where the concept of quasi-modes applies. Our
study is new in the following two aspects. Firstly, two-elliptic-tube systems
have yet to be explored in the context of collective waves, meaning that our
study will shed new light on how the dispersive properties of kink motions are
impacted by the joint effects of tube interactions and tube cross-sectional
shapes. When put to seismological contexts, our results can therefore help
assess the uncertainties in the key physical parameters that one deduces with
the customary practice where the joints effects are absent. Secondly, there
exist no IVP studies dedicated to the dispersive properties of damped kink
motions in two-tube systems to our knowledge. Our numerical analysis,
conducted with a self-developed code, will therefore help verify the T-matrix
results obtained by SL15.

The outline of this manuscript is as follows. Section 2 offers the
mathematical formulation of our IVP together with a description of our
numerical code. We then focus on circular tubes in Sect.3, testing our code
outputs against available results obtained with independent approaches.
Section 4 proceeds to present our numerical results for two-elliptic-tube
systems. The findings of this study are summarized in Sect.5, where some
concluding remarks are also given.

##  2 Problem Formulation and Solution Method

###  2.1 Governing Equations

We adopt pressureless, ideal MHD as our theoretical framework, in which the
primitive dependents are the mass density $\rho$, velocity $\@vec{v}$, and
magnetic field $\@vec{B}$. Let $(x,y,z)$ be a Cartesian system, and let the
equilibrium quantities be denoted with a subscript $0$. We take the
equilibrium magnetic field to be uniform and $z$-directed
($\@vec{B}_{0}=B_{0}\@vec{e}_{z}$). Only straight, static, field-aligned
configurations are of interest, meaning that $\@vec{v}_{0}=0$ and the
structuring is encapsulated in $\rho_{0}(x,y)$. The Alfv$\acute{\rm e}$n speed
$v_{\rm A}$ is defined by $v_{\rm A}^{2}=B_{0}^{2}/\mu_{0}\rho_{0}$, where
$\mu_{0}$ is the magnetic permeability of free space. With magnetically closed
structures in mind, we place two dense photospheres at $z=0$ and $z=L$.

We formulate an IVP to examine how our equilibrium responds to small-amplitude
perturbations (denoted by subscript $1$). It follows from linearized, ideal,
pressureless MHD equations that

$\displaystyle\rho_{0}\dfrac{\partial v_{1x}}{\partial t}$ $\displaystyle=$
$\displaystyle\dfrac{B_{0}}{\mu_{0}}\left(\dfrac{\partial B_{1x}}{\partial
z}-\dfrac{\partial B_{1z}}{\partial x}\right),$ (1)
$\displaystyle\rho_{0}\dfrac{\partial v_{1y}}{\partial t}$ $\displaystyle=$
$\displaystyle\dfrac{B_{0}}{\mu_{0}}\left(\dfrac{\partial B_{1y}}{\partial
z}-\dfrac{\partial B_{1z}}{\partial y}\right),$ (2)
$\displaystyle\dfrac{\partial B_{1x}}{\partial t}$ $\displaystyle=$
$\displaystyle B_{0}\dfrac{\partial v_{1x}}{\partial z},$ (3)
$\displaystyle\dfrac{\partial B_{1y}}{\partial t}$ $\displaystyle=$
$\displaystyle B_{0}\dfrac{\partial v_{1y}}{\partial z},$ (4)
$\displaystyle\dfrac{\partial B_{1z}}{\partial t}$ $\displaystyle=$
$\displaystyle-B_{0}\left(\dfrac{\partial v_{1x}}{\partial x}+\dfrac{\partial
v_{1y}}{\partial y}\right).$ (5)

We focus on axially standing motions by adopting the ansatz

$\begin{split}v_{1x}(x,y,z;t)&=\hat{v}_{x}(x,y;t)\sin(kz),\\\
v_{1y}(x,y,z;t)&=\hat{v}_{y}(x,y;t)\sin(kz),\\\
B_{1x}(x,y,z;t)&=\hat{B}_{x}(x,y;t)\cos(kz),\\\
B_{1y}(x,y,z;t)&=\hat{B}_{y}(x,y;t)\cos(kz),\\\
B_{1z}(x,y,z;t)&=\hat{B}_{z}(x,y;t)\sin(kz),\end{split}$ (6)

where $k=n\pi/L$ is the quantized axial wavenumber $(n=1,2,\cdots)$. Equations
(1) to (5) then become

$\displaystyle\dfrac{\partial\hat{v}_{x}}{\partial t}$ $\displaystyle=$
$\displaystyle-\dfrac{B_{0}}{\mu_{0}\rho_{0}}\left(k\hat{B}_{x}+\dfrac{\partial\hat{B}_{z}}{\partial
x}\right),$ (7) $\displaystyle\dfrac{\partial\hat{v}_{y}}{\partial t}$
$\displaystyle=$
$\displaystyle-\dfrac{B_{0}}{\mu_{0}\rho_{0}}\left(k\hat{B}_{y}+\dfrac{\partial\hat{B}_{z}}{\partial
y}\right),$ (8) $\displaystyle\dfrac{\partial\hat{B}_{x}}{\partial t}$
$\displaystyle=$ $\displaystyle kB_{0}\hat{v}_{x},$ (9)
$\displaystyle\dfrac{\partial\hat{B}_{y}}{\partial t}$ $\displaystyle=$
$\displaystyle kB_{0}\hat{v}_{y},$ (10)
$\displaystyle\dfrac{\partial\hat{B}_{z}}{\partial t}$ $\displaystyle=$
$\displaystyle-B_{0}\left(\dfrac{\partial\hat{v}_{x}}{\partial
x}+\dfrac{\partial\hat{v}_{y}}{\partial y}\right).$ (11)

Axial fundamentals ($k=\pi/L$) will be examined throughout, even though our
analysis can be readily adapted to any axial harmonic number $n$.

###  2.2 Energy Conservation Law

Energetics considerations turn out to be necessary. Consider a volume $V$ that
laterally occupies an arbitrary area $Q$ and is axially bounded by the planes
$z=0$ and $z=L$. Let the curve enclosing $Q$ be denoted by $\partial Q$. Let
$[\@vec{e}_{n},\@vec{e}_{t},\@vec{e}_{z}]$ further define a right-handed set
of orthonormal system at any point along $\partial Q$, with the normal
direction $\@vec{e}_{n}$ pointing away from $Q$. An energy conservation law
can be readily derived from Eqs. (1) to (5), reading

$\dfrac{{\rm d}E(t)}{{\rm d}t}=-F(t),$ (12)

where

$\displaystyle E(t)=\int_{Q}\epsilon dxdy,$ (13) $\displaystyle
F(t)=\oint_{\partial Q}dl_{t}\@vec{e}_{n}\cdot\hat{p}_{\rm
T}\left(\hat{v}_{x}\@vec{e}_{x}+\hat{v}_{y}\@vec{e}_{y}\right),$ (14)
$\displaystyle\epsilon(x,y;t)=\dfrac{1}{2}\rho_{0}\left(\hat{v}_{x}^{2}+\hat{v}_{y}^{2}\right)+\dfrac{1}{2\mu_{0}}\left(\hat{B}_{x}^{2}+\hat{B}_{y}^{2}+\hat{B}_{z}^{2}\right).$
(15)

Note that the ansatz (6) has been employed, and $\hat{p}_{\rm
T}=B_{0}\hat{B}_{z}/\mu_{0}$ is connected with the Eulerian perturbation of
total pressure. Note also that $dl_{t}$ denotes the elementary arclength in
the tangential direction defined by $\@vec{e}_{t}$. When multiplied by $L/2$,
the symbols $E$ and $F$ represent the instantaneous perturbation energy in $V$
and the net energy flux out of $V$, respectively. We nonetheless refer to $E$
($F$) as the total energy (the net energy flux) for brevity. Likewise, the
symbol $\epsilon$ in Eqs. (13) and (15) will be referred to as the energy
density. We also drop the hat from here onward.

###  2.3 Equilibrium Configuration and Initial Perturbation

Our equilibrium comprises a composite structure embedded in an otherwise
uniform external corona (denoted by subscript ${\rm e}$). By “composite” we
mean two identical tubes with elliptic cross-sections and separated by $d$.
Let $j$ label a tube, and let a 2D position vector
$\@vec{X}_{j}=(X_{j},Y_{j})$ denote the center of tube $j$. The tube centers
are placed on the $x$-axis without loss of generality, enabling a tube to be
referred to as either the left ($j={\rm L}$) or the right tube ($j={\rm R}$).
The tube centers are further taken to be symmetric about $x=0$, resulting in
$X_{\rm L,R}=\mp d/2$ and $Y_{\rm L,R}=0$. We consistently use subscript ${\rm
i}$ to denote the equilibrium quantities at either tube axis, meaning in
particular that the internal density and Alfv$\acute{\rm e}$n speed are
denoted by $\rho_{\rm i}$ and $v_{\rm Ai}$, respectively. Likewise, by $v_{\rm
Ae}$ we denote the Alfv$\acute{\rm e}$n speed evaluated with the external
density $\rho_{\rm e}$.

Our two-tube structure is characterized as follows. The tubes are taken to be
identical not only in cross-sectional shapes but in their orientations
relative to the $x$-axis. For simplicity, we consider only two orientations by
discriminating whether the major (hereafter “$x$-major”) or minor axis
(“$x$-minor”) of a tube aligns with the $x$-direction. Let $a_{x}$ and $a_{y}$
denote the spatial extent that a tube spans in the $x$\- and $y$-directions,
respectively. The $x$-major ($x$-minor) orientation then means that
$a_{x}>a_{y}$ ($a_{x}<a_{y}$). Both orientations are realized by the density
distribution ($j={\rm L}$ or ${\rm R}$),

$\displaystyle\rho_{j}(x,y)=\left\\{\begin{array}[]{lc}\rho_{\rm
i},&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
0\leq\bar{r}_{j}\leq 1-\bar{l};\\\\[5.69046pt] \dfrac{\rho_{\rm i}+\rho_{\rm
e}}{2}-\dfrac{\rho_{\rm i}-\rho_{\rm
e}}{2}\cos\dfrac{\pi(\bar{r}_{j}-1)}{\bar{l}},&\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ 1-\bar{l}\leq\bar{r}_{j}\leq
1,\end{array}\right.$ (18)

with the intermediate dimensionless variable $\bar{r}_{j}$ defined by

$\bar{r}_{j}(x,y)\coloneqq\sqrt{\left(\dfrac{x-X_{j}}{a_{x}}\right)^{2}+\left(\dfrac{y}{a_{y}}\right)^{2}}.$
(19)

We take $d\geq 2a_{x}$ such that tube $j$ can be unambiguously identified as
where $\bar{r}_{j}\leq 1$. We also use $a$ ($b$) to denote the semi-major
(semi-minor) axis, meaning that $[a_{x}=a,a_{y}=b]$ ($[a_{x}=b,a_{y}=a]$) for
the $x$-major ($x$-minor) orientation. The symbol $b$ will be favored when the
limiting case of circular tubes is examined ($a=b$). Equation (18) represents
a tube profile that continuously connects the internal density $\rho_{\rm i}$
to the external one $\rho_{\rm e}$ via some elliptic layer of width $\bar{l}a$
($\bar{l}b$) along the direction of the major (minor) axis. When $a=b$, this
profile is equivalent to the sinusoidal distribution introduced by Ruderman &
Roberts (2002, hereafter RR02) for modeling circular inhomogeneities (see
e.g., Van Doorsselaere et al., 2004; Soler et al., 2013; Chen et al., 2021,
for more applications). Somehow subtle is that the RR02 implementation for a
tube centered at $(X_{j},0)$ writes

$\displaystyle\rho_{j}(x,y)=\left\\{\begin{array}[]{lc}\rho_{\rm
i},&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
{r}_{j}\leq\mathcal{R}-\ell/2;\\\\[5.69046pt] \dfrac{\rho_{\rm i}+\rho_{\rm
e}}{2}-\dfrac{\rho_{\rm i}-\rho_{\rm
e}}{2}\sin\dfrac{\pi({r}_{j}-\mathcal{R})}{\ell},&\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\mathcal{R}-\ell/2\leq{r}_{j}\leq\mathcal{R}+\ell/2,\\\\[5.69046pt] \rho_{\rm
e},&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
{r}_{j}\geq\mathcal{R}+\ell/2.\end{array}\right.$ (23)

Here $r_{j}\coloneqq\sqrt{(x-X_{j})^{2}+y^{2}}$. Evidently, the geometrical
parameters $[\mathcal{R},\ell]$ characterizing the RR02 implementation are
connected to ours via

$\begin{split}&b=\mathcal{R}+\ell/2,\quad\\\
&\bar{l}=\dfrac{\ell}{b}=\dfrac{\ell/\mathcal{R}}{1+\ell/2\mathcal{R}}.\end{split}$
(24)

Our equilibrium configuration, illustrated in Fig. 1b, is a straightened
version of a composite system where two tubes are separated in the
$x$-direction and anchored in the photosphere (Fig. 1a). We see the $x$\-
($y$-) direction in Fig. 1b as horizontal (vertical) given that the $y-z$
plane can be identified as the tube plane. For two-tube systems, however, it
turns out to be also necessary to classify lower-order kink motions into the
$S_{x}$, $A_{x}$, $S_{y}$, and $A_{y}$ types by combining the symmetric and
polarization properties of the internal flow fields (see Fig. 2 in L08).
Proposed for two-identical-circular-tube systems, the subscript $x$ ($y$)
arises when the internal velocities are dominated by their $x$\-
($y$-)components, while the symbol $S$ ($A$) is such that the dominant
internal velocity component is symmetric (anti-symmetric) about $x=0$. This
classification is expected to carry over when elliptic tubes are examined.
Somehow complicated is that one needs to distinguish between the $x$-major and
$x$-minor arrangements (Figs. 1c and 1d), meaning for instance that the
frequencies of the $S_{x}$ motions are likely to be different for different
tube orientations.

Our equilibrium is perturbed via velocities. We consider only $S_{x}$ and
$A_{x}$ motions, given that kink motions are commonly observed to be
horizontally polarized (see Zhong et al. 2023 and references therein; see
e.g., Wang et al. 2008 and Aschwanden & Schrijver 2011 for observational
instances of vertical kink motions). Regardless of tube orientations, these
are excited by imposing

$\dfrac{\@vec{v}_{\rm ini}(x,y)}{v_{\rm Ae}}={\rm
e}^{4}\left\\{\exp\left[-\left(\dfrac{x+x_{\rm
ini}}{b}\right)^{2}\right]\pm\exp\left[-\left(\dfrac{x-x_{\rm
ini}}{b}\right)^{2}\right]\right\\}\@vec{e}_{x},$ (25)

where $x_{\rm ini}=X_{\rm R}+a_{x}+2b$, and we recall that $[X_{\rm
R}=d/2,b=\min(a_{x},a_{y})]$. The coefficient ${\rm e}^{4}$ is introduced only
for plotting purposes, the magnitude being immaterial for linear studies.
Equation (25) represents a pair of planar compressible perturbations
concentrated in the ambient corona, exciting $S_{x}$ ($A_{x}$) motions when
the plus (minus) sign applies.

###  2.4 Parameter Overview and Solution Method

The behavior of small-amplitude perturbations is determined by two sets of
parameters. We take the dimensional set to be $\\{\rho_{\rm e},b,v_{\rm
Ae}\\}$, which serves merely as normalizing constants in our context. Accept
that we are primarily interested in such timescales as the oscillation period
and some damping time. The importance of the dimensionless set then follows
from a straightforward dimensional analysis, which dictates that a timescale
$t_{\rm scale}$ is formally expressible as

$\dfrac{v_{\rm Ae}t_{\rm
scale}}{b}=\mathcal{F}\left[\mbox{orientation},\dfrac{\rho_{\rm i}}{\rho_{\rm
e}},\dfrac{a}{b},\bar{l},\dfrac{d}{b};\quad kb;\quad\mbox{perturbation
pattern}\right].$ (26)

The subgrouped parameters characterize the transverse structuring, axial
wavenumber, and initial perturbation, respectively. By “orientation” we refer
to either “$x$-major” or “$x$-minor”. By “perturbation pattern” we mean either
$S_{x}$ or $A_{x}$. We fix the axial wavenumber at $kb=\pi/30$ throughout,
meaning a tube length $L=30b$ for axial fundamentals. This $L/b$ suffices for
our purposes, lying toward the lower end of but within the accepted range for
AR loops imaged in the EUV (e.g., Aschwanden et al., 2004; Schrijver, 2007).
Table 1 briefly overviews our computations, which are grouped into four sets
to be detailed later. The density contrast $\rho_{\rm i}/\rho_{\rm e}$ is
fixed at either $3$ or $5$, both values being representative of AR loops
(e.g., Aschwanden et al., 2004, and references therein). It remains formidable
to exhaust the effects of the rest of the parameters, and we choose to fix
many.

We solve Eqs. (7) to (11) with the following procedure. A computational domain
$[-x_{\rm M},x_{\rm M}]\times[-y_{\rm M},y_{\rm M}]$ is discretized into a
uniform mesh with identical spacing in the $x$\- and $y$-directions
($\Delta_{x}=\Delta_{y}=\Delta$). The zero-gradient condition is implemented
for all unknowns at all boundaries. We evolve Eqs. (7) to (11) with the
classic MacCormack scheme (MacCormack, 1969), a popular finite-difference
algorithm second-order accurate in both space and time (see the textbooks by,
e.g., Anderson, 1995; Jardin, 2010, for more). The time step $\Delta_{t}$ is
set via some effective Courant number $c\coloneqq(v_{\rm
Ae}\Delta_{t}/\Delta_{x})^{2/3}+(v_{\rm Ae}\Delta_{t}/\Delta_{y})^{2/3}$ as
inspired by Hong (1996). A value of $c=0.99$ is employed for all the presented
results. We have verified that varying $c$ between $\sim 0.7$ and $\sim 1.3$
introduces no discernible difference, and numerical stability is consistently
maintained.

Some remarks on the grid setup are necessary. We start by noting that the
spatial spacing $\Delta$ restricts the timeframe in which numerical solutions
make physical sense, an aspect raised by Terradas et al. (2008b, hereafter
T08) who numerically solved an equivalent set of governing equations. Defining
the Alfv$\acute{\rm e}$n frequency $\omega_{\rm A}(x,y)=kv_{\rm A}(x,y)$, one
expects that some time-dependent phase-mixing length will emerge as (Mann et
al., 1995)

$L_{\rm ph}(x,y;t)=\dfrac{2\pi}{|\nabla\omega_{\rm A}(x,y)|t},$ (27)

which characterizes the transverse lengthscales of resonantly generated
Alfv$\acute{\rm e}$nic motions in the nonuniform portions in the system. It
suffices to consider only the nonuniform layer of one tube. Evidently, the
shortest phase-mixing length $L_{\rm ph}^{\rm min}(t)$ at a given instant
occurs at the strongest $|\nabla\omega_{\rm A}(x,y)|$. Equally evident is that
this strongest $|\nabla\omega_{\rm A}(x,y)|$ depends only on $\rho_{\rm
i}/\rho_{\rm e}$ and $\bar{l}$ even when elliptic tubes are examined. Let
$L_{\rm ph}^{\rm min}(500)$ denote $L_{\rm ph}^{\rm min}$ at $t=500b/v_{\rm
Ae}$, before which our computations are consistently terminated. It then
follows from the arguments by T08 that our signals are physically relevant
provided $L_{\rm ph}^{\rm min}(500)\geq 2\Delta$. As shown by Table 1, this
criterion is satisfied by the reference grid setup for any set of our
computations. Further justification of the spacing $\Delta$, together with
that for our reference domain size $[x_{\rm M},y_{\rm M}]$, will be presented
shortly.

Table 1: Summary of time-dependent solutions presented in the text

Set | $\rho_{\rm i}/\rho_{\rm e}$ | $\bar{l}$ | $L_{\rm ph}^{\rm min}(500)$ | Reference Grid | Section | Remarks
---|---|---|---|---|---|---
1 | $3$ | $\geq 0.2$ | $\gtrsim 0.032b$ | $[x_{\rm M},y_{\rm M}]=[12b,12b],\Delta=0.01b$ | §3.1 |  1-circular-tube
2 | $5$ | $\gtrsim 0.182$ | $\gtrsim 0.019b$ | $[x_{\rm M},y_{\rm M}]=[12b,10b],\Delta=0.005b$ | §3.2 |  2-circular-tube
3 | $3$ | $0.4$ | $0.063b$ | $[x_{\rm M},y_{\rm M}]=[30b,25b],\Delta=0.01b$ | §4.1 |  2-elliptic-tube $x$-major
4 | $3$ | $0.4$ | $0.063b$ | $[x_{\rm M},y_{\rm M}]=[30b,25b],\Delta=0.01b$ | §4.2 |  2-elliptic-tube $x$-minor

111All computations pertain to $kb=\pi/30$ or equivalently $L/b=30$ for axial
fundamentals.

##  3 Test Computations for Circular Tubes

The MacCormack scheme, while a textbook one, has not been applied to Eqs. (7)
to (11). Its applicability is therefore examined in this section via some test
computations for which the time-dependent behavior can be established or
expected with independent methods. Only circular tubes are of interest, and we
refer to $b$($=a$) as our tube radius. We start by examining the response of a
one-circular-tube system to axisymmetric (sausage) perturbations, following
our previous study (Li et al., 2022) to formulate both the equilibrium and the
initial perturbation. The Fourier-integral-based solutions, presented in Fig.
4 therein, agree closely with our MacCormack solutions. In particular, no
discernible numerical anisotropy shows up even though finite differences are
performed on a Cartesian grid to examine a non-planar equilibrium (see e.g.,
the review by Sescu, 2015, for more on numerical anisotropy). The rest of this
section focuses on kink perturbations.

###  3.1 Kink Motions in a One-Circular-Tube Setup

This subsection examines kink perturbations in a one-circular-tube
configuration, placing the tube center at the origin without loss of
generality. The equilibrium density is realized through Eqs. (18) and (19) by
taking $[a_{x}=a_{y}=b,X_{j}=0]$ and then dropping the subscript $j$. Kink
motions in such a configuration have been extensively studied (see the review
by Nakariakov et al., 2021, and references therein), readily allowing our
MacCormack computations to be compared with known results. By “known” we refer
to two sets of studies. Set one, presented in Sect. 4 of T08, adopts the RR02
implementation to examine the system response to an initial perturbation of
the form

$\@vec{v}_{\rm ini}(x,y)=v_{\rm Ae}\exp\left[-\dfrac{(y-y_{\rm
<EMAIL_ADDRESS>(28)

Figure 2 in T08 then presents the time sequence of $v_{y}$ sampled at the tube
center for a combination of parameters $[\rho_{\rm i}/\rho_{\rm
e}=3,k\mathcal{R}=\pi/20,\ell/\mathcal{R}=0.6]$ and $[y_{\rm
ini}=3\mathcal{R},\sigma=\mathcal{R}]$. We repeat the same experiment, using
Eq. (24) to address notational differences and adopting an identical
computational grid. Our finite-difference results are found to be consistent
with T08, despite that the numerical code therein adopts the finite-volume
methodology.

We proceed to perform an additional series of computations to examine whether
our code outputs agree with expectations for ideal kink quasi-modes. Let
$\omega=\Omega-i\gamma$ denote the complex-valued quasi-mode frequency, and
see $\Omega$ and $\gamma$ as positive. With the kink speed $c_{\rm kink}$
defined as

$c_{\rm kink}^{2}\coloneqq\dfrac{\rho_{\rm i}v_{\rm Ai}^{2}+\rho_{\rm e}v_{\rm
Ae}^{2}}{\rho_{\rm i}+\rho_{\rm e}},$ (29)

it is well established that (RR02; see also Goossens et al. 1992 and Soler et
al. 2013)

$\begin{split}&\Omega\approx kc_{\rm kink},\\\\[5.69046pt]
&\dfrac{\gamma}{\Omega}\approx\dfrac{1}{4}\dfrac{\bar{l}}{1-\bar{l}/2}\dfrac{\rho_{\rm
i}/\rho_{\rm e}-1}{\rho_{\rm i}/\rho_{\rm e}+1}\end{split}$ (30)

in the so-named thin-tube-thin-boundary limit (TTTB, $kb\ll 1$ and $\bar{l}\ll
1$). Note that the TTTB expressions are usually formulated in terms of the
RR02 notations $[\mathcal{R},\ell]$. Note also that the TTTB results may not
hold well beyond their range of applicability (e.g., Soler et al., 2014; Yu et
al., 2021). We therefore also evaluate the mode frequencies numerically with
the general-purpose finite-element code PDE2D (Sewell, 1988), computing ideal
kink quasi-modes as resistive eigenmodes (see Goossens et al. 2011 for
conceptual clarifications; see also Terradas et al. 2005 for the first
introduction of PDE2D to solar contexts). No restriction is necessary for $kb$
or $\bar{l}$.

The rest of this subsection is devoted to a fixed combination $[\rho_{\rm
i}/\rho_{\rm e}=3,kb=\pi/30]$, allowing only $\bar{l}$ to vary. Furthermore,
the reference grid (see Table 1) is consistently employed in our time-
dependent computations, where kink motions are excited by a perturbation of
the fixed form

$\dfrac{\@vec{v}_{\rm ini}(x,y)}{v_{\rm Ae}}={\rm
<EMAIL_ADDRESS>(31)

Equation (31) is essentially identical to Eq. (28) except that $v_{x}$ rather
than $v_{y}$ is perturbed. Figure 2a examines a representative case with
$\bar{l}=0.4$, showing the temporal evolution of the $x$-speed sampled at the
tube center ($v_{x}(0,0,t)$, the solid curve). This $v_{x}$ signal is seen to
feature some rapid variations for $t\lesssim 15b/v_{\rm Ae}$. By “rapid” we
mean timescales on the order of the transverse Alfv$\acute{\rm e}$n time
($b/v_{\rm Ai}$), which derives from the multiple reflections off the tube
boundaries of the perturbations that are imparted by the external driver to
the tube. The $v_{x}$ signal transitions toward a regular slower-varying
pattern afterwards, becoming monochromatic when $t\gtrsim 70b/v_{\rm Ae}$. We
choose to leave out the first extremum in this stage for safety, and fit the
segment encompassing the next six (the red asterisks) with an exponentially
damped cosine

$f_{\rm fit}(t)\propto\exp{\left(-\dfrac{t}{\tau}\right)}\cos\left(\dfrac{2\pi
t}{P}+\phi\right).$ (32)

The damping envelope from the best-fit is then plotted for the entire duration
in Fig. 2a by the dashed curves. One see that this best-fit envelope, while
deduced for some segment, offers a good description for the larger-time
behavior as well (say, $t\gtrsim 300b/v_{\rm Ae}$).

Figure 2b surveys a range of $\bar{l}$ by plotting the oscillation frequencies
($\Omega_{\rm fit}$, the black open circles) and damping rates ($\gamma_{\rm
fit}$, blue), which are translated from the best-fit periods and damping times
($\Omega_{\rm fit}=2\pi/P_{\rm fit}$ and $\gamma_{\rm fit}=1/\tau_{\rm fit}$).
Note that $\Omega$ is measured in units of the external Alfv$\acute{\rm e}$n
frequency $\omega_{\rm Ae}=kv_{\rm Ae}$, and $\gamma/\Omega$ is presented
rather than $\gamma$ itself. The curves in Fig. 2b further provide the quasi-
mode expectations with either the analytical TTTB expression (Eq. 30, the
solid lines) or the PDE2D computations (labeled “Resistive”, dashed). One sees
that the TTTB results provide a rather good approximation to the numerical
“Resistive” ones, the frequencies being practically the same and the damping
rates differing by $\gtrsim 10\%$ only when $\bar{l}\gtrsim 0.4$. One further
sees that the open circles agree well with the “Resistive” computations,
demonstrating that the internal flow fields practically evolve as an ideal
quasi-mode at least in the interval where the fitting is performed.

###  3.2 Kink Motions in a Two-Circular-Tube Setup

This subsection examines kink motions in a two-circular-tube configuration,
for which the distinction between $S_{x}$ and $A_{x}$ patterns becomes
necessary. The initial perturbation is therefore chosen to follow Eq. (25),
where $x_{\rm ini}$ evaluates to $d/2+3b$. Our code consistently adopts the
$[b,\bar{l}]$ notations. However, the mixed usage of the RR02 convention
$[\mathcal{R},\ell]$ turns out to be necessary. The reason is that this RR02
convention was adopted by SL15 to examine ideal quasi-modes from an
eigenvalue-problem perspective in a configuration physically identical to
ours. Briefly put, the study by SL15 is based on the TB-embedded T-matrix
formalism of scattering theory, the pertinent results being most
straightforward for our time-dependent computations to be compared with. By
“TB-embedded” we recall that the formalism therein addresses quasi-mode
damping with the Thin-Boundary connection formulae. By “pertinent” we
specifically refer to Fig. 4 therein, which presents the
$\ell/\mathcal{R}$-dependencies of the quasi-mode frequencies and damping
rates for a fixed combination $[\rho_{\rm i}/\rho_{\rm
e}=5,d/\mathcal{R}=2.5,k\mathcal{R}=\pi/100]$. Note that varying
$\ell/\mathcal{R}$ actually impacts our $\bar{l}$, $d/b$, and $kb$
simultaneously (see Eq. 24). We adopt a fixed $kb=\pi/30$ to save
computational time. This does not matter because we will measure our
timescales or frequencies in appropriate units (say, $\Omega$ in $\omega_{\rm
Ae}$) and the resulting readings do not depend on $k\mathcal{R}$ when
$k\mathcal{R}\lesssim\pi/20$ (see Fig. 3 in SL15, and note that
$k\mathcal{R}<kb$). Equation (24) readily converts the pair
$[\ell/\mathcal{R},d/\mathcal{R}]$ into our $[\bar{l},d/b]$.

####  3.2.1 Comparison with Expectations from T-matrix Theory

This subsection tests our time-dependent computations against the SL15
results, to be labeled “T-matrix” for clarity. It suffices to examine the time
sequences of the $x$-speed sampled at the left tube center (namely,
$v_{x}(-d/2,0;t)$); the right counterpart is either identical (for $S_{x}$) or
different only by sign (for $A_{x}$). Figure 3 focuses on the choice
$\ell/\mathcal{R}=0.4$, and plots the sampled $v_{x}$ for (a) the $S_{x}$ and
(b) the $A_{x}$ patterns by the solid curves. The same analysis as in Fig. 2a
is then repeated for both curves, whereby we single out the six extrema (the
red asterisks in Fig. 3) after the transitory phase to perform a fitting with
Eq. (32). The best-fit exponential envelope, plotted by the dashed curves for
the entire duration, is seen to well reflect the damping of kink motions for
much longer time. We nonetheless adopt a more conservative perspective to see
the internal flow field as evolving as an ideal quasi-mode in the restricted
segment, using the best-fit periods and damping times to evaluate the quasi-
mode frequencies and damping rates.

Figure 4 proceeds to examine the $\ell/\mathcal{R}$-dependencies of (a) the
oscillation frequencies ($\Omega$ in $\omega_{\rm Ae}$) and (b) damping-rate-
to-frequency ratios ($\gamma/\Omega$). The $S_{x}$ and $A_{x}$ patterns are
discriminated by the different colors. We present the best-fit values from our
time-dependent computations by the open circles, and overplot the T-matrix
results by the solid curves for comparison. Note that these T-matrix curves
are read from Fig. 4 of SL15. Our best-fit results are seen to compare
favorably with the T-matrix results, which is particularly true when
$\ell/\mathcal{R}\lesssim 0.3$. This statement holds despite the somewhat
visible difference in $\Omega$ at, say, $\ell/\mathcal{R}=0.2$, where the
best-fit value actually deviates from its T-matrix counterpart by only $\sim
2.7\%$. One further sees that the most significant departure occurs for the
$A_{x}$ damping rate when $\ell/\mathcal{R}=0.4$. However, this difference is
rather modest and reads $\sim 21.7\%$ in relative terms. We therefore conclude
that the dispersive properties of kink quasi-modes can be reasonably computed
by incorporating the shortcut TB formulae in the T-matrix framework for all
$\ell/\mathcal{R}$ examined here. Our time-dependent results, on the other
hand, further corroborate the SL15 conclusion that $A_{x}$ motions posses
higher frequencies and damp more rapidly than $S_{x}$ ones.

####  3.2.2 Subtleties with Domain Size and Grid Spacing

This subsection examines the influences that the numerical grid setup may have
on our time-dependent solutions. It suffices to present the representative
case where the $S_{x}$ pattern is examined for a specific
$\ell/\mathcal{R}=0.4$. We start with Fig. 5 by displaying the temporal
evolution of the $x$-speed at the left tube center ($v_{x}(-d/2,0;t)$, the
solid curves) for a number of domain sizes $[x_{\rm M},y_{\rm M}]$ as
discriminated by the different colors. The best-fit damping envelopes are
further plotted by the dashed curves for reference. All computations pertain
to a fixed grid spacing $\Delta=0.01b$, which is recalled to be different from
the reference value in Table 1. Some differences are seen among the solid
curves, the most prominent one being that the monochromatic cycles in the
$v_{x}$ signal are slightly delayed when the domain size increases. We find
that this domain size effect arises because the zero-gradient boundary
condition is not fully transparent to impinging motions initiated with the
pair of planar perturbations given by Eq. (25). Some weakly reflected
disturbances reach and interact with the composite structure, resulting in
some longer time for the internal flow fields to settle toward a regular
pattern for a larger domain. That the domain size effect is primarily
attributed to the form of the initial perturbation is further corroborated by
an additional series of computations where we multiply each exponential term
in Eq. (25) by a factor $\exp(-y^{2}/b^{2})$. The dependence on the domain
size is then found to disappear in the time-dependent solutions responding to
this localized perturbation. A monochromatic pattern remains clearly visible
after some transitory phase, the period being nearly identical to what we find
when Eq. (25) is employed. However, the damping envelope deviates considerably
from an exponential one, making it difficult to examine ideal quasi-modes from
an IVP perspective. The reasons for us to stick to Eq. (25) are twofold.
Firstly, what drives the internal field toward a quasi-mode behavior is beyond
our scope. Rather, we content ourselves with deriving the relevant periods and
damping times. Secondly, varying the domain size does not influence the best-
fit period ($P_{\rm fit}$) and only marginally impacts the best-fit damping
time ($\tau_{\rm fit}$). In fact, $P_{\rm fit}$ always evaluates to $\sim
112.5b/v_{\rm Ae}$ for the examined domains. The value of $\tau_{\rm
fit}v_{\rm Ae}/b$, on the other hand, changes from $320.1$ for $[x_{\rm
M}=12b,y_{\rm M}=10b]$ only modestly to $313.3$ for $[x_{\rm M}=18b,y_{\rm
M}=15b]$. An even weaker variation is found when the domain is further
extended to $[x_{\rm M}=24b,y_{\rm M}=20b]$, in which case $\tau_{\rm
fit}v_{\rm Ae}/b$ attains $314.3$.

One may ask why we chose a spacing $\Delta=0.01b$ different from the reference
value ($\Delta=0.005b$) when examining the domain size effect. The answer is
that the specific values of $\Delta$ do not impact the $v_{x}(-d/2,0;t)$
sequences for a given domain size if $\Delta$ is sufficiently small. By
“sufficiently” we mean some criterion related to the pair of equilibrium
quantities $[\rho_{\rm i}/\rho_{\rm e},\bar{l}]$ (see Sect. 2.4). We reach
this statement by varying only $\Delta$ for a substantial fraction of the
computations in Fig. 4. Consider only those with the reference domain size.
Figure 6 offers a representative case with $\ell/\mathcal{R}=0.4$ by
discriminating the $v_{x}(-d/2,0;t)$ sequences for different $\Delta$ with the
different linestyles. Relative to the computation with the reference $\Delta$
(the red dashed curve), no difference can be discerned when the spacing is
halved (blue dash-dotted) or doubled (black dotted). The reference $\Delta$ is
therefore overly demanding of computational resources for
$\ell/\mathcal{R}=0.4$, and the same can be said even for the most stringent
case where $\ell/\mathcal{R}=0.2$ provided that only $v_{x}(-d/2,0;t)$ is of
interest. That said, varying $\Delta$ necessarily entails differences in other
aspects of the flow fields. We illustrate this point by fixing
$\ell/\mathcal{R}$ again at $0.4$ and implementing the reference domain size.

Figure 7 presents the spatial profiles of the $x$-speed along the $y$-cut
through the left tube center at $t=340b/v_{\rm Ae}$. This instant is chosen
rather arbitrarily to be marginally inside the interval where our fitting
procedure is performed. A cut through a tube center, on the other hand, is
chosen to ease the description of some key dynamics. A number of grid spacings
are implemented and discriminated by the different colors. Only the lower half
($y\leq 0$) of a $y$-profile is displayed because the other half is symmetric
with respect to $y=0$. The vertical dash-dotted lines represent the borders of
the nonuniform layer. Focus for now on this layer, where one sees a series of
ripples in all computations. The most prominent effect associated with
$\Delta$ is that the ripples are better resolved when the numerical grid gets
finer, with the computed magnitudes of the local extrema being larger.
Physically speaking, however, all computations prove capable of capturing the
phase-mixing physics pertaining to the Alfv$\acute{\rm e}$n continuum. To
demonstrate this, we note that the Alfv$\acute{\rm e}$n frequency $\omega_{\rm
A}=kv_{\rm A}$ is locally symmetric about $x=-d/2$, meaning
$\partial\omega_{\rm A}/\partial x=0$ for the chosen cut. It then follows from
Mann et al. (1995) that the $v_{x}$ profile in the layer can be written as

$v_{x}\propto\mathcal{A}(y)\cos[\omega_{\rm A}(y)t+\varphi_{0}],$ (33)

where $\mathcal{A}(y)$ is some slow-varying envelope. Let $\omega_{\rm
Ai}=kv_{\rm Ai}$ be the internal Alfv$\acute{\rm e}$n frequency. Some
instantaneous phase variation $\phi_{\rm ac}(t)=(\omega_{\rm Ae}-\omega_{\rm
Ai})t$ is therefore accumulated over the continuum, leading to $N_{\rm
hc}=\lfloor\phi_{\rm ac}(t)/\pi\rfloor$ half-cycles with $\lfloor\cdot\rfloor$
being the floor function. Consequently, one expects $N_{\rm hc}$ or $N_{\rm
hc}+1$ extrema, where the uncertainty derives from the phase $\varphi_{0}$ in
Eq. (33). Plugging in the numbers yields a $\phi_{\rm ac}/\pi$ of $6.27$, and
therefore $6$ or $7$ extrema are expected. This expectation is faithfully
reproduced in Fig. 7. The same can actually be said for all instants when the
internal field behaves like a quasi-mode. Now move on to the portions outside
the nonuniform layer. The only difference that one can barely tell for
different $\Delta$ appears in the external flow fields immediately adjacent to
the layer. On the other hand, no difference can be discerned in the internal
motions.

Figure 8 examines the grid spacing effect from the energetics perspective,
differentiating $\Delta$ by the different colors. Equations (12) to (15) are
employed to evaluate all energetics-related quantities, and we take the area
$Q$ therein to be a box $[-(d/2+1.2b),(d/2+1.2b)]\times[-1.2b,1.2b]$. Figure
8a performs a gross energy balance analysis over $Q$ by presenting the
temporal profiles of the time-derivative of the instantaneous total energy
(${\rm d}E_{\rm tot}/{\rm d}t$, the solid curves) and those of the
instantaneous energy flux into $Q$ ($-F$, asterisks). The asterisks of
different colors cannot be told apart, the reason being that all computations
yield an identical set of dependents at the boundaries of $Q$. Some difference
nonetheless shows up among the solid curves, despite that ${\rm d}E_{\rm
tot}/{\rm d}t$ should balance $F$ (see Eq. 12). As intuitively expected, gross
energy conservation is maintained for a longer duration when a finer grid is
employed, with some $\sim 30\%$ difference between ${\rm d}E_{\rm tot}/{\rm
d}t$ and $-F$ staring to appear at $t\sim 90b/v_{\rm Ae}$ ($\sim 270b/v_{\rm
Ae}$) when $\Delta=0.01b$ ($0.0025b$). Note that the examined $Q$ is the union
of three mutually exclusive portions, namely the tube interiors (to be denoted
as “int”), nonuniform layers (“layer”), and exterior (“ext”). Figure 8b then
displays, against time, the energies in these individual portions. By far the
most significant effect that the grid spacing has on the energetics is that
the energy in the nonuniform layer is increasingly under-estimated when
$\Delta$ increases. This is understandable because small-scale variations in
the layer are less well resolved if our numerical grid goes coarser (see Fig.
7). However, there exists only some weak difference among the $E_{\rm ext}$
curves because only the flows in the immediate neighborhood of either tube are
somehow impacted by choices of $\Delta$. Furthermore, the internal energies
($E_{\rm int}$) consistently remain the same when $\Delta$ is varied,
reinforcing the notion that $\Delta$ has no influence on the internal flow
field.

##  4 Kink Motions in a Two-Elliptic-Tube Setup

Let us start by summarizing several numerical aspects of our computations.
Overall, we have verified that the MacCormack scheme is appropriate for our
purposes, with no issue arising from the application to nonplanar structures
of the finite-difference methodology on a Cartesian grid. The zero-gradient
boundary condition is somehow not fully transparent to perturbations excited
by our planar drivers, and the domain size somehow impacts how the internal
flows transition to a quasi-mode behavior. However, this domain size effect is
practically negligible on the best-fit periods ($P_{\rm fit}$) and damping
times ($\tau_{\rm fit}$) that we derive for the quasi-mode stage. Likewise,
the computed internal flows are not affected by the grid spacing $\Delta$
provided that $\Delta$ is sufficiently small.

This section is devoted to kink motions in a two-elliptic-tube configuration,
for which the $x$-major and $x$-minor orientations need to be discriminated.
The reference grid ($[x_{\rm M}=30b,y_{\rm M}=25b,\Delta=0.01b]$, see Table 1)
will be consistently adopted, enabling a meaningful comparison between the
tube orientations. This reference grid sufficies for our purposes. The
notations $[a,b,\bar{l}]$ are adopted throughout this section, where all
computations pertain to a fixed combination of physical parameters $[\rho_{\rm
i}/\rho_{\rm e}=3,\bar{l}=0.4,kb=\pi/30]$. We adjust only the dimensionless
tube separation ($d/b$) and the ratio of the semi-major to semi-minor axis
($a/b$) for a given orientation and a given perturbation pattern (see Eq. 26).
Tube overlapping is always avoided when the values of $a/b$ and $d/b$ are
experimented with. The equilibria are consistently perturbed with Eq. (25).

###  4.1 $S_{x}$ and $A_{x}$ Motions for the $x$-major Orientation

This subsection addresses the $x$-major orientation. We start with Fig. 9 to
present the temporal evolution of the $x$-speed sampled at the left tube
center ($v_{x}(-d/2,0;t)$, the solid lines) for both (a) the $S_{x}$ and (b)
the $A_{x}$ motions. The tube separation is fixed at $d=5b$, whereas a number
of values are examined for $a/b$ as discriminated by the different colors.
Only the variations after the rapid-varying phase are of interest for any
$v_{x}$ sequence. We somehow emphasize the segment exactly encompassing the
six extrema starting with the one labeled by the relevant asterisk, performing
a fitting procedure with Eq. (32). A best-fit exponential envelope results,
with an example plotted over the entire duration for the case where $a/b=2.5$.
With this example we demonstrate that an exponentially damped cosine is in
general adequate for describing the $v_{x}$ sequences after the rapid-varying
phase. The best-fit periods and damping times thus derived will be be
understood as pertaining to an ideal quasi-mode.

The temporal attenuation of the internal flow fields means energy
redistribution in ideal MHD. Let us take the case with $a/b=2.5$ in Fig. 9,
and note that the two tubes are actually in contact ($d=5b=2a$). Figure 10
presents the spatial distribution of the velocity fields
($\@vec{v}=v_{x}\@vec{e}_{x}+v_{y}\@vec{e}_{y}$, the blue arrows) and the
instantaneous energy density ($\epsilon$, filled contours) for (a) the $S_{x}$
and (b) $A_{x}$ motions at an arbitrarily chosen instant $t=260b/v_{\rm Ae}$.
A subarea of the second quadrant is singled out by each inset to emphasize the
$\epsilon$ distributions in the nonuniform layers of both tubes, with the red
curves delineating the layer boundaries. The white curve further indicates
where the local Alfv$\acute{\rm e}$n frequency equals the quasi-mode
frequency. Evidently, the white and red curves are contours of the equilibrium
density, and therefore part of the relevant ellipses (see Eq. 18). Figure 10
is actually taken from the attached animation. The arrows and filled contours
are plotted in a way that it makes sense to compare the $\@vec{v}$ or
$\epsilon$ strengths not only in one snapshot but between different instants
for a given orientation. In particular, the darker a portion is, the larger
the value of $\epsilon$ therein.

Consider first the ambient fields for both perturbation patterns, where by
“ambient” we refer to the flows some distance away from the tubes. Likewise,
by “tube flow” we mean the field excluding the ambient flow, and by “internal
flow” we refer specifically to the uniform portion inside either tube. For
both patterns, one sees from the animation that the ambient flows tend to vary
more rapidly than the internal ones. This feature is particularly clear for
the $S_{x}$ pattern because the ambient flows tend to be substantially
stronger than the internal flows as well. Regardless, a periodogram analysis
yields that the ambient flows follow primarily the external Alfv$\acute{\rm
e}$n frequency $\omega_{\rm Ae}$, whereas the variations of the internal flows
are by far dominated by the lower quasi-mode frequency. That the ambient
motions are not coordinated with the internal ones is a result of the form of
the initial exciter given by Eq. (25). The relevant physics is actually very
similar to what happens in the T08 study despite the configurational
differences, the key being the $y$-invariance of Eq. (25). It proves easier to
explain this by considering a uniform equilibrium with density $\rho_{\rm e}$,
for which Eqs. (7) to (11) can be combined to yield a Klein-Gordon equation if
the $y$-dependence is dropped (see Eq. 1 in Terradas et al., 2005). Suppose
that an individual component in Eq. (25) is applied. A dispersive fast wave
then results, whereby any fluid parcel in the system eventually end up in some
oscillatory wake at the Alfv$\acute{\rm e}$n frequency $\omega_{\rm Ae}$ (see
Fig. 1 in Terradas et al., 2005). Oscillatory wakes turn out to still exist in
the ambient flow when the two elliptic tubes are introduced, and when an
additional planar perturbation is implemented. Note that the initial
perturbations are external to our tubes. Note further that overall the large-
time behavior for the ambient flow actually comprises two component wakes,
each being the response to the corresponding component driver. The two
component wakes tend to interfere constructively (largely destructively) for
the $S_{x}$ ($A_{x}$) pattern, thereby explaining the strength of the ambient
flow relative to the internal one.

Now examine the tube flows. Focusing on the blue arrows, one sees from Fig. 10
that the gross patterns for both $S_{x}$ and $A_{x}$ are similar to Fig. 2 in
L08 despite that a two-circular-tube equilibrium with piece-wise constant
profile was addressed therein. Specifically, the internal velocity field is
more or less uniform, and a pair of vortical motions develop at the edge of
either tube. These two features, combined with the overall symmetric
properties of the internal field about $x=0$, demonstrate the robustness of
the classification scheme for lower-order kink motions. The primary difference
from L08, on the other hand, is that the vortical motions gradually evolve
from a simple dipolar behavior (see the first several instants in the
animation) into multiple shearing layers (Fig. 10). This evolution is well
known for kink motions in isolated tubes with either circular (e.g., Terradas
et al., 2008a; Pascoe et al., 2010; Antolin et al., 2014) or elliptic cross-
sections (e.g., Ruderman, 2003; Guo et al., 2020). What Fig. 10 demonstrates
is then that tube interactions do not compromise the phase-mixing physics
behind this evolution. That said, tube interactions do impact the detailed
development of the small-scale shears, as evidenced by the differences in the
velocity fields for the two different perturbation patterns. Suppose that tube
interactions are negligible. The tube flows are then expected to be identical
for the $S_{x}$ and $A_{x}$ patterns, meaning in particular a symmetric
distribution of the energy density ($\epsilon$) with respect to the minor axis
($x=-d/2=-2.5b$ here). However, this expectation roughly holds only for the
$S_{x}$ pattern (Fig. 10a inset), whereas the left portion in the left tube is
favored in the $\epsilon$ distribution for the $A_{x}$ pattern (Fig. 10b
inset). Besides phase-mixing, a closely related aspect that is not
fundamentally compromised by tube interactions is the resonant interplay
between the quasi-mode and the Alfv$\acute{\rm e}$n continuum (see Soler &
Terradas 2015 and references therein for the subtle distinction between
resonant absorption and phase-mixing). By this we refer to two features common
to the $S_{x}$ and $A_{x}$ patterns, and we concentrate on the insets in the
animation. Firstly, the attenuation of the internal field is accompanied by
the accumulation of perturbation energies in the nonuniform layer. Secondly,
this energy accumulation leads to a localized $\epsilon$ distribution, with
the strongest perturbations tending to the white resonance contour as time
proceeds for roughly one quasi-mode period after the rapid-varying phase.
Overall with Fig. 10 we conclude that the notions of phase-mixing and resonant
absorption remain applicable to kink motions in our two-elliptic-tube
configuration.

Figure 11 further examines the phase-mixing process by showing the $y$-cuts of
the $x$-speeds through the left tube center for both (a) the $S_{x}$ and (b)
the $A_{x}$ patterns. A number of instants are rather arbitrarily chosen and
discriminated by the different colors. The vertical dotted lines mark the
borders of the nonuniform layer, where shearing motions ($\partial
v_{x}/\partial y$ here) with increasingly small scales are seen to develop as
time proceeds. We choose to quantify this development by relating the number
of the $v_{x}$ extrema ($N_{\rm extrm}(t)$) to the instantaneous phase
variation accumulated over the Alfv$\acute{\rm e}$n continuum ($\phi_{\rm
ac}(t)$, see the discussion on Eq. 33). Note that $\phi_{\rm
ac}(t)=(\omega_{\rm Ae}-\omega_{\rm Ai})t$ is identical for both perturbation
patterns, corresponding specifically to $[2.82,4.23,5.64]$ half-cycles for the
examined instants. One then expects $[2^{+1}_{+0},4^{+1}_{+0},5^{+1}_{+0}]$
extrema in the $v_{x}$ profiles in Figs. 11a and 11b, which is indeed the
case. As can be readily verified, that $N_{\rm extrm}(t)=\lfloor[\phi_{\rm
ac}(t)/\pi\rfloor]$ at these instants actually constrains the phase
$\varphi_{0}$ in Eq. (33) to a very narrow range $[q\pi,(q+0.19)\pi]$ with
$q=0,1$.

Figure 12 proceeds to examine (a) the periods $P$ and (b) damping-time-to-
period ratios $\tau/P$ of kink quasi-modes pertaining to both the $S_{x}$ (the
solid curves) and the $A_{x}$ (dashed) patterns, showing how $P$ and $\tau$
depend on the ratio of the semi-major to semi-minor radius ($a/b$) for a given
tube separation ($d$). We also examine a number of $d$ as discriminated by the
different colors. Note that this survey adopts a fixed combination of physical
parameters $[\rho_{\rm i}/\rho_{\rm e}=3,\bar{l}=0.4,kb=\pi/30]$. The same set
of parameters is additionally adopted to evaluate $P$ and $\tau/P$ of the kink
quasi-mode for an isolated circular tube with the resistive eigenmode approach
(see Sect. 3.1), the results being plotted by the horizontal lines for
comparison. Note that the values for $P$ (in units of the axial
Alfv$\acute{\rm e}$n time $T_{\rm axial}=2\pi/\omega_{\rm Ae}=2L/v_{\rm Ae}$)
and $\tau/P$ vary little if $kb=\pi b/L$ is further reduced, meaning that
little will change if one examines axial fundamentals in AR loops with much
larger $L/b$ than adopted here. Note further that $a/b$ is required not to
exceed $d/2b$ for a given $d/b$ such that tube overlapping is avoided. For any
given pair $[a/b,d/b]$, Fig. 12 indicates that the $A_{x}$ motion possesses a
shorter period and damps more efficiently than the $S_{x}$ motion. This
behavior is identical to what SL15 found for two-circular-tube configurations.
The argument therein is that $A_{x}$ motions are more “forced”; the two tubes
move largely in a synchronous fashion for the $S_{x}$ pattern, whereas the
flows in between the two tubes periodically “collide” when the $A_{x}$ pattern
is examined. The same argument applies here although elliptic tubes are
addressed, with the insets of Fig. 10b already hinting at a more forced
behavior for $A_{x}$ motions. Also similar to SL15 is that, for a given $a/b$,
the difference between the values of $P$ (or $\tau/P$) for the $S_{x}$ and
$A_{x}$ motions tends to decrease monotonically with $d$. This agrees with the
intuitive expectation for a weaker tube interaction.

Let us pay more attention to the overall behavior for $P$ to increase
monotonically with $a/b$ for a given $d/b$. Somehow subtle is that tube
interactions also have some effect, because increasing $a/b$ actually brings
the two tubes effectively closer even if the distance between the tube centers
is fixed. Regardless, this tube interaction only plays a minor role, given
that the monotonical $a/b$-dependence takes place for different $d$ and for
the $S_{x}$ and $A_{x}$ motions alike. One may therefore expect the same
$a/b$-dependence for isolated elliptic tubes (or equivalently $d/b\to\infty$),
in which case it is no longer necessary to distinguish between the $S_{x}$ and
$A_{x}$ patterns. This was indeed seen in our previous IVP study (Guo et al.,
2020), where we offered some heuristic argument similar to the earlier one by
Ruderman (2003). For simplicity, suppose that the equilibrium configuration is
transversely structured in a piece-wise constant manner. The heuristic
argument then relies on the observation that the velocity normal to the tube
edge plays a central role for an isolated density-enhanced tube to communicate
with its surroundings. Evidently, a larger $a/b$ makes the tube edge more
elongated in the $x$-direction. Their velocities primarily $x$-directed, the
fluid parcels in the tube interior therefore become less aware of the
surroundings, meaning some enhanced effective inertia and hence a longer
period. The same argument can be invoked in this study, despite the subtlety
that the fluids surrounding one tube actually embed another tube.

Figure 12 allows us to say a few words on the seismological implications of
this subsection. Let us focus on axial fundamentals, and consider only the
$S_{x}$ motion for $[a/b=2.5,d/b=5]$ to ease our description. By
“seismological implications” we refer to the inference of the axial
Alfv$\acute{\rm e}$n time ($T_{\rm axial}=2L/v_{\rm Ae}$) and the
dimensionless nonuniform layer width ($\bar{l}$). Now assume that the relevant
oscillating tubes possess the “true” dimensionless values $[\rho_{\rm
i}/\rho_{\rm e}=3,\bar{l}=0.4]$ as implemented in our computations (note that
$L/b$ is immaterial when sufficiently large). Let $P_{\rm obs}$ and
$(\tau/P)_{\rm obs}$ be the measured values of the period and damping-time-to-
period ratio, and we take $(\tau/P)_{\rm obs}=6.39$ to comply with Fig. 12.
One therefore deduces a “true” value of $P_{\rm obs}/1.6$ for $T_{\rm axial}$,
which is recalled to be dimensional. Now suppose that neighboring elliptic
tubes with the $x$-major orientation are imaged in, say, the EUV when the line
of sight lies in the tube plane (see Fig. 1a), meaning that the tubes may be
readily mistaken as circular ones with a tube radius of $a$ (rather than $b$).
Suppose further that the density contrast is measured to be $\rho_{\rm
i}/\rho_{\rm e}=3$. Now let us follow the customary practice to infer $T_{\rm
axial}$ and $\bar{l}$ by seeing the $S_{x}$ motion as a kink quasi-mode in an
isolated circular tube. It matters little regarding whether one assigns $a$ or
$b$ to the tube radius. The dimensionless quasi-mode period $P/T_{\rm axial}$
depends essentially only on $\rho_{\rm i}/\rho_{\rm e}$, always attaining
$\sim 1.42$ when $\bar{l}\lesssim 0.6$ (see the green horizontal line in Fig.
12a; see also Fig. 2b). One therefore deduces a value of $P_{\rm obs}/1.42$
for $T_{\rm axial}$, which differs from the true value by $\sim 13\%$.
However, it is considerably more problematic when $(\tau/P)_{\rm obs}$ is
invoked to infer $\bar{l}$. Employing the resistive eigenmode approach, we
find that a $(\tau/P)_{\rm obs}=6.39$ results in a $\bar{l}\approx 0.18$,
which underestimates the true value by $\sim 55\%$. One may argue that this
difference is not that significant; after all the customary seismological
practice yields the correct qualitative picture that the oscillating tube
possesses a thin boundary. Our point, however, is that care needs to be
exercised when $(\tau/P)_{\rm obs}$ is put to quantitative use, with the
uncertainties in the deduced dimensionless layer width $\bar{l}$ readily
exceeding $\sim 50\%$ if one neglects the combined effect of tube cross-
sectional shapes and tube interactions.

###  4.2 $S_{x}$ and $A_{x}$ Motions for the $x$-minor Orientation

This subsection addresses the $x$-minor orientation. We start with Fig. 13
where the same set of physical quantities as in Fig. 9 is employed, and the
time sequences for the $x$-speed sampled at the left tube center are presented
in an identical format. A fitting procedure is once again performed for any
$v_{x}$ sequence over the duration exactly enclosing the six extrema starting
with the one indicated by the pertinent asterisk. Despite this, the best-fit
exponential envelopes are seen to describe the large-time behavior equally
well. We therefore take the view that the best-fit periods ($P$) and damping
times ($\tau$) pertain to the ideal quasi-modes supported by our two-elliptic-
tube configuration.

Figure 14 specializes to the case with $[a/b=2.5,d/b=5]$, following the same
format as in Fig. 10 to present the velocity fields (the blue arrows) and the
spatial distributions of the perturbation energy density ($\epsilon$, filled
contours) at a representative instant. Note that the tubes are quite some
distance apart in this case. Note also that Fig. 14 is a snapshot extracted
from the attached animation. Focusing on the velocity fields, one can safely
conclude from the animation that the notations of $S_{x}$ and $A_{x}$ make
physical sense as proposed by L08 for two-circular-tube systems. For both
perturbation patterns, the insets in the animation further indicate that the
attenuation of the internal flows is accompanied by the enhancement of
perturbation energy density $\epsilon$ in the nonuniform layers, with the
localization of $\epsilon$ strongly mediated by the Alfv$\acute{\rm e}$n
resonance. This latter point can be readily drawn from the close association
of the darkest portion in either $\epsilon$ distribution with the white curve
marking where the quasi-mode frequency matches the local Alfv$\acute{\rm e}$n
frequency.

Figure 15 presents, in a format identical to Fig. 11, the $y$-cuts of the
$x$-speed through the left tube center for a number of arbitrarily chosen
instants. The aim is also to show the key role that phase-mixing plays for the
gradual development of shearing motions with increasingly small scales in the
nonuniform layers. We follow Fig. 11 to quantify this with the aid of Eq.
(33), once again counting the instantaneous number $N_{\rm extrm}(t)$ of the
$v_{x}$ extrema. The pertinent instantaneous phase variations ($\phi_{\rm
ac}(t)$), on the other hand, are the same as in Fig. 11 and read
$[2.82,4.23,5.64]\pi$ for the examined instants. However, somehow different is
that the relevant values for $N_{\rm extrm}$ now may deviate from
$\lfloor\phi_{\rm ac}/\pi\rfloor$, attaining specifically $[3,4,5]$. This
deviation is readily attributable to the phase $\varphi_{0}$ in Eq. (33).
Conversely, this set of deviations can be readily verified to constrain
$\varphi_{0}$ to a rather narrow range of $[(q+0.19)\pi,(q+0.37)\pi]$ with
$q=0,1$.

Figure 16 moves on to survey (a) the quasi-mode periods $P$ and (b) damping-
time-to-period ratios $\tau/P$ for a substantial range of combinations
$[a/b,d/b]$. This figure is identical in format to Fig. 12. The values
examined for the tube separation $d$ consistently avoid tube overlapping, and
hence there is no need to constrain the ratio of the semi-major to semi-minor
axis ($a/b$) for a given $d/b$. Examine Fig. 16a first, from which one sees
the expected monotonic behavior for the differences between the $S_{x}$ and
$A_{x}$ periods to decrease with $d/b$ for a given $a/b$. For a given
separation $d/b$, one further sees that the values of $P$ tend to decrease
monotonically with $a/b$ for the $S_{x}$ and $A_{x}$ motions alike. However,
tube interactions remain possible to be partly responsible for this
$a/b$-dependence, because the area between the two tubes actually broadens
with $a/b$ even if $d/b$ is fixed. This tube-interaction effect can hardly be
dismissed given that the $S_{x}$ and $A_{x}$ values deviate more strongly when
$a/b$ increases for a given $d/b$. Having said that, it remains possible to
partly explain the monotonic $a/b$-dependence of $P$ by invoking the inertia
arguments initially offered for isolated density-enhanced elliptic tubes. Note
that the tube edges are now more elongated in the $y$-direction, whereas the
fluids parcels in either tube interior remain to oscillate primarily in the
$x$-direction. It therefore holds that an increase in $a/b$ makes the internal
fluid parcels better aware of the external medium, meaning some reduced
effective inertia and hence a shorter period.

Consider now Fig. 16b where the damping-time-to-period ratios ($\tau/P$) are
examined. Two features arise when $a/b\lesssim 1.5$, one being that $A_{x}$
motions attenuate more efficiently than $S_{x}$ ones, the other being that the
differences between the $S_{x}$ and $A_{x}$ values tend to decrease with tube
separation. These features resemble what happens for two-circular-tube
systems, and are expected because the cross-sectional shapes are not far from
circular. However, the behavior of $\tau/P$ at larger $a/b$ is considerably
more complicated than in Fig. 12b where the $x$-major orientation is
addressed. By “complicated” we specifically refer to the nonmonotoic
$d/b$-dependence of $\tau/P$ for either perturbation pattern. Take the case
with $a/b=2.5$. When $d/b$ increases, one sees an increase followed by some
decrease in $\tau/P$ for the $S_{x}$ motion, whereas $\tau/P$ tends to
decrease for the $A_{x}$ motion. This latter simplicity is only apparent,
because $\tau/P$ is bound to increase with $d/b$ again such that the $A_{x}$
values eventually join their $S_{x}$ counterparts when the tubes are
sufficiently separated. We have verified that the somehow involved behavior of
$\tau/P$ is not of numerical origin (see Sect. 3). Furthermore, it matters
little as to how to choose the segment of a $v_{x}$ sequence for fitting,
provided that this segment starts sufficiently later than the rapid-varying
phase (see Fig. 13). The most likely reason is then that the perturbations in
between the two tubes depend on $a/b$ and $d/b$ in an intricate fashion, and
hence an intricate dependence on $a/b$ and $d/b$ of the efficiency of tube
interactions. Regardless, the reason behind the $d/b$-dependence of $\tau/P$
is also behind, say, the nonmonotonic $a/b$-dependence of $\tau/P$ for $A_{x}$
motions for a given $d/b$ (see the dashed curves), and the occasional behavior
for $S_{x}$ motions to damp more efficiently (compare the black solid and
dashed curves).

Some remarks can be offered regarding the seismological implications of Fig.
16. Only axial fundamentals are of interest, and the tube plane is once again
assumed to be perpendicular to the plane-of-sky. We assess what uncertainties
the customary practice may introduce to the axial Alfv$\acute{\rm e}$n time
$T_{\rm axial}=2L/v_{\rm Ae}$ and dimensionless layer width $\bar{l}$ by
neglecting the joint effects of cross-sectional shapes and tube interactions.
The procedure is identical to what we performed for the $x$-major orientation;
we accordingly choose the $A_{x}$ values for $[a/b=2.5,d/b=3]$ given that they
deviate the most from the horizontal lines. Consider the period first. One
infers a true value of $T_{\rm axial}=P_{\rm obs}/1.13$ with Fig. 16a, whereas
the customary practice remains to yield a $T_{\rm axial}\approx P_{\rm
obs}/1.42$ or equivalently some underestimation by $\sim 20\%$. This
uncertainty in $T_{\rm axial}$ is larger than for the $x$-major orientation.
Now taking $(\tau/P)_{\rm obs}=3.44$ from Fig. 16b as measured, one deduces
with the customary practice a dimensionless layer width of $\bar{l}=0.25$ or
some underestimation of the true value ($\bar{l}=0.4$) by $\sim 38\%$. Somehow
smaller than for the $x$-major case, this uncertainty remains substantial
enough to corroborate our claim that the inhomogeneity scales deduced with the
customary practice need to be treated with caution.

##  5 Summary

With coordinated transverse displacements in neighboring active region (AR)
loops in mind, we have addressed damped kink motions in straight equilibria
where two identical parallel density-enhanced tubes with elliptic cross-
sections (elliptic tubes) are embedded in an ambient corona. Linear, ideal,
pressureless MHD was adopted throughout. We concentrated on axially standing
motions, formulating the ensuing two-dimensional (2D) initial value problem
(IVP) in the plane transverse to the equilibrium magnetic field. We identified
the direction connecting the tube centers as horizontal, discriminating two
tube orientations where either the major (dubbed “$x$-major”) or the minor
(“$x$-minor”) axis is horizontally placed. The system evolution was initiated
with external velocity drivers, implemented in such a way that the internal
flow fields are primarily horizontal and are either symmetric ($S_{x}$) or
anti-symmetric ($A_{x}$) with respect to the vertical axis about which our
equilibrium configuration is symmetric. The temporal evolution of the velocity
perturbations at tube centers was of particular interest, allowing us to
identify the quasi-mode stage as where the monochromatic time sequence follows
an exponentially damped envelope. We paid special attention to how the quasi-
mode periods and damping times depend on the tube separation and the cross-
sectional aspect ratio, thereby addressing the impact on the dispersive
properties of damped kink motions from the joint effects of tube interactions
and cross-sectional shapes.

Our numerical findings are summarized as follows. When two-circular-tube
equilibria are examined as a special case, the quasi-mode periods and damping
times found with our IVP approach are consistent with the modal analysis by
Soler & Luna (2015), thereby independently justifying both the analytical
T-matrix approach and the numerical results therein. We further find that the
notions of resonant absorption and phase-mixing are not undermined when
elliptic cross-sections are allowed for. Specifically, the nonuniform layers
around tube boundaries feature both the development of velocity shears with
increasingly small scales and the energy accumulation around where the local
Alfv$\acute{\rm e}$n frequency matches the quasi-mode frequency. The $A_{x}$
motions are found to possess shorter periods and damp more rapidly than the
$S_{x}$ motions for the $x$-major orientation, in which case the periods and
damping times for both perturbation patterns tend to increase with the major-
to-minor-axis ratio. When the $x$-minor orientation is addressed, the $A_{x}$
motions remain to oscillate more rapidly than the $S_{x}$ ones. However, the
periods for both patterns tend to decrease monotonically with the major-to-
minor-axis ratio. Somehow subtle are the joint effects of tube interactions
and cross-sectional shapes on the damping times, which may depend non-
monotonically on the tube separation or the major-to-minor-axis ratio.
Furthermore, the $A_{x}$ motions may occasionally damp less efficiently.
Neglecting these joint effects may introduce some uncertainty of $\sim 20\%$
($\sim 50\%$) to the axial Alfv$\acute{\rm e}$n time (the inhomogeneity
lengthscale) deduced with the period (damping time). Consequently, care may
need to be exercised when the damping times of kink motions are put to
standard seismological practice built on wave theories for isolated circular
tubes.

Some further remarks seem necessary before closing. Firstly, we note that our
IVP approach is not restricted to two-tube systems but capable of handling
rather generic 2D inhomogeneities. In particular, axially standing motions may
be examined for straight configurations comprising multiple tubular structures
with arbitrary cross-sections. Collective motions in multiple neighboring
resolved loops will then be better understood, together with the motions in
those apparently isolated loops that actually involve a multitude of
unresolved strands (see e.g., Sect.4.2 of the review by Reale 2014 for more on
this multi-stranded nature of AR loops). Secondly, that kink motions are
resonantly absorbed for configurations with continuous transverse structuring
is actually much-expected (e.g., the review by Goossens et al., 2011).
However, somehow largely unnoticed is that the decayless regime was also
observed for coordinated kink motions in adjacent AR loops (Wang et al.,
2012). Accept that the transverse structuring is likely to be continuous in
the relevant configurations. It then becomes necessary to address what
counteracts the ensuing resonant damping, thereby bringing the understanding
of decayless kink motions in multi-tube systems closer to the level of
sophistication reached for the decayless regime in isolated tubes (e.g.,
Nakariakov et al., 2016; Guo et al., 2019b; Ruderman et al., 2021). Thirdly,
coordinated motions in multiple tubular structures were identified not only in
AR loops but also in, say, groups of prominence threads (see e.g., Arregui et
al., 2018, for a review). To our knowledge, so far the only theoretical
investigation in the prominence context was the one by Soler et al. (2009)
with the T-matrix formalism. An IVP perspective is expected to be as fruitful,
provided that one further account for, say, the effects associated with
partial ionization (e.g., the review by Ballester et al., 2018).

###### Acknowledgements.

This research was supported by the National Natural Science Foundation of
China (41974200, 12273019, 12373055, 12203030, and 42230203). We gratefully
acknowledge ISSI-BJ for supporting the international team “Magnetohydrodynamic
wavetrains as a tool for probing the solar corona”, and ISSI-Bern for
supporting the international team “Magnetohydrodynamic Surface Waves at
Earth’s Magnetosphere and Beyond”.

## References

  * Aldhafeeri et al. (2021) Aldhafeeri, A. A., Verth, G., Brevis, W., et al. 2021, ApJ, 912, 50, doi: 10.3847/1538-4357/abec7a
  * Anderson (1995) Anderson, J. D. 1995, Computational Fluid Dynamics: the Basics with Applications (New York: McGraw-Hill, Inc.)
  * Anfinogentov et al. (2013) Anfinogentov, S., Nisticò, G., & Nakariakov, V. M. 2013, A&A, 560, A107, doi: 10.1051/0004-6361/201322094
  * Anfinogentov & Nakariakov (2019) Anfinogentov, S. A., & Nakariakov, V. M. 2019, ApJ, 884, L40, doi: 10.3847/2041-8213/ab4792
  * Anfinogentov et al. (2015) Anfinogentov, S. A., Nakariakov, V. M., & Nisticò, G. 2015, A&A, 583, A136, doi: 10.1051/0004-6361/201526195
  * Antolin et al. (2014) Antolin, P., Yokoyama, T., & Van Doorsselaere, T. 2014, ApJ, 787, L22, doi: 10.1088/2041-8205/787/2/L22
  * Arregui (2022) Arregui, I. 2022, Frontiers in Astronomy and Space Sciences, 9, 826947, doi: 10.3389/fspas.2022.826947
  * Arregui & Asensio Ramos (2011) Arregui, I., & Asensio Ramos, A. 2011, ApJ, 740, 44, doi: 10.1088/0004-637X/740/1/44
  * Arregui et al. (2018) Arregui, I., Oliver, R., & Ballester, J. L. 2018, Living Reviews in Solar Physics, 15, 3, doi: 10.1007/s41116-018-0012-6
  * Arregui et al. (2015) Arregui, I., Soler, R., & Asensio Ramos, A. 2015, ApJ, 811, 104, doi: 10.1088/0004-637X/811/2/104
  * Aschwanden et al. (1999) Aschwanden, M. J., Fletcher, L., Schrijver, C. J., & Alexander, D. 1999, ApJ, 520, 880, doi: 10.1086/307502
  * Aschwanden et al. (2004) Aschwanden, M. J., Nakariakov, V. M., & Melnikov, V. F. 2004, ApJ, 600, 458, doi: 10.1086/379789
  * Aschwanden et al. (2003) Aschwanden, M. J., Nightingale, R. W., Andries, J., Goossens, M., & Van Doorsselaere, T. 2003, ApJ, 598, 1375, doi: 10.1086/379104
  * Aschwanden & Schrijver (2011) Aschwanden, M. J., & Schrijver, C. J. 2011, ApJ, 736, 102, doi: 10.1088/0004-637X/736/2/102
  * Ballester et al. (2018) Ballester, J. L., Alexeev, I., Collados, M., et al. 2018, Space Sci. Rev., 214, 58, doi: 10.1007/s11214-018-0485-6
  * Bogdan & Zweibel (1985) Bogdan, T. J., & Zweibel, E. G. 1985, ApJ, 298, 867, doi: 10.1086/163670
  * Bogdan & Zweibel (1987) —. 1987, ApJ, 312, 444, doi: 10.1086/164889
  * Cally (1986) Cally, P. S. 1986, Sol. Phys., 103, 277, doi: 10.1007/BF00147830
  * Chen et al. (2021) Chen, S.-X., Li, B., Van Doorsselaere, T., et al. 2021, ApJ, 908, 230, doi: 10.3847/1538-4357/abd7f3
  * De Moortel & Nakariakov (2012) De Moortel, I., & Nakariakov, V. M. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 3193, doi: 10.1098/rsta.2011.0640
  * Edwin & Roberts (1983) Edwin, P. M., & Roberts, B. 1983, Sol. Phys., 88, 179, doi: 10.1007/BF00196186
  * Erdélyi & Morton (2009) Erdélyi, R., & Morton, R. J. 2009, A&A, 494, 295, doi: 10.1051/0004-6361:200810318
  * Erdélyi & Taroyan (2008) Erdélyi, R., & Taroyan, Y. 2008, A&A, 489, L49, doi: 10.1051/0004-6361:200810263
  * Gijsen & Van Doorsselaere (2014) Gijsen, S. E., & Van Doorsselaere, T. 2014, A&A, 562, A38, doi: 10.1051/0004-6361/201322755
  * Goddard et al. (2016) Goddard, C. R., Nisticò, G., Nakariakov, V. M., & Zimovets, I. V. 2016, A&A, 585, A137, doi: 10.1051/0004-6361/201527341
  * Goedbloed et al. (2019) Goedbloed, H., Keppens, R., & Poedts, S. 2019, Magnetohydrodynamics of Laboratory and Astrophysical Plasmas (Cambridge University Press), doi: 10.1017/9781316403679
  * Goossens et al. (2002) Goossens, M., Andries, J., & Aschwanden, M. J. 2002, A&A, 394, L39, doi: 10.1051/0004-6361:20021378
  * Goossens et al. (2008) Goossens, M., Arregui, I., Ballester, J. L., & Wang, T. J. 2008, A&A, 484, 851, doi: 10.1051/0004-6361:200809728
  * Goossens et al. (2011) Goossens, M., Erdélyi, R., & Ruderman, M. S. 2011, Space Sci. Rev., 158, 289, doi: 10.1007/s11214-010-9702-7
  * Goossens et al. (1992) Goossens, M., Hollweg, J. V., & Sakurai, T. 1992, Sol. Phys., 138, 233, doi: 10.1007/BF00151914
  * Guo et al. (2020) Guo, M., Li, B., & Van Doorsselaere, T. 2020, ApJ, 904, 116, doi: 10.3847/1538-4357/abc1df
  * Guo et al. (2019a) Guo, M., Van Doorsselaere, T., Karampelas, K., & Li, B. 2019a, ApJ, 883, 20, doi: 10.3847/1538-4357/ab338e
  * Guo et al. (2019b) Guo, M., Van Doorsselaere, T., Karampelas, K., et al. 2019b, ApJ, 870, 55, doi: 10.3847/1538-4357/aaf1d0
  * Hong (1996) Hong, H. 1996, Applied Mathematics Letters, 9, 99
  * Jardin (2010) Jardin, S. 2010, Computational Methods in Plasma Physics, Chapman & Hall/CRC Computational Science (CRC Press).  https://books.google.nl/books?id=gZzf_B56FDcC
  * Keppens et al. (1994) Keppens, R., Bogdan, T. J., & Goossens, M. 1994, ApJ, 436, 372, doi: 10.1086/174912
  * Klimchuk & DeForest (2020) Klimchuk, J. A., & DeForest, C. E. 2020, ApJ, 900, 167, doi: 10.3847/1538-4357/abab09
  * Klimchuk et al. (1992) Klimchuk, J. A., Lemen, J. R., Feldman, U., Tsuneta, S., & Uchida, Y. 1992, PASJ, 44, L181
  * Kucera et al. (2019) Kucera, T. A., Young, P. R., Klimchuk, J. A., & DeForest, C. E. 2019, ApJ, 885, 7, doi: 10.3847/1538-4357/ab449e
  * Li et al. (2020) Li, B., Antolin, P., Guo, M. Z., et al. 2020, Space Sci. Rev., 216, 136, doi: 10.1007/s11214-020-00761-z
  * Li et al. (2022) Li, B., Chen, S.-X., & Li, A.-L. 2022, ApJ, 928, 33, doi: 10.3847/1538-4357/ac5402
  * Li & Long (2023) Li, D., & Long, D. M. 2023, ApJ, 944, 8, doi: 10.3847/1538-4357/acacf4
  * Lin et al. (2003) Lin, Y., Engvold, O. r., & Wiik, J. E. 2003, Sol. Phys., 216, 109, doi: 10.1023/A:1026150809598
  * Luna et al. (2019) Luna, M., Oliver, R., Antolin, P., & Arregui, I. 2019, A&A, 629, A20, doi: 10.1051/0004-6361/201935850
  * Luna et al. (2008) Luna, M., Terradas, J., Oliver, R., & Ballester, J. L. 2008, ApJ, 676, 717, doi: 10.1086/528367
  * Luna et al. (2009) —. 2009, ApJ, 692, 1582, doi: 10.1088/0004-637X/692/2/1582
  * Luna et al. (2010) —. 2010, ApJ, 716, 1371, doi: 10.1088/0004-637X/716/2/1371
  * MacCormack (1969) MacCormack, R. W. 1969, AIAA Paper, 69, 354, doi: 10.2514/6.1969-354
  * Magyar & Van Doorsselaere (2016) Magyar, N., & Van Doorsselaere, T. 2016, ApJ, 823, 82, doi: 10.3847/0004-637X/823/2/82
  * Malanushenko & Schrijver (2013) Malanushenko, A., & Schrijver, C. J. 2013, ApJ, 775, 120, doi: 10.1088/0004-637X/775/2/120
  * Mann et al. (1995) Mann, I. R., Wright, A. N., & Cally, P. S. 1995, J. Geophys. Res., 100, 19441, doi: 10.1029/95JA00820
  * McCarthy et al. (2021) McCarthy, M. I., Longcope, D. W., & Malanushenko, A. 2021, ApJ, 913, 56, doi: 10.3847/1538-4357/abf4d5
  * Morton & Ruderman (2011) Morton, R. J., & Ruderman, M. S. 2011, A&A, 527, A53, doi: 10.1051/0004-6361/201016028
  * Nakariakov et al. (2016) Nakariakov, V. M., Anfinogentov, S. A., Nisticò, G., & Lee, D. H. 2016, A&A, 591, L5, doi: 10.1051/0004-6361/201628850
  * Nakariakov & Kolotkov (2020) Nakariakov, V. M., & Kolotkov, D. Y. 2020, ARA&A, 58, 441, doi: 10.1146/annurev-astro-032320-042940
  * Nakariakov & Ofman (2001) Nakariakov, V. M., & Ofman, L. 2001, A&A, 372, L53, doi: 10.1051/0004-6361:20010607
  * Nakariakov et al. (1999) Nakariakov, V. M., Ofman, L., Deluca, E. E., Roberts, B., & Davila, J. M. 1999, Science, 285, 862, doi: 10.1126/science.285.5429.862
  * Nakariakov & Verwichte (2005) Nakariakov, V. M., & Verwichte, E. 2005, Living Reviews in Solar Physics, 2, 3, doi: 10.12942/lrsp-2005-3
  * Nakariakov et al. (2021) Nakariakov, V. M., Anfinogentov, S. A., Antolin, P., et al. 2021, Space Sci. Rev., 217, 73, doi: 10.1007/s11214-021-00847-2
  * Nechaeva et al. (2019) Nechaeva, A., Zimovets, I. V., Nakariakov, V. M., & Goddard, C. R. 2019, ApJS, 241, 31, doi: 10.3847/1538-4365/ab0e86
  * Nisticò et al. (2013) Nisticò, G., Nakariakov, V. M., & Verwichte, E. 2013, A&A, 552, A57, doi: 10.1051/0004-6361/201220676
  * Ofman (2009) Ofman, L. 2009, ApJ, 694, 502, doi: 10.1088/0004-637X/694/1/502
  * Ofman & Wang (2008) Ofman, L., & Wang, T. J. 2008, A&A, 482, L9, doi: 10.1051/0004-6361:20079340
  * Okamoto et al. (2007) Okamoto, T. J., Tsuneta, S., Berger, T. E., et al. 2007, Science, 318, 1577, doi: 10.1126/science.1145447
  * Pascoe et al. (2010) Pascoe, D. J., Wright, A. N., & De Moortel, I. 2010, ApJ, 711, 990, doi: 10.1088/0004-637X/711/2/990
  * Petrova et al. (2023) Petrova, E., Magyar, N., Van Doorsselaere, T., & Berghmans, D. 2023, ApJ, 946, 36, doi: 10.3847/1538-4357/acb26a
  * Poedts & Kerner (1991) Poedts, S., & Kerner, W. 1991, Phys. Rev. Lett., 66, 2871, doi: 10.1103/PhysRevLett.66.2871
  * Reale (2014) Reale, F. 2014, Living Reviews in Solar Physics, 11, 4, doi: 10.12942/lrsp-2014-4
  * Roberts (2019) Roberts, B. 2019, MHD Waves in the Solar Atmosphere (Cambridge University Press), doi: 10.1017/9781108613774
  * Robertson & Ruderman (2011) Robertson, D., & Ruderman, M. S. 2011, A&A, 525, A4, doi: 10.1051/0004-6361/201015525
  * Robertson et al. (2010) Robertson, D., Ruderman, M. S., & Taroyan, Y. 2010, A&A, 515, A33, doi: 10.1051/0004-6361/201014055
  * Ruderman (2003) Ruderman, M. S. 2003, A&A, 409, 287, doi: 10.1051/0004-6361:20031079
  * Ruderman & Petrukhin (2023) Ruderman, M. S., & Petrukhin, N. S. 2023, MNRAS, 523, 2074, doi: 10.1093/mnras/stad1530
  * Ruderman et al. (2021) Ruderman, M. S., Petrukhin, N. S., & Pelinovsky, E. 2021, Sol. Phys., 296, 124, doi: 10.1007/s11207-021-01867-5
  * Ruderman & Roberts (2002) Ruderman, M. S., & Roberts, B. 2002, ApJ, 577, 475, doi: 10.1086/342130
  * Sakurai et al. (1991) Sakurai, T., Goossens, M., & Hollweg, J. V. 1991, Sol. Phys., 133, 227, doi: 10.1007/BF00149888
  * Schrijver (2007) Schrijver, C. J. 2007, ApJ, 662, L119, doi: 10.1086/519455
  * Schrijver et al. (2002) Schrijver, C. J., Aschwanden, M. J., & Title, A. M. 2002, Sol. Phys., 206, 69, doi: 10.1023/A:1014957715396
  * Schrijver & Brown (2000) Schrijver, C. J., & Brown, D. S. 2000, ApJ, 537, L69, doi: 10.1086/312753
  * Schrijver et al. (1999) Schrijver, C. J., Title, A. M., Berger, T. E., et al. 1999, Sol. Phys., 187, 261, doi: 10.1023/A:1005194519642
  * Sescu (2015) Sescu, A. 2015, Advances in Difference Equations, 2015, 1
  * Sewell (1988) Sewell, G. 1988, The Numerical Solution of Ordinary and Partial Differential Equations (San Diego: Academic Press)
  * Soler et al. (2013) Soler, R., Goossens, M., Terradas, J., & Oliver, R. 2013, ApJ, 777, 158, doi: 10.1088/0004-637X/777/2/158
  * Soler et al. (2014) —. 2014, ApJ, 781, 111, doi: 10.1088/0004-637X/781/2/111
  * Soler & Luna (2015) Soler, R., & Luna, M. 2015, A&A, 582, A120, doi: 10.1051/0004-6361/201526919
  * Soler et al. (2009) Soler, R., Oliver, R., & Ballester, J. L. 2009, ApJ, 693, 1601, doi: 10.1088/0004-637X/693/2/1601
  * Soler & Terradas (2015) Soler, R., & Terradas, J. 2015, ApJ, 803, 43, doi: 10.1088/0004-637X/803/1/43
  * Terradas et al. (2008a) Terradas, J., Andries, J., Goossens, M., et al. 2008a, ApJ, 687, L115, doi: 10.1086/593203
  * Terradas et al. (2008b) Terradas, J., Arregui, I., Oliver, R., et al. 2008b, ApJ, 679, 1611, doi: 10.1086/586733
  * Terradas et al. (2005) Terradas, J., Oliver, R., & Ballester, J. L. 2005, ApJ, 618, L149, doi: 10.1086/427844
  * Tian et al. (2012) Tian, H., McIntosh, S. W., Wang, T., et al. 2012, ApJ, 759, 144, doi: 10.1088/0004-637X/759/2/144
  * Tirry & Goossens (1996) Tirry, W. J., & Goossens, M. 1996, ApJ, 471, 501, doi: 10.1086/177986
  * Van Doorsselaere et al. (2004) Van Doorsselaere, T., Andries, J., Poedts, S., & Goossens, M. 2004, ApJ, 606, 1223, doi: 10.1086/383191
  * Van Doorsselaere et al. (2008a) Van Doorsselaere, T., Nakariakov, V. M., Young, P. R., & Verwichte, E. 2008a, A&A, 487, L17, doi: 10.1051/0004-6361:200810186
  * Van Doorsselaere et al. (2008b) Van Doorsselaere, T., Ruderman, M. S., & Robertson, D. 2008b, A&A, 485, 849, doi: 10.1051/0004-6361:200809841
  * Verwichte et al. (2009) Verwichte, E., Aschwanden, M. J., Van Doorsselaere, T., Foullon, C., & Nakariakov, V. M. 2009, ApJ, 698, 397, doi: 10.1088/0004-637X/698/1/397
  * Verwichte et al. (2004) Verwichte, E., Nakariakov, V. M., Ofman, L., & Deluca, E. E. 2004, Sol. Phys., 223, 77, doi: 10.1007/s11207-004-0807-6
  * Verwichte et al. (2013) Verwichte, E., Van Doorsselaere, T., White, R. S., & Antolin, P. 2013, A&A, 552, A138, doi: 10.1051/0004-6361/201220456
  * Wang & Sakurai (1998) Wang, H., & Sakurai, T. 1998, PASJ, 50, 111, doi: 10.1093/pasj/50.1.111
  * Wang et al. (2012) Wang, T., Ofman, L., Davila, J. M., & Su, Y. 2012, ApJ, 751, L27, doi: 10.1088/2041-8205/751/2/L27
  * Wang et al. (2008) Wang, T. J., Solanki, S. K., & Selwa, M. 2008, A&A, 489, 1307, doi: 10.1051/0004-6361:200810230
  * Watko & Klimchuk (2000) Watko, J. A., & Klimchuk, J. A. 2000, Sol. Phys., 193, 77, doi: 10.1023/A:1005209528612
  * White & Verwichte (2012) White, R. S., & Verwichte, E. 2012, A&A, 537, A49, doi: 10.1051/0004-6361/201118093
  * White et al. (2013) White, R. S., Verwichte, E., & Foullon, C. 2013, ApJ, 774, 104, doi: 10.1088/0004-637X/774/2/104
  * Williams et al. (2021) Williams, T., Walsh, R. W., & Morgan, H. 2021, ApJ, 919, 47, doi: 10.3847/1538-4357/ac0f76
  * Yi et al. (1991) Yi, Z., Engvold, O., & Keil, S. L. 1991, Sol. Phys., 132, 63, doi: 10.1007/BF00159130
  * Yu et al. (2021) Yu, H., Li, B., Chen, S., & Guo, M. 2021, Sol. Phys., 296, 95, doi: 10.1007/s11207-021-01839-9
  * Zajtsev & Stepanov (1975) Zajtsev, V. V., & Stepanov, A. V. 1975, Issledovaniia Geomagnetizmu Aeronomii i Fizike Solntsa, 37, 3
  * Zhong et al. (2023) Zhong, S., Nakariakov, V. M., Kolotkov, D. Y., et al. 2023, Nature Communications, 14, 5298, doi: 10.1038/s41467-023-41029-8
  * Zhong et al. (2022) Zhong, S., Nakariakov, V. M., Kolotkov, D. Y., Verbeeck, C., & Berghmans, D. 2022, MNRAS, 516, 5989, doi: 10.1093/mnras/stac2545
  * Zimovets & Nakariakov (2015) Zimovets, I. V., & Nakariakov, V. M. 2015, A&A, 577, A4, doi: 10.1051/0004-6361/201424960

Figure 1:  Schematics showing both the two-elliptic-tube configuration and
lower-order kink motions therein. Illustrated in (b) are two identical tubes
with elliptic cross-sections and separated in the $x$-direction. The $y-z$
plane is identified as the tube plane, the meaning of which is clearer in the
curved representation in panel a. This study examines only two tube
orientations (panels c and d), with the one labeled “$x$-major” (“$x$-minor”)
corresponding to the situation where the major- (minor-)axis of either tube is
aligned with the $x$-axis. The tube centers are placed symmetrically with
respect to the $y$-axis. Only axial fundamentals are considered, making it
possible to examine the three-dimensional (3D) wave fields with a 2D approach.
Lower-order kink motions are classified according to their internal flow
fields in the apex plane (panels c or d); hence such symbols as $S_{x}$ and
$A_{y}$. See text for more details.

Figure 2: Kink motions in isolated tubes with circular cross-sections for a
fixed pair of density contrast and dimensionless axial wavenumber $[\rho_{\rm
i}/\rho_{\rm e}=3,kb=\pi/30]$. A specific dimensionless layer width
$\bar{l}=0.4$ is examined in (a), the solid curve showing the $x$-speed
sampled at the tube center ($v_{x}(0,0,t)$). An exponentially damped cosine is
employed to fit the segment encompassing the extrema shown by the red
asterisks, with the resulting damping envelope plotted over the entire
interval (the dashed curves). Presented in (b) are the $\bar{l}$-dependencies
of the mode frequencies ($\Omega$, the black curves and symbols) and the
ratios of the decay rate to the frequency ($\gamma/\Omega$, blue) obtained
with various methods. The solid curves represent the analytical thin-tube-
thin-boundary (TTTB) expectations with Eq. (30). No assumption on tube length
or layer width is imposed in the “Resistive” computations (the dashed curves)
that evaluate ideal quasi-modes as resistive eigenmodes. Shown by the open
circles are the best-fit results from our time-dependent computations, where
kink motions are consistently excited with a perturbation given by Eq. (31).

Figure 3:  Kink motions with (a) $S_{x}$ and (b) $A_{x}$ patterns in a two-
circular-tube configuration. Plotted are the time sequences of the $x$-speed
sampled at the left tube center ($v_{x}(-d/2,0;t)$, the solid curves) together
with the best-fit damping envelopes (dashed). The time-dependent computations
pertain to a fixed combination $[\rho_{\rm i}/\rho_{\rm
e}=5,\ell/\mathcal{R}=0.4,d/\mathcal{R}=2.5,kb=\pi/30]$, with $\ell$ and
$\mathcal{R}$ being the nonuniform layer width and mean tube radius
implemented by Ruderman & Roberts (2002). A fitting procedure with an
exponentially damped cosine is performed over the duration exactly
encompassing the extrema indicated by the red asterisks, even though the best-
fit envelopes are plotted for the entire time interval. Kink motions are
excited with the perturbation described by Eq. (25) with appropriate signs.

Figure 4:  Comparison of (a) the oscillation frequencies $\Omega$ and (b)
damping rates $\gamma$ of ideal quasi-modes deduced from time-dependent
computations (the open circles) with those expected with the T-matrix
formalism (solid lines). With different colors we discriminate the results for
the $S_{x}$ and $A_{x}$ patterns. The T-matrix curves are read from Fig. 4 of
Soler & Luna (2015, SL15). Shown here is how $\Omega$ and $\gamma/\Omega$
depend on $\ell/\mathcal{R}$, while the rest of the parameters is fixed at
$[\rho_{\rm i}/\rho_{\rm e}=5,d/\mathcal{R}=2.5,kb=\pi/30]$. The
$[\ell/\mathcal{R},d/\mathcal{R}]$ notations follow from SL15 for the ease of
comparison.

Figure 5:  Temporal profiles of the $x$-speed sampled at the left tube center
($v_{x}(-d/2,0;t)$, the solid curves) for a number of domain sizes as labeled.
Overplotted by the dashed curves are the best-fit damping envelopes. All
computations pertain to the $S_{x}$ motion for a fixed combination of physical
parameters $[\rho_{\rm i}/\rho_{\rm
e}=5,\ell/\mathcal{R}=0.4,d/\mathcal{R}=2.5,kb=\pi/30]$. The grid spacing is
fixed at $\Delta=0.01b$.

Figure 6:  Temporal profiles of the $x$-speed sampled at the left tube center
for a number of grid spacings as labeled. All computations pertain to the
$S_{x}$ motion for a fixed combination of physical parameters $[\rho_{\rm
i}/\rho_{\rm e}=5,\ell/\mathcal{R}=0.4,d/\mathcal{R}=2.5,kb=\pi/30]$. The
computational domain size is fixed at $[x_{\rm M}=12b,y_{\rm M}=10b]$.

Figure 7:  Distributions of the $x$-speed at $t=340b/v_{\rm Ae}$ along the
$y$-cut through the left tube center for a number of grid spacings as labeled.
All computations pertain to the $S_{x}$ motion for a fixed combination of
physical parameters $[\rho_{\rm i}/\rho_{\rm
e}=5,\ell/\mathcal{R}=0.4,d/\mathcal{R}=2.5,kb=\pi/30]$. The computational
domain size is fixed at $[x_{\rm M}=12b,y_{\rm M}=10b]$. Only the lower half
($y\leq 0$) is shown for any $y$-profile, the other half being symmetric about
$y=0$. The vertical dash-dotted lines delineate the borders of the nonuniform
layer.

Figure 8:  Temporal profiles of some energetics-related quantities for several
grid spacings as differentiated by different colors. Energetics are
consistently examined for a box $[-(d/2+1.2b),(d/2+1.2b)]\times[-1.2b,1.2b]$,
in which the three mutually exclusive constituents are the tube interiors,
nonuniform layers, and exterior. Shown in (a) are the time derivative of the
total perturbation energy in the box (${\rm d}E_{\rm tot}/{\rm d}t$, the sold
curves), together with the instantaneous energy flux into this box ($-F$,
asterisks). Further plotted in (b) are the instantaneous total energy in the
interiors ($E_{\rm int}$), nonuniform layers ($E_{\rm layer}$), and the
exterior ($E_{\rm ext}$). All computations pertain to the $S_{x}$ pattern for
a fixed combination of physical parameters $[\rho_{\rm i}/\rho_{\rm
e}=5,\ell/\mathcal{R}=0.4,d/\mathcal{R}=2.5,kb=\pi/30]$. The computational
domain size is fixed at $[x_{\rm M}=12b,y_{\rm M}=10b]$.

Figure 9:  Kink motions with (a) $S_{x}$ and (b) $A_{x}$ patterns in a two-
elliptic-tube configuration with $x$-major orientation. Plotted are the time
sequences of the $x$-speed sampled at the left tube center ($v_{x}(-d/2,0;t)$,
the solid curves) for a number of ratios of the semi-major to semi-minor axis
($a/b$) as discriminated by the different colors. A fitting procedure is
performed, for each sequence, over the duration encompassing the six extrema
starting with the one denoted by an asterisk. The best-fit damping envelope is
plotted for the entire sequence, but only for $a/b=2.5$ to avoid over-crowding
the plots. The computations pertain to a fixed combination $[\rho_{\rm
i}/\rho_{\rm e}=3,\bar{l}=0.4,d/b=5,kb=\pi/30]$.

Figure 10:  Representative snapshots of the velocity field (the blue arrows)
and energy density distribution ($\epsilon$, filled contours) for (a) the
$S_{x}$ motion and (b) the $A_{x}$ one in a two-elliptic-tube configuration
with the $x$-major orientation. The inset in each panel emphasizes the
distribution of $\epsilon$ in a representative nonuniform layer, whose
boundaries are shown by the red solid curves. The white contour indicates
where the Alfv$\acute{\rm e}$n frequency $\omega_{\rm A}=kv_{\rm A}$ equals
the quasi-mode frequency deduced with the fitting procedure. Both computations
pertain to a fixed combination $[\rho_{\rm i}/\rho_{\rm
e}=3,a/b=2.5,\bar{l}=0.4,d/b=5,kb=\pi/30]$. This snapshot is extracted from
the animation attached to the current figure. The $\epsilon$ contours in all
snapshots are filled in such a way that darker portions correspond to larger
$\epsilon$. The blue arrows are also consistently scaled in all snapshots such
that longer arrows correspond to stronger velocities.

Figure 11:  Distributions of the $x$-speed along the $y$-cut through the left
tube center at a number of instants as labeled. Panels (a) and (b) correspond
to the $S_{x}$ and $A_{x}$ motions, respectively. The vertical dotted lines
represent the borders of the nonuniform layer. Only the lower half ($y<0$) of
a $y$-profile is displayed, the other half being symmetric with respect to
$y=0$. All computations are conducted for the $x$-major orientation with a
fixed combination of physical parameters $[\rho_{\rm i}/\rho_{\rm
e}=3,\bar{l}=0.4,d/b=5,kb=\pi/30]$.

Figure 12:  Dependencies of (a) the quasi-mode periods ($P$) and (b) damping-
time-to-period ratios ($\tau/P$) on the ratio of semi-major to semi-minor axis
($a/b$) for kink motions in a two-elliptic-tube configuration with the
$x$-major orientation. The $S_{x}$ and $A_{x}$ patterns are discriminated by
the linestyles. A number of dimensionless tube separations ($d/b$) are
examined as differentiated by the different colors, with $d/b$ constrained to
avoid tube overlapping. All results pertain to a fixed combination of physical
parameters $[\rho_{\rm i}/\rho_{\rm e}=3,\bar{l}=0.4,kb=\pi/30]$. The
horizontal lines represent the results for an isolated circular tube with the
same set of parameters.

Figure 13:  Similar to Fig. 9 but for the $x$-minor orientation.

Figure 14:  Similar to Fig. 10 but for the $x$-minor orientation. The
associated animation is also attached.

Figure 15:  Similar to Fig. 11 but for the $x$-minor orientation.

Figure 16:  Similar to Fig. 12 but for the $x$-minor orientation.