An innovative and automated method for vortex identification. I. Description of the SWIRL algorithm

Context. A universally accepted definition of what a vortex is has not yet been reached. Therefore, we lack an unambiguous and rigorous method for the identification of vortices in fluid flows. Such a method would be necessary to conduct robust statistical studies on vortices in highly dynamical and turbulent systems, such as the solar atmosphere. Aims. We aim to develop an innovative and robust automated methodology for the identification of vortices based on local and global characteristics of the flow. Moreover, the use of a threshold that could potentially prevent the detection of weak vortices in the identification process should be avoided. Methods. We present a new method that combines the rigor of mathematical criteria with the global perspective of morphological techniques. The core of the method consists in the estimation of the center of rotation for every point of the flow that presents some degree of curvature in its neighborhood. For that, we employ the Rortex criterion and combine it with morphological considerations of the velocity field. We then identify coherent vortical structures by clusters of estimated centers of rotation. Results. We demonstrate that the Rortex is a more reliable criterion than are the swirling strength and the vorticity for the extraction of physical information from vortical flows, because it measures the rigid-body rotational part of the flow alone and is not biased by the presence of pure or intrinsic shears. We show that the method performs well on a simplistic test case composed of two Lamb-Oseen vortices. We combine the proposed method with a state of the art clustering algorithm to build an automated vortex identification algorithm. (Abridged)


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 ⋆ An animation is available at https://www.aanda.org ⋆⋆ A Python implementation of the algorithm is publicly available. 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 λ 2 -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 1 (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 Q-criterion, the λ 2 -criterion, the swirling strengthfor 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 |Γ 1 |. 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. (2011Shelyag et al. ( , 2013, and Steiner & Rezaei (2012), to detect vortices and study their dynamics, while the enstrophy, defined as |ω| 2 , 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. (2020Silva 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. (2017Requerey et al. ( , 2018 employed Lagrangian tracers to visually identify vortices in super and meso-granular vertices, while Giagkiozis et al. (2018), andLiu 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 largescale 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.

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.

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, u, (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 of the rotation. For a rotational vortex, that is a flow rotating in a rigidbody 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.

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 = ∇u, which is the Jacobian matrix of the velocity field. If the flow is locally rotating, the velocity gradient tensor can be diagonalized as: where λ ± = λ cr ± iλ ci ∈ C and λ r ∈ 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: Moreover, the normalized eigenvector 2 associated with the real eigenvalue, u r , defines the rotation axis 3 . 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: where r 0 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.

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 = ∇u with one real and two complex conjugated eigenvalues, there exists a transposed quasi-triangular matrix V such that: where Q is an orthogonal matrix representing a rotation operator in three-dimensional space and V can be written as: The matrix V represents the velocity gradient tensor U in a rotated reference frame with the rotation axis parallel to the neŵ 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.
A118, page 3 of 16 A&A 668, A118 (2022) The various components of V represent rigid-body rotation (ϕ), stretching and compressing components (λ cr , λ r ), and shearing parts (ε, ξ, ν). 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: 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 , 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 , 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, 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.

Comparing criteria
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: 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: and compare it to the true period derived from the analytical formula Eq. (10), that is: 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.

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.

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 4 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 u = (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 0 = (x 0 , y 0 ), Rortex value R 0 0, and velocity u 0 = (v x,0 , v y,0 ). Moreover, we assume the velocity gradient tensor in x 0 , U| x 0 , to be known. Then, the EVC of the point x 0 can be computed based on two quantities: the curvature radius r, namely, the distance between the point x 0 and the center of curvature, and the radial direction, e r , that is the unit vector pointing from the point x 0 to the center of curvature.

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 0 and x 1 , 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 1 , u 1 , around x 0 as where ds = x 1 − x 0 and ∇u| x 0 = U| x 0 is the velocity gradient tensor computed in x 0 . If |ds| = |x 1 − x 0 | is sufficiently small, we can neglect the higher order terms. Moreover, the two points A118, page 5 of 16 A&A 668, A118 (2022) Fig. 3. Representation of a vortical velocity streamline. The velocity field at two points, x 0 and x 1 , is shown by green arrows and denoted as u 0 and u 1 , respectively. The two points are relatively close to each other and the vector joining them is labeled ds. The point x 0 is at distance r from the vortex center that is marked with a star. The purple unit vector, e r , represents the radial direction.
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 u 1 = 0. Moreover, if x 0 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: 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 0 . Therefore we also need to compute the radius of curvature to find the coordinates of the EVC.

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 R = 4π |R| . Therefore, the curvature radius, r, at the point x 0 can be estimated by combining Eq. (11) with the rotational velocity as follows: where v θ is the tangential component of the velocity field.
Assuming v θ = |u 0 |, where |u 0 | is the norm of the velocity field in x 0 , Eq. (15) can be expressed as and the curvature radius can be estimated based alone on the norm of the velocity field and the Rortex criterion in x 0 .
The assumption that v θ = |u 0 | 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 0 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 5 . 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.

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.

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 i j , 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: where χ is a kernel that depends on the Euclidean distance, d i j , 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: 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: 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 i j < d c , is around 1 to 2% of the total number of datapoints. Mehmood et al. (2016) proposed instead a nonparametric method 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 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 i j ).

Grid and vortex adaptation
The CFSFDP method is mainly limited by its computational cost, which scales as O(N 2 P ), 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 6 . 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 gridestimated 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: where ρ + i and ρ − 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 δ + i and δ − 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 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 A118, page 8 of 16 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 ρ ± i . The decision process is subject of the following subsection.

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: 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 δ = 1 2 ∆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 δ = 1 2 ∆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: where p γ > 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 p γ > 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 , 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.

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. 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.

Noise removal
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: 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 r r 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 A118, page 10 of 16 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 7 . In Appendix B, we review the structure of the code with a step-by-step flowchart.

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.

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 result is shown in the top panel of Fig. 10. For simplicity, 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 nonspiraling 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 7 https://github.com/jcanivete/swirl 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.

Orszag-Tang vortex test
The Orszag-Tang vortex problem (Orszag & Tang 1979) is one of the most widely known benchmark tests for MHD A118, page 11 of 16 A&A 668, A118 (2022) where B 0 = (4π) −1/2 , γ = 5/3, and p 0 = 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 A118, page 12 of 16 J. R. Canivete Cuissa and O. Steiner: An innovative and automated vortex identification method. I. the simulated flow, even in the presence of turbulence and magnetic fields. Applications to solar physical flows are planned for a follow-up paper.

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 A118, page 13 of 16 A&A 668, A118 (2022) map should be straightforward, as Eqs. (14) and (16) are valid in three-dimensions as well. In that case, we would obtain a threedimensional 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 thorough study on multi-dimensional clustering 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  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. 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. 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.