Issue 
A&A
Volume 639, July 2020



Article Number  A118  
Number of page(s)  12  
Section  The Sun and the Heliosphere  
DOI  https://doi.org/10.1051/00046361/202038060  
Published online  21 July 2020 
Vortices evolution in the solar atmosphere
A dynamical equation for the swirling strength
^{1}
Istituto Ricerche Solari Locarno (IRSOL), Via Patocchi 57, 6605 LocarnoMonti, Switzerland
^{2}
Center for Theoretical Astrophysics and Cosmology, Institute for Computational Science (ICS), University of Zurich, Winterthurerstrasse 190, 8057 Zürich, Switzerland
email: jcanivete@ics.uzh.ch
^{3}
LeibnizInstitut für Sonnenphysik (KIS), Schöneckstrasse 6, 79104 Freiburg i.Br., Germany
Received:
31
March
2020
Accepted:
20
May
2020
Aims. We study vortex dynamics in the solar atmosphere by employing and deriving the analytical evolution equations of two vortex identification criteria.
Methods. The two criteria used are vorticity and the swirling strength. Vorticity can be biased in the presence of shear flows, but its dynamical equation is well known; the swirling strength is a more precise criterion for the identification of vortical flows, but its evolution equation is not known yet. Therefore, we explore the possibility of deriving a dynamical equation for the swirling strength. We then apply the two equations to analyze radiative magnetohydrodynamical simulations of the solar atmosphere produced with the CO^{5}BOLD code.
Results. We present a detailed review of the swirling strength criterion and the mathematical derivation of its evolution equation. This equation did not exist in the literature before and it constitutes a novel tool that is suitable for the analysis of a wide range of problems in (magneto)hydrodynamics. By applying this equation to numerical models, we find that hydrodynamical and magnetic baroclinicities are the driving physical processes responsible for vortex generation in the convection zone and the photosphere. Higher up in the chromosphere, the magnetic terms alone dominate. Moreover, we find that the swirling strength is produced at small scales in a chaotic fashion, especially inside magnetic flux concentrations.
Conclusions. The swirling strength represents an appropriate criterion for the identification of vortices in turbulent flows, such as those in the solar atmosphere. Moreover, its evolution equation, which is derived in this paper, is pivotal for obtaining precise information about the dynamics of these vortices and the physical mechanisms responsible for their production and evolution. Since this equation is available, the swirling strength is now the ideal quantity to study the dynamics of vortices in (magneto)hydrodynamics.
Key words: magnetohydrodynamics / Sun: atmosphere / Sun: magnetic fields / turbulence
© ESO 2020
1. Introduction
Vortical motions in the solar atmosphere have received increasing attention in recent years. While a single vortical photospheric flow of multigranular spatial size was first observed by Brandt et al. (1988), Bonet et al. (2008) and Balmaceda et al. (2010) found small photospheric swirls of integranular scale by tracking magnetic bright points on their way to becoming engulfed by a downdraft. Subsequently, Bonet et al. (2010) and Vargas Domínguez et al. (2011) detected vortical motions of intergranular scale relying on the methods of local correlation tracking (LCT). Smallscale photospheric vortical flows were also reported, for example, by Attie et al. (2009) and Manso Sainz et al. (2011). In the chromosphere, WedemeyerBöhm & van der Rouppe Voort (2009) observed narrow rotating rings and ring fragments with the CRISP Fabry Pérot System of the Swedish Solar Telescope (SST) by recording narrow band filtergrams in the line core of Ca II 854.2 nm of a quite Sun region inside a coronal hole. The rings had a diameter of approximately 2 arcsec, corresponding to the size of a super large terrestrial hurricane. More recently, Tziotziou et al. (2018, 2019) found a persistent chromospheric swirl with a diameter of 6 arcsec and a lifetime of at least 1.7 h.
Subsequent observations and numerical simulations by WedemeyerBöhm et al. (2012) led to the conclusion that chromospheric swirls are the observational signatures of rotating coherent magnetic structures, which also bear imprints in the transition zone and corona. In particular, the simulations revealed that the chromospheric plasma is dragged by rotating “frozenin” magnetic field structures, which have their footpoints in the intergranular lanes of the photosphere and upper convection zone (Wedemeyer & Steiner 2014). Given that their shape resembles Earth’s atmosphere phenomena, they are usually referred to as (solar) “magnetic tornadoes” (Wedemeyer et al. 2013). In contrast to that scenario, Shelyag et al. (2013) argue that the photospheric vortexlike structures observed in numerical simulations are in reality torsional Alfvén waves propagating along the magnetic flux tubes which do not exhibit a tornadolike behavior.
One interesting aspect of magnetic tornadoes is that, by coupling the convection zone and the photosphere to the chromosphere and corona, they could act as channels for the transport of energy into the upper layers of the solar atmosphere, thus contributing to their heating. This idea has been proposed by WedemeyerBöhm et al. (2012), as they estimated through numerical simulations that magnetic tornadoes could carry a mean positive Poynting flux in the vertical direction of 440 W m^{−2} through torsional Alfvén waves.
It is not yet clear how magnetic tornadoes form. One possible solution is given by purely hydrodynamical vortices, which are observable as nonmagnetic bright points in simulations (Calvo et al. 2016), which form by the conservation of angular momentum as the plasma sinks in the downdrafts of photospheric intergranular lanes. This “bathtub” effect would capture magnetic fields that are frozen in the accreting plasma, making it rotate (WedemeyerBöhm et al. 2012; Wedemeyer & Steiner 2014).
One way of investigating the true nature of vortex flows is to study their dynamical generation in numerical simulations. Stein & Nordlund (1998) have shown with magnetic fieldfree hydrodynamical simulations that the generation of vorticity in the upper convection zone is principally due to baroclinicity, while Shelyag et al. (2011, 2012) used radiative magnetohydrodynamical (MHD) simulations to demonstrate that intergranular photospheric vortices are more efficiently produced by magnetic tension effects. All of these studies rely on vorticity as the criterion to identify the swirling motions in the simulations and on the vorticity equation to study their dynamics. However, it is well known that vorticity is not the optimal quantity to identify swirls in turbulent flows as it cannot differentiate between actual vortices and shear flows. This could lead to a misinterpretation of the numerical simulations and therefore to wrong conclusions.
An unambiguous vortex identification criterion is still an open problem in fluid mechanics (see, e.g., Kolář 2007; Liu et al. 2019), and numerous methods have been proposed and used. A promising candidate for solar atmospheric applications is the “swirling strength” (Zhou et al. 1999). This quantity, which is based on the eigenanalysis of the velocity gradient tensor, does not detect shear flows, is applicable to compressible hydrodynamics (Kolář 2009), and has already been used in the context of solar vortices by Moll et al. (2011) and Kato & Wedemeyer (2017). Nevertheless, a necessary ingredient for the study of the dynamics of vortical flows is the possibility to derive an analytical evolution equation for the criteria one seeks to use. The evolution equation for the vorticity is well known, easy to derive, and it is in fact used in the literature. Instead, there is no reference for an equation describing the time evolution of the swirling strength. Such an equation is necessary to understand whether the results obtained with the vorticity equation are biased by the ubiquitous shear flows in the solar atmosphere, or not.
In this paper, we review in detail the properties of the swirling strength criterion and we derive a new analytical equation for its evolution, the “swirling equation”. This equation is a new tool in the general context of vortices in (M)HD, which we apply to the particular case of numerical simulations of the solar atmosphere. To provide more details, we compare the vortex generation given by the vorticity equation and the swirling equation at different heights in the atmosphere, finding substantial differences, in particular in the upper convection zone. Moreover, we show a typical configuration of each one of the terms appearing in the swirling equation in the vicinity of a magnetic flux concentration, in order to gain insights into the physical processes at the origin of the vortices in these regions.
The paper is organized as follows: in Sect. 2 we present and discuss the two criteria, vorticity and swirling strength, while in Sect. 3 we derive in detail the swirling equation. In Sect. 4 we firstly study two simple toymodels in order to better understand the advantages of the swirling strength, then we use its evolution equation as a tool to analyze vortex generation in the numerical simulations. Finally, we summarize and conclude in Sect. 5.
2. Vortex identification criteria
In order to study the evolution of smallscale vortices in the solar atmosphere, we need a dynamical description of the criteria chosen for their identification. For this purpose, we start from the momentum equation of the ideal MHD system,
$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\mathit{v}=\frac{1}{\rho}\mathrm{\nabla}({p}_{\mathrm{g}}+{p}_{\mathrm{m}})+\frac{1}{\rho}(\mathit{B}\xb7\mathrm{\nabla})\mathit{B}\mathrm{\nabla}\mathrm{\Phi},\end{array}$$(1)
where v is the velocity vector, ρ is the plasma density, B is the magnetic field^{1} and Φ is a potential of conservative forces. We also have two pressure terms: the first one, p_{g}, is the usual atmospheric plasma pressure, while the latter represents the pressure due to the magnetic field, p_{m} = B^{2}/2. Finally, d/dt is the total material derivative, which is related to the partial time derivative ∂_{t} through
$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}={\partial}_{t}+\mathit{v}\xb7\mathrm{\nabla}.\end{array}$$(2)
From Eq. (1) we can derive dynamical equations of quantities which are sensible to the presence of rotational motions and can therefore be used as criteria for vortex identification.
Despite the very intuitive idea of vortex, its mathematical definition and, consequently, the choice of adequate criteria, reveals to be a complicated task, which is still under debate in fluid mechanics researches. Multiple criteria have been proposed in the past and new ones are yet to come. Between the most used in solar physics, we mention here the λ_{2}criterion (Jeong & Hussain 1995), the Γ detection functions (Graftieaux et al. 2001), and the Lagrangian Averaged Vorticity Deviation (LAVD) technique (Haller et al. 2016). However, for the purposes of this paper, we are only interested in two methods: the wellknown vorticity ω and the swirling strength λ (Zhou et al. 1999).
2.1. Vorticity
Vorticity is possibly the most common criterion in vortex identification. It is intuitively defined as the curl of the velocity field, ω ≔ ∇ × v. The direction of the vector indicates the rotation axis and the orientation of the vortex through the righthand rule. Moreover, the module is directly linked to the rotation period of a rigid body, T = 4π/ω. An interesting feature of vorticity is its evolution equation, which can be derived by taking the curl of Eq. (1), 3pt
$$\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\mathit{\omega}& =\stackrel{{\textstyle {T}_{\mathit{\omega}}^{0}}}{\stackrel{\u23de}{\mathit{\omega}(\mathrm{\nabla}\xb7\mathit{v})}}+\stackrel{{\textstyle {T}_{\mathit{\omega}}^{1}}}{\stackrel{\u23de}{(\mathit{\omega}\xb7\mathrm{\nabla})\mathit{v}}}\stackrel{{\textstyle {T}_{\mathit{\omega}}^{2}}}{\stackrel{\u23de}{\mathrm{\nabla}\frac{1}{\rho}\times \mathrm{\nabla}{p}_{\mathrm{g}}}}\hfill \\ \hfill & \underset{{\textstyle {T}_{\mathit{\omega}}^{3}}}{\underset{\u23df}{\mathrm{\nabla}\frac{1}{\rho}\times [\mathrm{\nabla}{p}_{\mathrm{m}}(\mathit{B}\xb7\mathrm{\nabla})\mathit{B}]}}+\underset{{\textstyle {T}_{\mathit{\omega}}^{4}}}{\underset{\u23df}{\frac{1}{\rho}\mathrm{\nabla}\times (\mathit{B}\xb7\mathrm{\nabla})\mathit{B}}},\hfill \end{array}$$(3)
and which is usually called the “vorticity equation”.
From Eq. (3), one can identify the physical mechanisms related to the generation or destruction of vorticity in MHD. Following Shelyag et al. (2011), we label different terms of the vorticity equation according to their physical meaning. The first two terms, ${T}_{\mathit{\omega}}^{0}$ and ${T}_{\mathit{\omega}}^{1}$, are directly related to vorticity and are called respectively vortex stretching and vortex tilting terms. We would like to point out that ${T}_{\mathit{\omega}}^{0}$ is not present in Shelyag et al. (2011) since we adopt a different prescription for the lefthand side of the equation. The term ${T}_{\mathit{\omega}}^{2}$ is the hydrodynamic baroclinic term, which is responsible for the generation of vorticity when the gradients of density and pressure are not parallel. Finally, we have two magnetic terms, ${T}_{\mathit{\omega}}^{3}$ and ${T}_{\mathit{\omega}}^{4}$. The first is the magnetic equivalent of ${T}_{\mathit{\omega}}^{2}$ and shall be denoted as magnetic baroclinic term, while the latter is related to magnetic tension.
This approach is used in Shelyag et al. (2011, 2012) to study the physical nature of magnetic photospheric intergranular vortices. However, we argue that these results could be biased by the fact that vorticity is not a reliable criterion for the identification of swirls in turbulent flows. It is in fact known that vorticity can assume large values also in the presence of shear flows (see e.g., Jeong & Hussain 1995), which can be ubiquitous in the solar atmosphere, in particular in the intergranular space and in the chromosphere. One possible solution to overcome this problem is to adopt a more sophisticated identification method which is not affected by shear flows. In what follows, we review in detail the criterion that we consider most appropriate for this study: the swirling strength.
2.2. Swirling strength
The swirling strength criterion λ (Zhou et al. 1999) is defined through the eigenanalysis of the velocity gradient tensor 𝒰,
$$\begin{array}{c}\hfill {\mathcal{U}}_{\mathit{ij}}\u2254{\partial}_{j}{v}_{i}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\leftrightarrow \phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\mathcal{U}\equiv \left[\begin{array}{ccc}{\partial}_{x}{v}_{x}& {\partial}_{y}{v}_{x}& {\partial}_{z}{v}_{x}\\ {\partial}_{x}{v}_{y}& {\partial}_{y}{v}_{y}& {\partial}_{z}{v}_{y}\\ {\partial}_{x}{v}_{z}& {\partial}_{y}{v}_{z}& {\partial}_{z}{v}_{z}\end{array}\right].\end{array}$$(4)
This matrix is invariant under Galilean transformations and in onetoone correspondence with vorticity via its antisymmetrized version: ω ↔ 𝒰 − 𝒰^{T}. In general, when 𝒰 has three discernible eigenvalues, they are either all real or one is real and there is a pair of complex conjugates (Haimes & Kenwright 1999).
Whenever a vortex is present in the flow, the velocity gradient tensor 𝒰 at that location is diagonalizable and can be decomposed in the following form,
$$\begin{array}{c}\hfill \mathcal{U}=\underset{{\textstyle \mathcal{P}}}{\underset{\u23df}{[{\mathit{u}}_{\mathrm{r}},{\mathit{u}}_{+},{\mathit{u}}_{}]}}\underset{{\textstyle \mathrm{\Lambda}}}{\underset{\u23df}{\left[\begin{array}{ccc}{\lambda}_{\mathrm{r}}& 0& 0\\ 0& {\lambda}_{+}& 0\\ 0& 0& {\lambda}_{}\end{array}\right]}}\underset{{\textstyle {\mathcal{P}}^{1}}}{\underset{\u23df}{{[{\mathit{u}}_{\mathrm{r}},{\mathit{u}}_{+},{\mathit{u}}_{}]}^{1}}},\end{array}$$(5)
where u_{r}, u_{+}, and u_{−} are the eigenvectors of 𝒰 associated with the eigenvalues λ_{r}, λ_{+}, and λ_{−}. They form respectively the matrix of eigenvectors 𝒫 (and its inverse 𝒫^{−1}) and the matrix of the eigenvalues Λ. To provide more details, λ_{r} is a real eigenvalue, while λ_{+} and λ_{−} are complex conjugates, that is, λ_{±} = λ_{cr} ± iλ_{ci}. In particular, their imaginary part, λ_{ci}, characterizes the strength of a rigid body swirling motion through T ≔ 2π/λ_{ci}. It is therefore natural to define the swirling strength as λ_{ci}, as done, for example, in Moll et al. (2011). Nevertheless, in the present paper, we use a different convention; in fact, since our goal is to derive a dynamical equation for the swirling strength and then compare it to the vorticity equation, Eq. (3), it is necessary that both criteria return the same value for a given vortex. Therefore, we decide to define the swirling strength as λ = 2λ_{ci}. In this way, whenever the flow follows a rigid body rotation of period T, both criteria return the same value, λ = ω = 4π/T. A more detailed description of the physical interpretation of λ_{r}, λ_{cr}, and λ_{ci} can be found in Appendix A.
Another important piece of information we can obtain from Eq. (5) is the vortex axis and its orientation. In fact, the eigenvector u_{r}, associated with the real eigenvalue λ_{r}, identifies the direction of the vortex axis. However, this is not sufficient to determine the orientation of the swirl, that is, whether the flow is turning in a clockwise or counterclockwise fashion around this axis. This arbitrariness occurs because eigenvectors are defined up to a multiplicative constant only; therefore their orientation in space can be reverted. Nevertheless, there exists a simple method to ensure the physicality of its orientation, which consists in checking the handedness of the ℝ^{3} basis which is formed by the eigenvectors of 𝒰. In fact, one can prove (see Appendix A) that if the basis is lefthanded, then the rotation is counterclockwise in the direction given by u_{r}, while if it is righthanded, it turns in a clockwise fashion. In practice, we force the eigenvectors to form a lefthanded basis by multiplying u_{r} by −1 where necessary. In this way, we ensure that all the rotations are counterclockwise in the direction defined by u_{r}, as prescribed by the righthand rule.
Once the orientation of u_{r} is fixed, we can define the “swirling vector” as λ ≔ λu_{r} and therefore be able to extract the same amount of information from a flow as we can do with vorticity. However, when it comes to evolution, vorticity still represents the only viable way since its dynamical equation is known. In fact, there is no reference in the literature to a dynamical equation describing the evolution of the swirling strength and, therefore, a research on the generation of vortices using this criterion was not possible. In what follows, we derive the swirling strength evolution equation, enabling therefore the study of vortex evolution with this more sophisticated criterion.
3. The swirling equation
In order to derive an evolution equation for the swirling strength λ or, more generally, for the swirling vector λ, we seek a tensor equation for 𝒰. Then, using the diagonalization properties given in Eq. (5), we isolate the equation relative to λ_{ci} obtaining therefore a “swirling equation”.
The starting point is the momentum equation, Eq. (1), to which we apply the gradient operator ∇ via the tensor product,
$$\begin{array}{c}\hfill \mathrm{\nabla}\left(\frac{\mathrm{d}\mathit{v}}{\mathrm{d}t}\right)=\mathrm{\nabla}(\frac{1}{\rho}\mathrm{\nabla}({p}_{\mathrm{g}}+{p}_{\mathrm{m}})+\frac{1}{\rho}(\mathit{B}\xb7\mathrm{\nabla})\mathit{B}\mathrm{\nabla}\mathrm{\Phi}),\end{array}$$(6)
where we define the tensor product between two vectors a and b, 𝒞 = ab, as the matrix 𝒞 whose components are given by 𝒞_{ij} = a_{i}b_{j}. The resulting Eq. (6) is therefore a tensor equation combining the different components of the gradient and the momentum equation.
As a first step, we develop the lefthand side of Eq. (6) by looking at a generic tensor component (i, j),
$$\begin{array}{cc}\hfill {\partial}_{i}\left(\frac{\mathrm{d}{v}_{j}}{\mathrm{d}t}\right)& ={\partial}_{i}({\partial}_{t}{v}_{j}+\sum _{k}({v}_{k}{\partial}_{k})\phantom{\rule{0.166667em}{0ex}}{v}_{j})\hfill \\ \hfill & ={\partial}_{t}\left({\partial}_{i}{v}_{j}\right)+\sum _{k}({v}_{k}{\partial}_{k})\phantom{\rule{0.166667em}{0ex}}{\partial}_{i}{v}_{j}+\sum _{k}({\partial}_{i}{v}_{k})\phantom{\rule{0.166667em}{0ex}}({\partial}_{k}{v}_{j})\hfill \\ \hfill & ={\partial}_{t}{\mathcal{U}}_{\mathit{ij}}^{\mathrm{T}}+\sum _{k}({v}_{k}{\partial}_{k})\phantom{\rule{0.166667em}{0ex}}{\mathcal{U}}_{\mathit{ij}}^{\mathrm{T}}+\sum _{k}{\mathcal{U}}_{\mathit{kj}}^{\mathrm{T}}{\mathcal{U}}_{\mathit{ik}}^{\mathrm{T}}\hfill \\ \hfill & =\frac{\mathrm{d}}{\mathrm{d}t}{\mathcal{U}}_{\mathit{ij}}^{\mathrm{T}}+{({\mathcal{U}}^{2})}_{\mathit{ij}}^{\mathrm{T}},\hfill \end{array}$$(7)
where we use the explicit form of the material derivative, Eq. (2), and the definition of the velocity gradient tensor 𝒰 given in Eq. (4). In this way, one can replace the velocity field in terms of 𝒰. Therefore, Eq. (6) can be rewritten in the following form,
$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}{\mathcal{U}}^{\mathrm{T}}={({\mathcal{U}}^{2})}^{\mathrm{T}}+{\mathcal{M}}^{\mathrm{T}},\end{array}$$
and by simply transposing,
$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\mathcal{U}={\mathcal{U}}^{2}+\mathcal{M},\end{array}$$(8)
where ℳ is the transpose of the tensor obtained in the righthand side of Eq. (6) by applying the gradient to the source terms. To be more specific, a generic component of ℳ is given by
$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathit{ij}}={\partial}_{j}(\frac{1}{\rho}{\partial}_{i}({p}_{\mathrm{g}}+{p}_{\mathrm{m}})+\frac{1}{\rho}\left(\sum _{k}{B}_{k}{\partial}_{k}\right){B}_{i}{\partial}_{i}\mathrm{\Phi}).\end{array}$$(9)
Successively, we use the diagonalization properties of 𝒰 in order to transform Eq. (8) into an evolution equation for Λ. In particular, by inverting Eq. (5) one finds the usual transformation law for matrices,
$$\begin{array}{c}\hfill \mathrm{\Lambda}={\mathcal{P}}^{1}\mathcal{U}\mathcal{P},\end{array}$$(10)
which leads to the following relations,
$$\begin{array}{c}\hfill \mathcal{P}\mathrm{\Lambda}=\mathcal{U}\mathcal{P},\phantom{\rule{1em}{0ex}}\mathrm{\Lambda}{\mathcal{P}}^{1}={\mathcal{P}}^{1}\mathcal{U}.\end{array}$$(11)
Moreover, given that 𝒫𝒫^{−1} = ℐ, where ℐ is the identity matrix, and that dℐ/dt = 0, one can easily prove the following relation,
$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}({\mathcal{P}}^{1})\mathcal{P}={\mathcal{P}}^{1}\frac{\mathrm{d}}{\mathrm{d}t}(\mathcal{P}).\end{array}$$(12)
Then, multiplying Eq. (8) by 𝒫^{−1} on the left and by 𝒫 on the right we obtain,
$$\begin{array}{c}\hfill {\mathcal{P}}^{1}(\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{U})\mathcal{P}={\mathcal{P}}^{1}{\mathcal{U}}^{2}\mathcal{P}+{\mathcal{P}}^{1}\mathcal{M}\mathcal{P}.\end{array}$$(13)
Using Eqs. (10)–(12), for the lefthand side of Eq. (13), we find,
$$\begin{array}{cc}\hfill {\mathcal{P}}^{1}(\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{U})\mathcal{P}\phantom{\rule{0.166667em}{0ex}}=& \phantom{\rule{0.166667em}{0ex}}\frac{\mathrm{d}}{\mathrm{d}t}\left({\mathcal{P}}^{1}\mathcal{U}\mathcal{P}\right)\frac{\mathrm{d}}{\mathrm{d}t}\left({\mathcal{P}}^{1}\right)\mathcal{U}\mathcal{P}{\mathcal{P}}^{1}\mathcal{U}\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{P},\hfill \\ \hfill =& \phantom{\rule{0.166667em}{0ex}}\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{\Lambda}\frac{\mathrm{d}}{\mathrm{d}t}\left({\mathcal{P}}^{1}\right)\mathcal{P}\mathrm{\Lambda}\mathrm{\Lambda}{\mathcal{P}}^{1}\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{P},\hfill \\ \hfill =& \phantom{\rule{0.166667em}{0ex}}\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{\Lambda}\mathrm{\Lambda}{\mathcal{P}}^{1}\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{P}+{\mathcal{P}}^{1}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathcal{P}\right)\mathrm{\Lambda},\hfill \\ \hfill =& \phantom{\rule{0.166667em}{0ex}}\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{\Lambda}[\mathrm{\Lambda},{\mathcal{P}}^{1}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathcal{P}\right)],\hfill \end{array}$$
where the square brackets indicate the commutator of two matrices, [𝒜,ℬ] ≔ 𝒜ℬ − ℬ𝒜. For the righthand side instead, we get,
$$\begin{array}{cc}\hfill {\mathcal{P}}^{1}{\mathcal{U}}^{2}\mathcal{P}+{\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\phantom{\rule{0.166667em}{0ex}}=& \phantom{\rule{0.166667em}{0ex}}{\mathcal{P}}^{1}\mathcal{U}\mathcal{I}\mathcal{U}\mathcal{P}+{\mathcal{P}}^{1}\mathcal{M}\mathcal{P},\hfill \\ \hfill =& \phantom{\rule{0.166667em}{0ex}}{\mathcal{P}}^{1}\mathcal{U}\mathcal{P}\phantom{\rule{0.166667em}{0ex}}{\mathcal{P}}^{1}\mathcal{U}\mathcal{P}+{\mathcal{P}}^{1}\mathcal{M}\mathcal{P},\hfill \\ \hfill =& \phantom{\rule{0.166667em}{0ex}}{\mathrm{\Lambda}}^{2}+{\mathcal{P}}^{1}\mathcal{M}\mathcal{P}.\hfill \end{array}$$
Combining back together the two sides of Eq. (13), we find an equation which couples the evolution of the eigenvalue matrix Λ with the eigenvector matrix 𝒫,
$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\mathrm{\Lambda}=[\mathrm{\Lambda},{\mathcal{P}}^{1}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathcal{P}\right)]{\mathrm{\Lambda}}^{2}+{\mathcal{P}}^{1}\mathcal{M}\mathcal{P}.\end{array}$$(14)
However, this coupling is just apparent since the commutator between a diagonal matrix and any other matrix results in an hollow matrix. The interested reader can refer to Appendix B for a proof. Consequently, the commutator in Eq. (14) between Λ and 𝒫^{−1}d𝒫/dt is hollow, thus it contributes to the tensor equations only in the offdiagonal terms. Since the total derivative of Λ is restricted to the diagonal terms of Eq. (14), the evolution of Λ and 𝒫 is decoupled and we can easily distinguish two sets of equations relative to each quantity. For the eigenvalues matrix Λ, we consider alone the equations relative to the diagonal terms,
$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{\Lambda}}_{\mathit{ii}}={\left({\mathrm{\Lambda}}^{2}\right)}_{\mathit{ii}}+{\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{\mathit{ii}},\end{array}$$(15)
where the index i refers to the ith eigenvalue of Λ, while for the eigenvector matrix 𝒫, we take into account the offdiagonal terms,
$$\begin{array}{c}\hfill {[\mathrm{\Lambda},{\mathcal{P}}^{1}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathcal{P}\right)]}_{\mathit{ij}}+{\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{\mathit{ij}}=0,\phantom{\rule{1em}{0ex}}i\ne j.\end{array}$$(16)
Since we are looking for a dynamical equation for the swirling vector, λ = λu_{r}, we are left with the extraction of the evolution equations for the swirling strength λ and the vortex axis u_{r} from Eqs. (15) and (16).
The evolution of the three eigenvalues present in Λ is dictated by the relative equation given in Eq. (15). Because of the definition of swirling strength, we can restrict ourselves to the one relative to the eigenvalue λ_{+} = λ_{cr} + iλ_{ci},
$$\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}{\lambda}_{+}& ={\lambda}_{+}^{2}+{\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{22}\phantom{\rule{0.166667em}{0ex}},\hfill \\ \hfill \frac{\mathrm{d}}{\mathrm{d}t}({\lambda}_{\mathrm{cr}}+\mathrm{i}{\lambda}_{\mathrm{ci}})& ={({\lambda}_{\mathrm{cr}}+\mathrm{i}{\lambda}_{\mathrm{ci}})}^{2}+{\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{22}\phantom{\rule{0.166667em}{0ex}}.\hfill \end{array}$$
By taking only the imaginary part, we recover an equation for λ_{ci}. Then, given our convention, we formulate the swirling equation as,
$$\begin{array}{ccc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\lambda =& 2\lambda {\lambda}_{\mathrm{cr}}+2\mathrm{Im}{\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{22}\phantom{\rule{0.166667em}{0ex}},\hfill & \hfill \\ \hfill =& 2\lambda {\lambda}_{\mathrm{cr}}\hfill & \hfill {T}_{\lambda}^{1}& \hfill \\ \hfill & 2\mathrm{Im}{\left\{{\mathcal{P}}^{1}\left[\mathrm{\nabla}(\frac{1}{\rho}\mathrm{\nabla}{\mathit{p}}_{\mathrm{g}})\right]\mathcal{P}\right\}}_{22}\hfill & \hfill {T}_{\lambda}^{2}& \hfill \\ \hfill & 2\mathrm{Im}{\left\{{\mathcal{P}}^{1}[\mathrm{\nabla}(\frac{1}{\rho}\mathrm{\nabla}{\mathit{p}}_{\mathrm{m}})(\mathrm{\nabla}\frac{1}{\rho})(\mathit{B}\xb7\mathrm{\nabla})\mathit{B}]\mathcal{P}\right\}}_{22}\hfill & \hfill {T}_{\lambda}^{3}& \hfill \\ \hfill & +2\mathrm{Im}{\left\{{\mathcal{P}}^{1}[\frac{1}{\rho}\mathrm{\nabla}((\mathit{B}\xb7\mathrm{\nabla})\mathit{B})]\mathcal{P}\right\}}_{22}\hfill & \hfill {T}_{\lambda}^{4}& \hfill \\ \hfill & 2\mathrm{Im}{\left\{{\mathcal{P}}^{1}[\mathrm{\nabla}(\mathrm{\nabla}\mathrm{\Phi})]\mathcal{P}\right\}}_{22}\phantom{\rule{0.166667em}{0ex}}.\hfill & \hfill {T}_{\lambda}^{5}\end{array}$$(17)
In order to give a physical interpretation to Eq. (17), we adopt a term labeling similar to the one in the vorticity equation, Eq. (3). In particular, we try to match terms that share the same physical mechanism for vortex generation. Therefore, ${T}_{\lambda}^{1}$ can be considered as a stretching term, since it involves λ_{cr}, which indicates compression or dilatation of the flow in the rotation plane, as we show in Appendix A. The second term, ${T}_{\lambda}^{2}$, is related to the hydrodynamic baroclinicity, while ${T}_{\lambda}^{3}$ is its magnetic correspondent. The generation of swirling strength by magnetic tension is instead represented by the term ${T}_{\lambda}^{4}$. Finally, we note that the last term, ${T}_{\lambda}^{5}$, which is associated to the potential of conservative forces, has no analog in the vorticity equation. When Φ represents the gravitational potential and we can assume a slowly varying unidirectional gravitation field, as it is the case in most simulations of the solar atmosphere, ${T}_{\lambda}^{5}$ can be safely neglected. However, this term could be relevant whenever the configuration of the gravitational potential is complex and inhomogeneous.
The equation derived for the swirling strength, Eq. (17), is similar in form to the evolution equation of vorticity, Eq. (3). The main difference consists in the multiplication of each source term present in ℳ by 𝒫^{−1} and 𝒫. One can physically interpret this difference recalling that the eigenvectors that compose 𝒫 form an eigenbasis of the velocity gradient tensor 𝒰. In such a basis, 𝒰 explicitly describes the characteristics of the flow through its eigenvalues. Therefore, 𝒫^{−1}ℳ𝒫 represents a change of basis of the source terms, and the basis defined by 𝒫 describes the proper coordinate system of the vortex, where the shears and the other components of the flow, which could mislead the vortex interpretation, have been implicitly taken into account in the basis vectors. Consequently, the terms appearing in Eq. (17) are related to the generation and the evolution of true vortices only.
Finally, we note that we could have chosen to pick the equation relative to the other complex conjugate eigenvalue, λ_{−} = λ_{cr} − iλ_{ci}, for the formulation of the swirling equation. This leads to an equation which, at first sight, appears different, but in reality is equivalent to Eq. (17). A proof is given in Appendix C. Moreover, one can also derive evolution equations for the λ_{cr} and λ_{r} parameters in a similar fashion.
Given the swirling equation, we are left with deriving an equation for the evolution of the eigenvector u_{r}. To do so, we turn to Eq. (16) using the property of commutators derived in Eq. (B.1),
$$\begin{array}{cc}& {[{\mathcal{P}}^{1}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathcal{P}\right),\mathrm{\Lambda}]}_{\mathit{ij}}={\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{\mathit{ij}},\hfill \\ \hfill & ({\lambda}^{j}{\lambda}^{i}){\left({\mathcal{P}}^{1}\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{P}\right)}_{\mathit{ij}}={\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{\mathit{ij}},\hfill \\ \hfill & {\left({\mathcal{P}}^{1}\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{P}\right)}_{\mathit{ij}}=\frac{1}{({\lambda}^{j}{\lambda}^{i})}{\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{\mathit{ij}},\hfill \\ \hfill & \sum _{k}{\mathcal{P}}_{\mathit{ik}}^{1}\frac{\mathrm{d}}{\mathrm{d}t}{P}_{\mathit{kj}}=\frac{1}{({\lambda}^{j}{\lambda}^{i})}{\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{\mathit{ij}},\phantom{\rule{0.166667em}{0ex}}i\ne j,\hfill \end{array}$$(18)
where in the second equation we use the fact that [𝒜,ℬ] = −[ℬ,𝒜] and in the last equation we explicitly multiply 𝒫^{−1} with d𝒫/dt and we recall that this equation is valid only for i ≠ j.
Equation (18) consists of six evolution equations, for all combinations of i and j such that i ≠ j, where j indexes the evolving eigenvector. Since we are interested in the evolution of u_{r}, we concentrate on the j = 1 components of Eq. (18),
$$\begin{array}{c}\hfill \sum _{k}{\mathcal{P}}_{\mathit{ik}}^{1}\frac{\mathrm{d}}{\mathrm{d}t}{\left({\mathit{u}}_{\mathrm{r}}\right)}_{k}=\frac{1}{({\lambda}_{\mathrm{r}}{\lambda}^{i})}{\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{i1}\phantom{\rule{0.166667em}{0ex}},\phantom{\rule{1em}{0ex}}i=2,3,\end{array}$$(19)
where we use the fact that λ^{1} = λ_{r} in our diagonalization convention of the velocity gradient tensor 𝒰. Equation (19) consists of two equations, coupling the evolution of the three components of the eigenvector u_{r}. We miss one equation in order to solve the system. One simple extra constraint comes from the normalization of u_{r}, which we fix to one, ${({\mathit{u}}_{\text{r}})}_{x}^{2}+{({\mathit{u}}_{\text{r}})}_{y}^{2}+{({\mathit{u}}_{\text{r}})}_{z}^{2}=1$, which can be expressed as
$$\begin{array}{c}\hfill {({\mathit{u}}_{\mathrm{r}})}_{x}\frac{\mathrm{d}}{\mathrm{d}t}{({\mathit{u}}_{\mathrm{r}})}_{x}+{({\mathit{u}}_{\mathrm{r}})}_{y}\frac{\mathrm{d}}{\mathrm{d}t}{({\mathit{u}}_{\mathrm{r}})}_{y}+{({\mathit{u}}_{\mathrm{r}})}_{z}\frac{\mathrm{d}}{\mathrm{d}t}{({\mathit{u}}_{\mathrm{r}})}_{z}=0.\end{array}$$
Thus, we have now a set of three equations that describe the dynamical evolution of the vortex axis.
Combining the results obtained, we can describe the evolution of the swirling vector λ by
$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\mathit{\lambda}={\mathit{u}}_{\mathrm{r}}\frac{\mathrm{d}}{\mathrm{d}t}\lambda +\lambda \frac{\mathrm{d}}{\mathrm{d}t}{\mathit{u}}_{\mathrm{r}},\end{array}$$(20)
where the first term represents the generation of swirling strength while the second characterizes the evolution of the vortex axis. We could therefore say that the term λdu_{r}/dt in Eq. (20) is the analog of the tilting term, ${T}_{\lambda}^{0}$, present in Eq. (3). However, in the present case, this term is not responsible for any vortex generation, but instead it describes the transition under constant swirling strength between components of λ when a vortex tilts its orientation.
4. Results and discussion
In this section, we compare the generation of vorticity and swirling strength, in particular with regard to their evolution equations, Eqs. (3) and (17). Firstly, we study two simple scenarios to better understand the limits of vorticity and to highlight the advantages of having an equation for the swirling strength available. Then, we apply the dynamical equations to highresolution radiative magnetohydrodynamic simulation data produced with the CO^{5}BOLD code (Freytag et al. 2012) to study the mechanisms of vortex generation. The simulation was carried out by starting from a previous, relaxed hydrodynamical model to which a uniform vertical magnetic field with a strength of 50 G was superimposed. This magnetic model was then further advanced with the MHD module of CO^{5}BOLD, which uses an HLL approximate Riemann solver (Harten et al. 1983) until relaxation of the system is achieved. The lateral boundary conditions are periodic, while at the top and at the bottom they are open under the condition that the net mass flux at the bottom boundary vanishes. The magnetic field is constrained to be vertical at both top and bottom boundaries but it can freely move in lateral directions.
The Cartesian computational domain encompasses a volume of a 9.6 × 9.6 × 2.8 Mm^{3} of the solar atmosphere, ranging from the surface layers of the convection zone to the chromosphere. The optical surface τ_{500} = 1 is approximately located in the middle of the height range. The size of the computational cells is 10 × 10 × 10 km^{3}, uniform in the entire box. The gravitational field is vertical and uniform with a constant value of log(g) = 4.44. More details on this simulation can be found in Calvo (2018). Here we use a sequence over a time period of about 10 min real time of solar evolution for which we have time instants stored every 10 s.
For the present study we consider vertical vortices only, that is, those for which the flow rotates around the zcoordinate axis, $\widehat{z}$. We are particularly interested in this kind of vortices since they potentially transport energy and mass across different layers of the solar atmosphere. Nevertheless, we would like to point out that the method is suitable to all kind of vortices.
4.1. Limits of vorticity
We start with a simple rotational flow with velocity v and corresponding vorticity ω given by,
$$\begin{array}{c}\hfill \mathit{v}=\left(\begin{array}{c}y\frac{\mathrm{\Omega}}{2}\\ x\frac{\mathrm{\Omega}}{2}\\ 0\end{array}\right),\phantom{\rule{2em}{0ex}}\mathit{\omega}=\mathrm{\nabla}\times \mathit{v}=\left(\begin{array}{c}0\\ 0\\ \mathrm{\Omega}\end{array}\right).\end{array}$$
In this example, the flow is performing a rigid, counterclockwise rotation of period T = 4π/Ω. Without loss of generality, we have chosen the axis of rotation to be along the $\widehat{z}$direction. Then, we can also compute the velocity gradient tensor 𝒰 and, by its diagonalization, find the swirling strength,
$$\begin{array}{cc}\hfill \mathcal{U}& =\left[\begin{array}{ccc}0& \frac{\mathrm{\Omega}}{2}& 0\\ \frac{\mathrm{\Omega}}{2}& 0& 0\\ 0& 0& 0\end{array}\right],\hfill \\ \hfill & =\left[\begin{array}{ccc}0& \mathrm{i}\frac{1}{\sqrt{2}}& \mathrm{i}\frac{1}{\sqrt{2}}\\ 0& \frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}\\ 1& 0& 0\end{array}\right]\left[\begin{array}{ccc}0& 0& 0\\ 0& \mathrm{i}\frac{\mathrm{\Omega}}{2}& 0\\ 0& 0& \mathrm{i}\frac{\mathrm{\Omega}}{2}\end{array}\right]\left[\begin{array}{ccc}0& 0& 1\\ \mathrm{i}\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}& 0\\ \mathrm{i}\frac{1}{\sqrt{2}}& \frac{1}{\sqrt{2}}& 0\end{array}\right],\hfill \\ \hfill & =\mathcal{P}\mathrm{\Lambda}{\mathcal{P}}^{1}.\hfill \end{array}$$(21)
Given our convention, we see that λ ≡ 2λ_{ci} = Ω. Moreover, by computing the real matrix $\stackrel{\sim}{\mathcal{P}}$ (defined in Appendix A) and its determinant, one can check that the ℝ^{3} basis in which the swirl is described is lefthanded. Thus, we conclude that the vortex axis is given by u_{r} = (0, 0, 1) and we find that, for this simple example, the swirling and the vorticity vectors are identical, λ = ω = (0, 0, Ω). This should not surprise us: indeed, for purely rotational flows, vorticity and swirling strength should be equivalent given the absence of shear flows.
At this point, we add a shear flow along the $\widehat{z}$axis proportional to the $\widehat{y}$ coordinate in such a way that the velocity and vorticity vectors are,
$$\begin{array}{c}\hfill \mathit{v}=\left(\begin{array}{c}y\frac{\mathrm{\Omega}}{2}\\ x\frac{\mathrm{\Omega}}{2}\\ y\xi \end{array}\right),\phantom{\rule{2em}{0ex}}\mathit{\omega}=\mathrm{\nabla}\times \mathit{v}=\left(\begin{array}{c}\xi \\ 0\\ \mathrm{\Omega}\end{array}\right),\end{array}$$
where ξ is the strength of the shear. We notice that the extra component in the velocity vector, v_{z} = yξ, which is not related to any kind of rotation in the flow, generates vorticity along the $\widehat{x}$axis. Computing once again the velocity gradient tensor and its diagonal decomposition we find,
$$\begin{array}{cc}\hfill \mathcal{U}& =\left[\begin{array}{ccc}0& \frac{\mathrm{\Omega}}{2}& 0\\ \frac{\mathrm{\Omega}}{2}& 0& 0\\ 0& \xi & 0\end{array}\right],\hfill \\ & =\left[\begin{array}{ccc}0& \frac{\mathrm{\Omega}}{2\xi c}& \frac{\mathrm{\Omega}}{2\xi c}\\ 0& \mathrm{i}\frac{\mathrm{\Omega}}{2\xi c}& \mathrm{i}\frac{\mathrm{\Omega}}{2\xi c}\\ 1& \frac{1}{c}& \frac{1}{c}\end{array}\right]\left[\begin{array}{ccc}0& 0& 0\\ 0& \mathrm{i}\frac{\mathrm{\Omega}}{2}& 0\\ 0& 0& \mathrm{i}\frac{\mathrm{\Omega}}{2}\end{array}\right]\left[\begin{array}{ccc}\phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\frac{2\xi}{\mathrm{\Omega}}& \phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}0& 1\\ \frac{\xi c}{\mathrm{\Omega}}& \mathrm{i}\frac{\xi c}{\mathrm{\Omega}}& 0\\ \frac{\xi c}{\mathrm{\Omega}}& \phantom{\rule{0.166667em}{0ex}}\phantom{\rule{0.166667em}{0ex}}\mathrm{i}\frac{\xi c}{\mathrm{\Omega}}& 0\end{array}\right],\hfill \end{array}$$(22)
where $c=\sqrt{1+{\mathrm{\Omega}}^{2}/(2{\xi}^{2})}$. Computing the matrix of real eigenvectors, $\stackrel{\sim}{\mathcal{P}}$, and checking the handedness of the corresponding basis according to Appendix A, one can conclude that the swirling strength vector λ = (0, 0, Ω) is not affected by the shear flow and correctly describes the properties of the vortex. One can therefore define the shear strength, ω_{sh} ≔ ω − λ, as the contribution to vorticity due to shear flows. In practice, the shear strength is a proxy for estimating how vorticity is biased by the shears present in the flow. Therefore, it can indicate the regions of the plasma where vorticity cannot be trusted for the identification of vortex flows.
These two simple examples are useful to get familiar with the criteria used to study vortices in MHD. However, when it comes to study realistic scenarios, the plasma dynamics are much more complex. Panels a and b of Fig. 1 illustrate the vertical components of vorticity ω_{z} and swirling strength λ_{z} for one particular vortex in an intergranular region with a strong vertical magnetic field. We also show in gray shades the optical surface of τ_{500} = 1. The presence of the vortex is highlighted by the streamlines (in green), which indicate rotation around the strong concentration of positive swirling strength and vorticity (red). However, we notice how the swirling strength is predominately positive and how, above the optical surface, it is mainly concentrated inside the vortex, while vorticity presents also negative filaments (blue) and looks slightly more dispersed. Considering vorticity alone, one could be lead to conclude that panel a of Fig. 1 shows two vortices counterrotating one with respect to the other. However, it is clear from panel b of Fig. 1 and from the streamlines that the negative vortex hinted by the vorticity is, in reality, a thin sheath of shear flow. As can be seen from panel c of Fig. 1, the shears are allover the atmospheric portion of interest, thus the importance of the swirling strength approach.
Fig. 1. Vertical component of (a) vorticity, (b) swirling strength, and (c) shear strength, computed in a small portion of the full simulation domain, measuring 0.6 × 0.6 × 1.8 Mm^{3} and encompassing an intergranular region with strong vertical magnetic field (B_{z} ≳ 500 G). The gray sheet displays the surface of optical depth τ_{500} = 1, while the plasma flow is rendered with green streamlines. These images were produced with the software ParaView (Ahrens et al. 2005). 
4.2. Comparing the generation of vorticity and swirling strength from their respective equations
This section investigates the generation of vortical flows in the solar atmosphere using the evolution equations of vorticity, Eq. (3), and of the swirling vector, Eq. (20). In particular, we want to study which physical process is dominant at different heights in the atmosphere, from the upper convection zone through the photosphere up to the chromosphere, and compare the results between the two approaches. For this purpose, we study the first term of Eq. (20) alone, since we are interested in the generation of the vertical component of vortices only and not in their tilting from one axis to another. In particular, u_{r}dλ/dt is related to the generation of swirling strength through Eq. (17). Therefore, in what follows, we define the vertical component of the swirling equation terms, ${T}_{i,z}^{\lambda}$, as being the $\widehat{z}$component of the generating part of Eq. (20),
$$\begin{array}{c}\hfill {({\mathit{u}}_{\mathrm{r}})}_{z}\frac{\mathrm{d}}{\mathrm{d}t}\lambda ={({\mathit{u}}_{\mathrm{r}})}_{z}\sum _{i}{T}_{i}^{\lambda}=\sum _{i}{T}_{i,z}^{\lambda}.\end{array}$$
Figure 2 shows the mean absolute values of the vertical components of the terms appearing in Eqs. (3) and (17) as a function of the height z. Moreover, we note that ${T}_{5}^{\lambda}$ is neglected in this study because the present simulation has a constant gravitational acceleration and, consequently, ∇(∇Φ) = 0.
Fig. 2. Mean absolute values of the vertical components of the terms appearing on the righthand side of the vorticity equation (left panel) and the swirling equation (right panel) as a function of z. The dotted curves represent the mean absolute value of the total production of vorticity and swirling strength in the simulation. The optical surface τ_{500} = 1 corresponds to z = 0. 
Concerning vorticity, we notice that the tilting term ${T}_{1}^{\omega}$ is the largest contributor from the convection zone up to the low photosphere, where it gets overtaken by the magnetic tension term ${T}_{4}^{\omega}$. From there on further up, magnetic effects are mainly responsible for the generation of vorticity. In particular, we see that in the chromospheric layers, the sum of all components is almost identical to the magnetic tension term alone, implying that the other terms are minor in this height range. We also notice a drop in vorticity generation in the photosphere, which can be ascertained to the steep drop in density at this height with the consequence that upflows expand. This results in a loss of rotational speed because of angular momentum conservation (Nordlund et al. 1997).
A similar study has already been worked out by Shelyag et al. (2011), who used the MURaM code (Vögler 2003; Vögler et al. 2005) to simulate the magnetized photosphere. We find very similar results to the ones presented there, apart from what concerns the terms ${T}_{1}^{\omega}$ and ${T}_{2}^{\omega}$, which seem to have been interchanged by Shelyag et al. (2011). This conclusion is also supported by the results of Rajaguru et al. (in prep.), which also show the presence of a strong tilting term in the upper convection zone and a subdominant hydrodynamical baroclinic term in simulations performed with the MURaM, CO^{5}BOLD, BIFROST (Gudiksen et al. 2011), and STAGGER (Nordlund et al. 1994; Nordlund & Galsgaard 1995) codes.
Yet, we argue that the results obtained with vorticity can be biased, since its equation also accounts for the generation of shears. Therefore, the right panel of Fig. 2 presents the same analysis using the swirling equation. In this case, we see that the stretching term, ${T}_{1}^{\lambda}$, is unimportant over the entire height range. In the convection zone, the two main contributors to the generation of swirling strength are the hydrodynamical baroclinic term ${T}_{2}^{\lambda}$ and the magnetic baroclinic term ${T}_{3}^{\lambda}$, while in the chromosphere we confirm that the magnetic terms, ${T}_{3}^{\lambda}$ and ${T}_{4}^{\lambda}$, are by far the most important ones. However, for the generation of vorticity, the effects of magnetic tension are stronger than those of magnetic pressure, while in the right panel of Fig. 2 we see that these two magnetic terms have almost identical mean contributions. In the photosphere, we observe a decrease in the total production of swirling strength for the same reason that leads to a drop in vorticity generation.
Interestingly, the mean total generation of swirling strength in the superficial layer of the convection zone and, in particular in the photosphere, is lower that the mean values of the plasma and magnetic baroclinic terms, implying that these two terms partially counterbalance their effects. This is intuitive if we think of a magnetic flux tube, where the external gas pressure must balance the sum of internal magnetic and gas pressures (Zwaan 1978). Thus, a positive gradient in magnetic pressure opposes a negative one in gas pressure, which results in opposite contributions to the baroclinic terms ${T}_{2}^{\lambda}$ and ${T}_{3}^{\lambda}$ in Eq. (17). On the other hand, this interpretation implies that the most important contributor to vortex generation in the photosphere and in the low chromosphere is the magnetic tension term, ${T}_{4}^{\lambda}$.
In conclusion, Fig. 2 tells us that there are substantial differences between the two approaches to identification of vortex generation, in particular regarding the convection zone. Obviously, these differences are caused by the generation of shears along with vortices, which is included in the righthand side terms of the vorticity equation and which can lead into erroneous analysis. Therefore, we think that the best tool to study vortex generation in the solar atmosphere is the swirling equation presented in Sect. 3. Furthermore, we can deduce that the main contributor to the production of vertical shear in the convection zone is the tilting term of vorticity, ${T}_{1}^{\omega}$. True vertical vortices, on the other hand, are mainly produced by baroclinic effects. In the photosphere and chromosphere, both methods suggest magnetic terms to be the main responsible for the generation of vortex flows. Nevertheless, the swirling equation shows that both the magnetic baroclinic and the magnetic tension terms are equally important, while according to vorticity the magnetic tension dominates alone. We can therefore infer that shears in these layers are mainly produced by magnetic tension effects.
4.3. Analysis of the swirling equation terms
We now focus on the swirling equation and see how the different terms locally contribute to the generation of vertical vortices. For this purpose, we study in more detail the photospheric vortex presented in Fig. 1 at three characteristic heights. We recall that, like in Sect. 4.2, we examine only the vertical component of the first term on the righthand side of the evolution equation for the swirling vector, Eq. (20), where we look at the various terms of dλ/dt as given by the swirling equation, Eq. (17). In what follows, positive values represent the generation of counterclockwise swirling strength along the $\widehat{z}$axis, while negative ones stand for clockwise perturbations.
Figure 3 shows the strength of the various swirling equation terms in a horizontal section in the convection zone of the subdomain under consideration, more precisely, at z = −1.0 Mm below the mean optical surface, τ_{500} = 1. It becomes clear that the two terms that show the strongest vertical values are ${T}_{2}^{\lambda}$ and ${T}_{3}^{\lambda}$, the hydrodynamical and magnetic baroclinic terms, as it was expected from the results of Sect. 4.2. The regions which exhibit strong generation of vertical swirling strength by these two terms harbor a strong magnetic field; in fact, the black contours show the boundaries of magnetic flux concentrations with B_{z} ≥ 500 G. This coexistence can be easily understood regarding the magnetic baroclinicity and the magnetic tension, since their terms are directly related to the magnetic field strength and its variations. Concerning ${T}_{2}^{\lambda}$ instead, we infer from this figure that the magnetic flux concentration is generating strong hydrodynamical baroclinicity. A close examination reveals that, inside the boundaries of the magnetic flux concentrations, the structure of ${T}_{2}^{\lambda}$ is essentially equivalent to the magnetic baroclinic one, ${T}_{3}^{\lambda}$, but with opposite sign. This picture confirms that the two baroclinic terms are partially counterbalancing each other and that this phenomenon is due to mechanical balance of pressures across the magnetic fluxtube boundary.
Fig. 3. Vertical component of the swirling equation terms in a horizontal section in the convection zone at z = −1.0 Mm. The section encompasses the photospheric vortex and surroundings of Fig. 1. The black contours show the boundaries of strong vertical magnetic flux concentrations with B_{z} ≥ 500 G. 
Outside the magnetic flux concentration, the two magnetic terms are weak compared to the hydrodynamical baroclinic one, which contributions are probably due to turbulent motions of the plasma. The stretching term ${T}_{1}^{\lambda}$ is generally weaker than the other terms, in particular inside the magnetic flux concentration, which is again in line with what has been previously found in Sect. 4.2.
The same analysis has been carried out at the mean optical surface (z = 0), which is shown in Fig. 4. In this layer, the dominant terms are still the baroclinic ones, ${T}_{2}^{\lambda}$ and ${T}_{3}^{\lambda}$, which are particularly strong inside the magnetic flux concentration. They show again a similar configuration but with opposite sign with respect to each other. At z = 0, the magnetic tension term starts to be comparable in strength to the baroclinic terms, in particular where the magnetic field is strong. Finally, the stretching term, ${T}_{1}^{\lambda}$, is subdominant also in the photosphere.
Fig. 4. Vertical component of the swirling equation terms in a horizontal section at the optical surface at z = 0 Mm. The section encompasses the photospheric vortex and surroundings of Fig. 1. The black contour shows the boundary of strong vertical magnetic flux concentration with B_{z} ≥ 500 G. 
The picture changes drastically as we move to the chromosphere, at z = 1.5 Mm above the optical surface, as shown in Fig. 5. Here, the swirling strength is principally produced by magnetic processes. In fact, the difference between the hydrodynamical terms, ${T}_{1}^{\lambda}$ and ${T}_{2}^{\lambda}$, and the magnetic ones, ${T}_{3}^{\lambda}$ and ${T}_{4}^{\lambda}$, is evident. At this height, in the solar atmosphere, the stretching and hydrodynamical baroclinic terms are comparable, as we can see from Fig. 2. Furthermore, the strength of the magnetic terms is not related to the magnetic field strength any more. In fact, as we can see in Fig. 5, the contributions of ${T}_{3}^{\lambda}$ and ${T}_{4}^{\lambda}$ to the production of swirling strength are not concentrated where the magnetic field is strong, but are rather homogeneously dispersed over the full cross section.
Fig. 5. Vertical component of the swirling equation terms in a horizontal section in the chromosphere at z = 1.5 Mm. The section encompasses the photospheric vortex and surroundings of Fig. 1. The black contour shows the boundary of strong vertical magnetic flux concentration with B_{z} ≥ 75 G. 
An interesting peculiarity of the swirling equation terms is their patchy structure, as can be seen in Figs. 3–5: they all share a pattern of sharp and chaotic change of sign from one neighboring smallscale region to another, and this structure appears at all characteristic layers. The size of the patches appears to be smaller in the convection zone (l ≲ 50 km) and, in general, in the concentrations of the magnetic field (see Fig. 4). In the chromosphere instead, the patches are on average slightly larger (l ≈ 100 km).
The swirling strength does not exhibit these patches (see, e.g., panel b of Fig. 1), which means that the contributions by the terms of the swirling equation generally do not accumulate but they do average to zero in time. This picture may indicate that a large part of these smallscale perturbations are chaotically produced. This situation is reminiscent of waves, which amplitudes are periodically enhanced and reduced by sinusoidal sources. Therefore, these patches could also be signatures of torsional Alfvén waves and Alfvéntype motions in the solar atmosphere. In fact, Figs. 3–5 show the instantaneous generation of both clockwise and counterclockwise swirling strength at a specific time of the simulation. In the convection zone and low photosphere, the production of a clockwise swirling motion in the plasma, followed by a counterclockwise one some time after, would induce a torsional wave in the “frozenin” magnetic field. Conversely, Alfvén waves in the chromosphere would induce successive positive and negative values in the swirling strength, produced by the magnetic terms of the swirling equation.
Torsional Alfvén waves have been observed by Jess et al. (2009) in a large bright point group with area 430 000 km^{2}. This supports our hypothesis because the perturbations in the swirling strength that we see are mainly confined within the boundaries of a vertically directed magnetic flux concentration. Moreover, Shelyag et al. (2013) found signatures of propagating Alfvén waves in numerical simulations of the photosphere, while Murawski et al. (2015) have shown that Alfvén waves can propagate within magnetic flux tubes from the photosphere to the corona. Battaglia (2020) found from magnetohydrodynamic simulation data unidirectional swirls propagating with Alfvén speed from the photosphere to the chromosphere, being generated by rotational motions at their magnetic footpoints. From the results presented in this paper we cannot confirm the presence of these waves in our simulations. However, the possibility that the smallscale perturbations we observe, or at least part of them, are signatures of these waves is exciting and deserves a more detailed study.
5. Conclusion
In this paper, we discussed the importance of “shearfree” detection criteria for a proper interpretation of vortices. Vorticity was the only criterion for which an evolution equation was known so far and, therefore, the only possible choice in order to study the dynamics of vortices. Since vorticity can lead to erroneous detections of vortices in turbulent flows, its evolution equation is in principle also not reliable for dynamical investigations. For this reason, we derived an evolution equation for the swirling strength, which is a more reliable choice for vortex identification and for which we provide a detailed description and physical interpretation.
The swirling equation represents a new theoretical result in (magneto)hydrodynamics. It can be easily generalized by modifying the underlying momentum equation, for example by adding viscosity terms or by removing the magnetic fields. With this equation, the dynamics of vortices can now be studied in more detail and confidence, since the presence of shear flows does not bias the results.
As a first application of the swirling equation, we investigated which one of its source terms is predominant in numerical simulations of the solar atmosphere. We found that, in the convection zone, the hydrodynamic and magnetic baroclinic terms dominate the production of the vertical component of the swirling strength. Concerning vertical vorticity, the tilting term dominates its production. Higher up, both the production of vertical swirling strength and vorticity is essentially due to magnetic effects. Nevertheless, while magnetic tension is the main factor contributing to the evolution of vorticity, for the swirling strength, both the magnetic baroclinic and the tension term are equally important.
We also showed the vertical component of the swirling equation terms in horizontal sections of a small portion of a simulation model at three characteristic heights. We found that, in the convection zone and in the photosphere, the production of vertical swirling strength is mainly located in strong concentrations of magnetic field. Furthermore, the hydrodynamic and magnetic baroclinic terms often have opposite effects canceling each other. In the chromosphere instead, the magnetic terms alone dominate and produce swirling strength virtually everywhere, in a rather isotropic fashion and quite independently of the location of the magnetic flux concentration in the photosphere. These figures also reveal that all the terms share a common smallscale patchy structure, which implies the production of patches of swirling strength of opposite orientation in almost contiguous regions. These swirls could possibly be caused by turbulent motions in the plasma, however, since their strength is increased within magnetic flux concentrations, they could indicate the presence of torsional Alfvén waves.
The swirling strength criterion and its equation are directly applicable to numerical simulations, since the velocity field and the other quantities are known at each time step and spatial point. From observations, the threedimensional velocity field could be derived by combining Doppler measurements with local correlation or feature tracking. However, a single wavelength measurement would be limited to a certain optical surface. The threedimensional velocity gradient tensor would be derived from observations in multiple spectral lines and wavelength regions, but it would be limited in accuracy and afflicted by noise. Nevertheless, in a first step, one could limit the application of the swirling strength criterion to a twodimensional optical surface. In this case, the detection of swirls in that surface would be based on the twodimensional velocity gradient tensor.
We expect that imaging and spectropolarimetry with the new Daniel K. Inouye Solar Telescope (DKIST) will reveal details of photospheric and chromospheric swirls of unprecedented high spatial and temporal resolution. Therefore, a deeper understanding of the physics related to vortices in the solar plasma is essential. In this sense, we believe that the swirling equation derived in the present paper is an important step forward to the understanding of the underlying mechanisms and processes of vortical motions in the solar atmosphere.
Acknowledgments
This work was supported by the Swiss National Science Foundation under grant ID 200020_182094. Special thanks are addressed to A. Battaglia, A. Bossart, G. Janett, P. Rajaguru, and G. Vigeesh for enriching comments and discussions. The numerical simulations were carried out by F. Calvo on Piz Daint at CSCS under project ID s560. This work has profited from discussions with the team of K. Tziotziou and E. Scullion (conveners) “The Nature and Physics of Vortex Flows in Solar Plasmas” at the International Space Science Institute (ISSI). O.S. wishes to acknowledge scientific discussions within the WaLSA team (Waves in the Lower Solar Atmosphere; https://www.WaLSA.team). The authors are grateful to the anonymous referee for valuable comments that helped improving the article.
References
 Ahrens, J., Geveci, B., & Law, C. 2005, in Visualization Handbook, eds. C. D. Hansen, & C. R. Johnson (Elsevier Inc.), 717 [CrossRef] [Google Scholar]
 Attie, R., Innes, D. E., & Potts, H. E. 2009, A&A, 493, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Balmaceda, L., Vargas Domínguez, S., Palacios, J., Cabello, I., & Domingo, V. 2010, A&A, 513, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Battaglia, A. 2020, Master’s thesis, ETHZürich, submitted [Google Scholar]
 Bonet, J. A., Márquez, I., Sánchez Almeida, J., Cabello, I., & Domingo, V. 2008, ApJ, 687, L131 [NASA ADS] [CrossRef] [Google Scholar]
 Bonet, J. A., Márquez, I., Sánchez Almeida, J., et al. 2010, ApJ, 723, L139 [NASA ADS] [CrossRef] [Google Scholar]
 Brandt, P. N., Scharmer, G. B., Ferguson, S., et al. 1988, Nature, 335, 238 [NASA ADS] [CrossRef] [Google Scholar]
 Calvo, F. 2018, Ph.D. Thesis, University of Geneva, sc. 5281 [Google Scholar]
 Calvo, F., Steiner, O., & Freytag, B. 2016, A&A, 596, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Freytag, B., Steffen, M., Ludwig, H. G., et al. 2012, J. Comput. Phys., 231, 919 [Google Scholar]
 Graftieaux, L., Michard, M., & Grosjean, N. 2001, Meas. Sci. Technol., 12, 1422 [NASA ADS] [CrossRef] [Google Scholar]
 Gudiksen, B. V., Carlsson, M., Hansteen, V. H., et al. 2011, A&A, 531, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Haimes, R., & Kenwright, D. 1999, On the Velocity Gradient Tensor and Fluid Feature Extraction (American Institute of Aeronautics & Astronautics), 315 [Google Scholar]
 Haller, G., Hadjighasem, A., Farazmand, M., & Huhn, F. 2016, J. Fluid Mech., 795, 136 [NASA ADS] [CrossRef] [Google Scholar]
 Harten, A., Lax, P. D., Leer, B., & v, 1983, SIAM Rev., 25, 35 [Google Scholar]
 Jeong, J., & Hussain, F. 1995, J. Fluid Mech., 285, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Jess, D. B., Mathioudakis, M., Erdélyi, R., et al. 2009, Science, 323, 1582 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kato, Y., & Wedemeyer, S. 2017, A&A, 601, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kolář, V. 2007, Int. J. Heat Fluid Flow, 28, 638 [CrossRef] [Google Scholar]
 Kolář, V. 2009, AIAA J., 47, 473 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, C., Gao, Y.S., Dong, X.R., et al. 2019, J. Hydrodyn., 31, 205 [Google Scholar]
 Manso Sainz, R., Martínez González, M. J., & Asensio Ramos, A. 2011, A&A, 531, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Moll, R., Cameron, R. H., & Schüssler, M. 2011, A&A, 533, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Murawski, K., Solov’ev, A., Musielak, Z. E., Srivastava, A. K., & Kraśkiewicz, J. 2015, A&A, 577, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nordlund, Å., & Galsgaard, K. 1995, A 3D MHD code for parallel computers, Tech. rep., Astron. Observ. (Copenhagen University) [Google Scholar]
 Nordlund, Å., Galsgaard, K., & Stein, R. F. 1994, in NATO Advanced Science Institutes (ASI) Series C, eds. R. J. Rutten, & C. J. Schrijver, NATO Adv. Sci. Inst. (ASI) Ser. C, 433, 471 [NASA ADS] [Google Scholar]
 Nordlund, A., Spruit, H. C., Ludwig, H. G., & Trampedach, R. 1997, A&A, 328, 229 [NASA ADS] [Google Scholar]
 Shelyag, S., Keys, P., Mathioudakis, M., & Keenan, F. P. 2011, A&A, 526, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shelyag, S., Fedun, V., Erdélyi, R., Keenan, F. P., & Mathioudakis, M. 2012, in Second ATSTEAST Meeting: Magnetic Fields from the Photosphere to the Corona, eds. T. R. Rimmele, A. Tritschler, F. Wöger, et al., ASP Conf. Ser., 463, 107 [NASA ADS] [Google Scholar]
 Shelyag, S., Cally, P. S., Reid, A., & Mathioudakis, M. 2013, ApJ, 776, L4 [NASA ADS] [CrossRef] [Google Scholar]
 Stein, R. F., & Nordlund, Å. 1998, ApJ, 499, 914 [Google Scholar]
 Tziotziou, K., Tsiropoula, G., Kontogiannis, I., Scullion, E., & Doyle, J. G. 2018, A&A, 618, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tziotziou, K., Tsiropoula, G., & Kontogiannis, I. 2019, A&A, 623, A160 [CrossRef] [EDP Sciences] [Google Scholar]
 Vargas Domínguez, S., Palacios, J., Balmaceda, L., Cabello, I., & Domingo, V. 2011, MNRAS, 416, 148 [NASA ADS] [Google Scholar]
 Vögler, A. 2003, PhD Thesis, Göttingen University [Google Scholar]
 Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wedemeyer, S., & Steiner, O. 2014, PASJ, 66, S10 [NASA ADS] [CrossRef] [Google Scholar]
 WedemeyerBöhm, S., & van der Rouppe Voort, L. 2009, A&A, 507, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 WedemeyerBöhm, S., Scullion, E., Steiner, O., et al. 2012, Nature, 486, 505 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Wedemeyer, S., Scullion, E., Steiner, O., de la Cruz Rodriguez, J., & Rouppe van der Voort, L. H. M. 2013, J. Phys. Conf. Ser., 440, 012005 [NASA ADS] [CrossRef] [Google Scholar]
 Zhou, J., Adrian, R. J., Balachandar, S., & Kendall, T. M. 1999, J. Fluid Mech., 387, 353 [NASA ADS] [CrossRef] [Google Scholar]
 Zwaan, C. 1978, Sol. Phys., 60, 213 [Google Scholar]
Appendix A: Physical interpretation of the swirling strength
The swirling strength is a precise criterion for vortex identification, however its convoluted definition causes its physical interpretation to be less straightforward than that of vorticity. This appendix is intended to shed some light on its mathematical definition and physical meaning.
We start from Eq. (5), which shows the explicit diagonalization of the velocity gradient tensor 𝒰. As becomes clear further down, whenever the flow is rotating, 𝒰 has one real eigenvalue, λ_{r}, and two complex conjugate ones, λ_{+} = (λ_{−})^{*}. Since λ_{+} and λ_{−} are complex conjugates, so are their relative eigenvectors, u_{+} = (u_{−})^{*}. To prove it, we take the defining relation for the u_{+} eigenvector, 𝒰u_{+} = λ_{+}u_{+}, and we complex conjugate it. We then obtain,
$$\begin{array}{cc}& {\left(\mathcal{U}{\mathit{u}}_{+}\right)}^{\ast}={\left({\lambda}_{+}{\mathit{u}}_{+}\right)}^{\ast},\hfill \\ \hfill & \mathcal{U}{\left({\mathit{u}}_{+}\right)}^{\ast}={\lambda}_{}{\left({\mathit{u}}_{+}\right)}^{\ast},\hfill \end{array}$$
where we use the fact that 𝒰 is a real matrix and λ_{+} = (λ_{−})^{*}. The last equation proves our proposition, since that relation is equivalent to the defining relation for the u_{−} eigenvector, 𝒰u_{−} = λ_{−}u_{−}.
Yet, in this formulation, both Λ and 𝒫 appearing in Eq. (5) are complex matrices, therefore, a physical interpretation is difficult to grasp. Consequently, we would like to transform them into real matrices; an easy way to do it is to use the following transformations,
$$\begin{array}{cc}& {\mathit{u}}_{\mathrm{cr}}=\frac{1}{\sqrt{2}}({\mathit{u}}_{+}+{\mathit{u}}_{}),\phantom{\rule{1em}{0ex}}{\mathit{u}}_{\mathrm{ci}}=\frac{1}{\sqrt{2}\mathrm{i}}({\mathit{u}}_{+}{\mathit{u}}_{}),\hfill \\ \hfill & {\lambda}_{\mathrm{cr}}=\frac{1}{2}({\lambda}_{+}+{\lambda}_{}),\phantom{\rule{1em}{0ex}}{\lambda}_{\mathrm{ci}}=\frac{1}{2\mathrm{i}}({\lambda}_{+}{\lambda}_{}).\hfill \end{array}$$
Then, we can prove that the velocity gradient tensor 𝒰 can be decomposed in the following form,
$$\begin{array}{c}\hfill \mathcal{U}=\underset{\stackrel{\sim}{\mathcal{P}}}{\underset{\u23df}{[{\mathit{u}}_{\mathrm{r}},{\mathit{u}}_{\mathrm{cr}},{\mathit{u}}_{\mathrm{ci}}]}}\underset{\stackrel{\sim}{\mathrm{\Lambda}}}{\underset{\u23df}{\left[\begin{array}{ccc}{\lambda}_{\mathrm{r}}& 0& 0\\ 0& {\lambda}_{\mathrm{cr}}& {\lambda}_{\mathrm{ci}}\\ 0& {\lambda}_{\mathrm{ci}}& {\lambda}_{\mathrm{cr}}\end{array}\right]}}\underset{{\stackrel{\sim}{\mathcal{P}}}^{1}}{\underset{\u23df}{{[{\mathit{u}}_{\mathrm{r}},{\mathit{u}}_{\mathrm{cr}},{\mathit{u}}_{\mathrm{ci}}]}^{1}}},\end{array}$$(A.1)
where $\stackrel{\sim}{\mathcal{P}}$, $\stackrel{\sim}{\mathrm{\Lambda}}$, and ${\stackrel{\sim}{\mathcal{P}}}^{1}$ are real matrices.
We analyze in more detail the new blockdiagonal matrix $\stackrel{\sim}{\mathrm{\Lambda}}$, where we can distinguish two blocks: the first one of dimensions 1 × 1 is associated with the eigenvector u_{r}. It can be seen as a dilatation operator, since when applied to a scalar it simply acts as a multiplicative constant. The eigenvalue λ_{r} can therefore be interpreted as a measure of how the flow is stretched or compressed along the direction given by u_{r}. The 2 × 2 block matrix instead can be further decomposed in two pieces,
$$\begin{array}{c}\hfill \left[\begin{array}{cc}{\lambda}_{\mathrm{cr}}& {\lambda}_{\mathrm{ci}}\\ {\lambda}_{\mathrm{ci}}& {\lambda}_{\mathrm{cr}}\end{array}\right]=\underset{\stackrel{\sim}{\mathcal{D}}}{\underset{\u23df}{\left[\begin{array}{cc}{\lambda}_{\mathrm{cr}}& 0\\ 0& {\lambda}_{\mathrm{cr}}\end{array}\right]}}+\underset{\stackrel{\sim}{\mathcal{R}}}{\underset{\u23df}{\left[\begin{array}{cc}0& {\lambda}_{\mathrm{ci}}\\ {\lambda}_{\mathrm{ci}}& 0\end{array}\right]}}.\end{array}$$
The $\stackrel{\sim}{\mathcal{D}}$ matrix can be interpreted as a 2D dilatation operator, for the same reasons we defined the 1 × 1 matrix in this way. The $\stackrel{\sim}{\mathcal{R}}$ matrix instead represents a infinitesimal rotation operator in the plane spanned by u_{cr} and u_{ci}, which orientation is clockwise in the direction defined by u_{r}. Therefore, we can say that λ_{cr} measures how stretched is the flow in the rotation plane, while λ_{ci} describes the strength of the swirling flow and u_{r} the vortex axis: hence the definition of swirling strength. For some visual intuition, the interested reader can refer to Haimes & Kenwright (1999).
The mathematical description here presented hides some arbitrariness in what concerns the orientation (clockwise or counterclockwise) of the swirl. In fact, eigenvectors are, in general, defined up to a multiplicative constant, which results in a arbitrary definition of their norm and orientation. This implies that, mathematically, we can multiply the eigenvector u_{r} by −1 and it would still be a correct solution of the eigenanalysis. However, physically, this would imply the inversion of the vortex orientation.
To be able to select the correct orientation of the eigenvector u_{r}, one has to check the handedness of the new basis. In fact, the matrix $\stackrel{\sim}{\mathcal{R}}$ describes a clockwise rotation only if the 3D basis in which it is defined is righthanded. The basis in question is not the standard one for a threedimensional space, but the one composed by the real vectors u_{r}, u_{ci}, and u_{cr} (in this order). Indeed, when we decompose the velocity gradient tensor 𝒰 as in Eq. (A.1), we are essentially performing a change of basis. The new basis, embodied in the matrix $\stackrel{\sim}{\mathcal{P}}$, is proper to the flow. Therefore, to find out if the orientation of the eigenvector is correct, one has to check the orientation of the basis itself. One easy way to do it, is to compute the determinant of the matrix formed by the three basis vectors, that is, $\stackrel{\sim}{\mathcal{P}}$: if $\mathrm{Det}(\stackrel{\sim}{\mathcal{P}})>0$ the basis is righthanded, while if $\mathrm{Det}(\stackrel{\sim}{\mathcal{P}})<0$ it is the opposite. In this way, one can fix the orientation of the real eigenvector and ensure that the swirling vector, λ ≡ λu_{r}, correctly describes the vortex flow.
Appendix B: Commutator with diagonal matrices
We prove that the commutator between a diagonal and a generic matrix results in a hollow matrix. In order to do that, here we consider the commutator between two matrices 𝒞 and 𝒟. If 𝒟 is diagonal, its terms can be expressed as 𝒟_{ij} ≡ d^{(i)}δ_{ij}, where δ_{ij} is the Kronecker’s delta and d^{(i)} are the diagonal terms,
$$\begin{array}{cc}\hfill {[\mathcal{C},\mathcal{D}]}_{\mathit{ij}}& ={\displaystyle \sum _{k}({\mathcal{C}}_{\mathit{ik}}{\mathcal{D}}_{\mathit{kj}}{\mathcal{D}}_{\mathit{ik}}{\mathcal{C}}_{\mathit{kj}}),}\hfill \\ \hfill & ={\displaystyle \sum _{k}({\mathcal{C}}_{\mathit{ik}}{d}^{(k)}{\delta}_{\mathit{kj}}{d}^{(i)}{\delta}_{\mathit{ik}}{\mathcal{C}}_{\mathit{kj}}),}\hfill \\ \hfill & =({d}^{(j)}{d}^{(i)}){\mathcal{C}}_{\mathit{ij}},\hfill \end{array}$$(B.1)
which implies that [𝒞,𝒟]_{ij} = 0 for i = j.
Appendix C: Alternative derivation of the swirling equation
In Sect. 3 we selected the (2, 2) component of Eq. (15) in order to derive the swirling equation. However, we also could have chosen the (3, 3) component. Here, we prove that the two approaches result in the same equation.
The (3, 3) component of Eq. (15) yields the dynamical equation of the λ_{−} eigenvalue,
$$\begin{array}{cc}& \frac{\mathrm{d}}{\mathrm{d}t}{\lambda}_{}={\lambda}_{}^{2}+{\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{33}\phantom{\rule{0.166667em}{0ex}},\hfill \\ \hfill & \frac{\mathrm{d}}{\mathrm{d}t}({\lambda}_{\mathrm{cr}}\mathrm{i}{\lambda}_{\mathrm{ci}})={({\lambda}_{\mathrm{cr}}\mathrm{i}{\lambda}_{\mathrm{ci}})}^{2}+{\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{33}\phantom{\rule{0.166667em}{0ex}}.\hfill \end{array}$$
Therefore, by taking the imaginary part, we get,
$$\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}t}\lambda =& 2\lambda {\lambda}_{\mathrm{cr}}2\mathrm{Im}{\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{33}\phantom{\rule{0.166667em}{0ex}},\hfill \end{array}$$
which is equivalent to Eq. (17) but for the last term. Therefore, to prove that the two approaches are equivalent we need to check that,
$$\begin{array}{c}\hfill \mathrm{Im}{\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{33}=\mathrm{Im}{\left({\mathcal{P}}^{1}\mathcal{M}\mathcal{P}\right)}_{22}\phantom{\rule{0.166667em}{0ex}}.\end{array}$$(C.1)
This can be done by considering the properties of the eigenvectors which compose the matrix 𝒫. In Appendix A we have shown that u_{+} = (u_{−})^{*}, then one can prove that a similar relation holds also for the vectors composing 𝒫^{−1}. In fact, if we describe 𝒫^{−1} with horizontal vectors,
$$\begin{array}{c}\hfill {[{\mathit{u}}_{\mathrm{r}},{\mathit{u}}_{+},{\mathit{u}}_{}]}^{1}\equiv \left[\begin{array}{c}{\mathit{t}}_{\mathrm{r}}^{\phantom{\rule{0.166667em}{0ex}}\mathrm{T}}\\ {\mathit{t}}_{+}^{\phantom{\rule{0.166667em}{0ex}}\mathrm{T}}\\ {\mathit{t}}_{}^{\phantom{\rule{0.166667em}{0ex}}\mathrm{T}}\end{array}\right]\phantom{\rule{0.166667em}{0ex}},\end{array}$$
the reader can check that t_{+} = (t_{−})^{*}.
Finally, using this decomposition, we can directly compute the two terms which we want to prove to be equivalent,
$$\begin{array}{c}\hfill {\mathcal{P}}^{1}\mathcal{M}\mathcal{P}=\left[\begin{array}{c}{\mathit{t}}_{\mathrm{r}}^{\phantom{\rule{0.166667em}{0ex}}\mathrm{T}}\\ {\mathit{t}}_{+}^{\phantom{\rule{0.166667em}{0ex}}\mathrm{T}}\\ {\mathit{t}}_{}^{\phantom{\rule{0.166667em}{0ex}}\mathrm{T}}\end{array}\right]\left[\begin{array}{ccc}{m}_{11}& {m}_{12}& {m}_{13}\\ {m}_{21}& {m}_{22}& {m}_{23}\\ {m}_{31}& {m}_{32}& {m}_{33}\\ \end{array}\right][{\mathit{u}}_{\mathrm{r}},{\mathit{u}}_{+},{\mathit{u}}_{}].\end{array}$$
By carrying out the computation for both (2, 2) and (3, 3) components, one can check that Eq. (C.1) holds true and that, therefore, the two approaches to the derivation of the swirling equation are equivalent.
All Figures
Fig. 1. Vertical component of (a) vorticity, (b) swirling strength, and (c) shear strength, computed in a small portion of the full simulation domain, measuring 0.6 × 0.6 × 1.8 Mm^{3} and encompassing an intergranular region with strong vertical magnetic field (B_{z} ≳ 500 G). The gray sheet displays the surface of optical depth τ_{500} = 1, while the plasma flow is rendered with green streamlines. These images were produced with the software ParaView (Ahrens et al. 2005). 

In the text 
Fig. 2. Mean absolute values of the vertical components of the terms appearing on the righthand side of the vorticity equation (left panel) and the swirling equation (right panel) as a function of z. The dotted curves represent the mean absolute value of the total production of vorticity and swirling strength in the simulation. The optical surface τ_{500} = 1 corresponds to z = 0. 

In the text 
Fig. 3. Vertical component of the swirling equation terms in a horizontal section in the convection zone at z = −1.0 Mm. The section encompasses the photospheric vortex and surroundings of Fig. 1. The black contours show the boundaries of strong vertical magnetic flux concentrations with B_{z} ≥ 500 G. 

In the text 
Fig. 4. Vertical component of the swirling equation terms in a horizontal section at the optical surface at z = 0 Mm. The section encompasses the photospheric vortex and surroundings of Fig. 1. The black contour shows the boundary of strong vertical magnetic flux concentration with B_{z} ≥ 500 G. 

In the text 
Fig. 5. Vertical component of the swirling equation terms in a horizontal section in the chromosphere at z = 1.5 Mm. The section encompasses the photospheric vortex and surroundings of Fig. 1. The black contour shows the boundary of strong vertical magnetic flux concentration with B_{z} ≥ 75 G. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.