© J. R. Canivete Cuissa and O. Steiner 2022

1 Introduction

Vortices are one of the fundamental features of fluid dynamics and turbulent flows. Although a vortex can conceptually be described as a fluid region rotating around a common axis, the task of compiling a mathematical definition of a vortex has proven to be incredibly challenging and it remains a subject of debate (see, e.g., Günther & Theisel 2018). A rigorous and objective definition is particularly needed for the identification of vortices in highly dynamical and turbulent flows, where detections based on the naive human intuition of swirling structures could bias the results.

Multiple detection methods and criteria have been proposed in the literature. Most of them can be categorized into two classes. The first and most widely known class consists of mathematical criteria based on local physical quantities related to the flow, such as the velocity field, the pressure, and their derivatives. Local criteria are defined for each point in space. In practice, on a discrete grid, the vortex identification involves a limited stencil of neighboring grid cells. Examples are the vorticity ω, the Q-criterion (Hunt et al. 1988), the λ₂-criterion (Jeong & Hussain 1995), the swirling strength λ (Zhou et al. 1999), the Γ functions, (Graftieaux et al. 2001), the instantaneous vorticity deviation (IVD), and the Lagrangian-averaged vorticity deviation (LAVD; Haller et al. 2016). A vortex is usually identified as a connected over-density region of one of these criteria.

Most of these methods are mathematically rigorous and physically consistent. For example, all the listed criteria are Galilean invariant and the IVD/LAVD is considered an objective measure¹ (Haller et al. 2016). However, a portion of these methods may be prone to miss vortices that are characterized by weak rotational velocities, since the mathematical criteria are related to the angular velocity of the flow and a threshold is usually employed to filter out low-magnitude, noisy signals. This problem is of particular concern with regard to the vorticity, the g-criterion, the λ₂-criterion, the swirling strength – for which weak rotational signals and the background turbulent noise are essentially indistinguishable, while, by definition, the Γ functions method requires a threshold on the value of |Γ₁|. The IVD/LAVD methods may be more robust in this regard as they are based on fluctuations of the vorticity with respect to the domain spatial mean, which should include, to a certain extent, the background noise. Moreover, false detections can happen when the flow is curved but does not perform a full rotation, since these methods only measure the local curvature at each grid point and do not take into account the global behavior of the flow.

The second class consists of morphological methods, which take properties of the flow in the entire domain (or a subset of it) into account. Morphological criteria are not necessarily defined in each point of space as vortices may be sparsely distributed. These methods rely on intuitive definitions of what a vortex is, such as the assertion that it consists of a “multitude of material particles rotating around a common center” (Lugt 1979), or resulting when “instantaneous streamlines mapped onto a plane normal to the vortex core exhibit a roughly circular or spiral pattern, when viewed from a reference frame moving with the center of the vortex core” (Robinson 1990).

Sadarjoen (1999) presented two techniques for the identification of vortices based on characterization of streamlines. The first, namely, the “curvature center method,” determines the center of curvature of streamlines and defines a vortex core as a region where curvature centers accumulate. The second, that is, the “winding-angle method,” detects vortices by clustering streamlines that appear to be curved around a close set of points. These methods are conceptually simple and consider the global features of the flow, but they usually are computationally expensive and are not easily extendable to three dimensions. Moreover, streamlines can give false impressions of vortical motions in unsteady and highly dynamical flows (see, e.g., Shelyag et al. 2013).

Many of the methods listed above have been employed to identify vortices in observations and numerical simulations of the solar atmosphere. Small-scale vortical motions appear to be ubiquitous in the quiet solar photosphere and chromosphere. Moreover, they are thought to play an important role in the transport of energy towards the upper layers of the solar atmosphere (see, e.g., Tziotziou et al. 2022, for a review).

Numerical simulations of the solar atmosphere are particularly suited for the application of mathematical criteria for the identification of vortices, since the velocity field and other quantities of interest are directly accessible at each point in the three-dimensional (3D) space. The vorticity was used by Stein & Nordlund (1998), Shelyag et al. (2011, 2013), and Steiner & Rezaei (2012), to detect vortices and study their dynamics, while the enstrophy, defined as |ω|², was used by Kitiashvili et al. (2012) to visualize three-dimensional vortical structures. Moll et al. (2011) first introduced the swirling strength in the context of photospheric swirls, a quantity that has also been adopted by Moll et al. (2012), Kato & Wedemeyer (2017), Canivete Cuissa & Steiner (2020), Yadav et al. (2020), and Battaglia et al. (2021). More recently, the IVD and LAVD criteria have also been employed in numerical simulations of the solar photosphere and chromosphere by Silva et al. (2020, Silva et al. 2021), and Aljohani et al. (2022).

Vortex identification methods can not be applied to observations in a straightforward manner since the necessary quantities are not immediately available from the data. However, it is possible to estimate the horizontal velocity field through the use of local correlation tracking (LCT) techniques. Therefore, identification methods that rely on the velocity field alone have also been employed on observational data. From horizontal velocity maps derived with LCT methods, Requerey et al. (2017, 2018) employed Lagrangian tracers to visually identify vortices in super and meso-granular vertices, while Giagkiozis et al. (2018), and Liu et al. (2019a,b) utilized the Γ functions. Silva et al. (2018) presented a comparison between the swirling strength, the Γ functions, and the LAVD criteria on horizontal velocity fields extracted from observations. An alternative approach has been put forward by Dakanalis et al. (2021), whereby the authors identify vortices directly from chromospheric filtergrams using the morphological winding-angle method.

A rigorous identification process is of fundamental importance to infer the statistical properties of vortical motions in the solar atmosphere and, consequently, their impact on chromospheric and coronal heating. The aforementioned detection methods greatly helped to shape our current understanding, but they are prone to errors and misidentifications. As discussed above, weakly rotating vortices can be problematic for mathematical criteria such as the vorticity, the swirling strength, and the Γ functions, especially in turbulent and dynamical flows such as those present in the solar atmosphere. Silva et al. (2018) showed that an extra criterion, the d-criterion, should be used with LAVD/IVD methods and with the algorithm proposed by Kato & Wedemeyer (2017) in the presence of strong shear flows, which are typical of integranular regions. Moreover, the identification of the vortex center and boundaries that is needed for a proper study of vortex dynamics and interaction can be obtained with the Γ functions and the LAVD/IVD methods only. In principle, the morphological method presented by Dakanalis et al. (2021) should be the preferred method because it takes the large-scale properties of the flow into account, but it is not well-suited for the analysis of numerical simulations.

Therefore, in this paper, we introduce a new vortex identification method based on a completely new technique. It takes into account the global features of the flow, as in the morphological methods, but also possesses the rigor of the mathematical criteria. Basically, it is a hybrid of the two classes presented above and it possesses the ability of recognizing coherent vortical structures from the velocity field. Moreover, it allows for a precise characterization of the vortex center and boundaries and it is robust against shear flows and noise.

The paper is organized as follows. In Sect. 2, we present the mathematical criteria adopted in this work. The method and the automated algorithm are described in detail in Sect. 3, while in Sect. 4 we test these on vortical and turbulent velocity fields and discuss the results. Finally, we present our summary and conclusions in Sect. 5. The application of the method to numerical simulations and observations of the solar atmosphere is deferred to a follow-up paper.

2 Theoretical background

In this section, we review some aspects of two of the most widely used mathematical criteria for vortical flows: the vorticity ω and the swirling strength λ. Moreover, we introduce a recently developed quantity known as Rortex (alternatively, Liutex) that we have adopted in the method we develop and present in this work.

2.1 Vorticity

Vorticity is the classical quantity to describe local rotational motions in fluid mechanics and it is defined as the curl of the velocity field, v,

$$\omega =\nabla \times v.$$(1)

The direction of the vorticity vector indicates the orientation of the rotation according to the right-hand rule, while its norm is proportional to the intensity ofthe rotation.

For a rotational vortex, that is a flow rotating in a rigid-body fashion around an axis, the norm of the vorticity vector is ω = |ω| = 2Ω, where Ω is the angular velocity of the fluid. However, this simple relation between the vorticity norm and the attributes of vortical flows does not hold for more realistic vortex models in fluid dynamics, such as Lamb–Oseen or Burgers vortices. Therefore, we cannot extract physical information about vortices directly from the vorticity. Moreover, the vorticity can assume large values also in the presence of shear flows and can therefore lead to false detections in non-rotating velocity fields.

2.2 Swirling strength

A suitable solution to the shear flow problem is given by the swirling strength criterion, λ, proposed by Zhou et al. (1999). It is computed through the eigen analysis of the velocity gradient tensor, U = ∇v, which is the Jacobian matrix of the velocity field. If the flow is locally rotating, the velocity gradient tensor can be diagonalized as:

$$U=\underset{P}{\underbrace{\begin{array}{c}\\ \left[u_{\text{r}},u_{+},u_{-}\right]\\ \end{array}}}\underset{\text{Λ}}{\underbrace{\left[\begin{array}{ccc}{\lambda }_{\text{r}}& 0& 0\\ 0& {\lambda }_{+}& 0\\ 0& 0& {\lambda }_{-}\end{array}\right]}}\underset{P^{-1}}{\underbrace{\begin{array}{c}\\ {\left[u_{\text{r}},u_{+},u_{-}\right]}^{-1},\\ \end{array}}}$$(2)

where λ_± = λ_cr ± iλ_ci ∈ ℂ and λ_r ∈ ℝ are the eigenvalues of U that form the eigenvalue matrix Λ, while u_± and u_r are the corresponding eigenvectors that form the change of basis matrix, P. We define the swirling strength λ as twice the imaginary part of the complex eigenvalues:

$$\lambda =2{\lambda }_{\text{ci}}.$$(3)

Moreover, the normalized eigenvector² associated with the real eigenvalue, u_r, defines the rotation axis³. We can then define the swirling strength vector as λ = λu_r, which carries the information about the strength and the direction of rotation. For a simple rotational vortex, the swirling strength vector is identical to the vorticity vector. For more information about the swirling strength vector, we refer to Canivete Cuissa & Steiner (2020).

The swirling strength criterion can be enhanced by considering the inverse spiraling compactness, ζ = λ_cr/λ_ci, where λ_cr is the real part of the complex eigenvalues of U (Chakraborty et al. 2005). The inverse spiraling compactness is a measure of the radial distance traveled by a test particle in the flow during one orbital period. It can be shown that the radial distance, r_n, of a particle spiraling in the vortex plane evolves as:

$$r_{\text{n}}=r_0\exp \left(2\pi n\zeta \right),$$(4)

where r₀ is the initial radius and n is the number of orbits (Chakraborty et al. 2005). If ζ = 0, then the flow is perfectly circular, while if ζ is positive (negative) the flow is spiraling outward (inward). Large (either positive or negative) values of ζ indicate very low spiraling compactness. Therefore, vortex regions identified by the swirling strength should be discarded if |ζ| is too large, since in these cases the orbital component of the velocity field is negligible compared to the radial one and the rotational motion is barely perceivable. In practice, we set λ = 0 wherever κ_ζ> ζ > δ_ζ, where κ_ζ and δ_ζ are free parameters that can be adjusted according to the compactness required. The typical values for these parameters are −κ_ζ ~ δ_ζ < 1.

The swirling strength solves the problem related to pure shear flows, but it still fails to predict the angular velocity of more complex and realistic vortices. Indeed, the swirling strength criterion cannot discern between the rigid-body rotational component of a flow and intrinsic shears, which are present in differentially rotating flows for example (see Sect. 2.3). This means that basic information about vortices, such as their period of rotation and their size, cannot be directly inferred from the simulation data using the swirling strength criterion, as we show in Sect. 2.4.

2.3 Rortex

The Rortex criterion, also called Liutex, is a state of the art mathematical quantity, proposed by Tian et al. (2018) and Liu et al. (2018), which allows for a precise quantification of the local rigid-body rotational part of the flow alone. Indeed, in complex but realistic vortex models, an intrinsic shear component is usually present besides the rigid-body rotational one. An example of such a model is the Lamb–Oseen vortex (see Sect. 2.4). The intrinsic shear component blends with the purely rotational one in the elements and eigenvalues of the velocity gradient tensor. Since the vorticity and the swirling strength cannot discern between these two components of the flow, both may become contaminated by the presence of an intrinsic shear.

To precisely represent the fluid curvature at a point, Liu et al. (2016) proposed to decompose the vorticity vector ω into a purely rotational part R and a non-rotational or shear component S = ω − R. The direction of R represents the rotation axis, while its norm is proportional to the angular velocity of the flow. For a purely rigid-body rotational flow, R = λ = ω at any point because there are no shears in such a flow. However, a series of tedious matrix transformations are in principle necessary to compute the vectors R and S and therefore extract the rotational component of the flow.

As shown by Xu et al. (2019), for any velocity gradient tensor U = ∇v with one real and two complex conjugated eigenvalues, there exists a transposed quasi-triangular matrix V such that:

$$U=QV{Q}^{\text{T}},$$(5)

where Q is an orthogonal matrix representing a rotation operator in three-dimensional space and V can be written as:

$$V=\left[\begin{array}{ccc}{\lambda }_{\text{cr}}& -\varphi & 0\\ \varphi +\epsilon & {\lambda }_{\text{cr}}& 0\\ \xi & v& {\lambda }_{\text{r}}\end{array}\right].$$(6)

The matrix V represents the velocity gradient tensor U in a rotated reference frame with the rotation axis parallel to the new $\widehat{z}$-axis. We recall that the rotation axis is given by the normalized eigenvector associated with the real eigenvalue λr of U, as described in Sect. 2.2.

The various components of V represent rigid-body rotation (ϕ), stretching and compressing components (λ_cr, λ_r), and shearing parts (ε, ξ, v). In particular, the ε coefficient is mixed with the purely rotational component ϕ and it is responsible for the intrinsic shear. From Eq. (5), it becomes clear why the swirling strength can be biased by the presence of intrinsic shear flows, as the eigenvalues of the matrix V are contaminated by the ε coefficient. Hence, the swirling strength directly depends on the strength of intrinsic shears and it cannot be used to measure the strength of the rigid-body rotational part of the flow.

To this end, the Rortex criterion is defined as:

$$R=2\varphi .$$(7)

By definition, it computes the strength of the rotational part of the flow alone. When there are no intrinsic shears (i.e., ϵ = 0), it is possible to show that the Rortex criterion is equivalent to the swirling strength by diagonalizing Eq. (6). Furthermore, we can compute the rotational part, R, of the vorticity vector, ω, by multiplying the Rortex criterion with the rotation axis given by the normalized eigenvector u_r,

$$R=R{u}_{\text{r}}.$$(8)

This vector is referred to as a Vortex or Rortex vector and it allows for a 3D characterization of the vortical flow.

To compute the Rortex criterion, we would have to find the orthogonal matrix, Q, such that the rotated local velocity gradient tensor, V, takes the form of Eq. (6), (see Gao & Liu 2018, for a guided procedure). Fortunately, Wang et al. (2019) and Xu et al. (2019) derived a simple formula to compute the Rortex criterion without having to dive into complex and computationally expensive matrix calculations. Accordingly,

$$R=2\varphi =\omega ·u_{\text{r}}-\sqrt{{\left(\omega ·u_{\text{r}}\right)}^2-{\lambda }^2},$$(9)

where ω is the vorticity vector (Sect. 2.1) and λ = 2λ_ci is the swirling strength criterion (Sect. 2.2). Here, the sign of the normalized eigenvector u_r is chosen such that ω · u_r > 0. Also, the Rortex criterion can be enhanced in the same way Chakraborty et al. (2005) did it for the swirling strength. Hence we set R = 0 wherever κ_ζ> ζ > δ_ζ, with κ_ζ and δ_ζ as the chosen thresholds (see Sect. 2.2).

The Rortex criterion allows for a precise quantification of the rotational strength of a vortex flow. Indeed, it is the only quantity not being contaminated by the intrinsic shear component, ε, in the velocity gradient tensor. Therefore, it is the only reliable quantity for the extraction of physical information about a vortex, such as the rotational period and the curvature radius.

2.4 Comparing entena

To prove the reliability of the Rortex criterion, we compare the three mathematical criteria presented in Sects. 2.1, 2.2, and 2.3 by applying them to a Lamb–Oseen vortex model. Lamb–Oseen vortices are analytical but realistic vortex models and are defined through the following velocity field in polar coordinates:

$$v_{\theta }=v_{\max }\left(1+\frac{1}{2\alpha }\right)\frac{r_{\max }}{r}\left[1-\exp \left(-\alpha \frac{r^2}{r_{\max }^2}\right)\right],v_r=0,$$(10)

where α = 1.256 establishes that r_max is the radius at which the rotational velocity v_θ is maximal and reaches the value v_max. We then define r_max to be the boundary of the vortex. In the following, we set v_max = 1 and r_max = 0.3.

The velocity field corresponding to the vortex model presented above and centered in (x, y) = (0.5, 0.5) is shown in the three panels of Fig. 1, with arrows proportional to the local strength of the velocity field, while the blue dashed lines represent the boundary of the vortex given by r = r_max = 0.3. The panels show (from left to right), the maps of the vorticity, ω, the swirling strength, λ, and the Rortex, R. As we can see, the three criteria identify the presence of local curvature in the flow. The swirling strength and the Rortex are non-zero only inside the boundary of the vortex, while vorticity shows positive values also for r > r_max. Moreover, the three criteria show different values in the vortical region.

To know which criteria is the most appropriate to measure the characteristics of such a vortical flow, we compare the rotational period, T, estimated with the vorticity, the swirling strength, and the Rortex at each point in the grid. Given the symmetry of the problem, the rotational period T depends only on the radius r. We assume that the angular velocity Ω and these quantities are related by c = 2 Ω _c, where c = ω, λ, R. Then, we can estimate the rotational period as:

$$T_{\text{c}}=\frac{2\pi }{\left|{\text{Ωc}}_{}\right|}=\frac{4\pi }{\left|c\right|},$$(11)

and compare it to the true period derived from the analytical formula Eq. (10), that is:

$$T_{\text{LO}}=\frac{2\pi }{\left|{\text{Ω}}_{\text{LO}}\right|}=\frac{2\pi r}{\left|v_{\theta }\right|},$$(12)

since Ω_LO = v_θ/r. We note that |c| in Eq. (11) can be a function of radius as is Ω_LO.

Figure 2 shows the radial profiles of such estimated rotational periods against the true value T_LO. The only criterion able to correctly estimate the true rotational period of the vortex flow as a function of the radius is the Rortex. This is because intrinsic shear contaminate the purely rotational component of Lamb–Oseen vortices. Therefore, the Rortex criterion is the optimal quantity to extract precise information about the rotational characteristics of this kind of vortices. This conclusion can be drawn also analytically and remains valid for many other vortex models, as we show in Appendix A.

In principle, we can define a vortical region as a collection of connected grid cells where R ≠ 0. However, this definition only takes into account the local properties of the flow and lacks of a more global perspective. Since a vortex is a coherent structure that has a spatial extension, it is important to consider the relation between all the fluid parcels that compose it. This is why, in the next section, we combine the Rortex criterion with a morphological method.

3 Method

In this section, we present a new vortex identification method. As stated above, this new method is based on both morphological and mathematical criteria. More precisely, it is an extension of the curvature center method proposed by Sadarjoen (1999), where we adopt the Rortex criterion to compute the coordinates of the curvature center. In this way, we solve the main flaw of the curvature center method, that is it being based on the curvature of streamlines of the flow. Our method takes the intuitive idea of the curvature center method and combines it with a state of the art mathematical criterion (the Rortex criterion) for a precise estimation of the curvature center position for all points where the velocity field is curved.

Fig. 1 Comparison between (a) the vorticity, (b) the swirling strength, and (c) the Rortex criteria for a Lamb–Oseen vortex. Arrows indicate the velocity vector field of the Lamb–Oseen vortex and their length is proportional to the magnitude of the flow. The blue-dashed circles denote the boundary of the vortex defined by r = rmax.

Fig. 2 Comparison of the estimated rotational period T derived from the vorticity, the swirling strength, and the Rortex criteria with the true rotational period for a Lamb–Oseen vortex. Blue vertical lines denote the boundary of the vortex defined by r = rmax. The thick gray curve indicates the true rotational period TLO computed from the analytical formula Eq. (10).

3.1 Estimated vortex center map

The basic idea of the method presented in this paper is to compute a map with all the “estimated vortex centers” (EVC), that is, the curvature center points from any parcel of flow presenting a rigid-body rotational component. Hereon, we refer to this map as the EVC map. Whenever a vortex with a rigid-body rotational component⁴ is present in the velocity field, a high density (or a cluster) of EVC points will show up in this map, since all points belonging to the vortex share a common curvature center. If instead the flow is characterized by incoherent random local curvatures only, the EVC points will be scattered in the domain and they will not form clusters. The reliability of such a map, and therefore of the method, depends on how precisely the vortex center can be estimated. In the following, we show how to mathematically estimate the center of curvature, that is, how to calculate the EVC map.

For simplicity, we consider only a two-dimensional flow, but this methodology can be extended to three dimensions. Given a generic velocity field v = (v_x, v_y) defined in a Cartesian grid with coordinates, x = (x, y), the vorticity, ω, the velocity gradient tensor, U, the real eigenvector, u_r, the swirling strength, λ, and the Rortex, R, can be computed for all the points in the grid according to Eqs.(1), (3), and (9). In the following we focus on a generic point with coordinates x₀ = (x₀, y₀), Rortex value R₀ ≠ 0, and velocity v₀ = (v_x,0, v_y,0). Moreover, we assume the velocity gradient tensor in x₀, $U{|}_{x_0}$, to be known. Then, the EVC of the point x₀ can be computed based on two quantities: the curvature radius r, namely, the distance between the point x0 and the center of curvature, and the radial direction, e_r, that is the unit vector pointing from the point x₀ to the center of curvature.

3.1.1 The radial direction

To estimate the radial direction in a generic point x, the naive idea would be to take the vector perpendicular to the velocity field in x. However, if the flow in this point has a radial velocity component or shear, the estimated vector will deviate from the true radial direction. Moreover, the perpendicular vector defines a plane in three dimensions, and we would need to rely on other quantities to determine the direction. Therefore, we adopted a method based on the Taylor expansion of the velocity field at any point x.

Let us consider two points, x₀ and x₁, which we assume to be along the radial direction of the vortex and relatively close to each other. The scenario is presented in Fig. 3. We can expand the velocity field in x₁, u₁, around x₀ as

$$v_1=v_0+\nabla v|{}_{x_0}ds+O\left({\left|ds\right|}^2\right),$$(13)

where ds = x₁ − x₀ and $\nabla v{|}_{x_0}=U{|}_{x_0}$ is the velocity gradient tensor computed in x₀. If |ds| = |x₁ − x₀| is sufficiently small, we can neglect the higher order terms. Moreover, the two points of interest are aligned along the radial direction by construction, hence, ds is proportional to e_r.

Both the velocity gradient tensor $U{|}_{x_0}$ and the distance vector ds are Galilean invariant, therefore we can safely set v₁ = 0. Moreover, if x₀ belongs to a vortical region, then $U{|}_{x_0}$ is diagonalizable (see, e.g., Canivete Cuissa & Steiner 2020) and therefore also invertible. Hence, we can write Eq. (13) as:

$$ds=U^{-1}{|}_{x_0}{v}_0,$$(14)

and, since ds ∝ e_r, we have a simple formula to compute the radial direction. It is important to notice that Eq. (14) alone only yields information about the straight line connecting the center of curvature and the point x₀. Therefore we also need to compute the radius of curvature to find the coordinates of the EVC.

Fig. 3 Representation of a vortical velocity streamline. The velocity field at two points, x0 and x1, is shown by green arrows and denoted as v0 and v1, respectively. The two points are relatively close to each other and the vector joining them is labeled ds. The point x0 is at distance r from the vortex center that is marked with a star. The purple unit vector, er, represents the radial direction.

3.1.2 The curvature radius

In Sect. 2.4, we have accurately estimated the rotation period of the Lamb–Oseen vortex using the Rortex criterion R, that is, $T_{\text{R}}=\frac{4\pi }{\left|R\right|}$. Therefore, the curvature radius, r, at the point x₀ can be estimated by combining Eq. (11) with the rotational velocity as follows:

$$\frac{4\pi }{\left|R_0\right|}=\frac{2\pi r}{\left|v_{\theta }\right|},$$(15)

where v_θ is the tangential component of the velocity field. Assuming v_θ = |v₀|, where |v₀| is the norm of the velocity field in x₀, Eq. (15) can be expressed as

$$r=\frac{2\left|v_0\right|}{\left|R_0\right|},$$(16)

and the curvature radius can be estimated based alone on the norm of the velocity field and the Rortex criterion in x₀.

The assumption that v_θ = |v₀| is valid if the flow has no radial component. In more realistic cases, Eq. (16) still holds as a first approximation since the tangential component of the velocity field must dominate over the radial one for the flow to be part of a vortex system. The estimation can probably be improved by considering the stretching, compression, and shearing components of the velocity gradient tensor, as seen in Sect. 2.3.

Finally, combining Eqs. (14) and (16), we are left with two points that are at a distance r from x₀ and on the straight line defined by e_r. To decide which one we ought to select as the EVC, we need to check the orientation of the local curvature given by the sign of the swirling strength⁵. If the swirling strength is positive, that is, the flow is curved in the anticlockwise direction, then the correct point is the one on the left with respect to the direction of the flow. If the sign of the criterion is negative, then the flow is is curved in the clockwise direction and the point on the right is to be selected.

In Fig. 4, we apply the method, step by step, to the same Lamb–Oseen vortex of Sect. 2.4. First, we compute the Rortex criterion on the whole grid (panel a). Then, wherever R ≠ 0, we estimate the radial direction vectors (panel b) and the curvature radii (panel c) according to Eqs. (14) and (16), respectively. Finally, we combine the results as discussed above and obtain the EVC map (panel d). As we can see, the method estimates very precisely the position of the vortex center as all the EVCs are clustered around it. The mean error of the estimation is of the order of the grid cell size, and it is mainly due to discretization errors in the numerical computation of the velocity field derivatives.

In more complicated scenarios, where turbulence or noise can distort the velocity field at the grid cell level, the accuracy of this method will not be as impressive as in Fig. 4. In such cases, a “multiple stencil” approach can be used to improve the method’s robustness. The idea is to compute and assemble together multiple EVC maps to help discern between true clusters and noisy points. To obtain different EVC maps, we compute the necessary derivatives using different stencils of grid cells but maintaining the same order of accuracy.

An ordinary, second order, central finite difference derivative uses three adjacent grid cells, which locations on the grid can be arbitrarily defined as −1, 0, and 1. We dub this derivative a stencil-1 derivative, because the distance between the used grid cells is 1 grid cell size. A stencil-2 derivative of the same type operates over grid cells separated by two grid cell sizes, that is at locations −2, 0, and 2. Both derivatives are second order in accuracy, but the stencil-1 derivative is better suited to capture variations in the velocity field at the grid cell size level, while larger stencils can be used to infer large-scale variations. By combining EVC maps computed with different stencils, we can discern small-scale perturbations from actual vortical structures and improve robustness against noise.

3.2 The algorithm

The method presented in Sect. 3.1 can be used as it is for a manual identification of vortices by searching for clusters of EVC points. However, when it comes to analyze large and turbulent velocity fields, as for example the ones resulting from numerical simulations of convection, it can be useful to have an algorithm that processes the EVC map and automatically finds the vortices in it. In the following, we present the automated detection algorithm based on the EVC map method.

3.2.1 Clustering method

Given the nature of the EVC map, we opt for a clustering algorithm. The idea is to classify EVCs into categories, or clusters, based on their distance relatively to each other. Each cluster represents a vortex, for which we can compute different properties. A similar approach has been adopted by Dakanalis et al. (2021) to automatically identify vortices from clusters of curved streamlines.

For this work, we chose to use the clustering by fast search and find of density peaks (CFSFDP) proposed by Rodriguez & Laio (2014). The basic assumption of the algorithm is that cluster centers can be defined as locations characterized by relatively high local densities of datapoints (here, EVCs) that are at relatively large distance from other datapoints with higher local density. Mathematically, these two criteria rely on the computation of two quantities: the local density of datapoints, ρ, and the spacing from other datapoints of higher local density, δ. Both of them depend only on the Euclidean distance, d_ij, between data points i and j, that is, the distance between two EVCs i and j. Therefore, no human interaction is required and the algorithm has shown high performance with different datasets and robustness against noise (Rodriguez & Laio 2014).

The local density, ρ_i, relative to each datapoint i is computed as:

$${\rho }_i={\displaystyle \sum _j\chi \left(d_{ij},d_{\text{c}}\right),}$$(17)

where χ is a kernel that depends on the Euclidean distance, d_ij, to neighboring points, j, and on an arbitrary cutoff distance, d_c. The diagonal terms, d_ii, of the distance matrix, that is, the self-distances, are 0 by definition. The role of the kernel χ is to evaluate the proximity of two data points. For this, Rodriguez & Laio (2014) used a cutoff kernel:

$$\chi \left(d_{ij},d_{\text{c}}\right)=\{\begin{array}{ll}1,\hfill & \text{if\hspace{0.17em}}{d}_{ij}-d_{\text{c}}<0,\hfill \\ 0,\hfill & \text{else,}\hfill \end{array}$$(18)

so that the density ρ is the number of datapoints (EVCs) within a hyper-sphere of radius, d_c. Another possibility is to use a Gaussian kernel of the type:

$$\chi \left(d_{ij},d_{\text{c}}\right)=\exp \left(-\frac{d_{ij}^2}{d_{\text{c}}^2}\right).$$(19)

The performance of the different kernels depend on the choice of the parameter d_c. As a rule of thumb, Rodriguez & Laio (2014) suggested choosing the cutoff parameter, d_c, such that the average number of neighbor data points, namely, those for which d_ij < d_c, is around 1 to 2% of the total number of datapoints. Mehmood et al. (2016) proposedinsteadanonparametricmethod to estimate the cutoff distance.

The spacing δ_i criterion is given by the distance between the generic datapoint i and the closest datapoint j with a higher local density, that is ρj > ρ_i. Mathematically, this statement can be written as

$${\delta }_i=\underset{j:{\rho }_j>{\rho }_i}{\min }\left(d_{ij}\right).$$(20)

This formulation ensures that only datapoints characterized by a local maximum in local density have a δ criterion considerably larger than the typical distance between neighboring datapoints. Therefore, datapoints with high density, ρ, and spacing, δ, are potential cluster centers. Equation (20) is applicable to all datapoints but the one with maximum local density. For this datapoint, the δ criterion is usually set to be δ = max (d_ij).

Fig. 4 Demonstration of the application of the EVC method to the Lamb–Oseen vortex of Fig. 1. Steps are as follows: (a) computation of the Rortex criterion, (b) estimation of the radial direction, (c) estimation of the curvature radius, and (d) computation of the EVCs. The blue-dashed circles denote the boundary of the vortex defined by r = rmax.

3.2.2 Grid and vortex adaptation

The CFSFDP method is mainly limited by its computational cost, which scales as $O\left(N_{\text{p}}^2\right)$, where N_P is the number of datapoints, mainly because of the local density computation (Sieranoja & Fränti 2019). Therefore, dealing with large datasets can be troublesome⁶. To work around this problem, we decided to transpose the algorithm to a grid-based approach. This allows us to considerably reduce the number of datapoints and to filter out noisy points in sparse areas (see, e.g., Xu et al. 2018). Moreover, our initial dataset is already grid-based, since we use velocity fields from numerical simulations or observations. Therefore, it is natural to force the EVCs and the resulting clusters to be defined on the same grid.

Another important point is that we seek to identify vortices, which can have two orientations with respect to an axis, whether clockwise or counter-clockwise. However, the CFSFDP method evaluates all datapoints equally and does not consider the orientation of the local curvature they stem from. This can lead to an incorrect single identification when two vortices with opposite orientation stand close to each other. To avoid this situation, we classify EVCs as clockwise or counter-clockwise based on the orientation of the flow curvature. This can be easily achieved by checking the sign of the swirling strength, λ, as done in Sect. 3.1.2. In this way, EVCs belonging to each one of the two classes are separately clustered.

The grid-vortex adaptation of the CFSFDP algorithm proceeds as follows. The EVCs coordinates are rounded to coincide with the grid of the original velocity field. Then we count how many EVCs there are in each grid cell, differentiating between clockwise and counter-clockwise EVCs. Every EVC associated with a clockwise curvature counts as −1, while counter-clockwise EVCs count as +1. The resulting sum for each grid cell is stored in the grid cardinality s. In this way, large amounts of clockwise EVCs rounded to the same grid cell generate a high negative cardinality s value, while a high positive value of s corresponds to a cluster of counter-clockwise EVCs. We ignore all grid cells that have |s| ≤ 1, since they do not contain any EVC or the one they contain can be considered as noise. From here on, all grid cells having |s| > 1 are labeled as grid-estimated vortex centers (G-EVCs) and their collection forms the G-EVC map.

The next step is to compute the local density ρ of G-EVCs instead of single EVCs. To do so, we modify Eq. (17) to take into account the grid cardinality value and the sign of each G-EVC, which for a generic point i is:

$$\begin{array}{lll}{\rho }_i^{+}\hfill & =\hfill & \left|s_i\right|+{\displaystyle \sum _{j,s_j>0}\left|s_j\right|\chi \left(d_{ij},d_{\text{c}}\right)}\hfill \\ {\rho }_i^{-}\hfill & =\hfill & \left|s_i\right|+{\displaystyle \sum _{j,s_j<0}\left|s_j\right|\chi \left(d_{ij},d_{\text{c}}\right),}\hfill \end{array}$$(21)

where ${\rho }_i^{+}$ and ${\rho }_i^{-}$ are the local densities of G-EVCs with positive (counter-clockwise) and negative (clockwise) cardinality, respectively. For the kernel, χ, we can use either Eq. (18) or Eq. (19). The first term, that is, |s_i|, accounts for the number of EVCs contained in each G-EVC, i, while the sum accounts for the local density due to neighboring G-EVCs. For the δ criterion, we use Eq. (20) but limit the calculation to the G-EVCs having the cardinality of the same sign. Hence, as for the density, we get two spacing criteria ${\delta }_i^{+}$ and ${\delta }_i^{-}$, which refer to the two classes (orientations) of G-EVCs.

This procedure greatly reduces the required computational cost while keeping a high level of precision. To speed-up the calculation even more or improve the accuracy, the grid can (in principle) be respectively coarsened or refined at will. Moreover, the clustering can be performed independently for clockwise and counter-clockwise vortices, therefore, this process can be parallelized.

Figure 5 demonstrates this procedure for an artificial flow composed of two Lamb–Oseen vortices of opposite orientation. Moreover, we added a random velocity field smoothed with a Gaussian filter, which serves as noise. The two vortices are defined through Eq. (10). In more detail, the clockwise vortex is smaller in size but stronger in rotational speed with Lamb–Oseen parameters α = 1.256, r_max = 0.07, and v_max = -0.9, while the counter-clockwise has α = 1.256, r_max = 0.2, and v_max = 0.6; hence it is larger but also slower. The average strength of the noise velocity field is |v_noise| = 0.226. The Cartesian grid has 200 × 200 grid points and covers an arbitrary domain of a size of 1.0 × 1.0.

The Rortex criterion captures the curvatures in the flow induced by the two vortices, as demonstrated in the top panel by the two large or intensive patches of positive (green) and negative (pink) values of R, respectively. The random motions caused by the Gaussian noise induce spurious signals (R ≠ 0) in the Rortex map. These signals are not related to the vortices, but to perturbations in the velocity field which cause the presence of local curvature in the flow. However, the flow in these places does not form a full circular or spiral pattern, and therefore these signals should not be regarded as vortices. If every connected region where R ≠ 0 is identified as a vortex structure, then all these spurious signals erroneously represent different vortical regions. Hence, the necessity to consider the global properties of the flow in a more elaborated identification algorithm.

The bottom panel of Fig. 5 shows the grid cardinality, s, of the G-EVC map. We notice how most of the EVCs cluster near the centers of the two vortices, resulting in G-EVCs with high positive values of s near the center of the counterclockwise vortex and high negative values around the clockwise one. The estimation is not as precise as the one shown in Fig. 4, as many points are scattered in the surroundings of the vortex centers. This is due to the presence of noise and to the interaction between the flows of the two vortices. Nevertheless, the peaks of |s| clearly indicate the presence of two vortex centers. The location of these centers correspond to the approximate location of the original centers of the two vortex models, which are indicated by two blue circles in Fig. 5.

In addition, there are small patches of G-EVCs scattered in the map. These points show uniformly low values of s, which is precisely the behavior expected for incoherent local curvature perturbations in the flow. Random curvatures in the flow show up in the mathematical criteria for the identification of vortices, as seen in the top panel of Fig. 5, but will not be selected as cluster centers since they are characterized by a relatively low density ${\rho }_i^{\pm }$. The decision process is subject of the following subsection.

Fig. 5 Application of the grid-adapted EVC method to an artificial flow composed of two Lamb–Oseen vortices of opposite orientation superposed by a Gaussian noise signal. Top panel: Rortex map. Bottom panel: G-EVC map. The cardinality s accounts for the number of EVC points which coordinates fall in the same grid cell. The blue circles show the approximate location of the two vortex centers. The regions where R ≠ 0 are outlined by a gray contour in the Rortex map.

3.2.3 Decision

The next step is to select the cluster centers based on the local density ρ = [ρ⁺, ρ⁻] and the spacing δ = [δ⁺, δ⁻] quantities. A universal selection criteria is still under debate. A common strategy is to define a couple of thresholds, ρ_th and δ_th, based on the distribution of ρ and δ within the dataset (Rodriguez & Laio 2014). Every point having local density and spacing larger than both thresholds is selected as a cluster center. However, the choice of the thresholds is arbitrary and depends on the type of dataset. For example, the local density and spacing thresholds can be defined as:

$$\begin{array}{cc}{\rho }_{\text{th}}=P_{\rho }{\mu }_{\rho },& {\delta }_{\text{th}}=p_{\delta }{\sigma }_{\delta },\end{array}$$(22)

where µ_ρ is the mean value of the local density ρ, σ_δ is the standard deviation of the spacing δ, and p_ρ, p_δ > 0 are free parameters. Moreover, to avoid biased thresholds ρ_th and δ_th when evaluations Eq. (22) in the grid-adapted version of the algorithm, we also need to take the number of datapoints (EVCs) hidden in each G-EVC into account. Therefore, we should add s_i − 1 mock datapoints with ρ = ρ_i and $\delta =\frac{1}{2}\text{Δ}x$ values for each G-EVC i, where ∆x is the grid cell size. In this way, each G-EVC i will be represented by one datapoint with ρ = ρ_i and δ = δ_i in the cluster center decision process, while the mock datapoints with ρ = ρ_i and $\delta =\frac{1}{2}\text{Δ}x$ account for the remaining EVCs, which are considered as close neighbors.

Alternatively, we can set a threshold for the γ criterion, where γ = ρδ. This quantity naturally highlights datapoints having both local density and spacing relatively high. The choice of the threshold is once more arbitrary and we opt for the following formulation:

$${\gamma }_{\text{th}}=p_{\gamma }{\delta }_{\min }{\rho }_{\max },$$(23)

where ρ_γ > 0 is again a free parameter, while δ_min and ρ_max are the minimum and the maximum of the spacing, δ, and local density, ρ, distributions, respectively. This choice ensures that, for ρ_γ > 1, datapoints with large local density but very small distance (that is close-neighbors) are not selected as cluster centers. Moreover, it relies on only one free parameter. Another possibility is to sort the values of γ in increasing order and identify a knee point in the curve, as proposed by Wang et al. (2016), but we have not yet tested this option so far.

Figure 6 shows two commonly used decision graphs for the visualization of the cluster centers selection (Rodriguez & Laio 2014). The top panel shows the normalized spacing criterion, δi/δmax, and the normalized density, ρ_i/ρ_max, for each datapoint, while in the bottom panel the γ criterion is plotted for each datapoint in increasing order of γ. The thresholds in ρ, δ, and γ, defined through Eqs. (22) and (23), are indicated as dotted and dashed orange lines, respectively. The chosen parameters are p_ρ = 1.0, p_δ = 0.5, and p_γ = 15.0. There are two datapoints, marked with a star, that satisfy both selection criteria and are therefore identified as cluster centers.

Fig. 6 Decision graphs for the artificial Lamb–Oseen vortex flow of Fig. 5 based on the ρ − δ thresholds of Eq. (22), (top) and on the γ threshold of Eq. (23), (bottom). In the bottom panel, n represents the datapoint number in increasing order of γ. Dotted lines represent the ρ and δ thresholds, while dashed lines show the γ threshold. Each point corresponds to a EVC. The green and pink starred EVCs satisfy both thresholds and are therefore selected as cluster centers.

3.2.4 Vortex identification

For each cluster center that fulfills the decision criteria set by the thresholds Eqs. (22) or (23), we proceed with the formation of an associated cluster of datapoints (in our case, G-EVCs). For this step, we assign each one of the remaining datapoints that have not already been identified as a cluster center to the same cluster of its nearest neighbor of highest density (Rodriguez & Laio 2014). The single EVCs are assigned to the same cluster of the G-EVCs they belong to.

At the end of this process, every G-EVC and consequently every EVC will belong to a cluster. Since there is a one-to-one correspondence between the EVC points and the grid cells they derive from, each cluster of EVCs corresponds to a group of grid cells showing local curvature in the velocity field. Therefore, each group of grid cells corresponds to a vortex candidate with its center given by the associated G-EVC cluster center. Figure 7 shows the two groups of grid cells (vortex candidates) found for the artificial flow of Fig. 5. The cluster centers found in Fig. 6 represent the centers of the identified candidate vortices and are marked with stars. The two main groups of grid cells belonging to the Lamb–Oseen vortical structures are identified by the method.

3.2.5 Noise removal

Figure 7 also shows patches of grid cells that have already been identified as noise. These are grid cells characterized by R ≠ 0, as shown in Fig. 5, that have been discarded because their EVCs are found either outside of the computational domain or alone in a grid cell of the G-EVC map. Therefore, we can assume that the curvature of the flow identified in these cells is not part of a coherent vortical structure. The associated EVCs are not considered in the clustering process.

There are also a few isolated patches of grid cells that have been associated to one of the candidate vortices by the clustering algorithm although they are clearly not related to them. These grid cells have not been discarded yet because their EVCs are accidentally associated with one of the two main clusters, or because they are part of other local clusters which however do not satisfy the decision criteria of the clustering algorithm. Indeed, a few more groups of G-EVCs with low grid cardinality |s| ≲ 20 are visible in the bottom panel of Fig. 5. Since these EVCs are considered in the cluster identification process, they are also inevitably associated to a cluster. Therefore, the corresponding grid cells belong to one of the candidate vortices. We seek to label these grid cells as noise as well and remove them from the corresponding candidate vortex.

In order to do so, we consider a grid cell that has been associated to a candidate vortex as noise if either its position or its EVC is “sufficiently distant” from the center of the cluster they belong to. In this way, grid cells can be removed if they or their EVCs are spatially unrelated to the vortex structure they supposedly belong to. We define “sufficiently distant” as a distance relative to the effective radius of the candidate vortex, r_eff, which is estimated based on the area covered by the grid cells belonging to that structure as:

$$r_{\text{eff}}=\sqrt{\frac{N_{\text{c}}}{\pi }}\text{\hspace{0.17em}Δ}x,$$(24)

where N_c is the number of grid cells that form the candidate vortex and ∆x is the grid cell size. Hence, grid cells for which the distance between the cluster center and their position or their corresponding EVC location are larger than p_r r_eff, where p_r is a free parameter, are labeled as noise and removed from the candidate vortex structure. We note that this is an iterative process because the removal of grid cells changes N_c and with it, the criterion p_rr_eff.

Finally, we try to identify and remove candidate vortices that do not show full circular patterns in their flow. For example, a sharp deviation in the flow induces some degree of local curvature (R ≠ 0) and, possibly, a cluster of EVCs, but it should not be classified as a vortex because the trajectories of particles in such a flow do not spiral. To differentiate non-spiraling coherent curvatures in the flow from actual vortices, we make use of the fact that grid cells belonging to a perfect vortex are isotropically distributed around the cluster center, while for non-spiraling coherent curvatures a strongly asymmetric distribution can be expected. Therefore, the effective radius computed with Eq. (24) results in a good estimate of the size of a vortical structure only, while for non-spiraling coherent curvatures, it underestimates their radius. Repeating the procedure for removing distant grid cells while updating with each iteration the number of grid cells, N_c, and the effective radius, r_eff, for a parameter, p_r, sufficiently close to 1, a cluster stemming from a non-spiraling coherent curvature in the flow should be depleted of all grid cells composing it, while a vortex should maintain approximately its total number of grid cells and size.

Figure 8 shows the final result, where the outlying patches of grid cells have been classified as noise with the procedure just described, while the two Lamb–Oseen vortices have been successfully identified. The effective radii of the two vortices, computed using Eq. (24), are 0.19 for the anti-clockwise vortex and 0.08 for the other one. Considering the interaction between the two vortices and the perturbations due to the random velocity field, the estimated radii are very close to the initial values, which are 0.20 and 0.07, respectively. The location of the vortex centers cannot be determined precisely because of the random Gaussian velocity field that perturbs the flow and therefore shifts, although only slightly only, their original position. Nevertheless, the identified vortex centers are found in the vicinity of the original positions of the Lamb–Oseen vortices, marked by blue circles in Fig. 7, which confirms the fitness of the method.

We implemented the new vortex identification method and the automated algorithm in a Python package named “SWirl Identification by Rotation-centers Localization” (SWIRL). The code is open source and can be found online⁷. In Appendix B, we review the structure of the code with a step-by-step flowchart.

Fig. 7 Groups of grid cells of candidate vortices (green and pink pixels) obtained by applying the grid-adapted CFSFDP algorithm to the G-EVC map of the artificial flow shown in Fig. 5. The cluster centers are marked with a star, while the blue circles show the approximate location of the two Lamb–Oseen vortex centers. Noisy grid cells are gray shaded.

Fig. 8 Vortices identified in the artificial flow shown in Fig. 5. The vortex centers are marked with a star, while the circles show their effective size according to Eq. (24). Grid cells determined to be noise are gray shaded.

4 Application and discussion

In this section, we test the vortex identification algorithm on more complex vortical flows than those of Sect. 3. First, we apply the SWIRL code to an artificial flow composed of nine vortices of different sizes and strengths and a Gaussian noisy signal. Then, in order to test the reliability of the code on a more turbulent flow, we apply it to a two-dimensional simulation of a magneto-hydrodynamical (MHD) Orszag-Tang vortex.

4.1 Artificial vortex flow

The artificial flow is composed of nine random Lamb–Oseen vortices and a Gaussian smoothed noisy velocity field as shown in Fig. 9. The parameters r_max, v_max, and α of the Lamb–Oseen model, as well as the coordinates of the vortex centers, are randomly generated. The size of the grid is 200 × 200 points covering a domain of 1.0 × 1.0 in arbitrary units. In the top panel, we show the Rortex map derived from the composed velocity field. After the computation of the EVC map, we obtain the grid cardinality s map shown in the bottom panel of the same figure. The approximate coordinates of the Lamb–Oseen vortex centers are marked with blue circles in both panels. We notice peaks of positive and negative grid cardinality s in correspondence to counter-clockwise and clockwise curved flows, respectively. We then proceed to an automated identification of vortices with the method presented in Sect. 3 and implemented in the SWIRL code.

The resultis shown in the top panel ofFig. 10. Forsimplicity, we show the identified vortices with a colored disk representing their effective circular area, instead of showing the individual grid cells that compose them. The size of the disk is estimated with Eq. (24), while the color indicates their orientation. Grid cells that have been identified as noise are instead colored in gray. Our algorithm identifies nine vortices in total, but only eight of those are related to the artificial Lamb–Oseen vortices inserted in the velocity field. The vortex placed at coordinates (x, y) ~ (0.20,0.75) has been discarded and recognized as noise. Inspecting the velocity field around that location with the help of the vector plot reveals that the flow there performs a non-spiraling coherent curvature. This is due to the presence of a nearby stronger vortex at coordinates (x, y) ~ (0.20,0.60), identified by the algorithm and labeled as v6. Since both vortices have the same orientation, the velocity field is dominated by the stronger one (vortex v6) and the flow is distorted around the weaker one.

The algorithm however identifies a vortex, labeled v5, around (x, y) ~ (0.05,0.95), that is not associated with any of the artificially generated vortices, but stems from the background random velocity field. The bottom panels of Fig. 10 show close-ups on the nine detected vortices with instantaneous streamlines of the velocity field instead of vector plots. The identified vortices are characterized by bounded streamlines around the vortex centers, which are marked by colored stars.

The identification of vortex v5 demonstrates an important feature of our algorithm. From the vector plot of Fig. 9, this vortex is not easily recognizable from visual inspection because the velocity field in that region is weak compared to the rest of the domain. Moreover, the values of the Rortex criterion are relatively small because of the weak angular velocity, as we can see from Fig. 9. Traditional algorithms based on mathematical criteria might miss this vortex, since for those a threshold in the criteria is likely put in place to filter out noise. In that case, the signal due to the vortex would be removed. However, the method presented in this paper does not require a threshold on the Rortex criterion; hence, it can also detect weak vortical flows.

Fig. 9 Artificial flow with nine different Lamb–Oseen vortices superposed by a Gaussian smoothed random velocity field. Top panel: Rortex map. Bottom panel: G-EVC map. Blue circles denote the position of the Lamb–Oseen vortex centers. The velocity field is represented by a vector field in both panels. The regions where R ≠ 0 are outlined with a gray contour in the Rortex map.

4.2 Orszag-Tang vortex test

The Orszag-Tang vortex problem (Orszag & Tang 1979) is one of the most widely known benchmark tests for MHD numerical schemes. It is a 2D problem and the initial conditions on a domain of size 1.0 × 1.0 are given as follows:

$$\begin{array}{lll}v\hfill & =\hfill & \left(-\sin \left(2\pi y\right),\sin \left(2\pi x\right)\right),\hfill \\ B\hfill & =\hfill & \left(-B_0\sin \left(2\pi y\right),B_0\sin \left(4\pi x\right)\right),\hfill \\ \rho \hfill & =\hfill & \gamma p_0,\hfill \end{array}$$(25)

where B₀ = (4π)^−1/2, γ = 5/3, and p₀ = 5/(12π). The initialized flow is unstable and breaks down into turbulence very quickly.

We ran the Orszag–Tang problem with the CO5BOLD code (Freytag et al. 2012) on a 512 × 512 grid. The pressure p and the large scale velocity field at t = 1.0 are shown in Fig. 11.

Multiple vortices are expected to form in the simulated turbulent flow. Given the complexity of the dynamics and the influence of magnetic fields, these vortices are most likely different from the Lamb–Oseen models we used in Sect. 4.1. Nevertheless, our approach is not limited to Lamb–Oseen vortex models since the Rortex criterion measures, by definition, the rigid-body rotational component of any velocity field. Therefore, we expect our methodology to be successful also in more complex flows, such as the Orszag–Tang vortex test or realistic simulations of the solar atmosphere. For example, Silva et al. (2020) showed that the profiles of the tangential velocity of vortices identified in numerical simulations of the solar atmosphere can be well fitted by a cubic polynomial.

We applied the SWIRL code to the simulated velocity field at t = 1 .0 and the results are shown in Fig. 12. To enhance the robustness of the identification process, we employed the “multiple stencil” approach (see Sect. 3.1.2) with stencils of 1, 3, 5, and 7 grid cells. In total, 37 swirls are identified in a mirror-like disposition, which evidences the symmetry of the problem. In the bottom panels of the same figure, we can appreciate the turbulent nature of the flow as shown by the instantaneous streamlines. Almost all the identified vortices correspond to locations in the flow where the instantaneous streamlines form circular or spiral patterns. The only exception is the one at the top of the B panel, where the flow seems to perform a sharp U turn but not a full vortical motion.

It is important to notice that the streamlines shown in the bottom panels of Fig. 12 do not reflect the trajectories that test particles would follow in such a flow, but, rather, they stem from the instantaneous velocity field. Therefore, we must be cautious when investigating the morphology of velocity fields from a single time instance, particularly when it changes on a timescale that is relatively short compared to the rotational period of the vortical flow (see, e.g., Shelyag et al. 2013). Nevertheless, we believe that the results demonstrate that the SWIRL code is able to correctly identify most of the vortical motions present in the simulated flow, even in the presence of turbulence and magnetic fields. Applications to solar physical flows are planned for a follow-up paper.

Fig. 10 Identified vortices in the artificial flow of Fig. 9. Top panel: position and the size of nine identified vortices. Green colored vortices rotate counter-clockwise, the pink ones rotate clockwise. The velocity field is displayed by a vector plot. Grid cells labeled as noise are shown in gray. Bottom panels: close-ups on the different vortices. The velocity field is represented by streamlines. Stars denote the center of the vortices and the radii are computed according to Eq. (24).

Fig. 11 Pressure, p, map of the two-dimensional MHD Orszag-Tang test at t = 1.0 carried out with the CO5BOLD code. The velocity field is represented by white arrows.

Fig. 12 Identified vortices in the Orszag-Tang test flow at t = 1.0. Top panel: position and the size of the identified vortices. Green colored vortices rotate counter-clockwise, while the pink ones rotate clockwise. Bottom panels: close-ups on four different regions of the flow with identified vortices. The velocity field is represented by instantaneous streamlines. Thesizeofeachvortexis computedaccording to Eq. (24).

5 Summary and conclusions

In this paper, we present a new method and an automated algorithm for the identification of vortical motions based on the velocity field alone. The core computation of this method is an estimation of curvature centers for all points in the velocity map exhibiting some degree of local curvature. For that, we combine an innovative and highly reliable mathematical criterion, the Rortex, with the global information of the flow derived from the velocity field. Hence, the method can be considered a hybrid between classical methods based on mathematical criteria and morphological methods. To our knowledge, this is the first time that such type of a vortex detection algorithm has been suggested, both in the domain of astrophysics and of computational fluid dynamics. We classify the estimated centers of rotation, or EVCs, with a grid-adapted version of the modern clustering algorithm CFSFDP. This automated procedure can piece together all the points that share a common center of rotation and therefore form a vortex. A first implementation of the algorithm, called SWIRL, is open source and can be found online.

We tested our code on a non-trivial artificial velocity field composed of nine different Lamb–Oseen vortices and a background Gaussian noise. The algorithm correctly detected all vortical motions of different angular velocities and sizes present in the flow, even one that was very weak and was generated by the random perturbations of the background noise. We also tested the SWIRL code with a MHD numerical simulation of a Orszag–Tang vortex system. The simulated flow is turbulent in nature and influenced by the magnetic field; hence, the identification process is much more complex than in the artificial case and misidentifications can be expected. Nevertheless, the results demonstrate the high level of accuracy of our methodology.

Therefore, the SWIRL code is a reliable tool for the study of vortical motions in complex and magneto-hydrodynamical astrophysical systems. It is robust against noise, turbulence, and shear flows; at the same time, weakly rotating vortices do not pose problems as thresholds in the Rortex values are not necessary. The vortex center is naturally given by the EVCs cluster center, while proper boundaries can be defined through the collection of grid cells composing the vortex structure in accordance with the definition of Lugt (1979) that a vortex is “a multitude of material particles (grid cells) rotating around a common center.” In this way, we can determine physical quantities of interest in any point of the vortical region, allowing for a trustworthy analysis of the vortex properties and dynamics. Finally, the SWIRL code can be directly applied to the velocity field of numerical simulations as well as to previously derived velocity fields of observations.

Regarding the computational performance, the method is very fast on small to medium grid ranges. For example, the full identification process carried out on a 200 × 200 grid for the artificial flow of Sect. 4.1 takes ~1 s on a single CPU. For the Orszag–Tang test instead, given the complexity of the flow and the larger size of the grid (512 × 512), the computation took ~60 s. It can become relatively slow on large domains, especially when numerous vortices are present, the main cause being the clustering step. However, we can always adopt a “divide and conquer” approach and run the code on smaller portions of the grid separately and combine the results in the end. In principle, it is also possible to parallelize the process. Finally, it is worth to mention that no cut-off value for the Rortex criterion needs to be imposed as done in many methods based on mathematical criteria, because EVCs due to noisy signals are scattered randomly in the domain and do not form clusters.

There are, however, two drawbacks to be noted. First, the method is not strictly Galilean invariant since the vortex structure is required to be at rest with respect to translations for the accurate estimation of its center of rotation. For solar applications, this is not much of a concern since vortices are usually anchored within intergranular lanes or vertices of intergranular lanes and move slowly relative to the vortical flow speed (see, e.g., Tziotziou et al. 2022). However, it could be a matter of concern in more dynamical scenarios. Second, the parameters of the clustering algorithm call for some fine tuning to perform as expected.

The algorithm can be further developed and improved. First of all, our implementation is limited to two dimensions. In principle, the computation of a fully three-dimensional EVC map should be straightforward, as Eqs. (14) and (16) are valid in three-dimensions as well. In that case, we would obtain a three-dimensional distribution of EVC points. However, the clustering process is likely to be computationally extremely expensive, because of the extra dimension and the consequently larger datasets. A thoroughstudy onmulti-dimensionalclustering algorithms is therefore required for a possible development in this direction. On the other hand, an improvement of the current status of the algorithm is also possible. The grid-based implementation of the CFSFDP clustering algorithm can be further refined to make it more accurate and robust. In particular, we would like to reduce the number of parameters required in order to make it as user-independent as possible. Techniques as the one proposed by Wang et al. (2016) could be implemented in the future to attain this goal.

Another area of possible improvement regards the accuracy of the computed EVCs. Indeed, Eqs. (14) and (16) are only approximate estimations of the radial direction and curvature radius, especially in complex and turbulent flows such as solar atmospheric ones. Improving the accuracy of these quantities, even slightly, would affect the clustering process and, ultimately, the performance of our method.

In conclusion, the method presented in this paper represents a new and reliable technique for the detection and analysis of vortices in numerical simulations of turbulent and (magneto)hydrodynamical flows. The publicly available implementation allows for a simple and quick usage of the algorithm. In a follow-up paper, we intend to apply our method to realistic numerical simulations of the solar and stellar atmospheres carried out with the CO5BOLD code. In the future, we would also like to compare our method to other commonly used automated vortex identification algorithms and to draw rigorous statistics of vortex populations from high-resolution observations of the solar atmosphere.

Movie

Movie (aa43740_22_animation) Access here

Acknowledgements

The authors acknowledge support by SNF under grant ID 200020_182094. 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” and with the team of P. Keys (convener) “WaLSA: Waves in the Lower Solar Atmosphere at High Resolution” (www.walsa.team) at the International Space Science Institute (ISSI). We would like to sincerely thank the anonymous referee for the constructive comments and F. Riva (IRSOL) for providing us with the Orszag-Tang test model.

Appendix A Analytical comparison between mathematical criteria

In Sect. 2.4, we show numerically that the Rortex is a mathematical criterion that accurately measures the rigid-body rotational part of the velocity field of a Lamb Oseen vortex model. However, we can prove it analytically as well, and generalize this results to different kinds of vortical motions by considering a few examples.

We start with the simplest vortex model, namely the rotational vortex. Without loss of generality, we center the vortex flow in (x, y) = (0,0), such that the associated velocity field can be defined as:

$$\begin{array}{lll}{v}_x\hfill & =\hfill & -\alpha y,\hfill \\ v_y\hfill & =\hfill & +\alpha x,\hfill \end{array}$$(A.1)

where α > 0 is a constant that defines the strength and the orientation of the rotation. We can compute the vorticity, the swirling strength, and the Rortex criteria with Eqs. (1), (2), and (9), respectively. The resulting values for these criteria are:

$$\omega =\lambda =R=2\alpha .$$(A.2)

Hence, for a simple rotational vortex, the three criteria give identical results. Moreover, the tangential component of the velocity is v_θ = αr, where $r=\sqrt{x^2+y^2}$ is the radius in polar coordinates. Employing Eqs. (11) and (12), we can check that the three criteria yield the correct rotational period, that is T = 2π/α. This is not a surprise since a rotational vortex behaves as a rotating rigid-body and therefore consists of pure rotation.

Next, we use a generalized version of the rotational vortex model where the tangential velocity scales as a power of the radius, that is,

$$v_{\theta }=\alpha r^{\beta },$$(A.3)

where β ∈ ℝ. The rotational vortex model is a special case of Eq. (A.3) with β = 1, while the irrotational vortex model can be obtained with β = −1. For the case of Eq. (A.3), the velocity field in Cartesian coordinates is given by:

$$\begin{array}{lll}{v}_x\hfill & =\hfill & -\alpha r^{\beta -1}y,\hfill \\ v_y\hfill & =\hfill & +\alpha r^{\beta -1}x.\hfill \end{array}$$(A.4)

Using again Eqs. (1), (2), and (9), we obtain the following values for the three criteria:

$$\begin{array}{lll}\omega \hfill & =\hfill & \alpha r^{\beta -1}\left(1+\beta \right),\hfill \\ \lambda \hfill & =\hfill & 2\alpha r^{\beta -1}\sqrt{\beta },\hfill \\ R\hfill & =\hfill & 2\alpha r^{\beta -1}.\hfill \end{array}$$(A.5)

In this case, each criterion yields a different value and therefore predicts a different period of rotation. The true period of rotation is given by Eq. (12) with vθ given by Eq. (A.3), that is, T = 2π/(αr^β−1). Using Eq. (11) to estimate the period of rotation from the different criteria, we obtain:

$$\begin{array}{lll}{T}_{\omega }\hfill & =\hfill & \frac{4\pi }{\alpha r^{\beta -1}\left(1+\beta \right)},\hfill \\ T_{\lambda }\hfill & =\hfill & \frac{4\pi }{\alpha r^{\beta -1}\sqrt{\beta }},\hfill \\ T_{\text{R}}\hfill & =\hfill & \frac{2\pi }{\alpha r^{\beta -1}}\hfill \end{array}$$(A.6)

which demonstrates that the Rortex correctly measures the rotational part of this flow, while the other two criteria are both biased by the presence of intrinsic shears.

It is also interesting to notice that the Rortex criterion can in principle also measure the rotationalpartofanirrotationalvortex (i.e., when β = −1). It is well known that both the vorticity and the swirling strength values are always zero in the presence ofan irrotational vortex, as also shown in Eq. (A.5). The vorticity is null because of the 1 + β term, while the $\sqrt{\beta }$ term in the swirling strength would render the eigenvalue of the velocity gradient tensor purely real, hence, λ = 0. The Rortex criterion does not suffer from these problems. However, in practical applications, we use Eq. (9) to numerically compute its value. Therefore, if ω = 0 and λ = 0, then also R = 0.

Finally, we test the mathematical criteria on a Lamb–Oseen vortex defined through Eq. (10). Without loss of generality, we set α = 1.0, r_max = 1.0, and v_max = 2/3, so that the velocity field in Cartesian coordinates is expressed as:

$$\begin{array}{lll}{v}_x\hfill & =\hfill & -\frac{y}{r^2}\left(1-\exp \left(-r^2\right)\right),\hfill \\ v_y\hfill & =\hfill & +\frac{x}{r^2}\left(1-\exp \left(-r^2\right)\right),\hfill \end{array}$$(A.7)

and the analytical rotational period is:

$$T=\frac{2\pi r^2}{1-\exp \left(-r^2\right)}.$$(A.8)

The calculations needed for the vorticity, swirling strength, and Rortex are slightly more involved, finally result in:

$$\begin{array}{lll}\omega \hfill & =\hfill & 2\exp \left(-r^2\right),\hfill \\ \lambda \hfill & =\hfill & \frac{2\exp \left(-r^2\right)}{r^2}\sqrt{\left(\exp \left(r^2\right)-1\right)\left(1+2r^2-\exp \left(r^2\right)\right)},\hfill \\ R\hfill & =\hfill & \frac{2}{r^2}\left(1-\exp \left(-r^2\right)\right).\hfill \end{array}$$(A.9)

Similarly to the generalized rotational vortex model, each criterion yields a different value and therefore also a different period of rotation. We may check with Eq. (11) that the only quantity predicting the correct value for the rotation period is the Rortex criterion. Therefore, the Rortex is the only reliably criterion for the extraction of physical information on the Lamb–Oseen vortex from the velocity field.

Appendix B Algorithm pseudo-code

The algorithm presented in Sect. 3 is shown in a diagrammatic fashion in Fig. B.1. The identification process is performed over a two-dimensional velocity field, which is the input of the algorithm. We also provide an animation of the step-by-step process on a low resolution test case⁸.

The first goal is to compute the mathematical criterion that is used in the method, which is the Rortex. The Rortex can be computed on the entire grid of the two-dimensional velocity field⁹ according to Eq. (9) for different parameters and with different stencils of grid cells. This step corresponds to the top panel of Fig. 5.

Once the Rortex map is obtained, we proceed with the computation of the (G-)EVC map. For that, one selects the grid cells where the velocity field is characterized by some degree of curvature, that is R ≠ 0, and computes the associated estimated radial direction (Eq. 14) and curvature radius (Eq. 16). The EVCs are computed by combining these two estimated quantities, and the collection of EVC points forms the EVC map. The G-EVC map can then be obtained following the procedure explained in Sect. 3.2.2. The underlying idea is that the G-EVC map should show large positive (negative) values close to the center of counter-clockwise (clockwise) vortices, as shown in the bottom panel of Fig. 5.

Fig. B.1 Flowchart algorithm of the SWIRL code.

At this point, we can employ the grid version of the CFSFDP algorithm to automatically cluster EVC points and find candidate vortices. The first step is to characterize each G-EVC by its local density, ρ, and spacing, δ, defined in Eqs. (21) and (20), respectively. Then, we find the cluster centers using the thresholds for δ, ρ, or γ = ρδ given as parameters. These parameters must be tweaked according to the characteristics of the flow for optimal results. An example of this selection process is given in Fig. 6.

Once the cluster centers have been identified, the clustering of the remaining EVCs follows the process described in Sect. 3.2.4. At the end, every cluster of EVCs corresponds to a groups of grid cells characterized by R ≠ 0 which form the different candidate vortices. Figure 7 shows two candidate vortices obtained from the two main clusters of (G-)EVCs shown in the bottom panel of Fig. 5.

Given a set of candidate vortices, the last step consists of an iterative noise removal procedure. In this way, we can discard misidentified candidate vortices or remove noisy grid cells that were wrongly associated to a vortex. The iterative procedure is described in Sect. 3.2.5. The final outputs of the algorithm are the identified vortices and the grid cells classified as noise, as shown, for example, in Fig. 8. Each vortex consists of a set of grid cells which EVCs are sufficiently close together, meaning that the input velocity field defined on these grid cells is performing a global rotation around a common axis.

References

All Figures

	Fig. 1 Comparison between (a) the vorticity, (b) the swirling strength, and (c) the Rortex criteria for a Lamb–Oseen vortex. Arrows indicate the velocity vector field of the Lamb–Oseen vortex and their length is proportional to the magnitude of the flow. The blue-dashed circles denote the boundary of the vortex defined by r = rmax.
In the text

	Fig. 2 Comparison of the estimated rotational period T derived from the vorticity, the swirling strength, and the Rortex criteria with the true rotational period for a Lamb–Oseen vortex. Blue vertical lines denote the boundary of the vortex defined by r = rmax. The thick gray curve indicates the true rotational period TLO computed from the analytical formula Eq. (10).
In the text

Fig. 3 Representation of a vortical velocity streamline. The velocity field at two points, x0 and x1, is shown by green arrows and denoted as v0 and v1, respectively. The two points are relatively close to each other and the vector joining them is labeled ds. The point x0 is at distance r from the vortex center that is marked with a star. The purple unit vector, er, represents the radial direction.

In the text

	Fig. 4 Demonstration of the application of the EVC method to the Lamb–Oseen vortex of Fig. 1. Steps are as follows: (a) computation of the Rortex criterion, (b) estimation of the radial direction, (c) estimation of the curvature radius, and (d) computation of the EVCs. The blue-dashed circles denote the boundary of the vortex defined by r = rmax.
In the text

Fig. 5 Application of the grid-adapted EVC method to an artificial flow composed of two Lamb–Oseen vortices of opposite orientation superposed by a Gaussian noise signal. Top panel: Rortex map. Bottom panel: G-EVC map. The cardinality s accounts for the number of EVC points which coordinates fall in the same grid cell. The blue circles show the approximate location of the two vortex centers. The regions where R ≠ 0 are outlined by a gray contour in the Rortex map.

In the text

Fig. 6 Decision graphs for the artificial Lamb–Oseen vortex flow of Fig. 5 based on the ρ − δ thresholds of Eq. (22), (top) and on the γ threshold of Eq. (23), (bottom). In the bottom panel, n represents the datapoint number in increasing order of γ. Dotted lines represent the ρ and δ thresholds, while dashed lines show the γ threshold. Each point corresponds to a EVC. The green and pink starred EVCs satisfy both thresholds and are therefore selected as cluster centers.

In the text

	Fig. 7 Groups of grid cells of candidate vortices (green and pink pixels) obtained by applying the grid-adapted CFSFDP algorithm to the G-EVC map of the artificial flow shown in Fig. 5. The cluster centers are marked with a star, while the blue circles show the approximate location of the two Lamb–Oseen vortex centers. Noisy grid cells are gray shaded.
In the text

	Fig. 8 Vortices identified in the artificial flow shown in Fig. 5. The vortex centers are marked with a star, while the circles show their effective size according to Eq. (24). Grid cells determined to be noise are gray shaded.
In the text

	Fig. 9 Artificial flow with nine different Lamb–Oseen vortices superposed by a Gaussian smoothed random velocity field. Top panel: Rortex map. Bottom panel: G-EVC map. Blue circles denote the position of the Lamb–Oseen vortex centers. The velocity field is represented by a vector field in both panels. The regions where R ≠ 0 are outlined with a gray contour in the Rortex map.
In the text

Fig. 10 Identified vortices in the artificial flow of Fig. 9. Top panel: position and the size of nine identified vortices. Green colored vortices rotate counter-clockwise, the pink ones rotate clockwise. The velocity field is displayed by a vector plot. Grid cells labeled as noise are shown in gray. Bottom panels: close-ups on the different vortices. The velocity field is represented by streamlines. Stars denote the center of the vortices and the radii are computed according to Eq. (24).

In the text

	Fig. 11 Pressure, p, map of the two-dimensional MHD Orszag-Tang test at t = 1.0 carried out with the CO5BOLD code. The velocity field is represented by white arrows.
In the text

Fig. 12 Identified vortices in the Orszag-Tang test flow at t = 1.0. Top panel: position and the size of the identified vortices. Green colored vortices rotate counter-clockwise, while the pink ones rotate clockwise. Bottom panels: close-ups on four different regions of the flow with identified vortices. The velocity field is represented by instantaneous streamlines. Thesizeofeachvortexis computedaccording to Eq. (24).

In the text

	Fig. B.1 Flowchart algorithm of the SWIRL code.
In the text