Vortex Flows in the Solar Atmoshpere: Automated Identification and Statistical Analysis

Aims. Due to the fundamental importance of vortices on the photosphere, in this work we aim to fully automate the process of intensity vortex identification to facilitate a more robust statistical analysis of their properties. Methods. Using on-disk observational data of the Fe I continuum, the process of vortex identification is fully automated, for the first time in solar physics, with the help of an established method from hydrodynamics initially employed for the study of eddies in turbulent flows (Graftieaux et al. 2001). Results. We find that the expected lifetime of intensity vortices is much shorter (~ 17s) compared with previously observed magnetic bright point swirls. Our findings suggest that at any time there are 1.4e6 such small-scale intensity vortices covering about 2.8% of the total surface of the solar photosphere.


Introduction
Traditionally photospheric intensity flow fields have been traced using local correlation tracking of (magnetic) bright points and the revealed vortex flows have been identified by eye. This manual approach has two major shortcomings, i) it introduces observational bias into the statistical analysis and ii) a large number of vortex flow fields are most likely missed simply due to the sheer scale of the task, which also has adverse effects on the variance of the statistical analysis. Small-scale vortices in the quiet Sun regions are widely accepted to form due to turbulent convection and the bathtub effect (e.g. Shelyag et al. 2011;Kitiashvili et al. 2012a;Shelyag et al. 2012). Solar photospheric vortex flows have drawn the attention of researchers as they have the potential to excite a wide range of MHD waves, e.g. slow and fast magneto-acoustic as well as Alfvén (Fedun et al. 2011;. Vortex flows also appear to have a prominent role in both direct and alternating current models of solar atmospheric heating. In direct current models, neighboring magnetic flux tubes (or strands) can become magnetically twisted under the influence of photospheric vortices. This, in turn, implies that current sheets may develop at the interface between such strands allowing the possibility of magnetic reconnection (Parker 1972(Parker , 1983bKlimchuk 2015). In alternating current models, photospheric vortices can be seen as MHD wave drivers (Fedun et al. 2011; and as precursors to large scale solar tornadoes (Wedemeyer-Böhm et al. 2012;Wedemeyer et al. 2013;Amari et al. 2015). These tornadoes have an estimated net positive Poynting flux of 440W m −2 that is more than adequate to heat the quiet solar atmosphere whose energy flux requirement is estimated to be approximately 300W m −2 (Withbroe & Noyes 1977). Unfortunately, despite the increasing interest in these coherent flows in the solar photosphere, the number of observations reported in the literature is still based on small sample sizes, with reports often associated with only a few detected events (see for example Steiner et al. 2010;Palacios et al. 2012;Park et al. 2016), or tens of observations (e.g. Bonet et al. 2008Bonet et al. , 2010Vargas Domínguez et al. 2011).
In this paper we present a fully automated method to identify vortex flows, namely the center of circulation and their flow boundary that is based on local correlation tracking (Fisher & Welsch 2008) applied to photospheric intensity observations, combined with an established method for identifying vortices used in the study of turbulence (Graftieaux et al. 2001). Subsequently, we estimate characteristic vortex parameters, such as lifetime diameter, mean perpendicular velocity and area. The main results of this paper are the following. There is an abundance of smallscale intensity vortices in the quiet Sun and their typical lifetimes are approximately 17 seconds. We estimate that at any given time, the expected number of vortices in the photosphere is 1.4 × 10 6 and that they occupy 2.8% of the photosphere. Although the area of these vortices may appear small in the photosphere, even if only a tenth of these vortex flows reaches the lower corona they may occupy more than 17% of its total area.

Observations
The observations investigated here were carried out between 08:07:24-09:05:46 UT on the 21st June 2012, with the CRisp Imaging SpectroPolarimeter (CRISP) at the Swedish 1-m Solar Telescope (SST: Scharmer et al. 2003Scharmer et al. , 2008 on La Palma. The image resolution of the CRISP observations is 0.059 ′′ per pixel. A quiet Sun region very close to disk center was observed with an effective Field-Of-View (FOV) of 55×55 arcsec, centered on solar-x =−3.1 ′′ and solar-y=69.9 ′′ . The required accuracy of the pointing of the CRISP FOV, in heliocentric coordinates, was achieved through co-alignment of bright-points observed with the CRISP wideband images, together with co-temporal continuum images in 170.0nm from the Solar Dynamics Observatory / Atmospheric Imaging Assembly (SDO/AIA: Lemen et al. 2012

Vortex Identification Process
The automated vortex identification methodology we present splits into four stages: i) pre-processing, ii) velocity field estimation, iii) vortex identification and, iv) vortex lifetime estimation. The intensity maps obtained from observations have varying intensity at different times that appears to be due to atmospheric effects. These intensity variations are a few standard deviations from the mean and the effect is global. To counter these effects image histogram equalization (e.g. Pizer et al. 1987) was used in the following way: -First, the expected distribution of intensities is estimated by means of averaging the histogram distributions across all frames. The rationale for this is that the Sun is not expected to change its general power emission spectrum during the time of the observation. -Once the expected intensity distribution has been obtained, histogram equalization is applied to all frames using that distribution as a reference.
This procedure is fast and efficiently removes inter-frame flickering, and, improves the numerical stability of the LCT method.
The pre-processing stage removes rapid intensity fluctuations, while preserving the relative counts as much as possible. Subsequently, we remove the mild seeing effects in our observations using a Gaussian filter. Although it is well known that atmospheric seeing is a nonlinear effect (see e.g. November & Simon 1988), given that the seeing conditions were good, this simple averaging method produced similar results to destreching algorithms and is computationally more efficient. To avoid the reduction of the temporal resolution a moving-average Gaussian filter is employed with a 3dB attenuation at quarter the Nyquist frequency 1/(8T ) where T is the cadence. The velocity field estimation is performed using fast local correlation tracking (FLCT) (Fisher & Welsch 2008) with a Gaussian apodizing window with a width of σ = 10 pixels (Louis et al. 2015).
Subsequently, for the vortex identification, we implement a proven and established method from the study of turbulence in fluid dynamics. Once the velocity field estimates are found we implement the same approach as Graftieaux et al. (2001) to identify the vortex centers and boundaries. Graftieaux et al. (2001) defined two functions Γ 1 and Γ 2 , for the identification of the vortex centers and boundaries, respectively. The function Γ 1 used in this work is, Here, S = {x m : ||x m − x p || 2 ≤ R} is a disk of radius R about the point x p , || · || 2 is the Euclidean norm, 1 z is a unit vector normal to the plane and |S| is the cardinality of S. Γ 1 defines a scalar field and its magnitude achieves a maximum at unity. Graftieaux et al. (2001) shows that this function achieves this maximum when x p is at the center of an axisymmetric vortex. However, given that ideal axisymmetric vortices are quite uncommon, the threshold for classifying a point in S as a potential vortex center is reduced to 0.9, and, the local maximum of these points is classified as the vortex center. For the identification of the vortex boundary, we use the discrete version of Γ 2 , defined as follows, wherev p is the mean velocity in the neighborhood of the point x p . It is shown in Graftieaux et al. (2001) that in the inner core of a vortex the magnitude of Γ 2 is larger than 2/π. Flows with values of Γ 2 < 2/π are dominated by strain and when Γ 2 = 2/π we have a pure shear.
Let us now calculate the vortex centers and their boundaries at every time instance, however, we still need to estimate the lifespan of these vortices. For this purpose we assume that the vortex center can move at approximately the sound speed of the photosphere, about 10km s −1 (Nordlund et al. 2009). If the speed of the vortex centers is comparable to the sound speed, this would suggest that the maximum distance a center could traverse from one frame to the next would be 82.5km, which is almost 2 pixels at the spatial resolution of our data. However, at present the vortex formation mechanism has not been clearly established and if such flows in the photosphere are formed as shown in Figure 1, the speed of their center may be much larger than the sound speed. What we suggest in Figure 1 is the following, the edges of the granules are represented as line segments (red and blue line segments in Figure 1). We define the points where the vertical component of the velocity transits from being mostly positive, as is on granules, to being negative, as is the case in the inter-granular lanes. Due to the dynamic nature of the granulation pattern on the photosphere, their edges are in constant relative motion with respect to the edges of neighboring granules. This relative motion, when combined with counter streaming flows of two neighboring granules, can drive vortex flows whose centers can move (v ) at much larger speed compared with the relative speed that generated them (see Figure 1). Therefore, using a conservative estimate we assume that vortices that are within a 4 pixel radius in two consecutive times, are in fact the same vortex. 1. Cartoon of the proposed physical mechanism modelling the high velocity of vortex centers. The line segments yL and yR, shown in blue and red color, respectively, represent the edges of two neighboring granules. In this instance, the two edges are moving towards each other with speed |v|. The streamlines in the plane outline the velocity field near the edges of the granules, with vL and vR denoting the velocity field in the left and right granule, respectively. The velocity of the vortex center is labeled v . The blue streamlines in the z-direction depict magnetic field lines above the vortex center. The black arrow shows the evolution.

Results and Statistical Analysis
A representative example of the results obtained from the vortex identification process is shown in Figure 2. The grayscale denotes intensity, normalized in the range 0 to 1 corresponding to black and white, respectively. Overplotted is the LCT estimate of the surface velocity field. The red-filled circles mark counter-clockwise vortex flows (positive), blue circles correspond to clockwise flows (negative) and the orange line delimits the vortex flow boundary. Figure 3 and Figure 4 show statistical results based on a sample size of N = 26, 988 vortices. As the only source of information that is used here is based on LCT applied to intensity observations, we refer to the identified vortex flows as intensity vortices. This is to acknowledge line-ofsight integration effects and temperature variations that, from a practical standpoint, lead to the estimated velocity field being a weighted average of the plasma motions at different heights within the spectral line formation height (Nordlund et al. 2009).
We find that the expected lifetime of such vortices is independent of orientation (see a) and b) in Figure 3). In fact, we have found no statistically significant deviations in the distributions of positively or negatively oriented vortices for any of the measured parameters, i.e. lifetime, space and time density, diameter, area or perpendicular speed (see Figure 3 and 4). What is intriguing however, is that for the majority of vortices (approximately 85%) their expected lifetimes are less or equal to three times the cadence (24.75s). This is much shorter when compared with similar features identified by tracking magnetic bright points (BPs) (e.g. Bonet et al. 2008). The apparent discrepancy could be attributed to errors in LCT, where very short lived structures are the result of errors in the identified velocity field. Notwithstanding this limitation, LCT velocity maps have been shown to be a good first order approximation to the velocity field (Verma et al. 2013;Louis et al. 2015).  Figure 3), we estimate that there will be τ · d = 0.244 Mm −2 vortices at any time. In turn, this implies that there are continuously 1.48 × 10 6 vortices over the whole photosphere. Under the assumption that the expected lifetime, as well as, space and time density in other regions of the photosphere are similar to our observations. If intensity vortices are indeed closely correlated with the actual velocity field, then based on the simulation results reported by Kitiashvili et al. (2012a,b) we anticipate that their expected size will decrease with the advent of higher spatial resolution observations. A prime example of near future expected capability is the Daniel K. Inouye Solar Telescope (DKIST) whose visible broadband imager is planned to have spatial resolution of 16km to 25km per pixel, at 430.4nm and 656.3nm, respectively, and cadence of 3.2s (Berger & ATST Science Team 2013). Figure 4, panels e) and f ) show the distribution of the average perpendicular speed within the vortex boundary. This is calculated by projecting the velocity vector at every point within the vortex to a vector perpendicular to the ray emanating from the vortex center. Lastly, panel g) in Figure 4 provides an estimate of the percent of the area of the photosphere covered by intensity vortices at any time.

Discussion and Conclusions
Vortex flows in the solar atmosphere may contribute significantly to the energy flux requirements for heating the quiet Sun atmosphere. However, for that connection to be established strong evidence is required: i) vortex flows motions are ubiquitous in the solar atmosphere, ii) that these motions appear at different heights, e.g. photosphere, chromosphere and corona. We have shown, that the automated identification approach described in this work results in a significantly larger number of identified vortices compared with previous observational studies. This is evidence consolidating the fact that small-scale vortices are prevalent in the solar photosphere. Most interestingly, an overwhelming majority of these vortices have lifetimes that are often much shorter than previously believed, which suggests that these flows are highly dynamic in nature.
Due to the episodic nature of the formation of these small-scale vortices, any magnetic field through them will be supplied with a broadband impulse comprised of both torsional and radial components which will generate propagating MHD waves. The presence of a magnetic field in vortices is consistent if we recall that their location is in the inter-granular lanes where the magnetic field concentrations are highest. Both observational and numerical simulations (e.g. Fedun et al. 2011; support the idea that MHD waves with a broad frequency range can be generated by vortex flows. However this is to be expected on more fundamental grounds due to a particular duality in frequency space. Namely, localization in time, leads to spread (broadening) in the frequency domain and vice versa. In physical terms this implies that the rapidity of vortex formation alongside with deviation from axi-symmetry offer a wave driver that results in waves of different frequencies, albeit with different amplitudes. Regarding energy transport to the upper layers of the atmosphere, numerical simulations suggest that vortex driven MHD waves (Amari et al. 2015) are a feasible mechanism.
The most compelling differences compared with previous reports (e.g. Bonet et al. 2008Bonet et al. , 2010;  #)). In the last column we also provide the number of vortices on which the statistical analysis is based on. The values in the last row in the table are in parentheses since they do not correspond to observations in the photosphere, however, these are included for reference.
Vargas Domínguez et al. 2011) are in the expected lifetime and space and time density. Table 1 shows summary statistics comparing the main results in this work with previous studies that are based on more than 3-4 observed vortices. In our view, there are at least two explanations for this mismatch. First, vortex flows and any type of feature tracking in observations, is time consuming and error prone when performed manually. This increases the likelihood of bias and increased variance. Also, our estimate of lifetimes relies on the accuracy of LCT for the surface velocity field identification, which although has been shown to have reasonable correlation with the true velocity field (Louis et al. 2015), is only a first approximation to small-scale motions in inter-granular lanes. Notwithstanding this uncertainty, given that our chosen automated technique is straightforward to implement, the results can be cross-validated by other studies, which, in our view, is extremely important. The similarities in scale and location, of vortex flows identified in this work with small-scale whirlpools (or swirls) reported by Bonet et al. (2008) and Vargas Domínguez et al. (2011, lead us to conjecture that these are, in fact, the same flow features in the quiet Sun. If that is indeed the case, then remarkably the number of vortices in the photosphere would be an order of magnitude larger than previous estimates. The results in this work could suggest that previously identified magnetic tornadoes by Wedemeyer-Böhm et al. (2012), may also be more numerous by an order of magnitude. If intensity vortices prove to be the root cause of solar tornadoes (Amari et al. 2015), then this would suggest that 10% of the photospheric vortices reach the lower corona forming magnetic tornadoes. This has the extraordinary implication that at least 17% of the area of the lower corona is constantly supplied with a positive Poynting flux of 440W m −2 , as opposed to 1.2% implied from Wedemeyer-Böhm et al. (2012). The assumptions in this estimate are that the photospheric vortices that do extend up to the lower corona have a mean radius of 1.5Mm, their corrected expected number is 1.48 × 10 5 vortices at every time, i.e. 10% of our intensity vortices on the photosphere, instead of 1.1 × 10 4 (see Table 1) and that the average net positive Poynting flux in each magnetic tornado is 440W m −2 . These implications may be exciting, however, at present we are far from establishing a link between intensity vortices and magnetic tornadoes. Nevertheless, this would be an interesting direction of future enquiry as this could be key in resolving the quiet Sun heating problem. Fig. 3. Estimates of (a and b) vortex lifetime mass function, and (c,d and e) the number of vortices per M m 2 · minute. The red circles denote the best fit of a parametric mass density function (PMF). In this case, the Geometric distribution was a best fit for the lifetimes of the vortices. The orange line, and the white font E on its right, is the expected value calculated from the empirical distribution of the data. Values with a hat indicate best fit parameter estimates for the particular distribution, and, E(·) is the expected value.