© ESO 2021

1. Introduction

Mechanisms for heating the upper solar atmosphere and corona, be they magnetohydrodynamical (MHD) waves or magnetic reconnection, generally rely on the transverse motion of magnetic footpoints in the photosphere. Excitations of kink tube waves or compressive fast modes are often considered to arise from “granular buffeting”, by which photospheric magnetic flux concentrations are rapidly and stochastically moved in a sideward direction. Another option is given by turbulent motions in the convection zone, which generate a whole spectrum of transversal excitations. Also the heating by magnetic reconnection is commonly considered to be induced by transverse motions of magnetic footpoints. These can lead to the twisting and braiding of magnetic flux tubes forming dissipative electric current sheets.

Different from transverse excitations, vortical motions of the magnetic footpoints have been less in the focus of research, mainly for lack of observational support. In the past, vortical motions and torsional Alfvén waves were of rather hypothetical nature, often serving theoretical models for the driving of spicules (e.g., Stenflo 1975; Hollweg 1981; Fedun et al. 2011a). This state of affairs has drastically changed over the past decade.

In the photosphere, swirls of granular scale were detected by tracking magnetic bright points (BPs) in the intergranular space (Bonet et al. 2008; Balmaceda et al. 2010; Manso Sainz et al. 2011) and by the means of local correlation tracking (LCT) (Bonet et al. 2010; Vargas Domínguez et al. 2011; Requerey et al. 2017, 2018). In the chromosphere, vortical motions with a typical size of about 2 arcsec were observed by Wedemeyer-Böhm & Rouppe van der Voort (2009) and Wedemeyer-Böhm et al. (2012) using Ca II 854 nm narrowband filtergrams obtained with the CRisp Imaging SpectroPolarimeter (CRISP) instrument of the Swedish 1-m Solar Telescope (SST). Moreover, these chromospheric swirls were colocated with magnetic BPs in the photosphere, suggesting a magnetic origin. Higher up, Wedemeyer-Böhm et al. (2012) detected imprints of chromospheric swirls in the transition region and the corona in the spectral lines of He II 30.4 nm and Fe IX 17.1 nm, respectively, using recordings of the Atmospheric Imaging Assembly (AIA) instrument of NASA’s space-based Solar Dynamics Observatory (SDO). These observations suggested that chromospheric swirls are the observational signatures of rotating coherent magnetic structures, which reach from the convection zone to the outer layers of the atmosphere.

Morton et al. (2013) found the chromospheric counterpart of photospheric swirls in the form of a quasi-periodic torsional motion, in time series of Hα narrowband filtergrams obtained with the ROSA (Rapid Oscillations in the Solar Atmosphere) instrument attached to the Dunn Solar Telescope. Park et al. (2016) detected swirls consisting of spiral arms in Hα filtergrams taken with CRISP at SST, which coexisted with strong upflows in the upper chromosphere as seen from Mg II k line Dopplergrams obtained with the space-based Interface Region Imaging Spectrograph (IRIS). More recently, a 1.7 h persistent vortex flow of 6 arcsec diameter was observed by Tziotziou et al. (2018, 2019) with CRISP at SST in the cores of the Hα and Ca II 854.2 nm lines, while no BPs were observed in the line wings of Hα and Ca II. Finally, Shetye et al. (2019) found, from spectral imaging in the lines of Hα and Ca II 854 nm along with polarimetry in Fe I 630.2 nm, that the rotation of magnetic flux concentrations in the photosphere matches the chromospheric swirl. However, they reported that they could not determine whether a swirl is a gradual response to the photospheric motion or an actual propagating Alfvénic wave.

Alfvén waves were theoretically predicted by Alfvén (1942), but their detection in the Sun, which provides preferential conditions for their existence, proved to be a difficult task. Nevertheless, over the past two decades, multiple observations have revealed their presence and propagation in the solar atmosphere and corona (e.g., Tomczyk et al. 2007; Okamoto & De Pontieu 2011; Okamoto et al. 2007). Jess et al. (2009) found signatures of propagating torsional Alfvén waves in the photosphere in oscillatory phenomena associated with a conglomeration of magnetic BPs. De Pontieu et al. (2012) reported rotational motion in spicules, which they interpreted to be torsional Alfvén waves, and Srivastava et al. (2017) detected ubiquitous high frequency (∼12 − 42 mHz) torsional motions in thin spicular-type structures in the chromosphere that resemble torsional Alfvén waves. Liu et al. (2019a) called attention to ubiquitous torsional Alfvén pulses by correlating photospheric and chromospheric swirls. Moreover, Grant et al. (2018) provided observational evidence of Alfvén waves heating chromospheric plasma in a sunspot umbra through the formation of shocks.

Vortex flows and their connection with Alfvén waves have also been extensively studied with MHD numerical simulations. Small-scale swirls appear regularly in simulations of the solar surface convection and atmosphere (see, e.g., Stein & Nordlund 1998; Vögler 2004; Shelyag et al. 2011; Moll et al. 2011, 2012; Steiner & Rezaei 2012; Kitiashvili et al. 2013; Calvo et al. 2016; Liu et al. 2019b; Yadav et al. 2020). Shelyag et al. (2013) identified apparent vortex-like motions in magnetic flux tubes of the photosphere as torsional Alfvén waves by analyzing time–distance diagrams and plotting particle tracks. Wedemeyer & Steiner (2014) put forward a model consisting of two vortex systems, where the upper chromospheric and photospheric swirling plasma is tightly coupled to the “frozen-in” magnetic field lines, which have their footpoints within the lower intergranular vortex flow. This system of vortical structures was dubbed a “magnetic tornado” (Wedemeyer-Böhm et al. 2012). Fedun et al. (2011b), Vigeesh et al. (2012), and Murawski et al. (2015) studied the generation, propagation, and energy transfer of torsional Alfvén waves in modeled solar magnetic flux tubes, and Liu et al. (2019a) showed that azimuthal perturbations of the magnetic field in the upper photosphere can generate Alfvén pulses, which carry the information of the rotational motion into the chromosphere. However, the exact mechanism at the origin of the swirls in realistic numerical simulations is still unclear.

A particularly interesting aspect of these swirls is that they possibly generate torsional Alfvén waves that propagate upward and provide an efficient mechanism for chromospheric and coronal heating. From numerical simulations, Wedemeyer-Böhm et al. (2012) estimated a net vertical Poynting flux associated with swirling motions of 440 W m⁻² at the interface between the chromosphere and the corona, while Liu et al. (2019a) estimate a minimal nonthermal energy flux of 33 − 131 W m⁻² in the middle chromosphere.

In this paper, we study in detail “magnetic swirls”, that is, the relation between swirls, magnetic field perturbations, and Alfvén waves. In particular, we investigate the MHD processes that lead to the propagation of vortical motions from the photosphere to the chromosphere, where they appear as chromospheric swirls. For that purpose, we analyze realistic radiative MHD numerical simulations carried out with the CO5BOLD code (Freytag et al. 2012) by using the swirling strength criterion and its evolution equation (Zhou et al. 1999; Canivete Cuissa & Steiner 2020).

The paper is organized as follows: Sect. 2 reviews a number of theoretical concepts that are used for the analysis. Section 3 gives details on the numerical simulations and Sect. 4 presents the analysis and discusses the results. A summary and conclusions are given in Sect. 5.

2. Theoretical background

This section reviews a few theoretical concepts that are used for the analysis of the simulations. These concern properties of Alfvén waves, the computation of energy fluxes, and the concept of swirling strength. Moreover, we introduce a new quantity, the “magnetic swirling strength”, which is used for detecting twists in magnetic fields.

2.1. Alfvén waves

Alfvén waves are perturbations in the plasma with the magnetic tension acting as a restoring force. The perturbations are transverse to the propagation direction and magnetic field lines; therefore, one speaks of shear or torsional Alfvén waves. When the perturbation propagates across the magnetic field, causing compression and rarefaction of the magnetic field and plasma, one speaks of compressional Alfvén waves or fast mode waves (see, e.g., Priest 2014). For the purposes of this paper, we focus solely on the former.

We consider with Fig. 1 a static, incompressible plasma in the ideal MHD approximation, with constant density ρ and uniform and vertical magnetic field B₀ = B₀e_z. Moreover, we suppose that the magnetic field dominates the equilibrium so that the hydrodynamical pressure and gravity can be neglected. Although these conditions are not satisfied in the solar atmosphere, this abstraction still serves to describe some basic behavior of the Alfvén wave propagation in the solar photosphere and chromosphere. Introducing the magnetic field perturbation with B = B₀ + b, where |b| ≪ |B₀|, the velocity perturbation v, and linearizing the system of MHD equations using incompressibility (see, e.g., Priest 2014, Chap. 4), yields two important characteristics of Alfvén wave propagation.

Fig. 1. Torsional Alfvén wave propagating in the direction k. Perturbations in the magnetic field, b, and in the plasma velocity field, v, are normal to the plane spanned by k and B0.

The first one relates the velocity perturbation to the magnetic field perturbation

$$\begin{array}{c}\hfill \mathit{v}=-{\displaystyle \frac{\omega \mathit{b}}{\mathit{k}·{\mathit{B}}_0}}=\{\begin{array}{cc}\hfill -(v_{\mathrm{A}}/B_0)\mathit{b}& \mathrm{if}0\le \vartheta <\pi /2,\hfill \\ \hfill (v_{\mathrm{A}}/B_0)\mathit{b}& \mathrm{if}\pi /2<\vartheta \le \pi ,\hfill \end{array}\end{array}$$(1)

where ω is the angular frequency of the plane wave, k the wave-vector indicating the propagation direction, ${\mathit{v}}_{\mathrm{A}}=B_0/\sqrt{4\pi \rho }$ is the Alfvén speed, and ϑ is the angle between e_z and k. We note that, if ϑ < π/2, the perturbations v and b are anti-parallel, while they are parallel if ϑ > π/2. In the following, for simplicity, we only consider the case in which the wave travels in the direction of the magnetic field, that is, 0 ≤ ϑ < π/2. Furthermore, one can prove that both the magnetic field and the velocity perturbation are azimuthal with respect to the static magnetic field, B₀, and the wave-vector k can point in any direction as depicted in Fig. 1. This implies that the linearized version of the magnetic pressure is null because p_m = b ⋅ B₀/8π = 0. Since plasma pressure and gravity are not considered, the only driving force of torsional Alfvén waves is the magnetic tension.

The second result is the dispersion relation

$$\begin{array}{c}\hfill {\omega }^2=v_{\mathrm{A}}^2{(\mathit{k}·{\mathit{e}}_z)}^2.\end{array}$$

From this equation, it is possible to obtain the group velocity of an Alfvén wave packet, which corresponds to the velocity at which the energy is transmitted,

$$\begin{array}{c}\hfill {\mathit{v}}_{\text{group}}={\displaystyle \frac{\partial \omega }{\partial \mathit{k}}}=v_{\mathrm{A}}{\mathit{e}}_z\mathrm{if}0\le \vartheta <\pi /2.\end{array}$$(2)

It is equivalent to the Alfvén speed and energy is carried along the equilibrium magnetic field.

2.2. Energy fluxes

Alfvén waves can contribute to the energy transport in the solar atmosphere. To study this aspect, the MHD Poynting flux vector is employed,

$$\begin{array}{c}\hfill \mathit{S}={\displaystyle \frac{1}{4\pi }}\mathit{B}\times (\mathit{v}\times \mathit{B})={\displaystyle \frac{1}{4\pi }}(\mathit{B}·\mathit{B})\mathit{v}-{\displaystyle \frac{1}{4\pi }}(\mathit{B}·\mathit{v})\mathit{B}.\end{array}$$(3)

Since we are interested in the vertical energy flux, the z-component of Eq. (3) is expanded following Shelyag et al. (2012),

$$\begin{array}{c}\hfill S_z=\underset{S_z^{\mathrm{v}}}{\underbrace{{\displaystyle \frac{1}{4\pi }}{v}_z(B_x^2+B_y^2)}}-\underset{S_z^{\mathrm{h}}}{\underbrace{{\displaystyle \frac{1}{4\pi }}{B}_z(v_x{B}_x+v_y{B}_y)}},\end{array}$$(4)

where the first term, $S_z^{\text{v}}$, is related to vertical motions of horizontal magnetic field and the second term, $S_z^{\text{h}}$, corresponds to the vertical flux generated by horizontal motions of magnetized plasma.

The contribution of the swirls to the mean Poynting flux over a given field of view of area A_FOV is then given by

$$\begin{array}{c}\hfill {\overline{S}}_z=S_z{\overline{N}}_{\mathrm{s}}\frac{\pi {\overline{r}}_{\mathrm{s}}^2}{A_{\mathrm{FOV}}},\end{array}$$(5)

where S_z is the mean contribution from a single swirl, ${\overline{N}}_{\mathrm{s}}$ is the average number of swirls that can be observed at any time in the given field of view, be it of a simulation or an observation, and ${\overline{r}}_{\mathrm{s}}$ is the average swirl radius.

A swirl can also produce a mechanical energy flux, which is given by

$$\begin{array}{c}\hfill F_z^{\mathrm{r}}=\frac{1}{A_{\mathrm{s}}}{\displaystyle \int {\int }_{A_{\mathrm{s}}}}\mathrm{d}x\mathrm{d}y\frac{1}{2}\rho v_{\mathrm{r}}^2{v}_{\mathrm{p}},\end{array}$$(6)

where the integration is taken over the assumed circular area A_s associated with the swirl, ρ = ρ(x, y) is the density of the plasma, and v_r = v_r(x, y) is the rotational velocity. Furthermore, v_p = v_p(x, y) is the propagation speed of the swirl in the vertical direction. For simplicity, we assume A_s to lie in a horizontal plane and v_r to be the velocity projected into this plane. Assuming that the flow rotates around a common axis, we can express the rotational velocity in terms of the vertical component of the swirling strength vector as ${\mathit{v}}_{\mathrm{r}}=\frac{1}{2}r{\lambda }_z$ (see Sect. 2.3), hence

$$\begin{array}{c}\hfill F_z^{\lambda }=\frac{1}{8A_{\mathrm{s}}}{\displaystyle \int {\int }_{A_{\mathrm{s}}}}\mathrm{d}x\mathrm{d}y\rho r^2{\lambda }_z^2{v}_{\mathrm{p}},\end{array}$$(7)

where $r=\sqrt{x^2+{\mathit{y}}^2}$ is the distance from the center of the swirl and λ_z = λ_z(x, y) is the associated vertical component of the swirling strength vector.

The contribution of swirling motions to the mean mechanical energy flux ${\overline{F}}_z$ is then

$$\begin{array}{c}\hfill {\overline{F}}_z=F_z^{\mathrm{r},\lambda }{\overline{N}}_{\mathrm{s}}\frac{\pi {\overline{r}}_{\mathrm{s}}^2}{A_{\mathrm{FOV}}},\end{array}$$(8)

where $F_z^{\text{r},\lambda }$ is the mean mechanical flux density of a single swirl computed with either Eq. (6) or Eq. (7). We notice that S_z, F_z, ${\overline{F}}_z$ and ${\overline{S}}_z$ are energy flux densities.

2.3. Swirling strength

A vortex or swirl can be intuitively described as the rotation of fluid parcels around a common axis. Despite this simple concept, a rigorous mathematical definition is still an open issue in fluid mechanics (see, e.g., Kolář 2007; Liu et al. 2019c). The vorticity is the classical quantity to describe rotational flows and it has been adopted to investigate vortex flows in numerical simulations of the solar convection zone (see, e.g., Nordlund et al. 1997) and atmosphere (see, e.g., Shelyag et al. 2011, 2013). However, Jeong & Hussain (1995) showed that the vorticity is not a suitable tool for vortex identification since it cannot distinguish between non-rotational shear flows and actual vortices. This can lead to misidentifications in a dynamical and turbulent system such as the solar atmosphere. Therefore, we decided to employ the swirling strength criterion, which is not affected by shears (Zhou et al. 1999).

The swirling strength is based on the eigenanalysis of the velocity gradient tensor 𝒰_ij ≔ ∂_jv_i, that is, the Jacobian matrix of the local velocity field. If a vortex is present in the flow, the velocity gradient tensor can be locally diagonalized and it will exhibit one real and a pair of complex conjugated eigenvalues

$$\begin{array}{c}\hfill \mathcal{U}=\underset{\mathcal{P}}{\underbrace{[{\mathit{u}}_{\mathrm{r}},{\mathit{u}}_{+},{\mathit{u}}_{-}]}}\underset{\mathrm{\Lambda }}{\underbrace{\left[\begin{array}{ccc}{\lambda }_{\mathrm{r}}& 0& 0\\ 0& {\lambda }_{+}& 0\\ 0& 0& {\lambda }_{-}\end{array}\right]}}\underset{{\mathcal{P}}^{-1}}{\underbrace{{[{\mathit{u}}_{\mathrm{r}},{\mathit{u}}_{+},{\mathit{u}}_{-}]}^{-1}}},\end{array}$$(9)

where λ_± = λ_cr ± iλ_ci and λ_r ∈ ℝ. The swirling strength λ is then defined through the imaginary part of the complex eigenvalues, λ_ci. In this paper, we adopt the same convention as in Canivete Cuissa & Steiner (2020) and define the swirling strength as λ ≔ 2λ_ci. For a rigid body rotation, the period of revolution T of the flow is then computed as T = 4π/λ and the vorticity ω ≔ |∇ × v| equals the swirling strength λ. Furthermore, the rotation axis and the orientation of the vortex can be inferred from the normalized eigenvector associated with the real eigenvalue, u_r. Hence, it is possible to define a swirling vector, λ ≔ λu_r, which carries the necessary information to characterize the vortex in three dimensions. For a detailed review of the swirling strength, the reader can refer to Canivete Cuissa & Steiner (2020).

The swirling strength is invariant under Galilean transformations and can be used in the context of compressible (magneto-)hydrodynamics (Kolář 2009). This criterion has already been successfully applied in studies regarding vortex flows in the solar atmosphere (see Moll et al. 2011; Kato & Wedemeyer 2017; Yadav et al. 2020). Furthermore, Canivete Cuissa & Steiner (2020) derived the evolution equation for the swirling strength. A similar equation was already known for the vorticity (see, e.g., Shelyag et al. 2011) but, in this case, it can also be biased by shear flows. Therefore, the swirling equation is a valuable tool for the investigation of the physical processes responsible for the dynamics of vortices. It can be derived from a general momentum equation of (magneto-) hydrodynamics: For the specific case of ideal MHD, it can be formulated as

$$\begin{array}{cccc}\hfill \frac{\mathrm{D}\lambda }{\mathrm{D}t}=& -2\lambda {\lambda }_{\mathrm{cr}}\hfill & \hfill T^1& \hfill \\ \hfill & -2\mathrm{Im}{\left\{{\mathcal{P}}^{-1}[\mathrm{\nabla }\otimes (\frac{1}{\rho }\mathrm{\nabla }{p}_{\mathrm{g}})]\mathcal{P}\right\}}_{22}\hfill & \hfill T^2& \hfill \\ \hfill & -2\mathrm{Im}{\left\{{\mathcal{P}}^{-1}[\mathrm{\nabla }\otimes (\frac{1}{\rho }\mathrm{\nabla }{p}_{\mathrm{m}})-(\mathrm{\nabla }\frac{1}{\rho })\otimes \frac{(\mathit{B}·\mathrm{\nabla })\mathit{B}}{4\pi }]\mathcal{P}\right\}}_{22}\hfill & \hfill T^3& \hfill \\ \hfill & +2\mathrm{Im}{\left\{{\mathcal{P}}^{-1}[\frac{1}{\rho }\mathrm{\nabla }\otimes \frac{(\mathit{B}·\mathrm{\nabla })\mathit{B}}{4\pi }]\mathcal{P}\right\}}_{22}\hfill & \hfill T^4& \hfill \\ \hfill & -2\mathrm{Im}{\{{\mathcal{P}}^{-1}[\mathrm{\nabla }\otimes (\mathrm{\nabla }\mathrm{\Phi })]\mathcal{P}\}}_{22},\hfill & \hfill T^5\end{array}$$(10)

where λ_cr is the real component of the complex eigenvalues of 𝒰, 𝒫 and 𝒫⁻¹ are, respectively, the matrix composed of the eigenvectors of 𝒰 and its inverse as shown in Eq. (9), p_g is the atmospheric gas pressure, p_m = B²/8π represents the pressure owing to the magnetic field, and Φ corresponds to a potential of conservative forces. We notice that the terms between curly brackets are matrices since the symbol ⊗ denotes the tensor product between two vectors. We are only interested in the (2, 2) component of the resulting matricial multiplication. Finally, the material derivative on the left-hand side of Eq. (10) is defined as

$$\begin{array}{c}\hfill \frac{\mathrm{D}\lambda }{\mathrm{D}t}=\frac{\partial \lambda }{\partial t}+(\mathit{v}·\mathrm{\nabla })\lambda ,\end{array}$$

where ∂λ/∂t defines the local production of swirling strength, while (v⋅∇)λ accounts for the swirling strength advected with the plasma.

Following Canivete Cuissa & Steiner (2020), we give a physical interpretation for each one of the terms appearing in Eq. (10). The first term, T¹, is a stretching term, while T² and T³ are related to the baroclinic effects: The former is a pure hydrodynamical process and the latter is due to magnetic fields. Then, T⁴ is associated with magnetic tension effects and T⁵ describes the generation of swirling strength by conservative forces.

For the statistical analyses, plots in Sect. 4, and the appendices, a swirling strength threshold is applied to ease the interpretation and reduce noise (Zhou et al. 1999). It is also motivated by the request that swirls perform a significant fraction of a full rotation over the time span of the simulation. Therefore, swirls having a period longer than 10 min, that is, λ < 2.09 × 10⁻² Hz, are discarded.

Fig. 2. Time instant of the swirling strength λ in a slice of dimension 9.6 × 2.6 × 2.8 Mm3 of the whole physical domain of the numerical simulation. The red contour indicates the surface of optical depth τ500 = 1. The yellow contour corresponds to the isosurface of plasma-β = 1. Swirls with a period larger than 10 min (λ < 2.09 × 10−2 Hz) are not shown and the values λ > 1 Hz are saturated. This rendering has been produced with ParaView (Ahrens et al. 2005).

2.4. Magnetic swirling strength

Torsional Alfvén waves are characterized by azimuthal perturbations in the vertically directed magnetic field lines. In order to identify this torsion in numerical simulations, we define the criterion of magnetic swirling strength, λ^B. This quantity is derived in the same way as the swirling strength, but the matrix to be diagonalized is the magnetic gradient tensor, ℬ_ij ≔ ∂_jB_i. Hence, the magnetic swirling strength is defined through the imaginary component of the complex eigenvalues of ℬ, that is, ${\lambda }^B:=2{\lambda }_{\mathrm{ci}}^B$. Following the same reasoning as for the swirling strength, the magnetic swirling vector ${\mathit{\lambda }}^B={\lambda }^B{\mathit{u}}_{\mathrm{r}}^B$ can be defined, which describes the three-dimensional twisting of the magnetic field. It is worth noticing that the units of the magnetic swirling strength are [G cm⁻¹] and that therefore this quantity does not describe a rotational flow of the magnetic field lines, but rather their twisting around the axis represented by the real eigenvector, ${\mathit{u}}_{\mathrm{r}}^B$. Thus, a magnetic field line rotating together with the plasma flow in a rigid body fashion will produce no magnetic swirling strength, while a magnetic field line with an helical structure will present a finite value of magnetic swirling strength. Consequently, this criterion is well suited to detect azimuthal torsions in the magnetic field and, therefore, imprints of torsional Alfvén waves.

For torsional Alfvén waves in a static, incompressible plasma dominated by the magnetic field, we can find a simple relation between the swirling strength λ and the magnetic swirling strength λ^B. By taking the tensor product between ∇ and Eq. (1) it is possible to obtain

$$\begin{array}{c}\hfill \mathcal{U}=-{\displaystyle \frac{1}{\sqrt{4\pi \rho }}}\mathcal{B}\mathrm{if}0\le \vartheta <\pi /2,\end{array}$$(11)

which can be diagonalized on both sides since $1/\sqrt{4\pi \rho }$ is a real scalar. Given the definitions of λ and λ^B, we get

$$\begin{array}{c}\hfill \lambda =-{\displaystyle \frac{1}{\sqrt{4\pi \rho }}}{\lambda }^B\mathrm{if}0\le \vartheta <\pi /2,\end{array}$$(12)

which simply states that, for a torsional Alfvén wave propagating in direction of the magnetic field (ϑ < π/2), the swirling strength and the magnetic swirling strength are proportional to each other and have opposite sign (see also Fig. 1).

3. Numerical simulation

The high-resolution and high-cadence simulations analyzed in this study have been carried out with the radiative MHD code CO5BOLD (Freytag et al. 2012), which solves the coupled system of compressible MHD equations in an external gravitational field and includes nonlocal, frequency-dependent radiative transfer.

The simulations are performed on a three-dimensional Cartesian grid of size 9.6 × 9.6 × 2.8 Mm³, with a grid cell size of 10 km in each spatial dimension. The vertical component ranges from about 1300 km below the optical surface τ₅₀₀ = 1 to 1500 km above it. Hence, the simulation domain encompasses a small volume near the solar surface, ranging from the surface layers of the convection zone to the middle chromosphere. The gravitational field is uniform and vertical with a value of log(g) = 4.44.

The initial condition of the simulation is given by a previously relaxed hydrodynamical model, to which a uniform and vertical magnetic field of 50 G was superimposed. The system is then advanced with the MHD module of CO5BOLD (Schaffenberger et al. 2005) until relaxation. The lateral boundary conditions are periodic for both plasma and magnetic fields, while the top and bottom ones are open under the condition that the net mass flux at the bottom is zero. An entropy inflow is supplied to maintain an average surface effective temperature of T_eff = 5770 K. The magnetic field is constrained to be vertical at the top and bottom boundaries. More details on the simulation setup are given in Calvo (2018, Sect. 2), in particular run d3gt57g44v50fc.

For the present study, a time series of 441 three-dimensional data cubes with a cadence of 2 s, amounting to about 5 TB of data, is analyzed. It starts from t = 5520 s of d3gt57g44v50fc. This relatively high-cadence has proven to be necessary to be able to capture the evolution and the dynamics of the detected vortices. The series spans a period of approximately 15 min.

Fig. 3. Statistical properties of swirling motions and mean physical quantities as a function of height z. (a) Spearman’s correlation coefficient r between the swirling strength λ and the magnetic field strength |B|, r(λ, |B|) in blue, and between the swirling strength λ and the magnetic swirling strength λB, r(λ, λB) in red. The green curve represents r(λ, |B|) calculated in regions in which β < 1 and above z = 0 only. The curves are means over eleven different time instants with a cadence of one minute and the shaded areas correspond to the standard deviations obtained from the temporal variations. Data points with λ < 2.09 × 10−2 Hz are excluded. (b) Mean magnetic field strength ⟨|B|⟩ (red) and mean swirling strength ⟨λ⟩. The dark blue and red curves refer to the same MHD time instants as in (a), while the light blue curve refers to 14 time instants of a purely hydrodynamical simulation.

4. Results and discussion

This section starts by taking a global view on the simulation, evidencing the correlation between vortices and magnetic field. It then concentrates on a single swirl event, which is shown to be a propagating torsional Alfvén pulse, and concludes with an estimate of the energy carried by the pulse and by swirling motions in the entire simulation domain. Both global and local analyses are supplemented by appendices.

4.1. Relation between swirls and magnetic fields

Figure 2 shows a volume rendering of the distribution of swirling strength seen across a vertical section through the simulation domain at an arbitrary location and time instant. Three regions can be distinguished. The first one extends from the bottom of the box to z ≈ 0 km, the height that roughly corresponds to the surface of optical depth τ₅₀₀ = 1 indicated in red color. This region coincides with the convection zone and it is characterized by an almost isotropic distribution of swirls, mostly produced by the convective and turbulent motions of the plasma.

At the photospheric and chromospheric level, there exist two distinct regions: The first one shows tubes of intense swirling strength aligned with the magnetic field, while the other one is characterized by a severe depletion of swirls. The yellow, funnel-shaped surface, which sharply separates these two regions, is given by the isosurface of plasma-β = 1, where β = p_g/p_m is the ratio between gas pressure p_g and magnetic pressure p_m. This clear separation shows evidence of strong coupling between swirling motions and magnetic fields. Above the yellow surface and within the funnels β < 1, therefore, the dynamics are governed by the magnetic fields, while the region of depleted swirling strength is characterized by β > 1, meaning that the gas pressure dominates.

We notice, however, that the region in between the magnetic flux concentrations in the photosphere is not completely void of swirls. There exist arches of swirling strength that are low lying and arches that reach farther out, as well as vertically directed swirls, which rise above the solar surface as it can be seen for example in the ranges 3000 km < x < 5000 km and 7000 km < x < 9000 km in Fig. 2. This kind of swirls and swirl arches were also noted and described by Muthsam et al. (2010) and Moll et al. (2011) and can probably to a great extent be identified with horizontal vortex tubes that form at the edges of granules (Steiner et al. 2010).

Appendix A compares the time instant of Fig. 2 with a magnetic field-free simulation and with a simulation of initial magnetic field strength of 200 G. It also quantifies the orientation of swirls as a function of height.

Next, a statistical study is performed, with the purpose of evaluating the degree of correlation between the swirling strength λ, the magnetic field strength |B|, and the magnetic swirling strength λ^B. To do so, Spearman’s r rank coefficient is employed, which assesses how well two quantities are monotonically correlated: r = 1 indicates perfect rank correlation, r = −1 stands for rank anti-correlation, while r = 0 means that the two quantities are uncorrelated. Figure 3a shows this coefficient, as a function of height z, for the correlation between λ and |B|, r(λ, |B|), and between λ and λ^B, r(λ, λ^B).

From Fig. 3a, we see that a strong rank correlation exists between swirling strength and magnetic swirling strength, r(λ, λ^B), in particular in the near surface convection zone and in the photosphere (r ≈ 0.8). Into the chromosphere, the correlation coefficient decreases to r ≈ 0.6. Therefore, we conclude that a strong correlation exists between λ and λ^B in the solar atmosphere and, consequently, also between plasma vortices and torsional perturbations in the magnetic field. The lower values in the chromosphere can be explained by the upper boundary condition of vanishing horizontal magnetic field component, which prevents torsional perturbations at the very top of the computational domain. Given Eq. (12), one can tentatively link this high correlation to torsional Alfvén waves, which is the topic of Sect. 4.2.

Concerning the correlation between swirling strength and magnetic field strength, r(λ, |B|), Spearman’s coefficient is r ≈ 0.35 in the convection zone, drops to r ≈ 0.3 in the photosphere and rises again to r ≈ 0.4 around the classical temperature minimum near z = 500 km before decreasing to r ≈ 0.0 in the middle chromosphere. The behavior of the curve is similar to r(λ, λ^B), but the overall value is always smaller and the variation is more pronounced. The low value in the convection zone can be explained by the turbulent motion of the plasma. Indeed, the source of vortices is not related to the strength of the magnetic field, which is almost everywhere much smaller than the kinetic equipartition value in the convection zone. Swirls that are carried from the convection zone to the photosphere rapidly lose angular velocity owing to the steep decrease in mass density and correspondingly strong expansion of the plasma (Nordlund et al. 1997). This explains the rapid decrease in the average swirling strength above z = 0 km shown in Fig. 3b. Together with the many nonmagnetic or weakly magnetized swirl arches in between the magnetic flux concentrations, this reduces the correlation between swirling strength and magnetic field in the photosphere.

However, the correlation stays high for swirls within the magnetic funnels (β < 1), which is demonstrated by the green curve in Fig. 3a. While weak-field swirl-arches become less abundant with height, the correlation rises throughout the photosphere, reaching a maximum of 0.43 at z = 500 km. Higher up in the atmosphere, the correlation r(λ, |B|) rapidly drops because the magnetic field becomes increasingly homogeneous and its strength drops, while the swirling strength steeply increases and remains highly structured as it can be seen from Figs. 2 and 3b¹. This behavior is distinctly different from a magnetic field-free simulation in which the swirling strength does not show this steep increase (light blue curve in Fig. 3b).

Fig. 4. Time sequence of a single swirl event from t = 5760 s to t = 6160 s. From top to bottom row: vertical component of the magnetic field Bz at z = 0 km, continuum intensity I, vertical velocity vz at z = 700 km, vertical component of the swirling vector λz at z = 700 km, vertical component of the magnetic swirling vector ${\lambda }_z^B$ at z = 700 km, and the bin-5 intensity I5. Maps of λz and ${\lambda }_z^B$ also show the streamlines of the velocity field and the magnetic field projected into the horizontal plane at z = 700 km, respectively.

Fig. 5. Time sequence of the swirling strength distribution. Each panel represents a bi-dimensional histogram normalized to the maximum density fraction at each height level z. The bin sizes are Δz = 10 km and Δλ = 3.62 × 10−3 Hz. These histograms refer to stacks of quadratic, plane-parallel cross sections of 1000 km side length centered on the swirl event of Fig. 4. Data points with λ < 2.09 × 10−2 Hz are excluded. The black arrows point to the upward propagation of a local peak of swirling strength.

4.2. Evidence of Alfvén pulses

In this subsection, we study in detail one of the swirl events that occurred in the course of the simulation with the 50 G initial field strength. A list of eight supplementary events is presented in Appendix C. The six rows of Fig. 4 show close-ups on various observable quantities that characterize the event. The cadence of the close-ups is 50 s. The top row corresponds to a proxy of a magnetogram. It exhibits a positive polarity (upwardly pointing) magnetic flux concentration with B_z ≈ 1500 G at z = 0 km and around t = 5910 s, while it has merely half that strength before and after this maximum. The flux concentration is responsible for the bright knot located within the integranular vertex visible in the continuum-intensity maps of the second row. Next to this flux concentration is a weaker one of inverse polarity.

The third, fourth, and fifth rows show, respectively, the vertical component of the velocity field, v_z, the vertical component of the swirling vector, λ_z, and the vertical component of the magnetic swirling vector, ${\lambda }_z^B$, all in a cross section at a height of z = 700 km. We see the development of a clockwise rotating swirl (negative λ_z) from the start of the time series, best visible from t ≈ 5860 s onward in both v_z and λ_z (third and fourth rows). The positive polarity magnetic flux concentration in the top row and the BP in the second row appear to be rotating clockwise too, that is, in the same direction as the plasma does. The rotation is unidirectional (clockwise) over the full time period from 5760 s to 6160 s, that is, over 400 s. The fifth row exhibits that there also exists an azimuthal perturbation in the magnetic field but opposite to the velocity field. This aspect will be discussed in the paragraph after next.

Finally, the last row shows the maps of the synthetic bin-5 intensity, I₅. It corresponds to the fifth opacity band of the non-gray radiative transfer used in the simulation and represents an average of intensities from strong spectral lines of large opacities (Ludwig 1992). The τ₅₀₀ = 1 level computed with this opacity bin is located in the upper photosphere to lower chromosphere, and, therefore, I₅ can be taken as a proxy for chromospheric line core intensities. In these maps, we recognize a chromospheric swirl with a maximal diameter of ≈1.4 arcsec, which is in the range of sizes of the ones reported by Wedemeyer-Böhm & Rouppe van der Voort (2009).

From Fig. 4 we can deduce an important aspect of this event: The vertical component of the swirling vector λ_z and of the magnetic swirling vector ${\lambda }_z^B$ have opposite signs, that is, opposite orientations. In fact, it can be seen from the streamlines in the corresponding panels that the plasma is rotating in a clockwise fashion (negative swirling strength), while the magnetic field lines are counterclockwise twisted (positive magnetic swirling strength). This is a characteristic of upwardly propagating torsional Alfvén waves in a positive polarity (upwardly directed) magnetic field, as was pointed out in Sect. 2. Therefore, the fourth and fifth rows of Fig. 4 represent the first piece of evidence of a perturbation similar to a torsional Alfvén wave in this swirl event.

The time sequence of Fig. 4 shows the evolution of the swirl in a horizontal section at z = 700 km only. In order to gain a perspective of its evolution in the vertical direction, Fig. 5 gives the temporal evolution of the histograms of the swirling strength as a function of z. The histograms refer to a region of interest of 1.0 × 1.0 Mm² centered on the swirl and are taken with a cadence of 10 s. This region of interest is almost as large as the individual maps of Fig. 4: It is large enough to ensure that the swirl stays within this region at all heights 0 km < z < 1500 km and small enough to avoid disturbances from neighboring swirling motions.

From Fig. 5 one recognizes a decline of swirls (both in density and in strength) at photospheric levels due to reasons already mentioned in Sect. 4.1. However, there is a small bump of large swirling strengths appearing at (t, z) = (5880 s, 400 km), marked with an arrow. It grows and moves upward in time, proceeding ≈100 km every 10 s. The bump that started in the photosphere then seems to merge with other local over-densities at t = 5930 s to form a large peak around z = 1000 km, which continues its ascent after t = 5940 s. We interpret the ascension of this over-density of large swirling strength to be the signature of an upwardly propagating swirl, which starts in the photosphere at z ≈ 400 km and becomes manifest at chromospheric levels at around t = 5910 s, as was already visible in Fig. 4.

To further investigate the upward propagation of the swirling motion, Fig. 6 presents three time–distance diagrams. In these diagrams, the quantities at each height and instant (time–distance point) are averages over a horizontal plane cross-section of side length 150 km centered on the magnetic swirl event shown in Fig. 4. This side length corresponds to approximately 0.2 arcsec, which is distinctly smaller than the side length of the maps shown in Fig. 4. Panel a shows the evolution of the swirling strength λ and b the vertical component of the swirling vector λ_z.

Fig. 6. Time–distance diagrams of (a) the swirling strength λ, (b) the vertical component of the swirling vector λz, and (c) the vertical component of the magnetic swirling vector ${\lambda }_z^B$ for the swirl event shown in Fig. 4. Values at each time step and height are averaged over a finite horizontal plane of 150 km side length. The paths of two test particles moving with Alfvén plus bulk speed (green) and with sound plus bulk speed (red) along the swirl are shown in panel b.

One easily recognizes a crest of large swirling strength in both panels, starting at (t, z) = (5830 s, 0 km) and continuing up to at least (t, z) = (5970 s, 1050 km). The similarity between these two diagrams tells us that the swirl axis is to a large degree vertically directed. The crest in λ_z fades above z ≈ 900 km, while that of λ can be followed further up to z ≈ 950 km. A more careful analysis reveals that this disappearance is not due to a fading of the swirl itself but to the bending of the swirl out of the strictly vertical detection slit of 150 × 150 km² cross section used to establish these time–distance plots. In fact, from Fig. 5, which is based on a wider cross-section, it can be readily seen that the swirl propagates up to the top of the computational domain, showing still a peak in swirling strength at (t, z) = (5960 s, 1200 km).

From panels a and b of Fig. 6 it can firstly be seen that this swirl is not a rigidly rotating magnetic structure but it propagates wave-like throughout the atmosphere in the upward direction. However, at any time, there is swirling motion distributed over a large height range and only the local maximum (the crest) of swirling strength is confined to a narrower height range. Second, we notice from panel b that the swirling strength is always negative, thus the rotation of the swirl is clockwise. Only in chromospheric layers above z = 800 km there are traces of counterclockwise motion, but these do not arise from a precursory wave train starting from the photosphere. Thus, we do not deal with an oscillatory wave but with an unidirectional pulse. In fact, there exists a sequence of equally directed pulses as it can be best seen from panel b at around z = 600 km with a sequence of local maxima in swirling strength at times t = 5920 s, 5970 s, and 6020 s.

Panel c of Fig. 6 shows the evolution of the vertical component of the magnetic swirling vector ${\lambda }_z^B$. From it, it is possible to observe that the propagation of negative λ_z occurs in parallel with the propagation of positive ${\lambda }_z^B$. This behavior confirms that the vortex flow in the plasma is linked to a twist of opposite orientation in the magnetic field, which is a characteristic of torsional Alfvén waves propagating in direction of the magnetic field, as was already pointed out in the discussion of Fig. 4.

The value of ${\lambda }_z^B$ is large in the photosphere and then fades away in the chromosphere for three reasons. First, the magnetic field is weaker and therefore the value of ${\lambda }_z^B$ is smaller higher up in the atmosphere than in the photosphere. Second, the magnetic swirl becomes curved in the chromosphere bending outside of the detection slit used for establishing this time–distance diagram. Finally, the boundary conditions impose the magnetic field to be vertically directed at the top boundary, hence twists in the upper layers are suppressed. We also note that the twist angle of the magnetic field (inclination with respect to the vertical direction) does not have to be large but can be imperceptibly small where β ≪ 1.

Next, the propagation speed of the swirl is investigated. The green dotted curve of Fig. 6b shows the path of an imaginary test particle moving at speed v(z) = ⟨v_A(z)⟩ ⋅ cos(⟨θ⟩) + ⟨v_z(z)⟩, that is, with the mean local Alfvén speed plus the mean local, vertical plasma bulk-speed, starting in the photosphere at z ≈ 100 km from the ridge that leads highest up. Here, the mean is taken over the horizontal plane cross-section of side length 150 km of the detection slit and ⟨θ⟩ is the mean inclination of the swirl, that is, the angle of the swirling strength vector with respect to the z-axis. The red particle moves with speed v(z) = ⟨c_s(z)⟩ ⋅ cos(⟨θ⟩) + ⟨v_z(z)⟩, which corresponds to the mean sound plus plasma bulk-speed. It starts at the same position as the green curve.

We notice how the green trajectory perfectly follows the propagation of the vertical component of the swirling vector, which means that the swirl is moving with the local Alfvén plus bulk speed. Different from that, the red curve strongly deviates from the ridge of maximal vertical swirling strength starting from approximately z = 400 km, indicating that the swirl does not propagate with sound speed. This is the second piece of evidence of the Alfvénic nature of the swirl event.

It is essential to notice that the vertical component of both the swirling vector and the magnetic swirling vector do not change sign in time. This means that the swirling motion, as well as the perturbation of the magnetic field, are unidirectional. Therefore, one must think of a torsional Alfvén pulse instead of a generic, oscillatory, torsional Alfvén wave. The existence of such torsional pulses in the solar atmosphere has been suggested by Liu et al. (2019a) using numerical simulations with a torsional perturbative driver acting on a preexisting magnetic flux tube. Here, it has been demonstrated that such torsional Alfvén pulses also emerge from a fully self-consistent, realistic, numerical simulation of the solar atmosphere.

The present section focused on a single, isolated swirl event. More often, simulations show swirling motions in the chromosphere to appear as a result of the superposition of two or more individual swirls. Examples of the superposition of swirls are displayed and discussed in Appendices B and C.

Fig. 7. Time–distance diagram of (a) the swirling strength λ and of the vertical component of the following quantities: (b) swirling vector λz, (c) magnetic swirling vector ${\lambda }_z^B$, (d) partial derivative of the swirling vector ∂λz/∂t, (e) advection term $T_z^{\text{adv}}={[(\mathit{v}\cdot \nabla )\mathit{\lambda }]}_z$, (f) stretching term $T_z^1$, (g) hydrodynamical baroclinic term $T_z^2$, (h) magnetic baroclinic term $T_z^3$, (i) sum of baroclinic terms $T_z^{\text{bar}}=T_z^2+T_z^3$, and (j) magnetic tension term $T_z^4$. The black dashed contours correspond to the λz = −0.1 Hz contours, which highlight the propagation of the Alfvén pulse. All plots are relative to the swirl event shown in Fig. 4 and are derived from a fixed, one pixel wide vertical slit through the center of that swirl.

4.3. Dynamics of the pulse

This section focuses on the dynamics of the magnetic swirl shown in Fig. 4, in particular, on the physical processes responsible for the production and propagation of the swirling motion and Alfvén pulse. For this purpose, we employ the swirling equation, Eq. (10), to reveal the local production of the vertical component of the swirling vector by the various source terms. Each one of these source terms is related to a different physical mechanism, as described in Sect. 2.3.

Figure 7 shows time–distance diagrams of various quantities. Different from Fig. 6, the time–distance slit in this case is given by a single, vertical column of computational cells (single line-of-sight), located in the center of the magnetic swirl. Panel a shows the time–distance plot for the swirling strength λ, b the vertical component of the swirling vector λ_z, and c the vertical component of the magnetic swirling vector ${\lambda }_z^B$, similar to Fig. 6. However, with the single line-of-sight slit, we see two separate pulses of swirling strength, both with the same sense of rotation. They propagate upward at local Alfvén speed: the first starting at (t, z) = (5800 s, 0 km), the second at (t, z) = (5850 s, −50 km). These pulses are well paired with perturbations in the magnetic field, as can be seen from panel c. These plots reproduce in more detail what was already observed from Fig. 6.

Figures 7d,e show the time–distance diagrams of the local production and the advection of the vertical component of the swirling vector, ∂λ_z/∂t and [(v ⋅ ∇)λ]_z, respectively. Together, these two quantities constitute the material derivative that appears on the left-hand side of Eq. (10). We notice that the local production of negative swirling strength between −200 km ≲ z ≲ 200 km and 5760 s ≲ t ≲ 5860 s in panel d corresponds to the appearance of negative swirling strength in panel b. The black dashed lines are contours of λ_z = −0.1 Hz. They are intended to visualize the boundaries of the two pulses and guide the eye in identifying the relevant source terms. Also from panel d, we notice that the two upwardly propagating pulses are characterized by the production of negative swirling strength on the left (leading) boundary of the two stripes, which is compensated for, approximately 20 s later, by the production of positive swirling strength on the right (trailing) boundary. We conclude that these two sources cause the upward evolution of the swirl.

Contributions from the advection term to the local production of swirling strength are mostly subdominant in the photosphere and into the low chromosphere, as one can see from panel e. Therefore, we conclude that the swirling strength is mainly locally produced in these regions. Only in the upper part of the chromospheric layers, the advection term produces some non-negligible positive vertical swirling strength². This production can only contribute to the reduction of the strong pulse in the height range 600 km < z < 800 km around t = 5920 s and to the fading of the tip of swirling strength around (t, z) = (5900 s, 800 km), but it is not the main cause. Therefore, if the dynamics of the swirling strength alone cannot explain the fading of the swirl beyond 850 km, we conclude that it must be the tilted swirl axis to cause it, as already hypothesized in Sect. 4.2. This tilt also leads to the reduction of the vertical propagation speed of the swirl, seen to set in already beyond z ≈ 600 km.

The second row of Fig. 7 shows the time–distance diagrams of the various terms of Eq. (10). It serves to infer which physical processes are responsible for the production of negative swirling strength at the photospheric level and for the propagation of the pulses. The gravitational term $T_z^5$ is omitted because the simulations employ a constant gravitational field and thus $T_z^5$ = 0. First of all, analyzing the different source terms, we notice that the stretching term, $T_z^1$, is much weaker than the other terms. Second, the two baroclinic terms, $T_z^2$ and $T_z^3$, have similar configurations but are of opposite sign. Therefore, panel i shows the sum $T_z^2+T_z^3$ from which it can be seen that most of their contributions cancel out. This result is in accordance with the findings of Canivete Cuissa & Steiner (2020).

The counter-balancing of the two baroclinic effects inside and across the boundary of a stationary magnetic flux tube can be explained by the required equilibrium between magnetic pressure, magnetic tension forces, and gas pressure. However, the magnetic flux tube hosting the swirl is not stationary as can be seen from Fig. 4. The dynamics of the magnetic structure and of the surrounding plasma can break the balance between the pressures, which can result in a net production of swirling strength, as visible from panel i. Thus, it becomes clear that the local production of negative swirling strength seen in panel d between −200 km ≲ z ≲ 100 km and 5780 s ≲ t ≲ 5860 s is caused by the imbalance between the two baroclinic terms. More specifically, around z = 0 km and between 5790 s  ≲ t ≲ 5850 s the baroclinic magnetic term, $T_z^3$, which is producing negative swirling strength, overcomes the hydrodynamical one, $T_z^2$, which instead produces positive swirling strength.

Because the time–distance slit is located in the center of the magnetic flux concentration, where the β = 1 surface dips below z = 0 km into the convection zone (see Fig. 2), the magnetic field dominates the gas pressure and the dynamics in this location. This suggests that the origin of the swirl is due to magnetic effects alone that could arise from an MHD instability or from magnetic annihilation, possibly owing to nearby inverse polarity fields such as in the case of quiet Sun Ellerman bombs (QSEBs, Joshi et al. 2020). Alternatively, the dominance of the baroclinic magnetic term may also arise from a weakening of the hydrodynamic baroclinic term caused by, for example, a sudden low pressure region or the bath tub effect (see Wedemeyer-Böhm et al. 2012) in the convective motion surrounding the magnetic flux concentration. Revealing the true origin of the swirls calls for a multidimensional analysis of the different source terms of swirling strength, which goes beyond the scope of the present paper but shall become subject of a subsequent study.

Having dealt with the triggers of such an event, we next pay attention to the propagation of the swirl. As already mentioned further above, the upward propagation is caused by the production of negative swirling strength and its successive destruction by the production of a positive counter-part, approximately 20 s later. The inspection of the source-term diagrams of Fig. 7 reveals that the magnetic tension term, $T_z^4$, is the main responsible for the production and destruction of λ_z along the two pulses of swirling strength, as can be seen from Fig. 7j. The stretching term, $T_z^1$, plays no role in the propagation of the swirl, the hydrodynamical baroclinic term $T_z^2$ is overcome by its magnetic counterpart $T_z^3$, which in turn only partially contributes to the production of the first pulse of swirling strength and to the destruction of the second. Therefore, we can conclude that magnetic tension forces are responsible for the propagation of the swirl in the photosphere and low chromosphere. This analysis provides further evidence of the Alfvénic nature of the vortex flow since the magnetic tension is the driving force of Alfvén waves.

Fig. 8. Time–distance diagrams of the vertical component of the Poynting flux vector Sz (a) and of the two terms that compose it, $S_z^{\text{v}}$ (b) and $S_z^{\text{h}}$ (c), for the swirl event shown in Fig. 4. Values at each time step and height level are averages over a finite horizontal plane of 150 km side length.

4.4. Energetic considerations

The purpose of the present section is to investigate how much energy swirling events and associated Alfvén pulses can channel up from the photosphere to the upper atmospheric layers. To this end, we analyze the vertical component of the Poynting flux vector and the mechanical energy flux relative to the swirling motions, as given in Sect. 2.2. The first subsection concentrates on the single swirl event of Fig. 4 and the second one evaluates the mean Poynting flux over the full computational domain.

4.4.1. Energy transported by the swirl of Fig. 4

Figure 8a shows the time–distance diagram of the vertical component of the Poynting flux, S_z, associated with the swirl event of Fig. 4. Like in Fig. 6, the distance slit stretches from z = −200 km to z = 1200 km. It has a finite horizontal, quadratic extension of 150 km side length, and S_z is taken to be the average over this finite cross section. At the location of origin of the Alfvén pulse at around z = 0 (as was seen from Fig. 7), the upwardly directed Poynting flux is on the order of 10⁶ W m⁻². The flux decreases with height but it is still on the order of 10⁵ W m⁻² at chromospheric levels. Panels b and c of Fig. 8 show the corresponding time–distance diagrams of the terms $S_z^{\text{v}}$ and $S_z^{\text{h}}$ respectively, which are defined in Eq. (4). We notice that the upward flux is entirely due to the term $S_z^{\text{h}}$, which means that the energy is carried by the horizontal motions of the magnetized plasma. Hence, the magnetic swirl and associated Alfvén pulse are responsible for the upwardly directed energy transfer.

Fig. 9. Vertical component of the Poynting flux vector Sz of the swirl event presented in Fig. 4 in a horizontal section at z = 700 km. Arrows indicate the velocity vector field projected into the horizontal plane. Their length is proportional to the magnitude of the horizontal flow.

Figure 9 reveals the spiral pattern outlined by the vertical component of the Poynting flux, S_z, in a horizontal section centered on the swirl and at (t, z) = (5980 s,  700 km). It roughly follows that of the swirling plasma flow, as it can be seen from the overlaid velocity field. This reinforces the idea that the swirl is carrying magnetic energy through the solar atmosphere. The maximum positive Poynting flux of the swirl is S_z, max = 3.5 × 10⁵ W m⁻², while the mean vertical Poynting flux is S_z = 34.5 kW m⁻², where the mean was taken over a circular area of diameter 1.2 arcsec centered on the swirl and over the time period from t = 5760 s to t = 6050 s.

From these numbers, the contribution of the swirl alone to the mean vertical Poynting flux through the cross section of the full computational domain of 9.6 × 9.6 Mm² is 236 W m⁻². This value is smaller but on the order of the 440 W m⁻² reported by Wedemeyer-Böhm et al. (2012) for the contribution to a cross section of 8 × 8 Mm², however, at a height of 2000 km, that is, at the base of the transition zone. We attribute this weak Poynting flux of the swirl of Fig. 4 to its relatively weak magnetic field and to its bending above z ≈ 700 km.

Liu et al. (2019a) count from observations with the Ca II H Broadband Filter Imager (BFI) of the Solar Optical Telescope (SOT) of Hinode ${\overline{N}}_{\mathrm{s}}=48.2$ swirls at any time in a field of view of A = 800 Mm². Assuming that these swirls are similar to the one of Fig. 4 with a mean Poynting flux of S_z = 34.5 kW m⁻² over a circular cross section of diameter 1.2 arcsec, we obtain from Eq. (5) an average energy flux of ${\overline{S}}_z=1.2\mathrm{kW}{\mathrm{m}}^{-2}$ at a height of z = 700 km. This flux is not enough to compensate the radiative losses of approximately 4.3 kW m⁻² (excluding Lyman α) in the semiempirical model of the chromosphere by Vernazza et al. (1981) but could still account for an important source of energy in the upper chromosphere and the corona.

In the same way as for the Poynting flux, it is possible to estimate the mean mechanical energy flux owing to swirls. We take again the single swirl of Fig. 4 as a model and average again over a fixed circular area of diameter 1.2 arcsec and the time period from t = 5760 s to t = 6050 s. Furthermore, we assume as propagation speed the local Alfvén plus bulk speed in the vertical direction, v_p = v_Acos(θ) + v_z, as done in Sect. 4.2, to compute the mechanical energy fluxes according to Eqs. (6) and (7).

Using Eq. (6), the estimate for the mean vertical mechanical flux owing to horizontal velocities is $F_z^{\text{r}}=6.1\text{kW}{\text{m}}^{-2}$, while Eq. (7) yields $F_z^{\lambda }=15.2\text{kW}{\text{m}}^{-2}$. These formulas give only an approximate upper limit of the true mechanical energy flux because of a radial velocity component that was not removed in the computation of $F_z^{\text{r}}$, while $F_z^{\lambda }$ is overestimated because of contribution from peripheral swirling strength not pertaining to the swirl under consideration. Nevertheless, the computed mechanical flux is only a fraction of the Poynting flux, which reflects the circumstances that the plasma-β ≲ 1 in the considered region.

Assuming, once again, that the swirl shown in Fig. 4 is representative of the ${\overline{N}}_{\mathrm{s}}=48.2$ swirls detected at any time in a field-of-view of A = 800 Mm², and using Eq. (8) to compute the upper limit of the average mechanical energy flux owing to swirls, we find, at z = 700 km, ${\overline{F}}_z^{\mathrm{r}}=218\mathrm{W}{\mathrm{m}}^{-2}$ and ${\overline{F}}_z^{\lambda }=544\mathrm{W}{\mathrm{m}}^{-2}$, using Eq. (6) and Eq. (7), respectively. For comparison, Liu et al. (2019a) derive, from observations alone, a lower limit of the average mechanical energy flux from swirls of 33–131 W m⁻² at a height of approximately z = 1000 km.

Fig. 10. Vertical component of the Poynting flux vector Sz over the entire computational domain at z = 700 km and t = 5800 s. Black circles indicate the location of some of the swirl events listed in Table C.1 also shown in Fig. C.1.

Fig. 11. Mean vertical component of the Poynting flux vector Sz (blue) and of the terms $S_z^{\text{h}}$ (green) and $S_z^{\text{v}}$ (orange) as a function of height z. The spatial mean is obtained by averaging over the horizontal cross sections of the full computational domain of 9.6 × 9.6 Mm2, while the temporal mean is computed from 21 different time instants spread over the time period of the high-cadence series of 15 min. The shaded areas correspond to the standard deviation obtained from the temporal variation.

4.4.2. Mean energy transport

Looking at the mean net Poynting flux across the full computational domain of 9.6 × 9.6 Mm², it turns out that this flux is much larger than the extrapolation from the single, relatively weak, isolated swirl event of Fig. 4, and that the bulk of energy transport is delivered by the more complex events of superposition of swirls (see Appendix B). This situation is illustrated with Fig. 10, which shows the vertical component of the Poynting flux on the horizontal cross section at z = 700 km for the same time instant as Fig. C.1, that is, t = 5800 s. In this figure, the largest single contribution to the mean vertical Poynting flux of 18.4 kW m⁻² comes from event b near the center of the field of view with a peak Poynting flux of 2123.7 kW m⁻². Comparing it with Fig. C.1, it can be seen that this event starts from a relatively large magnetic footpoint patch with a complex but clearly swirling motion in the bin-5 intensity (see the animation of Fig. C.1). For comparison, event f corresponds to the single, isolated swirl of Fig. 4.

Figure 11 shows the temporal and spatial mean of the vertical component of the Poynting flux S_z together with the terms $S_z^{\text{v}}$ and $S_z^{\text{h}}$ as a function of height z. The shaded area surrounding each curve indicates the standard deviation of the temporal variation over the time period of 15 min of the high-cadence simulation³. First we observe that all curves converge to zero flux at the upper boundary, which is a consequence of the top boundary condition for the magnetic field. Therefore, magnetic disturbances are reflected at the top boundary with the consequence that the net vertical Poynting flux vanishes in the upper part of the atmosphere⁴. Second, the contribution due to horizontal motions (including swirls) is always positive, while vertical motions contribute negatively to the Poynting flux. Above z ≈ 1000 km, the negative contribution is even slightly dominating, possibly due to the boundary condition of vertical magnetic field at the top, which favors vertical motions in the upper part of the atmosphere. In this height range, β < 0.1 in regions of significant Poynting flux so that the top boundary condition strongly influences the magnetic field down to z ≈ 1000 km because of the stiffness of this field. However, at the base of the chromosphere, which we take here again to be located at z = 700 km, the mean flux is 12.8 ± 6.5 kW m⁻², which, for a comparison, would be sufficient to compensate the radiative losses of 4.3 kW m⁻² in the chromosphere.

We also notice that at this height almost all the vertical Poynting flux is carried by the horizontal motions of the magnetized plasma (S_z ≈ $S_z^{\text{h}}$), that is, largely by swirling motions. This result agrees with the findings reported by Yadav et al. (2020). Moreover, they demonstrate that the largest contribution to the Poynting flux comes from small-scale vortices that cluster to form the larger, observable swirls. This is consistent with our finding that superposition of swirls can carry large amounts of magnetic energy as in the case of event b of Fig. 10. This large magnetic energy flux is tantalizing, considering the findings of Fossum & Carlsson (2005) that high-frequency acoustic waves are not sufficient to heat the solar chromosphere and that the emission from the middle and upper chromosphere must be balanced by processes related to the magnetic field: a problem that is still unsolved (Carlsson et al. 2019). Moll et al. (2012) showed with the help of simulations similar to the present ones that Ohmic dissipation of magnetized swirling motions induce substantial local heating, although mostly in the photosphere.

However, a word of caution is in order here. The initial condition of a homogeneous vertical magnetic field together with the periodic lateral boundary conditions, of which the simulation keeps a memory for all times, favor the formation of vertically extended vortices. These simulations may adequately represent regions with a predominantly vertically directed magnetic field such as magnetic network patches or plage regions. With a more heterogeneous magnetic field, such as that generated by a turbulent dynamo in very quiet regions, magnetically induced swirls are probably less numerous and hence the Poynting flux less vigorous.

5. Summary and conclusions

We carried out numerical radiation MHD simulations of the solar atmosphere with the purpose of studying the origin and propagation of vortical motions that occur in conjunction with magnetic fields and bear great resemblance to observed chromospheric swirls. The simulations are of fairly high spatial and temporal resolution and reach from the convectively unstable subsurface layers to the chromosphere. For the analysis, we mostly use a run that started with a homogeneous vertical magnetic field of a strength of 50 G, but use for comparison also runs that started with a magnetic field of 200 G and without magnetic field.

Looking first at the distribution of swirling strength, which is a measure for the vigor of swirling motions, we find a drastic difference between the simulations with and those without a magnetic field. The simulations with magnetic field, which show the usual funnel-shaped magnetic flux concentrations in the photosphere, have the swirling motion strongly concentrated within these low β funnels and aligned with the magnetic field. Different from that, the simulation without magnetic field possesses more arch-shaped vortex tubes in the photosphere and almost isotropically distributed vortices in the chromosphere. These results agree with and confirm earlier findings of Moll et al. (2011) and Steiner & Rezaei (2012).

Beyond this relation between magnetic field and swirling motion, we find a tight relationship between swirling motions and perturbations in the magnetic field of intense magnetic flux tubes. In particular, we show with the help of a statistical analysis that the swirling strength and the magnetic swirling strength are highly correlated in the solar atmosphere. While swirling strength indicates a vortex in the plasma, the magnetic swirling strength indicates a twist in the magnetic lines of force. The latter typically occurs in the presence of torsional Alfvén waves. Therefore, the tight correlation between these two quantities is a first indication of the relationship between swirling motions and torsional Alfvén waves in the simulated solar atmosphere.

We then analyze one particular swirl event that was identified in the numerical simulations (see Fig. 4). The well developed chromospheric swirl is colocated with a positive polarity magnetic BP in an intergranular lane of the deep photosphere and presents a twist in the magnetic field. Therefore, it can be identified as a magnetic swirl. Using time–distance diagrams to study its evolution in time and altitude, we find three pieces of evidence of the Alfvénic nature of this magnetic swirl. First, the vortical motions and the magnetic perturbations, in the upwardly directed (positive polarity) magnetic field, have opposite orientation; second, the swirl propagates upward with the local Alfvén speed and, third, the driving force that is responsible for such propagation is the magnetic tension. These are all characteristics of torsional Alfvén waves. However, we do not observe an oscillatory wave, but instead pulses of unidirectional vortical motions and twisted magnetic field lines. From all that, we conclude that the swirl subsists in an upwardly propagating torsional Alfvén pulse. The other eight swirls listed in Table C.1 also share these characteristics and are therefore also identified as torsional Alfvén pulses. A similar conclusion was also reached by Liu et al. (2019a) from observations and modeling. The novelty of the present work is that these Alfvén pulses are not deliberately excited but naturally arise from a self-consistent numerical solution of the MHD equations and radiative transfer.

Regarding the origin of the Alfvén pulse, we cannot say with certainty if it is predominantly hydrodynamic or magnetic. From the analysis based on the swirling equation, we find that magnetic baroclinicity in the top convection zone and the deep photosphere is responsible for the creation of the swirl. Hence, magnetic dynamics alone would launch the pulse. However, the present time–distance analysis is based on a single, vertical slit through the center of the swirl, which does not follow the temporal evolution of the swirl nor its geometrical shape. This can lead to misinterpretations of the diagrams, especially in the chromosphere where the plasma flows are fast. Therefore, it is necessary to consider the full three-dimensional structure of the individual terms of the swirling equation as a function of time in order to make progress. Such a three-dimensional analysis was beyond the scope of the present paper, but it can be expected to reveal whether the magnetic baroclinicity in the center of the swirl is caused by magnetic effects alone (e.g., magnetic diffusion or an MHD instability such as that responsible for quiet Sun Ellerman bombs) or by hydrodynamic effects in the surroundings of the magnetic flux concentration (e.g., by the bathtub effect or the granular buffeting).

The upwardly propagating torsional Alfvén pulses carry a substantial amount of magnetic and mechanical energy from the photosphere to the outer layers of the atmosphere. At the top boundary of the computational domain, the Poynting flux vanishes because of the adopted boundary condition that forces the field to become strictly vertical. Therefore, we limited our analysis to the base of the chromosphere for which we took here the level z = 700 km.

For the swirl of Fig. 4, we evaluate at this height level a mean Poynting flux in the upward direction of S_z = 34.5 kWm⁻². This flux is generated by the horizontal motions of the magnetized plasma alone and is therefore an effect of the horizontally swirling motion. Taking the observed frequency of chromospheric swirls from Liu et al. (2019a) and assuming they were all similar to the swirl of Fig. 4, we obtain a mean Poynting flux of ${\overline{S}}_z=1.2\mathrm{kW}{\mathrm{m}}^{-2}$ and for the corresponding mean mechanical energy flux we find an upper limit of ${\overline{F}}_z=218-544\mathrm{W}{\mathrm{m}}^{-2}$. This energy flux could source the energy budget of the higher layers, but it is not enough to compensate alone for radiative losses at chromospheric level.

However, the swirl of Fig. 4 is relatively weak; when evaluating the average Poynting flux across the full computational domain, still at a height of z = 700 km, we obtain 12.8 ± 6.5 kW m⁻². This large amount mainly stems from more complex events of superposition of swirls, which originate in relatively large, complex, and intense magnetic footpoints. Assessing these numbers one should bear in mind that the initial condition of a homogeneous, vertical magnetic field favors the formation of vertically extended vortices and may be rather applicable to regions with a predominantly vertically directed magnetic field such as network and plage regions than to quiet Sun regions.

Would it be possible to directly observe the Alfvénic nature of chromospheric swirls? We believe that it should be feasible with high resolution and high sensitivity spectropolarimetry in two spectral lines, preferably sampling the photosphere and the chromosphere. The challenge would consist in measuring the transverse magnetic field accurately enough to prove the twist of the magnetic field relative to the vortex motion that would be detected with the help of LCT. The two height levels would inform about the propagation speed of the swirl. High photon flux is needed for accurately measuring the transverse Zeeman effect, which requires a large aperture solar telescope such as the D. K. Inouye Solar Telescope (DKIST).

Movie

Movie 1 associated with Fig. C.1.(movie_IcoBin5) Access here

Acknowledgments

AFB wishes to acknowledge support by the Swiss National Science Foundation (SNF) through grant 200021_189180 (Solar Orbiter STIX), JRCC support by SNF under grant ID 200020_182094, and FC support through the CHROMATIC project (2016.0019) of the Knut and Alice Wallenberg foundation. Special thanks are addressed to L. Harra who directed AFB’s Master’s thesis from which the present work derives. This work has profited from discussions with the team of K. Tziotiou and E. Scullion (conveners) “The Nature and Physics of Vortex Flows in Solar Plasma” at the International Space Science Institute (ISSI) and with the Waves in the Lower Solar Atmosphere (WaLSA; www.WaLSA.team) team (S. Jafarzadeh, convener), which is supported by the Research Council of Norway (project number 262622). Simulations were carried out at the Swiss National Supercomputing Centre (CSCS) under project ID s560 with support from SNF under grant ID 200020_157103. Thanks are also extended for valuable input by G. Vigeesh and P. Rajaguru, and for very helpful and encouraging comments by an anonymous referee.

References

Appendix A: Strength and distribution of vortices

This appendix compares three different CO5BOLD simulations: one purely hydrodynamic without magnetic field and two MHD ones with an initial homogeneous, vertical magnetic field of 50 G, and 200 G. We shall refer to them as Hydro, MHD 50 G, and MHD 200 G hereafter. The MHD 200 G run corresponds to d3gt57g44v200fc of Calvo (2018), while the Hydro run to d3gt57g44roefc of Calvo et al. (2016) and Calvo (2018).

Figure A.1 shows instants of the swirling strength λ for each of the three simulations: from top to bottom, MHD 200 G, MHD 50 G, and Hydro. The τ₅₀₀ = 1 surface is displayed in red. Swirls with λ < 2.09 × 10⁻² Hz, that is, with a period larger than 10 min, are neglected, while strong swirls with λ > 1.0 Hz are saturated. Comparing the three panels of Fig. A.1, one recognizes that below the surface of τ₅₀₀ = 1, the three simulations share a nearly isotropic distribution of swirls. This behavior was also found by Moll et al. (2011). These swirls are obviously caused by the turbulent plasma motions in the surface layers of the convection zone. However, the Hydro simulation shows a higher density of swirls than the MHD simulations. This is because the magnetic fields hamper the formation of vortices in the convective layers of the simulation.

The photosphere of the MHD 200 G time instant harbors an almost homogeneous structure of vertically oriented swirls and only a few swirl arches are present near the τ₅₀₀ = 1 surface. In the MHD 50 G case, the vertical swirls are clustered in funnel-shaped structures and the arches are more prominent than in MHD 200 G. Different from these two cases, vertical swirls are essentially absent in the Hydro simulation, while swirl arches are omnipresent in the photosphere. Again, these results well agree with and confirm earlier results of Moll et al. (2011, 2012).

A similar behavior is observed in the chromosphere: the stronger the initial magnetic field, the more abundant the vertical swirls are. However, we observe that the strongest swirls appear at the top of the MHD 50 G simulation. This suggests that too strong magnetic fields can possibly inhibit the formation of strong vertical vortices in the chromosphere. In the Hydro simulation, we see a quasi-isotropic distribution of swirls at chromospheric levels, which are most probably related to hydrodynamic shocks taking place in this region of the field-free solar atmosphere (see, e.g., Moll et al. 2012). Therefore, not all swirls found at chromospheric levels are necessarily related to the magnetic field, even though the magnetism clearly plays a cardinal role in the generation of vertical vortices.

The above more qualitative assessment is confirmed by the quantitative analysis of Fig. A.2, which shows for each simulation a bi-dimensional histogram of the inclination angle cos(θ) of the swirling strength vector at each height z. The angle θ is defined through

$$\begin{array}{c}\hfill cos(\theta )=\frac{\mathit{\lambda }·{\widehat{\mathit{e}}}_z}{|\mathit{\lambda }|}=\frac{{\lambda }_z}{\lambda },\end{array}$$

being the angle between the swirling strength vector λ and the unit vector in the vertical direction ${\widehat{\mathit{e}}}_z$. Therefore, horizontally oriented swirls are binned close to cos(θ) = 0, while vertically oriented swirls are close to cos(θ) = 1. We choose to plot the histogram of cos(θ) instead of θ in order to easily recognize isotropy. In this way, an isotropic angle distribution at a certain height z is characterized by a flat distribution in the histograms of Fig. A.2.

In the convection zone, that is, for z ≲ 0 km, all three simulations show a flat distribution, which confirms the nearly isotropic distribution of swirls in that region. Also, as estimated before, there are on average more swirls in the Hydro simulation (panel c) than in the MHD 50 G (panel b) and the MHD 200 G (panel a) simulations. Near the optical surface, that is, around z = 0 km, there is a sudden break of isotropy. Two families of vortices can be spotted. First, a cluster of horizontally oriented swirls (cos(θ)≈0), visible in all three panels: They stem from the arches of swirling strength that are seen to populate the near surface layers. Second, a cluster of vertically oriented swirls (cos(θ)≈1) present only in the MHD simulations. Both groups of swirls were already seen in Fig. A.1. We notice that the number of horizontally oriented swirls slightly increases as the initial magnetic field decreases. The opposite happens for the vertically oriented swirls, which are more numerous in the MHD 200 G simulation than in MHD 50 G. Higher up, the vortices are essentially uniquely vertical in the case of MHD 200 G (panel a), still predominantly vertical in MHD 50 G (panel b) and almost isotropically distributed in the Hydro simulation (panel c). This confirms once again the visual perception given by Fig. A.1.

Caution is indicated concerning the regions close to the top and to the bottom boundaries of the computational domain. For the MHD simulations, the magnetic field is forced to be vertical at the top and bottom boundaries. In the purely hydrodynamic case, this condition is not present. Nevertheless, a large number of weak, vertical swirls are encountered at the top and bottom boundaries of all three simulations. Most of them have been removed by the application of the threshold, λ = 2.09 × 10⁻² Hz, in Figs. A.1 and A.2. However, an over-density of vertical swirls is still present at the bottom boundary of panels a, b, c, and at the top boundary of the Hydro simulation.

Panel d of Fig. A.2 shows the maximum value of the swirling strength, λ_max, as a function of z for the three simulations. From it, we notice that in the convection zone the stronger the initial magnetic field, the weaker the swirls. However, the overall behavior is similar and the local peak is reached in all three cases around z ≈ −500 km. Thus, relatively strong vortices can develop in the convection zone, independently of the initial magnetic configuration of the simulation.

In the photosphere, a local minimum in swirling strength is reached. It is a result of the slowdown of vortical motions in the rapidly expanding convective overshoots as a consequence of angular momentum conservation (Nordlund et al. 1997). For the two MHD simulations, the minimum is reached in the lower half of the photosphere, while the Hydro simulation reaches its minimum at z ≈ 500 km. This difference can be attributed to strong arches of swirling strength in the hydrodynamical simulation, which permeate the photosphere and extend precisely up to z ≈ 500 km, as it can be seen in Fig. A.1. These arches are less developed in the MHD simulations, probably as a consequence of the magnetic canopy.

Fig. A.1. Volume rendering of the swirling strength seen from the side for three different simulations at arbitrary time instants: (top) MHD simulation with an initial vertical, homogeneous magnetic field of 200 G; (middle) MHD simulation with an initial field of 50 G; (bottom) purely hydrodynamical simulation. The red surface indicates the τ500 = 1 surface. Swirls with a period larger than 10 min (λ < 2.09 × 10−2 Hz) have been neglected and values above 1 Hz are saturated.

Finally, looking at the photosphere and chromosphere, we observe that λ_max increases exponentially for the MHD 50 G simulation, while this growth is only linear for the MHD 200 G simulation. In the Hydro case instead, there is a growth only in the low chromosphere up to approximately z ≈ 1200 km. Above this height, the swirls become weaker again. Clearly, magnetic fields favor the generation of swirls. However, it seems that too strong magnetic fields can also suppress the formation of strong vortices. We surmise that there exists an optimal value of the initial magnetic field for which the strength of the vortices is maximal.

Fig. A.2. Distribution of the inclination cos(θ) of the swirling strength vector as a function of altitude z for three different simulations: (a) MHD simulation with an initial vertical, homogeneous magnetic field of 200 G, (b) MHD simulation with an initial field of 50 G, and (c) purely hydrodynamical simulation. The inclination θ refers to the angle between the swirl axis and the vertical direction. Horizontally oriented swirls are binned close to cos(θ) = 0, while vertically oriented swirls are close to cos(θ) = 1. The maximum value of swirling strength λmax at each height z, for each simulation, and with a minimum of 10 occurrences per bin, is shown in panel d. The occurrences are means over 5 different time instants of 960 × 960 × 280 data points for each simulation. The bin sizes are Δz = 10 km and Δcos(θ) = 3.62 × 10−3. Swirls with λ < 2.09 × 10−2 Hz are neglected.

Appendix B: Superposition of swirls

Section 4.2 considers the torsional Alfvén pulses of a single swirl event. An important characteristic of that event is that its magnetic footpoint is small compared with others that can be found in the simulation and that it is quite isolated from the rest of magnetic flux concentrations. In this sense, the swirl of Fig. 4 is an ideal, clear-cut case, which we have intentionally chosen for ease of analysis. Whenever the magnetic footpoint is larger or when it consists of several footpoints close together, as may be representative of a magnetic network patch, a filigree, or a magnetic knot, the corresponding intensity pattern in the chromosphere becomes more involved. These sorts of complex magnetic structures can be energetically more important since multiple swirls may coexist in there. In the light of this phenomenon, we introduce the concept of “superposition of swirls”.

Figures B.1 and C.2 both show the time sequence of a superposition of swirls in which the magnetic footpoint is distinctly larger and has a more complex structure than that of Fig. 4. The figures are organized in the same way as Fig. 4. Especially in the case of Fig. B.1, the bin-5 intensity, I₅, does no longer show a clear swirling motion, as was the case in Fig. 4. However, this does not imply that no swirling motions are present. Indeed, the vertical component of the swirling vector clearly shows two distinct swirls rotating in opposite direction: a clockwise vortex in the top left corner of the field of view and a counterclockwise one in the bottom right corner.

In principle, the same analysis done in Sect. 4 also applies in the case of these two swirls. The more complex pattern of the bin-5 intensity can be explained by the geometry and dynamics of the large footpoint. To see this, Fig. B.2 depicts the evolution of the magnetic field strength at z = 0 km for the same field of view as it is shown in Fig. B.1. Each map is overlaid with the vertical component of the swirling vector in regions where λ > 10⁻³ Hz and λ^B > 10⁻⁶ G cm⁻¹. Here, red indicates positive (counterclockwise) swirling strength and blue stands for negative (clockwise) swirling strength. The distribution of the swirling strength in this region is fragmented since within the same magnetic structure several spots of swirling strength with different orientations occur. The reason for this is that the whole magnetic footpoint is not rigidly rotating, but possesses substructure that can have different relative motions and rotations. Consequently, different rotation directions within the same magnetic footpoint can be expected to occur. As the magnetic fields expand with height and start to fill the entire available volume at chromospheric levels, these individual rotational motions result in a superposition of different chromospheric swirls. (Provided that the vertical propagation effectively takes place, as is the case here). However, a swirling motion detected at z = 0 km does not necessarily propagate up to the chromosphere. Some of these swirls may be the footpoints of vortex arches, they may be too weak to reach the chromosphere, or may become dominated by neighboring upwardly propagating perturbations.

Fig. B.1. Time sequence of a superposition event from t = 5750 s to t = 5830 s. From top to bottom row: vertical component of the magnetic field Bz at z = 0 km, continuum intensity I, vertical velocity vz at z = 700 km, vertical component of the swirling vector λz at z = 700 km, vertical component of the magnetic swirling vector ${\lambda }_z^B$ at z = 700 km, and the bin-5 intensity I5. Maps of λz and ${\lambda }_z^B$ also show the streamlines of the velocity field and the magnetic field projected into the horizontal plane at z = 700 km, respectively.

Fig. B.2. Time sequence of the absolute magnetic field strength at z = 0 (gray scale) in the same field of view as shown in Fig. B.1. Plotted on top of the field strength is the vertical component of the swirling vector λz in regions where λ > 10−3 Hz and λB > 10−6 G cm−1. Red indicates counterclockwise and blue clockwise rotation, that is, positive and negative values of λz, respectively.

Figure B.2 shows different patches of clockwise (negative) swirling strength within the left part of the magnetic footpoint. This then becomes manifest as a clockwise vortical motion in the plasma visible in the top left corner of λ_z at z = 700 km from t = 5780 s onward (fourth row of Fig. B.1). We also notice that the horizontal component of the positive polarity magnetic field has opposite direction with respect to the velocity field, indicating the torsional Alfénic nature of the pulse. Within the right part of the magnetic footpoint, we observe a dominance of patches of counterclockwise (positive) swirling strength. These become manifest as a counterclockwise vortical motion in the upper photosphere to lower chromosphere, which can be seen in the bottom right corner of λ_z in Fig. B.1.

The consequences for chromospheric observations can be guessed from the bin-5 intensity maps, which are distinctly different from the ones shown in Fig. 4. They still show traces of swirling motions, but they are less obvious than in the case of an unidirectional swirl event. The superposition of swirls is a ubiquitous phenomenon in our simulations since large footpoints are more frequent than small, isolated ones. The present interpretation is consistent with the findings of Yadav et al. (2020). Using MHD simulations, they studied vortex flows at different spatial scales and concluded that the observed large vortices likely consists of clusters of smaller swirls that are not resolved by the observations.

In the case of superposed swirls, caution is indicated when analyzing time–distance diagrams. Because of the horizontal motion of the footpoints, neighboring swirls of opposite rotation may enter and leave the time–distance slit, giving the impression of an oscillatory torsional motion. Thus, when performing this type of analysis, care must be taken that the very same perturbation is traced along the distance slit in time. This effect could possibly explain the oscillatory, torsional motion detected by Shelyag et al. (2013) from time–distance diagrams of the vorticity in the top region of their simulation domain.

In conclusion, we interpret complex rotational motions above large footpoints as the superposition of multiple, unidirectional swirls, which are caused by separate, upwardly propagating, torsional Alfvén pulses. They, in turn, originate in the dynamics of substructure elements of the magnetic footpoint.

Appendix C: Supplementary swirls

Fig. C.1. Time instant t = 5800 s showing (left) the continuum intensity with the typical granulation pattern of the solar surface and (right) the bin-5 intensity I5 as a proxy for a chromospheric line core intensity. White circles indicate the location of some of the swirl events listed in Table C.1. An animation of this figure is available online.

Fig. C.2. Time sequence of a superposition event from t = 5680 s to t = 5840 s. From top to bottom row: vertical component of the magnetic field Bz at z = 0 km, continuum intensity I, vertical velocity vz at z = 700 km, vertical component of the swirling vector λz at z = 700 km, vertical component of the magnetic swirling vector ${\lambda }_z^B$ at z = 700 km, and the bin-5 intensity I5. Maps of λz and ${\lambda }_z^B$ also show the streamlines of the velocity field and the magnetic field projected into the horizontal plane at z = 700 km, respectively.

Fig. C.3. Time sequence of a single swirl event from t = 5780 s to t = 5940 s. From the top to the bottom row: vertical component of the magnetic field Bz at z = 0 km, continuum intensity I, vertical velocity vz at z = 700 km, vertical component of the swirling vector λz at z = 700 km, vertical component of the magnetic swirling vector ${\lambda }_z^B$ at z = 700 km, and the bin-5 intensity I5. Maps of λz and ${\lambda }_z^B$ also show the streamlines of the velocity field and the magnetic field projected into the horizontal plane at z = 700 km, respectively.

Table C.1 lists nine upwardly propagating swirls found in the high-cadence simulation run described in Sect. 3. This is a non-exhaustive list of events in this time series: We have not attempted to create nor to use an automatized algorithm for detecting swirls and we discard swirls that do not propagate close to the vertical direction. The detection has been carried out manually by looking for chromospheric swirls in the bin-5 intensity I₅, over-densities in swirling strength, concentric streamlines in the horizontal velocity field, or vertical propagation of disturbances in magnetic flux concentrations near the τ₅₀₀ = 1 surface. All the events listed in Table C.1 are unidirectional, upwardly propagating swirling motions, originating from magnetic footpoints near the τ₅₀₀ = 1 surface and possessing the characteristics of torsional Alfvén pulses.

Of the nine events of Table C.1, f and c were already presented in Figs. 4 and B.1, respectively. Another two events, b and e, are depicted in Figs. C.2 and C.3, respectively. For the rest of events coordinates, time interval, and approximate size are given. The time interval does not necessarily reflect the lifetime of a swirl because it can be affected by a superposition event. Instead, it corresponds to the interval in which the event can be traced in the bin-5 intensity, hence, the interval in which it can be seen in chromospheric layers. The events of Table C.1 are marked with a circle in the animation of Fig. C.1 available online. Figure C.1 represents the instant t = 5800 s of that animation. The left panel shows the continuum intensity, the right panel the bin-5 intensity I₅. The chosen time instant contains five of the nine events of Table C.1.

All Tables

All Figures

	Fig. 1. Torsional Alfvén wave propagating in the direction k. Perturbations in the magnetic field, b, and in the plasma velocity field, v, are normal to the plane spanned by k and B0.
In the text

Fig. 2. Time instant of the swirling strength λ in a slice of dimension 9.6 × 2.6 × 2.8 Mm3 of the whole physical domain of the numerical simulation. The red contour indicates the surface of optical depth τ500 = 1. The yellow contour corresponds to the isosurface of plasma-β = 1. Swirls with a period larger than 10 min (λ < 2.09 × 10−2 Hz) are not shown and the values λ > 1 Hz are saturated. This rendering has been produced with ParaView (Ahrens et al. 2005).

In the text

Fig. 3. Statistical properties of swirling motions and mean physical quantities as a function of height z. (a) Spearman’s correlation coefficient r between the swirling strength λ and the magnetic field strength |B|, r(λ, |B|) in blue, and between the swirling strength λ and the magnetic swirling strength λB, r(λ, λB) in red. The green curve represents r(λ, |B|) calculated in regions in which β < 1 and above z = 0 only. The curves are means over eleven different time instants with a cadence of one minute and the shaded areas correspond to the standard deviations obtained from the temporal variations. Data points with λ < 2.09 × 10−2 Hz are excluded. (b) Mean magnetic field strength ⟨|B|⟩ (red) and mean swirling strength ⟨λ⟩. The dark blue and red curves refer to the same MHD time instants as in (a), while the light blue curve refers to 14 time instants of a purely hydrodynamical simulation.

In the text

Fig. 4. Time sequence of a single swirl event from t = 5760 s to t = 6160 s. From top to bottom row: vertical component of the magnetic field Bz at z = 0 km, continuum intensity I, vertical velocity vz at z = 700 km, vertical component of the swirling vector λz at z = 700 km, vertical component of the magnetic swirling vector ${\lambda }_z^B$ at z = 700 km, and the bin-5 intensity I5. Maps of λz and ${\lambda }_z^B$ also show the streamlines of the velocity field and the magnetic field projected into the horizontal plane at z = 700 km, respectively.

In the text

Fig. 5. Time sequence of the swirling strength distribution. Each panel represents a bi-dimensional histogram normalized to the maximum density fraction at each height level z. The bin sizes are Δz = 10 km and Δλ = 3.62 × 10−3 Hz. These histograms refer to stacks of quadratic, plane-parallel cross sections of 1000 km side length centered on the swirl event of Fig. 4. Data points with λ < 2.09 × 10−2 Hz are excluded. The black arrows point to the upward propagation of a local peak of swirling strength.

In the text

Fig. 6. Time–distance diagrams of (a) the swirling strength λ, (b) the vertical component of the swirling vector λz, and (c) the vertical component of the magnetic swirling vector ${\lambda }_z^B$ for the swirl event shown in Fig. 4. Values at each time step and height are averaged over a finite horizontal plane of 150 km side length. The paths of two test particles moving with Alfvén plus bulk speed (green) and with sound plus bulk speed (red) along the swirl are shown in panel b.

In the text

Fig. 7. Time–distance diagram of (a) the swirling strength λ and of the vertical component of the following quantities: (b) swirling vector λz, (c) magnetic swirling vector ${\lambda }_z^B$, (d) partial derivative of the swirling vector ∂λz/∂t, (e) advection term $T_z^{\text{adv}}={[(\mathit{v}\cdot \nabla )\mathit{\lambda }]}_z$, (f) stretching term $T_z^1$, (g) hydrodynamical baroclinic term $T_z^2$, (h) magnetic baroclinic term $T_z^3$, (i) sum of baroclinic terms $T_z^{\text{bar}}=T_z^2+T_z^3$, and (j) magnetic tension term $T_z^4$. The black dashed contours correspond to the λz = −0.1 Hz contours, which highlight the propagation of the Alfvén pulse. All plots are relative to the swirl event shown in Fig. 4 and are derived from a fixed, one pixel wide vertical slit through the center of that swirl.

In the text

	Fig. 8. Time–distance diagrams of the vertical component of the Poynting flux vector Sz (a) and of the two terms that compose it, $S_z^{\text{v}}$ (b) and $S_z^{\text{h}}$ (c), for the swirl event shown in Fig. 4. Values at each time step and height level are averages over a finite horizontal plane of 150 km side length.
In the text

	Fig. 9. Vertical component of the Poynting flux vector Sz of the swirl event presented in Fig. 4 in a horizontal section at z = 700 km. Arrows indicate the velocity vector field projected into the horizontal plane. Their length is proportional to the magnitude of the horizontal flow.
In the text

	Fig. 10. Vertical component of the Poynting flux vector Sz over the entire computational domain at z = 700 km and t = 5800 s. Black circles indicate the location of some of the swirl events listed in Table C.1 also shown in Fig. C.1.
In the text

Fig. 11. Mean vertical component of the Poynting flux vector Sz (blue) and of the terms $S_z^{\text{h}}$ (green) and $S_z^{\text{v}}$ (orange) as a function of height z. The spatial mean is obtained by averaging over the horizontal cross sections of the full computational domain of 9.6 × 9.6 Mm2, while the temporal mean is computed from 21 different time instants spread over the time period of the high-cadence series of 15 min. The shaded areas correspond to the standard deviation obtained from the temporal variation.

In the text

Fig. A.1. Volume rendering of the swirling strength seen from the side for three different simulations at arbitrary time instants: (top) MHD simulation with an initial vertical, homogeneous magnetic field of 200 G; (middle) MHD simulation with an initial field of 50 G; (bottom) purely hydrodynamical simulation. The red surface indicates the τ500 = 1 surface. Swirls with a period larger than 10 min (λ < 2.09 × 10−2 Hz) have been neglected and values above 1 Hz are saturated.

In the text

Fig. A.2. Distribution of the inclination cos(θ) of the swirling strength vector as a function of altitude z for three different simulations: (a) MHD simulation with an initial vertical, homogeneous magnetic field of 200 G, (b) MHD simulation with an initial field of 50 G, and (c) purely hydrodynamical simulation. The inclination θ refers to the angle between the swirl axis and the vertical direction. Horizontally oriented swirls are binned close to cos(θ) = 0, while vertically oriented swirls are close to cos(θ) = 1. The maximum value of swirling strength λmax at each height z, for each simulation, and with a minimum of 10 occurrences per bin, is shown in panel d. The occurrences are means over 5 different time instants of 960 × 960 × 280 data points for each simulation. The bin sizes are Δz = 10 km and Δcos(θ) = 3.62 × 10−3. Swirls with λ < 2.09 × 10−2 Hz are neglected.

In the text

Fig. B.1. Time sequence of a superposition event from t = 5750 s to t = 5830 s. From top to bottom row: vertical component of the magnetic field Bz at z = 0 km, continuum intensity I, vertical velocity vz at z = 700 km, vertical component of the swirling vector λz at z = 700 km, vertical component of the magnetic swirling vector ${\lambda }_z^B$ at z = 700 km, and the bin-5 intensity I5. Maps of λz and ${\lambda }_z^B$ also show the streamlines of the velocity field and the magnetic field projected into the horizontal plane at z = 700 km, respectively.

In the text

	Fig. B.2. Time sequence of the absolute magnetic field strength at z = 0 (gray scale) in the same field of view as shown in Fig. B.1. Plotted on top of the field strength is the vertical component of the swirling vector λz in regions where λ > 10−3 Hz and λB > 10−6 G cm−1. Red indicates counterclockwise and blue clockwise rotation, that is, positive and negative values of λz, respectively.
In the text

	Fig. C.1. Time instant t = 5800 s showing (left) the continuum intensity with the typical granulation pattern of the solar surface and (right) the bin-5 intensity I5 as a proxy for a chromospheric line core intensity. White circles indicate the location of some of the swirl events listed in Table C.1. An animation of this figure is available online.
In the text

Fig. C.2. Time sequence of a superposition event from t = 5680 s to t = 5840 s. From top to bottom row: vertical component of the magnetic field Bz at z = 0 km, continuum intensity I, vertical velocity vz at z = 700 km, vertical component of the swirling vector λz at z = 700 km, vertical component of the magnetic swirling vector ${\lambda }_z^B$ at z = 700 km, and the bin-5 intensity I5. Maps of λz and ${\lambda }_z^B$ also show the streamlines of the velocity field and the magnetic field projected into the horizontal plane at z = 700 km, respectively.

In the text

Fig. C.3. Time sequence of a single swirl event from t = 5780 s to t = 5940 s. From the top to the bottom row: vertical component of the magnetic field Bz at z = 0 km, continuum intensity I, vertical velocity vz at z = 700 km, vertical component of the swirling vector λz at z = 700 km, vertical component of the magnetic swirling vector ${\lambda }_z^B$ at z = 700 km, and the bin-5 intensity I5. Maps of λz and ${\lambda }_z^B$ also show the streamlines of the velocity field and the magnetic field projected into the horizontal plane at z = 700 km, respectively.

In the text