# The effect of linear background rotational flows on magnetoacoustic modes of
a photospheric magnetic flux tube

S. J. Skirvin1,2, V. Fedun2, S. S. A. Silva2, T. Van Doorsselaere1, N. Claes1,
M. Goossens1 and G. Verth3
1Centre for mathematical Plasma Astrophysics, Mathematics Department, KU
Leuven, Celestijnenlaan 200B bus 2400, B-3001 Leuven, Belgium.
2Plasma Dynamics Group, Department of Automatic Control & Systems Engineering,
The University of Sheffield, Sheffield, S3 7RH, UK
3Plasma Dynamics Group, Department of Mathematics and Statistics, The
University of Sheffield, Sheffield, S3 7RH, UK
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Magnetoacoustic waves in solar magnetic flux tubes may be affected by the
presence of background rotational flows. Here, we investigate the behaviour of
$m=0$ and $m=\pm 1$ modes of a magnetic flux tube in the presence of linear
background rotational flows embedded in a photospheric environment. We show
that the inclusion of a background rotational flow is found to have little
effect on the obtained eigensolutions for the axisymmetric $m=0$ sausage mode.
However, solutions for the kink mode are dependent on the location of the flow
resonance modified by the slow frequency. A background rotational flow causes
the modified flow resonances to possess faster phase speeds in the thin-tube
(TT) limit for the case $m=1$. This results in solutions for the slow body and
slow surface kink modes to follow this trajectory, changing their dispersive
behaviour. For a photospheric flux tube in the TT limit, we show that it
becomes difficult to distinguish between the slow surface and fast surface
kink ($m=1$) modes upon comparison of their eigenfunctions. 2D velocity field
plots demonstrate how these waves, in the presence of background rotational
flows, may appear in observational data. For slow body kink modes, a swirling
pattern can be seen in the total pressure perturbation. Furthermore, the tube
boundary undergoes a helical motion from the breaking of azimuthal symmetry,
where the $m=1$ and $m=-1$ modes become out of phase, suggesting the resulting
kink wave is circularly polarised. These results may have implications for
seismology of magnetohydrodynamic waves in solar magnetic vortices.

###### keywords:

magnetohydrodynamics (MHD) – waves

††pubyear: 2022††pagerange: The effect of linear background rotational flows
on magnetoacoustic modes of a photospheric magnetic flux tube–The effect of
linear background rotational flows on magnetoacoustic modes of a photospheric
magnetic flux tube

##  1 Introduction

High resolution observations have revealed, over the last few decades, that
magnetohydrodynamic (MHD) waves are ubiquitous throughout the Sun’s atmosphere
(Nakariakov et al., 1999; Aschwanden et al., 1999; De Pontieu et al., 2007;
Tomczyk et al., 2007; Morton et al., 2015; Grant et al., 2015; Keys et al.,
2018; Stangalini et al., 2022; Bate et al., 2022). Understanding MHD wave
properties from a theoretical point of view is of the utmost importance in
solar physics as these waves could contribute to the heating of local plasma
and may also be used as a proxy to determine sub-resolution plasma properties.
Additionally, observed MHD wave behaviour can be used to provide an estimate
of the properties of local plasma that cannot be measured directly, for
example the magnetic field in the corona. Furthermore, MHD waves may also be
responsible in the formation of some observed phenomena, e.g. jets, in the
solar atmosphere (De Pontieu et al., 2004; Rouppe van der Voort et al., 2007;
Scullion et al., 2011).

The uniform cylindrical waveguide model (see e.g. Wilson, 1979; Spruit &
Zweibel, 1979; Edwin & Roberts, 1983) has provided a foundation for analytical
studies into MHD wave investigations under solar atmospheric conditions. These
studies have shown that the trapped wave modes of a magnetic flux tube can be
described by the number of nodes present in each geometrical direction. The
axisymmetric sausage mode has zero nodes in the azimuthal direction $m=0$,
whereas the (linearly polarised) non-axisymmetric kink mode has a single
azimuthal node $m=1$ or $m=-1$. Whilst this uniform model is simplistic in
nature, due to the uniform plasma considered inside and outside the waveguide,
it has been modified over recent decades to model specific configurations
which better match observational data (see e.g. Van Doorsselaere et al., 2004;
Verth et al., 2007; Erdélyi & Fedun, 2010; Aldhafeeri et al., 2021; Ruderman &
Petrukhin, 2022). The cylinder model can be further extended to incorporate
additional physical environments that may be common configurations in magnetic
flux tubes, such as including a background rotational flow. Magnetic flux
tubes with background rotational flows are a common configuration observed in,
e.g., flux tubes rooted in intergranular lanes, solar tornadoes and spicules
(Bonet et al., 2010; Wedemeyer-Böhm et al., 2012; Tziotziou et al., 2018;
Shetye et al., 2019). Such structures display dynamic characteristics and
excite a wide range of MHD waves which couple different layers of the solar
atmosphere. As a result, these structures act as a natural conduit for the
transfer of mass, momentum and energy throughout the solar atmosphere.
Furthermore, rotational flows naturally appear in numerical MHD simulations of
regions in the solar atmosphere with vortex drivers (see, e.g. Fedun et al.,
2011a; Fedun et al., 2011b; Shelyag et al., 2011, 2012, 2013; González-Avilés
et al., 2017; González-Avilés et al., 2018; Snow et al., 2018) and also in
magnetoconvection simulations (Yadav et al., 2020, 2021; Silva et al., 2021).

The presence of a background rotational flow would manifest itself as a non-
zero azimuthal component of the background velocity field vector. Considering
$m=\pm 1$ transverse kink modes, the effect of introducing a background
rotational flow into the model would break the symmetry of the system with
respect to the direction of azimuthal wave propagation. Similarly, it has been
shown that the inclusion of a steady vertical background plasma flow, aligned
with the magnetic field, can obtain solutions which also break the symmetry of
the forward and backward propagating waves (Nakariakov & Roberts, 1995; Terra-
Homem et al., 2003; Soler et al., 2009; Skirvin et al., 2022). The presence of
the longitudinal flow results in a Doppler-shifted wave frequency, depending
upon the amplitude of the flow and the vertical wavenumber, which may shift
certain wave modes into different physical regimes. Furthermore, a background
longitudinal flow may also affect the wavelength and damping length of
resonant absorption for kink modes, causing more efficient damping of backward
propagating modes compared to forward propagating modes (Soler et al., 2011).

Circularly polarised kink modes, previously observed in chromospheric magnetic
elements (Stangalini et al., 2017) and sunspots (Jess et al., 2017), have also
been studied recently by Magyar et al. (2022) where the authors discuss the
differences that resonant absorption and phase mixing have on linearly and
circularly polarised kink waves. Magyar et al. (2022) also conclude, upon
analysis of the Doppler signatures for both polarisation states, that there is
very little difference between the two polarisation’s in Doppler observations.
Circularly polarised kink waves have been shown to develop in coronal loops
with twisted magnetic fields (Terradas & Goossens, 2012; Ruderman & Terradas,
2015).

Previous analytical studies have investigated the stability status of rotating
flux tubes, as an azimuthal velocity shear across the waveguide boundary may
be susceptible to the Kelvin-Helmholtz instability (KHI) (Soler et al., 2010;
Zaqarashvili et al., 2015; Zhelyazkov & Chandra, 2019). However, these studies
assume zero plasma-$\beta$ and focus on coronal conditions only, ignoring the
slow magnetoacoustic modes entirely. The stability of a magnetic flux tube
with a linear background magnetic twist and rotational flow component was
studied by Cheremnykh et al. (2018) who found that the $m=0$ sausage mode
becomes unstable for azimuthal flow speeds that create a centrifugal force
which can overcome the magnetic tension, whereas, the $m=1$ kink mode can only
become unstable for sufficiently large values of longitudinal (axial) flow
speed. Whilst the stability studies mentioned above investigate the
susceptibility of the KHI to various wave modes in magnetic flux tubes with
background flows, we stress that it is not the goal of this study. The primary
aim of the present work, however, is to provide a description of how $m=0$ and
$m=\pm 1$ wave modes would manifest themselves in solar magnetic flux tubes,
under photospheric conditions, in the presence of background rotational flows
and to aid mode interpretations of observational and numerical data.

This paper is presented as follows, in Section 2 the governing equations
describing a magnetic flux tube in the presence of a background rotational
flow is presented, along with a brief description of the numerical eigensolver
implemented in this work. Section 3 presents the results of an investigation
into the effect of a linear profile of background rotational flow on the
properties of magnetoacoutsic waves under photospheric conditions. Section 3
looks closely at the obtained eigenvalues, one dimensional eigenfunctions and
two/three dimensional visualisations of the kink mode in the presence of a
background rotational flow. Finally, in Section 4, we summarise the findings
presented in this work and discuss avenues of future research.

##  2 Method

The ideal MHD equations adopted in this study are:

$\displaystyle\frac{d\rho}{dt}+\rho\nabla\cdot\textbf{v}$ $\displaystyle=0,$
(1) $\displaystyle\rho\left(\frac{d\textbf{v}}{dt}\right)$
$\displaystyle=-\nabla
p+\frac{1}{\mu_{0}}\left(\nabla\times\textbf{B}\right)\times\textbf{B},$ (2)
$\displaystyle\frac{d}{dt}\left(\frac{p}{\rho^{\gamma}}\right)$
$\displaystyle=0,$ (3) $\displaystyle\frac{\partial\textbf{B}}{\partial t}$
$\displaystyle=\nabla\times\left(\textbf{v}\times\textbf{B}\right),$ (4)
$\displaystyle\nabla\cdot\textbf{B}$ $\displaystyle=0,$ (5)

where $\rho$, v, $p$, B, $\gamma$ and $\mu$ denote plasma density, plasma
velocity, plasma pressure, magnetic field, ratio of specific heats (taken
$\gamma=5/3$) and the magnetic permeability respectively. We consider a
cylindrical geometry $(r,\varphi,z)$ where the initial equilibrium has a
radially spatially dependent form of the velocity field vector
$\textbf{v}_{0}=(0,v_{\varphi}(r),0)$. In our model the magnetic field is
taken to be straight and uniform with no background axial plasma flow. In the
lower solar atmosphere it should be noted that magnetic flux tubes possess
significant non-vertical magnetic field as a result of flux tube expansion to
maintain pressure balance, however, this effect can be considered to be
negligible in the current study due to the analysis in the local plasma
environment. Since the equilibrium quantities depend on $r$ only, the
perturbed quantities can be Fourier-analysed with respect to the ignorable
coordinates $\varphi,z$ and time $t$ and put proportional to:

$\text{exp}\left[i\left(m\varphi+kz-\omega t\right)\right],$

where $m$ is the azimuthal wave number, $k$ the vertical wavenumber and
$\omega$ the wave frequency.

After linearising Equations (1)-(5), we arrive at a system of two differential
equations containing the total pressure perturbation $\hat{P}_{T}$ and the
radial displacement perturbation $r\hat{\xi}_{r}$ (see e.g. Sakurai et al.,
1991; Goossens et al., 1992), which can be written as:

$\displaystyle D\frac{d}{dr}\left(r\hat{\xi}_{r}\right)$
$\displaystyle=C_{1}r\hat{\xi}_{r}-C_{2}r\hat{P}_{T},$ (6) $\displaystyle
D\frac{d\hat{P}_{T}}{dr}$ $\displaystyle=C_{3}\hat{\xi}_{r}-C_{1}\hat{P}_{T},$
(7)

where,

$\displaystyle D$
$\displaystyle=\rho\left(c^{2}+v_{A}^{2}\right)\left(\Omega^{2}-k^{2}v_{A}^{2}\right)\left(\Omega^{2}-k^{2}c_{T}^{2}\right),$
(8) $\displaystyle\Omega$ $\displaystyle=\omega-\frac{m}{r}v_{\varphi},$ (9)
$\displaystyle c^{2}$ $\displaystyle=\frac{\gamma p}{\rho},\ \ \ \
v_{A}^{2}=\frac{B_{z}^{2}}{\mu\rho},\ \ \ \
c_{T}^{2}=\frac{v_{A}^{2}c^{2}}{\left(c^{2}+v_{A}^{2}\right)},$ (10)
$\displaystyle C_{1}$
$\displaystyle=Q\Omega^{2}-2m\left(c^{2}+v_{A}^{2}\right)\left(\Omega^{2}-k^{2}c_{T}^{2}\right)\frac{T^{2}}{r^{2}},$
(11) $\displaystyle C_{2}$
$\displaystyle=\Omega^{4}-\left(c^{2}+v_{A}^{2}\right)\left(\frac{m^{2}}{r^{2}}+k^{2}\right)\left(\Omega^{2}-k^{2}c_{T}^{2}\right),$
(12) $\displaystyle C_{3}$
$\displaystyle=D\left\\{\rho\left(\Omega^{2}-k^{2}v_{A}^{2}\right)+r\frac{d}{dr}\left[-\rho\left(\frac{v_{\varphi}}{r}\right)^{2}\right]\right\\}+$
(13)
$\displaystyle+Q^{2}-4\left(c^{2}+v_{A}^{2}\right)\left(\Omega^{2}-k^{2}c_{T}^{2}\right)\frac{T^{2}}{r^{2}},$
$\displaystyle Q$
$\displaystyle=-\left(\Omega^{2}-k^{2}v_{A}^{2}\right)\frac{\rho
v_{\varphi}^{2}}{r},$ (14) $\displaystyle T$ $\displaystyle=\rho\Omega
v_{\varphi}.$ (15)

It should be noted here that we assume $B_{\varphi}=v_{z}=0$, which simplifies
the fully inclusive
($B_{\varphi}=B_{\varphi}(r),B_{z}=B_{z}(r),v_{z}=v_{z}(r),\rho=\rho(r)$) set
of equations previously noted in literature by Goossens et al. (1992).
Quantities $c^{2}$, $v_{A}^{2}$, and $c_{T}^{2}$ define the squares of the
local sound speed, Alfvén speed and cusp (tube) speed respectively. The
quantity $\Omega$ represents the Doppler shifted frequency as a result of the
background rotational plasma flow.

Equations (6)-(7) can be combined to create a single differential equation in
either $r\hat{\xi}_{r}$:

$\frac{d}{dr}\left[f(r)\frac{d}{dr}\left(r\hat{\xi}_{r}\right)\right]-g(r)\left(r\hat{\xi}_{r}\right)=0,$
(16)

where,

$f(r)=\frac{D}{rC_{2}},$ (17)
$g(r)=\frac{d}{dr}\left(\frac{C_{1}}{rC_{2}}\right)-\frac{1}{rD}\left(C_{3}-\frac{C_{1}^{2}}{C_{2}}\right),$
(18)

or $\hat{P}_{T}$:

$\frac{d}{dr}\left[\tilde{f}(r)\frac{d\hat{P}_{T}}{dr}\right]-\tilde{g}(r)\hat{P}_{T}=0,$
(19)

where,

$\tilde{f}(r)=\frac{rD}{C_{3}},$ (20)
$\tilde{g}(r)=-\frac{d}{dr}\left(\frac{rC_{1}}{C_{3}}\right)-\frac{r}{D}\left(C_{2}-\frac{C_{1}^{2}}{C_{3}}\right).$
(21)

For a non-uniform plasma, the governing Equations (6)-(7) possess regular
singularities where the wave frequency matches the local characteristic
frequencies at:

$\omega=\frac{m}{r}v_{\varphi}(r)\pm kv_{A},$ (22)
$\omega=\frac{m}{r}v_{\varphi}(r)\pm kc_{T}.$ (23)

Equations (22) and (23) define the flow continua modified by the local Alfvén
($kv_{A}$) and slow ($kc_{T}$) frequencies, respectively. In ideal MHD, the
wave solutions existing inside the continua, with positions given by Equations
(22) and (23), are known as ‘quasi-modes’ where the wave frequency becomes a
complex quantity (De Groof & Goossens, 2000; Goedbloed & Poedts, 2004;
Geeraerts et al., 2022). The nature of the solutions lying inside the continua
is not discussed in the present study and instead will be the focus of future
work.

Both Equations (16) and (19) have no known closed form analytical solutions,
without making assumptions that somehow reduce the mathematical complexity.
Therefore, investigating the properties of wave modes propagating within an
equilibrium which is non-uniform must be done numerically. The numerical
approach used in this study is based on the eigensolver applied in Skirvin et
al. (2021); Skirvin et al. (2022) for non-uniform magnetic slabs and non-
uniform flux tubes, respectively. The eigensolver implements the numerical
shooting and bisection methods whilst also relying on fundamental properties
of the sausage and kink modes. The numerical shooting method solves Equations
(16) and (19) ensuring continuity of $\hat{P}_{T}$ and $\hat{\xi}_{r}$ across
the boundary of the flux tube. This technique has been applied before in solar
physics by, e.g. Tirry & Goossens (1996); Pinter et al. (1998); Andries et al.
(2000); Taroyan & Erdélyi (2002, 2003). These studies also utilise the jump
conditions introduced by Sakurai et al. (1991); Goossens et al. (1992) to deal
with the regular singularities that appear in the governing equations,
allowing the authors to study the resulting quasi-modes. However, following
the primary objective of the present study, we only consider eigenmodes with
real valued wave frequency and wavenumber.

##  3 Magnetic flux tube in the presence of a linear background rotational
flow

In this section, a magnetic flux tube in the presence of a linear rotational
background flow is investigated. For all cases the magnetic flux tube is
otherwise uniform such that the equilibrium plasma density and magnetic field
is constant across the flux tube. A profile comparable to the magnetic twist
profile incorporated by Erdélyi & Fedun (2007) and Erdélyi & Fedun (2010) is
chosen but applied to the azimuthal velocity field component $v_{\varphi}$
instead. A rotational flow can be either clockwise or counter-clockwise in the
reference frame relative to the observer. The only difference between a
clockwise rotational flow and an anti-clockwise rotational flow will be the
sign in front of $v_{\varphi}$ and the direction of shifted wave frequency
relative to the flow. In the following sections, a magnetic flux tube is
presented with an equilibrium azimuthal flow component, acting in a counter-
clockwise direction, which takes the form:

$v_{\varphi}(r)=A\left(\frac{r}{a}\right)^{\alpha},$ (24)

where $A$ is the amplitude of the rotational flow, $\alpha$ is the parameter
(exponent) dictating the radial profile of the rotational flow and $a$
indicates the location of the boundary of the flux tube, taken to be $a=1$ in
our study. The case when $\alpha=1$ for example corresponds to a linear
rotational flow, which will be the focus of this study, (see e.g. Figure 1).
It can be shown that when $v_{\varphi}$ is linear, Equations (6)-(15)
simplify, and in many cases the dependence on the radial coordinate is
removed. The rotational flow is constant with height $z$ in all cases
considered in this work. Obtaining an equilibrium in a magnetic cylinder with
a background rotational velocity component is not as mathematically simple as
the scenario of a uniform magnetic cylinder. In order to maintain total
pressure balance across the waveguide the following expression must be
satisfied (Goossens et al., 2011):

$\frac{d}{dr}\left(p+\frac{B_{0z}^{2}}{2\mu}\right)=\frac{\rho
v_{\varphi}^{2}(r)}{r}=\rho A^{2}r.$ (25)

Integration of Equation (25) yields:

$p+\frac{B_{0z}^{2}}{2\mu}=\rho A^{2}\frac{r^{2}}{2},$ (26)

where the constant of integration is absorbed into the gas pressure term, $p$,
and corresponds to the plasma pressure on the axis of the cylinder where the
amplitude of the flow is zero (see e.g. Cheremnykh et al., 2018). Under the
photospheric conditions considered in this work, the total pressure balance is
achieved by an increase in temperature to balance the increase in azimuthal
flow amplitude towards the boundary of the flux tube. For configurations where
the amplitude of the rotational flow is weak (e.g. $A<0.5c_{i}$), then the
change in spatial behaviour of the plasma pressure and temperature is small,
but must not be dismissed.

The presence of a background rotational flow not only modifies the equilibrium
pressure balance relationship, but also affects the continuity conditions on
the boundary of the waveguide. Considering a magnetic flux tube in the
presence of a background rotational flow, the resulting boundary continuity
conditions state:

$\displaystyle\hat{\xi}_{re}\Bigr{\rvert}_{r=a}$
$\displaystyle=\hat{\xi}_{ri}\Bigr{\rvert}_{r=a},$ (27)
$\displaystyle\hat{P}_{Te}\Bigr{\rvert}_{r=a}$
$\displaystyle=\left(\hat{P}_{Ti}+\frac{\rho_{0i}v_{\varphi}^{2}}{a}\hat{\xi}_{ri}\right)\Bigr{\rvert}_{r=a}.$
(28)

The change in boundary conditions are accounted for in the numerical
eigensolver, and a pair of eigenvalues will only be retrieved for values
satisfying the above conditions for each respective case study.

Finally, when the background rotational flow is taken to be linear with
respect to the radial direction, Equations (22) and (23) become:

$\displaystyle\omega$ $\displaystyle=Am\pm kv_{A},$ (29) $\displaystyle\omega$
$\displaystyle=Am\pm kc_{T}.$ (30)

Equations (29) and (30) no longer define continuum regions as the resonant
positions no longer cover a range of frequencies, rather they define singular
resonant locations which correspond to the flow resonance position modified by
the Alfvén and slow frequencies respectively. The locations of these resonant
positions depend heavily on the amplitude of the rotational flow in the setup
presented in this work.

###  3.1 Linear rotating magnetic flux tube under photospheric conditions

Figure 1: Equilibrium background rotational flow profiles for cases with
increasing amplitude for a photospheric cylinder. In all cases the profiles
are linear with respect to spatial coordinate $r$ (i.e. $\alpha=1$). The
amplitude of the rotational flow increases linearly up to a boundary value of
$A=0.01$ (green line), $A=0.05$ (black line), $A=0.1$ (yellow line), $A=0.15$
(red line) and $A=0.25$ (blue line). The boundary of the flux tube is located
at $r=1$.

In this section, a magnetic flux tube under photospheric conditions
($v_{Ae}<c_{i}<c_{e}<v_{Ai}$) in the presence of a linear background
rotational flow is investigated. For all photospheric cases in this work, the
numerical values correspond to $c_{i}=1$, $c_{e}=1.5c_{i}$, $v_{Ai}=2c_{i}$,
$v_{Ae}=0.5c_{i}$ and $\rho_{i}=1$. This choice of equilibrium parameters
results in a density contrast between the internal and external plasma to be
roughly $\rho_{i}/\rho_{e}=0.567$. The background velocity vector inside the
waveguide can be written as $\mathbf{v}_{0i}=(0,Ar,0)$. The flow outside the
cylinder is zero which results in a velocity shear across the cylinder
boundary at $r=a$, however the value of $A$ is chosen to be small and both
sub-Alfvénic, sub-sonic and below the threshold for the Kelvin-Helmholtz
instability (Soler et al., 2010). This choice of amplitude $A$ also agrees
with observed values of photospheric flows when compared with the local
sound/Alfvén speed (Bonet et al., 2008). Shown in Figure 1 are the linear
profiles of background rotational flow considered in this section. In all
cases the flow amplitude is proportional to the radial distance from the
center of the flux tube, at $r=0$, up to the boundary at $r=1$, however the
amplitude is allowed to vary.

####  3.1.1 Eigenvalues

(a)

(b)

(c)

(d)

Figure 2: Dispersion diagrams for a photospheric cylinder with a linear
background rotational flow of varying flow amplitudes. The different cases
with varying amplitude, displayed on the top of each panel, are shown
corresponding to those in Figure 1. The red curves indicate solutions for the
$m=0$ sausage mode and the blue curves show the $m=+1$ kink mode solutions.
The green curve in all cases shows the flow resonance locations modified by
the slow frequency as given by Equation (30). Dash-dotted lines indicate the
equilibrium characteristic kink speed, $c_{k}$ sound speed, $c$, and tube
speed, $c_{T}$, internal (subscript ‘i’) and external (subscript ‘e’) to the
flux tube. The dashed box region shown in Figure 2(c) is investigated in more
detail in Figure 5.

Figure 2 highlights the change in eigenvalues for the different cases of flow
profiles considered in Figure 1. For now, we focus on the forward propagating
wave modes ($\omega/k>0$). The axisymmetric $m=0$ sausage mode appears
unaffected, or at least not significantly affected, by the presence of the
background flow, which is a result of the $m=0$ mode being the axisymmetric
mode in the azimuthal direction. As a result of this, the rotational flow does
not break the symmetry of this mode. This can also be seen analytically, by
setting $m=0$ in Equations (6)-(15), the majority of terms containing the
presence of the (linear) background flow are removed when $m=0$. Furthermore,
when $m=0$, the resonant locations described by Equations (29) and (30) reduce
to the resonant positions corresponding to a uniform magnetic flux tube (Edwin
& Roberts, 1983). Whilst Equations (6)-(15) do contain some terms, for example
in the variables $Q$ and $T$, which are dependent on the background rotational
flow and that remain for the $m=0$ solution, we do not observe any significant
modification to the obtained eigenvalues for the sausage mode in Figure 2 when
compared to the analytical solutions for the uniform magnetic flux tube,
although very minor effects can be seen when the amplitude of the rotational
flow is sufficiently large ($A>0.25$). Investigating this result both
analytically and numerically may be a focus of future work.

Conversely, there is a considerable effect on the $m=1$ kink mode solutions
due to the presence of a background rotational flow. In the long wavelength
(thin-tube) limit, the phase speeds of the slow body and slow surface kink
modes tend to an infinite phase speed, similar to the case study of a linear
background magnetic twist (Erdélyi & Fedun, 2010), where they may even enter
the leaky regime. As the amplitude of the azimuthal flow is increased, the
corresponding phase speeds of the slow kink modes also increases for all
values of wavenumber. The flow resonance locations modified by the slow
frequency given by Equation (30) is shown by the green line in Figure 2. The
slow body and slow surface kink modes follow this resonance curve in the long
wavelength limit and even undergo an avoided crossing (Abdelatif, 1990; Mather
& Erdélyi, 2016; Allcock & Erdélyi, 2017) where the slow surface mode
approaches the fast surface mode. This avoided crossing implies a transfer of
properties between the fast and slow surface modes. In the long wavelength
limit it may not be appropriate to refer to the slow surface mode as a slow
mode anymore as, in the reference frame of the observer, it possesses phase
speeds similar to that of the fast surface mode. Therefore, using the observed
phase speeds alone may introduce some difficulty when distinguishing between
the fast and slow surface kink modes in rotating photospheric flux tubes. It
may be more appropriate to differentiate between the two modes by comparing
their eigenfunctions both parallel and perpendicular to the magnetic field,
which is discussed in more detail in Section 3.1.2. In the case of a uniform
photospheric cylinder, the slow body modes tend toward $c_{Ti}$ in the long
wavelength limit (Edwin & Roberts, 1983; Priest, 2014), however with the
inclusion of a background rotational flow, these modes now encounter the
modified flow resonance point where they may become resonantly damped.
Furthermore, the dispersive nature of the fast surface kink mode is also
modified by the presence of a background rotational flow. In a uniform
photospheric flux tube, the phase speed of the fast surface mode tends to the
kink speed in the long wavelength limit (Edwin & Roberts, 1983). However,
similar to the slow modes, the fast surface mode encounters the modified flow
resonance in the long-wavelength limit, where the phase speed of this mode
also appears to increase sharply depending on the amplitude of the flow. The
fast mode may also enter the leaky regime in the long-wavelength limit under
the assumed equilibrium configuration, which may be important in the context
of solar observations.

Of course, the kink mode possesses an azimuthal wavenumber which can be either
positive or negative. Granular buffeting in the lower solar photosphere due to
convective motions beneath the solar surface, may excite both $m=1$ and $m=-1$
kink modes. In the uniform cylinder model of a solar waveguide, the
‘traditional’ kink mode is considered to display signatures that resemble a
periodic transverse displacement of the waveguide. This is because, due to the
symmetry of the model, the opposite rates of rotation set up a standing wave
in the azimuthal direction. However, introducing a background rotational flow
breaks this azimuthal symmetry and, as a result, the $m=1$ and $m=-1$ modes
are expected to behave differently and the resulting observed mode will no
longer be a standing wave in the azimuthal direction. In Figure 3, we show the
obtained wave solutions for the $m=1$ and $m=-1$ modes in a photospheric flux
tube with a linear background rotational flow of amplitude $A=0.1$. For the
case when the kink wave is rotating with the flow, the phase speed of the fast
surface kink mode increases in the long-wavelength limit and enters the leaky
regime for small $ka$. This suggests that the fast kink surface mode, when
propagating in a thin waveguide with the background rotational flow, should
not be seen in observations of rotating photospheric structures, for example
in magnetic bright points or flux bundles rooted within intergranular lanes.
However, the obtained solutions for the fast surface mode are different for
the two cases when $m=1$ and $m=-1$. This is due to the $m=1$ mode rotating in
the same direction as the background flow, in our specific configuration.
Therefore, the $m=1$ mode constructively interferes with the background
rotational flow, increasing its phase velocity, whereas the $m=-1$ mode,
rotating against the flow, destructively superimposes with the flow, hence
lowering its phase velocity. The most notable difference between the case when
$m=1$ and $m=-1$ can be seen by the dispersive behaviour of the slow surface
and body modes. In both cases, the slow modes follow the trajectory of the
flow resonance locations modified by the slow frequency as given by Equation
(30), however for the case when $m=-1$, this curve decreases with decreasing
$ka$. As a result, the difference in phase speed for a given wavenumber
between $m=1$ and $m=-1$ is greater for the slow modes, with a more dramatic
difference seen as $ka$ approaches zero. This example of the phase speed
difference between the forward propagating slow body and surface modes, for
$m=1$ and $m=-1$, demonstrates how the background flow breaks the symmetry of
the kink mode.

Figure 3: The dispersion diagrams showing the obtained solutions for the $m=1$
and $m=-1$ modes denoted by the blue curves. The green curve shows the
behaviour of Equation (30) for both cases of $m=\pm 1$, respectively. The
respective zoom in plots highlight the difference in the behaviour of the
various category of modes for different values of $m$.

####  3.1.2 Eigenfunctions

Figure 4: The resulting eigenfunctions for the slow surface kink mode ($m=1$)
for all linear cases of rotational flow with the colour scheme consistent with
Figure 1. A wavenumber value of $k=0.6$ was chosen for all plots.

Shown in Figure 4 are the spatial eigenfunctions for the slow surface kink
mode at a fixed wavenumber $ka=0.6$ for varying rotational flow amplitudes
with $m=1$. The colour scheme shown in the eigenfunctions is consistent with
that for the rotational flow profiles shown in Figure 1. It can be seen that
increasing the amplitude of the equilibrium linear rotational flow, changes
the spatial behaviour of the observable eigenfunctions. For the case of
$A=0.01$ which corresponds to a very small rotational flow parameter, the
eigenfunctions still obey a ‘surface-like’ structure, that is, the amplitude
of the radial displacement and velocity perturbations possesses a maximum at
the boundary of the flux tube and decays away from the boundary. However,
increasing the amplitude of the background rotational flow causes the radial
displacement perturbation to increase towards the centre of the flux tube,
such that the maximum displacement perturbation is no longer at the point
where $r=a$. This results in an eigenfunction that shares striking
similarities to that of the fundamental kink mode, and may therefore be
misinterpreted in observational data.

To further emphasise this point, it is possible to plot the eigenfunctions of
$\hat{P}_{T}$ and $\hat{\xi}_{r}$ for eigenvalues of a similar phase speed on
either side of the modified flow resonance point. One of these solutions
corresponds to the slow surface kink mode and the other is the fast surface
kink mode.

(a)

(b)

Figure 5: Panels showing the obtained wave solutions in a zoomed region
denoted by the dashed box in Figure 2(c). The eigenfunctions $\hat{P}_{T}$,
$\hat{\xi}_{r}$, $\hat{\xi}_{z}$ and $\hat{\xi}_{\varphi}$ for the (a) fast
and (b) slow magnetoacoustic surface kink mode solutions in a photospheric
flux tube with a background rotational flow given by $v_{\varphi}=0.1r$. Both
panels (a) and (b) are for the same $v_{ph}=\omega/k=1.3$ indicated by the
blue dot on the dispersion diagram. All plots are normalised such that the
external value of each eigenfunction is equal to unity at the boundary.

Shown in Figure 5 are these eigenfunctions for the slow and fast
magnetoacoustic kink modes at a similar phase speed. It can be seen that the
normalised eigenfunctions for $\hat{P}_{T}$, $\hat{\xi}_{r}$ and
$\hat{\xi}_{\varphi}$ are difficult to distinguish between the slow and the
fast surface modes. Furthermore, the modes no longer display the typical
characteristics of the surface mode anymore. In particular, the main
characteristic of a surface mode from the uniform cylinder model is that it
possesses a maximum amplitude of radial displacement perturbation at the
boundary of the waveguide, which is no longer the case when a rotational
background plasma flow is present. Investigating the nature of perturbations
parallel and perpendicular to the magnetic field may aid in distinguishing
between the slow and fast surface kink modes. Slow modes tend to propagate
(mainly) along the magnetic field lines, which in our study are straight and
vertical, so we should expect $\hat{\xi}_{z}$ to dominate for the slow mode
when compared to the fast mode, which may propagate at an angle across the
magnetic field. In Figure 5, we also show the normalised $z$ component of the
displacement perturbation for both the fast and slow surface kink modes. As
expected, $\hat{\xi}_{z}$ is dominant over $\hat{\xi}_{r}$ for the slow
surface mode (see e.g. Moreels & Van Doorsselaere, 2013), however, the same is
also true for the fast surface mode. Although, the absolute magnitude of the
normalised $z$ component is greater for the slow surface mode compared to the
fast surface modes, by roughly a factor of $2$, suggesting that the component
of displacement is still more dominant for the slow surface kink mode in a
rotating photospheric flux tube. It should be noted that, displayed in some
plots of the eigenfunctions, the amplitude of the resulting eigenfunctions
increases as we approach the center of the cylinder. This is a feature which
relates to the fact that there is a singularity in the set of Equations (16)
and (19) at $r=0$. However, this is purely an artifact of the plotting
technique and does not affect the eigensolver obtaining the solutions.

####  3.1.3 2D and 3D velocity fields

Following the analysis from Section 2, the wave solutions are put proportional
to $\text{exp}\left[i\left(m\varphi+kz-\omega t\right)\right]$, therefore, it
is possible to convert the one dimensional radial eigenfunctions, for example
those shown in Figures 4 and 5, into two and three dimensional plots, to
better represent an observer’s perspective.

Figure 6: Snapshots of the $2$D velocity field for the scenario of a slow body
mode in a photospheric flux in the presence of a background flow given by
$v_{\varphi}=0.1r$. The plots are arranged in a $3\times 2$ configuration
where the left hand column shows the velocity perturbation only, whereas the
right hand column shows the perturbed velocity field plus the background
velocity field. The top row corresponds to the solution for the $m=1$ mode,
the middle row shows the solution for the $m=-1$ mode and the bottom row shows
the resulting velocity field and total pressure perturbation for the sum of
the $m=1$ and $m=-1$ modes. In all panels the velocity vectors are normalised
by their maximum values. The colour contour denotes the normalised total
pressure perturbation, $\hat{P}_{T}$, which is the same for both the left and
right columns. The boundary of the flux tube is highlighted by the solid blue
line. An animated movie of this Figure is available online.

In Figure 6, we show a snapshot of the two dimensional (converted into
Cartesian $x$ and $y$) velocity field for the perturbations alongside the
addition of the background flow onto the perturbations. This snapshot is
chosen specifically at the time when the $m=1$ and the $m=-1$ perturbations
effectively cancel each other such that the sum of the perturbed velocity
field is zero. Of course, due to the breaking of symmetry from the presence of
the background rotational flow, the boundary of the flux tube is no longer at
equilibrium, as the $m=1$ and $m=-1$ modes are slightly out of phase with one
another. In this snapshot, the perturbed total pressure, displayed by the
colour contour, can be seen to display a swirling characteristic, albeit with
a small magnitude, and is clearly visible in the bottom panels of Figure 6.
This swirling behaviour is present purely in the perturbations, signifying the
effect that a rotating waveguide has on the perturbed eigenfunctions. This
swirling behaviour of the total pressure perturbation may represent an
observational signature of the slow body kink mode propagating in a rotating
photospheric flux tube, although detecting this may be a challenge for
observers due to the small absolute value of magnitude. In the plots in the
right hand column of Figure 6, we add the background velocity field to the
perturbed velocity field. The background velocity field is given by
$v_{\varphi}=0.1r$ and acts counter-clockwise in the $xy$ plane. The addition
of the background velocity field completely changes the observed distribution
of the velocity field for the slow body kink mode. In the case of a uniform
cylinder, the velocity field for the slow body kink mode can be seen emanating
from two ‘islands’, which correspond to two anti-nodes in the total pressure
eigenfunction. This can be seen in the left hand column showing the
perturbations only for the case of $m=1$ and $m=-1$, however, when the
background flow is added, the velocity field no longer displays these typical
characteristics. It should be noted that both the perturbed velocity field and
the background velocity field are normalised separately to their respective
maximum values to aid these visualisations, as it is to be expected that the
perturbed components will be significantly smaller in absolute value than the
background quantities.

Figure 7: Same as Figure 6 but a snapshot at a later time. An animated movie
of this Figure is available online.

In Figure 7, we show the same quantities as those displayed in Figure 6 at a
later time. As we have shown in Figure 3, there is a significant difference in
the phase speed for the $m=1$ and $m=-1$ forward propagating slow body waves.
Therefore, we would expect these perturbations to become increasingly out of
phase with one another as time progresses, a behaviour which is displayed in
Figure 7. The two modes becoming out of phase with one another can be seen by
comparing the top row and the middle rows of Figure 7 (which display the $m=1$
and $m=-1$ modes respectively). For the case of a uniform cylinder with no
background rotational flow, we would expect that the $m=1$ and $m=-1$ modes
would rotate in opposite directions but perfectly in phase with one another,
however, this is no longer true when a rotational background flow is added. As
a result of one mode rotating with the background flow (in our case this is
the $m=1$ mode), and the other mode propagating against the background flow
($m=-1$ mode), the sum of the two modes results in a total pressure
perturbation which is also propagating (rotating) in the azimuthal direction.
Furthermore, the boundary of the flux tube can be seen to rotate (see
associated online movie). The observed rotating motion when combining the
individual $m=1$ and $m=-1$ wave modes presents another observational
signature that can be sought to identify the kink mode in the presence of a
background rotational flow.

(a)

(b)

Figure 8: This figure shows the three dimensional structure of the slow body
kink mode for (a) a uniform magnetic flux tube without the presence of a
background flow and (b) a magnetic flux tube in the presence of a background
rotational flow, by visualising the normalised total pressure perturbation
$\hat{P}_{T}$. On the left panels, the magnetic flux tube is immersed in the
volume rendering of $\hat{P}_{T}$. The three cross-sectional cuts (shown as
coloured rings at three different heights, $z=0.0,2.5$ and $5.0$) correspond
to the corresponding right subplots. These three subplots show the LIC
visualisation at the same heights. The white rings represent the boundary of
the flux tube. An animated version of this Figure is available online.

Figure 8 shows a three dimensional visualisation of how the sum of the $m=\pm
1$ slow body kink mode manifests itself in a uniform magnetic flux tube and
also for the case of a magnetic flux tube in the presence of a background
rotational flow. Figure 8 also displays a Line Integral Contour (LIC) plot at
three different heights in the flux tube which are denoted by white rings in
the 3D plot. A LIC visualisation aids imaging the velocity field of the plasma
(see, e.g. Cabral & Leedom, 1993). In other words, LIC visualisations display
the path lines that an object would follow if it were placed into the fluid
with a stationary velocity field. It can be seen that the LIC for the slow
body kink mode in a uniform magnetic flux tube (Figure 8(a)) remains unchanged
with height, as the wave propagates and the structure oscillates in one plane,
in other words it is linearly polarised. In addition, the LIC also displays a
‘double X-point’ close to $x=0$ and $y=\pm 1$, near the boundary in the plane
in which the structure does not oscillate. This behaviour is due to the
rotational motions that are inherently associated with the non-axisymmetric
kink mode (Goossens et al., 2014). However, once a rotational flow is present,
the observed evolution of the flux tube motion is affected. Firstly, we can
observe that the three dimensional behaviour of the total pressure (density)
perturbation begins to represent a helical structure, in contrast to the
sinusoidal oscillation in the uniform case, as the perturbation now also
displays a rotational behaviour. This suggests that the resulting kink wave is
circularly polarised, as a result of the $\pm 1$ modes becoming out of phase
with one another. Therefore, we show that, in the presence of a background
rotational flow, the kink mode may manifest itself as a circularly polarised
mode, similar to a twisted magnetic flux tube (Terradas & Goossens, 2012;
Ruderman & Terradas, 2015).

Furthermore, the LIC visualisations show that the velocity field evolution is
modified in the presence of a background rotational flow. Firstly, the ‘double
X-point’ feature seen in the LIC for the uniform flux tube no longer displays
the same characteristics in the scenario including a background rotational
flow. Instead, the ‘X-points’ can be seen at different locations at certain
times as the wave propagates upwards. Secondly, the evolution of the velocity
field is no longer constant and does not appear to oscillate in a plane on a
single axis anymore. Similar to the behaviour seen in Figures 6 and 7, the LIC
visualisation shows that the velocity field also rotates as it evolves, and
the traditional signatures of the kink mode in a uniform flux tube are no
longer present. These features could possibly be retrieved from numerical
simulations of kink waves in rotating structures.

##  4 Conclusions

In the present study, we have extended previous studies investigating the
properties of MHD waves in photospheric waveguides, by introducing a linear
azimuthal component to the background velocity field. In order to conduct this
investigation, we utilised a previously developed numerical eigensolver
(Skirvin et al., 2021; Skirvin et al., 2022) to model a rotating flux tube in
a photospheric environment with a background $v_{\varphi}$ component. For the
inclusion of a linear rotational flow, very similar results to those of the
linear magnetic twist were recovered (Erdélyi & Fedun, 2010). We find that the
obtained kink mode solutions possess phase speeds along the magnetic field
which tend to infinite values in the long wavelength limit where they may
become leaky. We have shown this to be a result of the modes being guided by
the trajectory of the flow resonance location modified by the slow frequency.
As a result of considering a linear background flow, the equations describing
the continua regions reduce to single point locations. Therefore, there
becomes a point where, under photospheric conditions, the slow surface kink
mode and the slow body kink mode approach the same solution at the resonant
point and undergo an avoided crossing where their properties become mixed.
Comparison of the eigenfunctions for the slow surface kink mode and fast
surface kink mode, in the presence of a linear background rotational flow at
similar phase speeds, indicated that identification of the two modes becomes
extremely difficult in the long wavelength limit using the total pressure and
radial displacement perturbations. However, comparison of $\hat{\xi}_{z}$
shows that this component is still dominant for slow modes, comparable to
uniform theory. Furthermore, we find that the axisymmetric $m=0$ sausage mode
remains unaffected by any background azimuthal component. Analytically this
can be understood by examining the governing set of Equations (6)-(15) when
setting $m=0$ and noticing many terms either simplifying or disappearing
altogether. In addition, it suggests that sausage mode observations in the
lower solar atmosphere in e.g. pores and sunspots, may not be a suitable wave
mode to conduct atmospheric-seismology, if the structure is in the presence of
any magnetic twist or background rotational flows.

We have also presented 2D plots showing the velocity field with and without
the inclusion of the background flow (Figures 6-7) and 3D plots showing the
perturbation of the normalised total pressure $\hat{P}_{T}$ (Figure 8). These
plots clearly show the breaking of symmetry between the $m=1$ and $m=-1$ kink
modes when a flux tube is in the presence of a rotational background flow.
Furthermore, inspection of the total pressure perturbations shows that it
displays signatures of a swirling motion when the flux tube experiences its
maximum displacement. It is hoped that these visualisations will aid future
interpretations of observational data from high resolution instruments on
state of the art telescopes such as DKIST. Finally, we have also produced
three dimensional visualisations of the kink mode propagating in a
photospheric flux tube with a linear background rotational flow incorporated.
We have also shown, in Figure 8, the LIC visualisations at different heights
displaying the evolution of the velocity field as the wave propagates along
the magnetic field, which may be retrieved from numerical simulations.

In this study, we have suggested potential observational signatures of
magnetoacoustic kink modes in photospheric flux tubes with background
rotational flows. These signatures share striking characteristics with
previous observations of circularly polarised kink modes in the lower solar
atmosphere (Jess et al., 2017; Stangalini et al., 2017). However, naturally,
the magnitude of the background flows are much greater than those of the
perturbations arising from the waves, which then raises the important question
of how to actually observe this phenomena in the solar atmosphere. Both the
rotational motion of the waveguide combined with the swirling pattern seen in
the total pressure perturbation are observational signatures that may be
detected with current and future telescopes, however distinguishing between
the contribution of the background and perturbed components may be
challenging.

Furthermore, we have focused only on a magnetic flux tube in a photospheric
environment. Under the photospheric conditions adopted in this work, the
Alfvén continuum given by Equation (22) exists in the leaky regime at phase
speeds above the cut-off speed at $v_{ph}=c_{e}$. Therefore this resonant
region is not discussed in detail in this work, however, it will become
important for studies under coronal conditions, where contributions from both
Equations (22) and (23) become important. In the solar atmosphere, vortices
can be buffeted from multiple directions, for example convective motion due to
convection cells. In this study, we have considered just one fixed driver that
oscillates from side to side, however in reality the picture will be more
complicated due to drivers in additional directions. These perturbations along
different axis will manifest themselves as apparent rotation, regardless of a
background rotational flow, this is work to be investigated in the future.
Additional future work includes a closer inspection on the effect that
rotational flows have on the sausage mode. We have shown that the sausage mode
may be affected by the presence of background rotational flow in magnetic flux
tubes, however, more insight is required to quantify exactly how the sausage
mode is affected and how it manifests itself in rotating flux tubes. Finally,
future work to be investigated includes modelling a background rotational flow
with a realistic radial profile. For example Silva et al. (2020) have shown,
using information from MURAM simulations, that the radial profile of the
azimuthal velocity component in solar vortex tubes can be accurately modelled
with a cubic polynomial. In this case, the background rotational flow becomes
non-linear, such that the value of parameter $\alpha$ in Equation (24) no
longer equals $1$. We expect that in this example, the MHD spectrum is densely
occupied by the (modified) flow continua, consequently, trapped modes may find
it difficult to exist. Obtaining the eigenvalues that lie inside the continua
can be achieved by modifying the numerical eigensolver employed in this work
or by using an alternative 1D eigensolver such as Legolas (Claes et al.,
2020). Modelling a non-linear background rotational flow such as this is an
objective for future studies.

## Acknowledgements

This work has been supported by STFC (UK). SJS is grateful to STFC for the PhD
studentship project reference (2135820). VF, GV, and SSAS are grateful to
Science and Technology Facilities Council (STFC) grant ST/V000977/1, and The
Royal Society, International Exchanges Scheme, collaboration with Brazil
(IES191114) and Chile (IE170301). TVD was supported by the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant agreement No 724326) and the C1 grant TRACEspace of Internal
Funds KU Leuven. The results received support from the FWO senior research
project with number G088021N. SJS and VF would like to thank the International
Space Science Institute (ISSI) in Bern, Switzerland, for the hospitality
provided to the members of the team on ‘The Nature and Physics of Vortex Flows
in Solar Plasmas’. SJS and GV wish to acknowledge scientific discussions with
the Waves in the Lower Solar Atmosphere (WaLSA; https://www.WaLSA.team) team,
which is supported by the Research Council of Norway (project no. 262622).
This research has also received financial support from the ISEE, International
Joint Research Program (Nagoya University, Japan) and the European Union’s
Horizon 2020 research and innovation program under grant agreement No. 824135
(SOLARNET).

Software: Numpy (Harris et al., 2020), matplotlib (Hunter, 2007), SciPy
(Virtanen et al., 2020). The numerical code, Sheffield Dispersion Diagram Code
(SDDC) introduced and applied in this work is available on the Plasma Dynamics
Group (PDG) website 111https://sites.google.com/sheffield.ac.uk/pdg/solar-
codes?authuser=0 along with the user manual which explains some cases shown in
this work. This code and the accompanying tools have been developed using
Python an open-source and community-developed programming language.

## Data Availability

Accessibility of data used in this research is available upon request from the
authors.

## References

  * Abdelatif (1990) Abdelatif T. E., 1990, Sol. Phys., 129, 201
  * Aldhafeeri et al. (2021) Aldhafeeri A. A., Verth G., Brevis W., Jess D. B., McMurdo M., Fedun V., 2021, ApJ, 912, 50
  * Allcock & Erdélyi (2017) Allcock M., Erdélyi R., 2017, Sol. Phys., 292, 35
  * Andries et al. (2000) Andries J., Tirry W. J., Goossens M., 2000, ApJ, 531, 561
  * Aschwanden et al. (1999) Aschwanden M. J., Fletcher L., Schrijver C. J., Alexander D., 1999, ApJ, 520, 880
  * Bate et al. (2022) Bate W., et al., 2022, ApJ, 930, 129
  * Bonet et al. (2008) Bonet J. A., Márquez I., Sánchez Almeida J., Cabello I., Domingo V., 2008, ApJ, 687, L131
  * Bonet et al. (2010) Bonet J. A., et al., 2010, ApJ, 723, L139
  * Cabral & Leedom (1993) Cabral B., Leedom L. C., 1993, in Proceedings of the 20th Annual Conference on Computer Graphics and Interactive Techniques. SIGGRAPH ’93. Association for Computing Machinery, New York, NY, USA, p. 263–270, doi:10.1145/166117.166151, https://doi.org/10.1145/166117.166151
  * Cheremnykh et al. (2018) Cheremnykh O., Fedun V., Ladikov-Roev Y., Verth G., 2018, ApJ, 866, 86
  * Claes et al. (2020) Claes N., De Jonghe J., Keppens R., 2020, ApJS, 251, 25
  * De Groof & Goossens (2000) De Groof A., Goossens M., 2000, A&A, 356, 724
  * De Pontieu et al. (2004) De Pontieu B., Erdélyi R., James S. P., 2004, Nature, 430, 536
  * De Pontieu et al. (2007) De Pontieu B., et al., 2007, Science, 318, 1574
  * Edwin & Roberts (1983) Edwin P. M., Roberts B., 1983, Sol. Phys., 88, 179
  * Erdélyi & Fedun (2007) Erdélyi R., Fedun V., 2007, Sol. Phys., 246, 101
  * Erdélyi & Fedun (2010) Erdélyi R., Fedun V., 2010, Sol. Phys., 263, 63
  * Fedun et al. (2011a) Fedun V., Shelyag S., Verth G., Mathioudakis M., Erdélyi R., 2011a, Annales Geophysicae, 29, 1029
  * Fedun et al. (2011b) Fedun V., Shelyag S., Erdélyi R., 2011b, ApJ, 727, 17
  * Geeraerts et al. (2022) Geeraerts M., Vanmechelen P., Van Doorsselaere T., Soler R., 2022, arXiv e-prints, p. arXiv:2203.02359
  * Goedbloed & Poedts (2004) Goedbloed J. P. H., Poedts S., 2004, Principles of Magnetohydrodynamics
  * González-Avilés et al. (2017) González-Avilés J. J., Guzmán F. S., Fedun V., 2017, ApJ, 836, 24
  * González-Avilés et al. (2018) González-Avilés J. J., Guzmán F. S., Fedun V., Verth G., Shelyag S., Regnier S., 2018, ApJ, 856, 176
  * Goossens et al. (1992) Goossens M., Hollweg J. V., Sakurai T., 1992, Sol. Phys., 138, 233
  * Goossens et al. (2011) Goossens M., Erdélyi R., Ruderman M. S., 2011, Space Sci. Rev., 158, 289
  * Goossens et al. (2014) Goossens M., Soler R., Terradas J., Van Doorsselaere T., Verth G., 2014, ApJ, 788, 9
  * Grant et al. (2015) Grant S. D. T., et al., 2015, ApJ, 806, 132
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hunter (2007) Hunter J. D., 2007, Computing in Science and Engineering, 9, 90
  * Jess et al. (2017) Jess D. B., et al., 2017, ApJ, 842, 59
  * Keys et al. (2018) Keys P. H., et al., 2018, ApJ, 857, 28
  * Magyar et al. (2022) Magyar N., Duckenfield T., Van Doorsselaere T., Nakariakov V. M., 2022, A&A, 659, A73
  * Mather & Erdélyi (2016) Mather J. F., Erdélyi R., 2016, ApJ, 822, 116
  * Moreels & Van Doorsselaere (2013) Moreels M. G., Van Doorsselaere T., 2013, A&A, 551, A137
  * Morton et al. (2015) Morton R. J., Tomczyk S., Pinto R., 2015, Nature Communications, 6, 7813
  * Nakariakov & Roberts (1995) Nakariakov V. M., Roberts B., 1995, Sol. Phys., 159, 213
  * Nakariakov et al. (1999) Nakariakov V. M., Ofman L., Deluca E. E., Roberts B., Davila J. M., 1999, Science, 285, 862
  * Pinter et al. (1998) Pinter B., Cadez V. M., Goossens M., 1998, A&A, 332, 775
  * Priest (2014) Priest E. R., 2014, Magnetohydrodynamics of the Sun
  * Rouppe van der Voort et al. (2007) Rouppe van der Voort L. H. M., De Pontieu B., Hansteen V. H., Carlsson M., van Noort M., 2007, ApJ, 660, L169
  * Ruderman & Petrukhin (2022) Ruderman M. S., Petrukhin N. S., 2022, Sol. Phys., 297, 72
  * Ruderman & Terradas (2015) Ruderman M. S., Terradas J., 2015, A&A, 580, A57
  * Sakurai et al. (1991) Sakurai T., Goossens M., Hollweg J. V., 1991, Sol. Phys., 133, 227
  * Scullion et al. (2011) Scullion E., Erdélyi R., Fedun V., Doyle J. G., 2011, ApJ, 743, 14
  * Shelyag et al. (2011) Shelyag S., Fedun V., Keenan F. P., Erdélyi R., Mathioudakis M., 2011, Annales Geophysicae, 29, 883
  * Shelyag et al. (2012) Shelyag S., Mathioudakis M., Keenan F. P., 2012, ApJ, 753, L22
  * Shelyag et al. (2013) Shelyag S., Cally P. S., Reid A., Mathioudakis M., 2013, ApJ, 776, L4
  * Shetye et al. (2019) Shetye J., Verwichte E., Stangalini M., Judge P. G., Doyle J. G., Arber T., Scullion E., Wedemeyer S., 2019, ApJ, 881, 83
  * Silva et al. (2020) Silva S. S. A., Fedun V., Verth G., Rempel E. L., Shelyag S., 2020, ApJ, 898, 137
  * Silva et al. (2021) Silva S. S. A., Verth G., Rempel E. L., Shelyag S., Schiavo L. A. C. A., Fedun V., 2021, ApJ, 915, 24
  * Skirvin et al. (2021) Skirvin S. J., Fedun V., Verth G., 2021, MNRAS, 504, 4077
  * Skirvin et al. (2022) Skirvin S. J., Fedun V., Silva S. S. A., Verth G., 2022, MNRAS, 510, 2689
  * Snow et al. (2018) Snow B., Fedun V., Gent F. A., Verth G., Erdélyi R., 2018, ApJ, 857, 125
  * Soler et al. (2009) Soler R., Oliver R., Ballester J. L., 2009, ApJ, 693, 1601
  * Soler et al. (2010) Soler R., Terradas J., Oliver R., Ballester J. L., Goossens M., 2010, ApJ, 712, 875
  * Soler et al. (2011) Soler R., Terradas J., Goossens M., 2011, ApJ, 734, 80
  * Spruit & Zweibel (1979) Spruit H. C., Zweibel E. G., 1979, Sol. Phys., 62, 15
  * Stangalini et al. (2017) Stangalini M., et al., 2017, ApJ, 840, 19
  * Stangalini et al. (2022) Stangalini M., et al., 2022, Nature Communications, 13, 479
  * Taroyan & Erdélyi (2002) Taroyan Y., Erdélyi R., 2002, Physics of Plasmas, 9, 3121
  * Taroyan & Erdélyi (2003) Taroyan Y., Erdélyi R., 2003, Physics of Plasmas, 10, 266
  * Terra-Homem et al. (2003) Terra-Homem M., Erdélyi R., Ballai I., 2003, Sol. Phys., 217, 199
  * Terradas & Goossens (2012) Terradas J., Goossens M., 2012, A&A, 548, A112
  * Tirry & Goossens (1996) Tirry W. J., Goossens M., 1996, ApJ, 471, 501
  * Tomczyk et al. (2007) Tomczyk S., McIntosh S. W., Keil S. L., Judge P. G., Schad T., Seeley D. H., Edmondson J., 2007, Science, 317, 1192
  * Tziotziou et al. (2018) Tziotziou K., Tsiropoula G., Kontogiannis I., Scullion E., Doyle J. G., 2018, A&A, 618, A51
  * Van Doorsselaere et al. (2004) Van Doorsselaere T., Debosscher A., Andries J., Poedts S., 2004, A&A, 424, 1065
  * Verth et al. (2007) Verth G., Van Doorsselaere T., Erdélyi R., Goossens M., 2007, A&A, 475, 341
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Wedemeyer-Böhm et al. (2012) Wedemeyer-Böhm S., Scullion E., Steiner O., Rouppe van der Voort L., de La Cruz Rodriguez J., Fedun V., Erdélyi R., 2012, Nature, 486, 505
  * Wilson (1979) Wilson P. R., 1979, A&A, 71, 9
  * Yadav et al. (2020) Yadav N., Cameron R. H., Solanki S. K., 2020, ApJ, 894, L17
  * Yadav et al. (2021) Yadav N., Cameron R. H., Solanki S. K., 2021, A&A, 645, A3
  * Zaqarashvili et al. (2015) Zaqarashvili T. V., Zhelyazkov I., Ofman L., 2015, ApJ, 813, 123
  * Zhelyazkov & Chandra (2019) Zhelyazkov I., Chandra R., 2019, Sol. Phys., 294, 20