Observations

The SOT data employed consist of blue continuum (FG-blue) images sampled with the wide-band imager with a central wavelength of 4504.5 Å and a band width of 4 Å, and chromospheric images at the Ca II H line core with a central wavelength of 3968.5 Å. The FG-blue data were collected between 05:48:03 UT and 08:29:59 UT on the 5 March 2007, and the Ca II H data between 05:48:06 UT and 07:09:28. Images were targeted at a quiet region centred at x c = 5.3″, y c = 4.1″. The cadence is ~6.42 s. Each of the images has a FOV of ~56″ × 28″ (40.1 Mm × 20.1 Mm), with a image size of 1024 × 512 px2 and a pixel size of 0.0545″ (39.2 km). The level-1 fits files were processed with the fg_prep.pro program available in the SolarSoft IDL packages.

Swirl Detection Using ASDA

For every pixel P in a preprocessed, scientifically ready image, two dimensionless parameters24 are defined as:

$$\begin{array}{*{20}{l}} {\Gamma _1(P)} \hfill & = \hfill & {\widehat {\mathbf{z}} \cdot \frac{1}{N}\mathop {\sum}\limits_S {\frac{{{\mathbf{n}}_{PM} \times {\mathbf{v}}_M}}{{|{\mathbf{v}}_M|}}} ,} \hfill \\ {\Gamma _2(P)} \hfill & = \hfill & {\widehat {\mathbf{z}} \cdot \frac{1}{N}\mathop {\sum}\limits_S {\frac{{{\mathbf{n}}_{PM} \times ({\mathbf{v}}_M - \overline {\mathbf{v}} )}}{{|{\mathbf{v}}_M - \overline {\mathbf{v}} |}}} .} \hfill \end{array}$$ (1)

Here, S is a two-dimensional region containing the target point P, with a size of N pixels. M is a point within S. n PM denotes the normal vector pointing from point P to M. \(\overline {\mathbf{v}}\) is the average velocity vector within the region S and v M is the velocity vector at point M. Symbols | | and × are for the mode of vectors and the cross product, respectively. The velocity field is estimated using the Fourier Local Correlation Tracking (FLCT) method25,26 on successive images from the observations. \(\overline {\mathbf{z}}\) denotes the normal vector perpendicular to the observation surface pointing towards the observer.

It has been demonstrated that, |Γ 1 | peaks at the centre and |Γ 2 | is larger than 2/π within the edge of a swirl24. To find all swirls in a given frame of observations, the levels of ±2/π of Γ 2 are contoured to find out all candidates of swirls; and then candidates with peak |Γ 1 | values greater than a given threshold (0.89)14,24 are confirmed as swirls. The given threshold removes all candidates with expanding/shrinking speeds larger than half of their rotating speeds14. Detailed tests on a series of synthetic data and realistic numerical simulation data with varying spatial resolutions, have been performed and proved ASDA as a proficient and astute method in accurately detecting solar atmospheric swirls14.

Parameters of swirls

Given the velocity field and edge of a swirl, its effective radius (R) and average rotating speed (v r ) are determined by:

$$\begin{array}{*{20}{l}} R \hfill & = \hfill & {\sqrt {\frac{A}{\pi }} ,} \hfill \\ {v_r} \hfill & = \hfill & {\frac{1}{k}\mathop {\sum}\limits_{i = 1}^k {{\mathbf{v}}_i \cdot {\mathbf{n}}_i} .} \hfill \end{array}$$ (2)

Here, A is the area of the swirl. k is the total number of points at the edge of the swirl. v i and n i are the velocity and the normal vector perpendicular to the local radial direction (from the centre of the swirl) of the ith point at the edge, respectively.

Lifetimes of swirls at a single height of osbervations are estimated following the method proposed in the ASDA paper14. Suppose, there are two swirls (S 1 and S 2 ) detected in two successive frames. S 1 is detected at time t 0 and S 2 at time t 0 + Δt, where Δt is the cadence of the observation. S 1 and S 2 will be considered as the same swirl, if:

$${\mathbf{c}}_{\mathbf{1}} + {\mathbf{v}}_{{\mathbf{c1}}} \cdot {\mathrm{\Delta }}t \subset S_2.$$ (3)

Here, c 1 is the location of the centre of swirl S 1 . v c1 is the speed of the centre of S 1 . The symbol ⊂ means belonging to. Taking into account that a swirl may experience changes to its rotational motion through time, we then allow swirls to be missing from one frame when evaluating their lifetimes, i.e., S 1 and S 3 are still considered as the same swirl if c 1 + v c1 ⋅ 2Δt⊂S 3 is true where S 3 is a swirl detected at time t 0 + 2Δt.

Estimation of correlation indices and overlaps

Many difficulties exist in directly comparing swirls at two different layers L 1 and L 2 (for example L 1 in the photosphere and L 2 in the chromosphere) to find their overlaps: first of all, swirls confirmed by ASDA are only part of all the candidates. A photospheric candidate might be confirmed as a swirl, but its chromospheric correspondence (if there is any) could be not due to reasons such as too large expanding speed. This case would not be rare because of the common flux tube expansion from the photosphere to the chromosphere. Secondly, magnetic field in the low atmosphere is more or less inclined27, it is unlikely that a photospheric swirl and its chromospheric correspondence stay at exactly the same horizontal location. And, thirdly the shapes of swirls are found to be irregular, meaning that a photospheric swirl and its chromospheric correspondence are unlikely to be 100% overlapped, even if they stay at exactly the same horizontal location.

To bypass the above difficulties, we have developed the following method to estimate the correlation and overlaps of two given layers (L 1 and L 2 ):

Firstly, binarize the Γ 2 maps of L 1 and L 2 , set all points with values larger than 2/π as 1 and less than −2/π as −1. All other points are set to 0. The resulted binarized maps are \(\Gamma _2^1\) and \(\Gamma _2^2\) for layer L 1 and L 2 , respectively.

Secondly, sum the absolute value of all points in \(\Gamma _2^1\) as T 1 , \(\Gamma _2^2\) as T 2 respectively. Define T as the smaller one between T 1 and T 2 .

Thirdly, multiply \(\Gamma _2^1\) and \(\Gamma _2^2\), and obtain the correlation map C. Then, we define (∑C)/T as the correlation index (CI) between layers L 1 and L 2 . The above procedures determine that CI can only range between −1 to 1, with a higher CI value implying a higher correlation between layers L 1 and L 2 . However, the CI of two given layers cannot give us information of how many swirls are overlapped between layers L 1 and L 2 . For example, a CI of 0.01 could either mean only 1% of the swirls in one layer are overlapped with those in the other layer, or all swirls are overlapped but each swirl has only 1% points overlapped, in extreme cases.

And finally, for a given swirl S detected in layer L 1 , we map all of its points in to the C map obtained in the previous step and calculate the percentage of points which correspond to positive C values as its own correlation index (CI S ). Swirl S is then labelled to have its correspondence (to be overlapped) in layer L 2 if CI > 0 and CI S > t h , where CI is the correlation index between layers L 1 and L 2 . t h is a positive value defined by \(\Delta ^2/\overline A\), where Δ is the pixel size (39.2 km for the utilized SOT observations) and \(\overline A\) is the average area of all swirls (\(= \pi \overline R ^2\), \(\overline R\) is the average effective radius ~290 km), meaning that on average there is at least one point within swirl S corresponding to a positive correlation.

However, we are aware of that the above processes cannot fully overcome all difficulties raised above, especially the influence of the inclined magnetic fields. We suggest that the number of photospheric swirls which were found to have their correspondences in the chromosphere using the above method should have been under-estimated.

Numerical simulations

The 3D MHD simulations have been performed using the Sheffield Advanced Code (SAC)19, which solves the ideal MHD equations with the presence of arbitrary perturbations in a gravitationally stratified and magnetised atmosphere. SAC separates background and perturbation variables in order to accurately resolve perturbations in a stratified atmosphere. This approach has been developed with the capability to perform simulations even in the non-linear regime. SAC has been tested and used to study wave propagations along a flux tube previously10. The governing equations of SAC are briefly recapped in the Supplementary Note 2.

The initial density and temperature profiles of the simulation performed in this work have been constructed using the VAL IIIC model28. An axisymmetric and self-similar expanding flux tube, with a magnetic field strength of 800 G at its footpoint constructed following previous literatures10,29, is then embedded into the ambient atmosphere. The analytic equations and parameters used for the construction of the magnetic flux tube are available in the Supplementary Note 2. The computational domain of the particular simulation is −1.0 ≤ x ≤ 1.0, −1.0 ≤ y ≤ 1.0 and 0.4 ≤ z ≤ 1.6 Mm, simulating from the upper photosphere to the upper chromosphere. The above domain is resolved with a grid size of (129, 129, 259) points in the x−, y− and z-direction, respectively. Thus the spatial resolutions in the x−, y− and z-direction are 15.5, 15.5 and 4.6 km, respectively.

The rotational driver is introduced at the bottom of the flux tube at the simulation time t = 0, generated from the following formula10:

$$v_r = v_0 \cdot {\mathrm{exp}}\left( { - \frac{{r^2}}{{\delta r^2}}} \right) \cdot {\mathrm{exp}}\left[ { - \frac{{(z - z_0)^2}}{{\delta z^2}}} \right] \cdot {\mathrm{sin}}\left( {\frac{{\pi t}}{P}} \right).$$ (4)