© The Authors 2024

1. Introduction

Observations carried out over the past two decades indicate that small-scale vortical motions are ubiquitous in the quiet solar atmosphere. Many of the vortex detections have been obtained by individually following the trajectories of bright points (BPs) and small-scale magnetic structures (Bonet et al. 2008; Balmaceda et al. 2010; Manso Sainz et al. 2011) or by visual tracking of swirling photospheric and chromospheric features, such as rings, filaments, and arcs (Wedemeyer-Böhm & Rouppe van der Voort 2009; Wedemeyer-Böhm et al. 2012; Park et al. 2016; Tziotziou et al. 2018, 2019; Shetye et al. 2019). Another approach is to use local correlation tracking (LCT) techniques to identify spiraling motions in the morphology of the estimated horizontal velocity fields (Bonet et al. 2010; Vargas Domínguez et al. 2011; Requerey et al. 2017, 2018). A review on vortical motions in the solar atmosphere is presented in Tziotziou et al. (2023).

These methods yielded precious results on the characteristic sizes, lifetimes, and rotational periods of photospheric and chromospheric swirls. However, they are possibly biased by the human conception of the definition of a vortex, because the detection processes rely on the visual identification of swirling motions from images, time sequences, and velocity field maps derived from observations. A more pragmatic approach consists in using mathematical criteria and geometrical methods to limit the effects of human subjectivity in the identification process. A number of such vortex identification methods can be found in the literature (see, e.g., Günther & Theisel 2018, for a review).

Many of these methods have been employed to study small-scale vortical motions in the solar atmosphere, especially in the context of numerical simulations where the necessary physical quantities are directly accessible. For example, studies using the Γ functions (Graftieaux et al. 2001) to study photospheric vortices in observations yielded average diameters of ∼0.5 Mm and lifetimes of ∼0.3 min (Giagkiozis et al. 2018; Liu et al. 2019a). In the chromosphere, Dakanalis et al. (2022) reported swirls with a mean diameter of 2.6 ± 0.6 Mm and a lifetime of about 10.3 min using a new morphological approach (Dakanalis et al. 2021). From simulations, chromospheric swirls with typical diameters of 0.7 ± 0.3 Mm and average lifetimes of ∼1.0 min have been found by Kato & Wedemeyer (2017) using the swirling strength criterion (Zhou et al. 1999). We further refer the reader to Stein & Nordlund (1998), Moll et al. (2011, 2012), Shelyag et al. (2011, 2013), Kitiashvili et al. (2012), Steiner & Rezaei (2012), Silva et al. (2018, 2020, 2021), Liu et al. (2019b), Canivete Cuissa & Steiner (2020), Yadav et al. (2020), Battaglia et al. (2021), and Aljohani et al. (2022) for a nonexhaustive list of studies employing mathematical and geometrical methods to analyze swirling motions in observations and numerical simulations of the solar atmosphere.

However, a universally accepted and rigorous method for vortex identification has not yet been found. Indeed, all the proposed methods present shortcomings when applied to the magnetized, turbulent, and highly dynamical flows of the solar atmosphere (see, e.g., Canivete Cuissa & Steiner 2022, for a discussion). In Canivete Cuissa & Steiner (2022, hereafter Paper I), we presented a new method for the automated identification of vortices, called the estimated vortex center (EVC) method. It combines the accuracy and quantitative aspects of mathematical criteria with the global and morphological perspective of the curvature center method proposed by Sadarjoen & Post (1999). We implemented the method in a Python package called SWirl Identification by Rotation-centers Localization (SWIRL; Canivete Cuissa 2022), which is open source on GitHub¹.

The SWIRL algorithm was tested on an artificial velocity field composed of nine Lamb-Oseen vortex models with a random Gaussian noise and on the turbulent flow resulting from a magneto-hydrodynamical (MHD) Orszag-Tang vortex system. In particular, the MHD Orszag-Tang test yields a flow with a diverse spectrum of MHD modes, shocks, and turbulence (see, e.g., Londrillo & Del Zanna 2000). Consequently, accurate vortex identification in such a complex flow poses a significant challenge to any dedicated algorithm. The results showed the reliability and robustness of the algorithm in the presence of noise, turbulence, and magnetic fields. Moreover, as the EVC method does not require the use of a threshold, vortices with rotational velocities that are comparable to the noisy background velocity field are not precluded from being identified. Therefore, the SWIRL algorithm proved to be suitable for identifying vortices in astrophysical velocity fields.

There are many open questions regarding small-scale swirls in the solar atmosphere. For example, the typical size, number density, strength, and lifetime of these events have been the subject of multiple observational and numerical studies. However, the obtained results most certainly depend on the spatial resolution of the instrumentation or simulation (Yadav et al. 2020), and on the methods employed (see, e.g., Silva et al. 2018; Tremblay et al. 2018). It is also not yet clear whether coherent vortical structures observed in the lower solar atmosphere can extend into the corona (Breu et al. 2022).

Moreover, being tightly coupled to the small-scale magnetic field of the Sun, small-scale vortical motions could be associated with torsional Alfvén waves (Wedemeyer-Böhm et al. 2012; Shelyag et al. 2013). Signatures of torsional Alfvén waves in the solar atmosphere have been found by, for example, Jess et al. (2009), Okamoto & De Pontieu (2011), De Pontieu et al. (2012), and Srivastava et al. (2017), while Liu et al. (2019a) and Battaglia et al. (2021) reported upwardly propagating torsional Alfvénic pulses related to chromospheric swirls in observations and numerical simulations, respectively. These studies, among others, indicate that the energy flux associated with vortical events can sustain the radiative losses in the chromosphere, and therefore these events can contribute to local heating.

The role that small-scale swirls may play in the dynamics of the solar atmosphere calls for a rigorous method for their identification. In particular, a robust statistical analysis of their properties is required to asses their real impact on chromospheric and coronal heating. In this paper, we demonstrate that the method presented in Paper I also reliably identifies swirls in the turbulent and highly dynamical flow of a three-dimensional, MHD numerical simulation of the solar atmosphere. Moreover, we carry out a statistical analysis on the properties and characteristics of small-scale swirling motions in those simulations.

The paper is organized as follows. In Sect. 2, we briefly describe the numerical simulations and the vortex identification method used in this work. In Sect. 3, we present and discuss the performance of the method when applied to numerical simulations of the solar atmosphere and a statistical analysis of the identified swirls. Finally, we summarize our findings and present our conclusions in Sect. 4.

2. Methods

2.1. Numerical simulations

We employed realistic numerical simulations of the solar atmosphere obtained with the radiative MHD code CO5BOLD (Freytag et al. 2012). The size of the Cartesian simulation box is 9.6  ×  9.6  ×  2.8 Mm³ and the cell size is 10 km in each spatial direction. The number of grid cells is therefore 960 × 960 × 280. The average optical surface τ₅₀₀ = 1, which we label as z = 0 km, is found at ∼1300 km from the bottom of the box. Therefore, the simulation domain represents a small volume near the solar surface, which includes the surface layers of the convection zone, the photosphere, and up to the middle chromosphere. The average stratifications of density, temperature, and the root mean square (rms) of the vertical component of the velocity field are shown in Fig. 1.

Fig. 1. Average stratifications of density, ρ, temperature, T, and rms of the vertical component of the velocity field, vz,  rms. The profiles represent averages, both temporally across the 30 time instances of the CO5BOLD simulation, and spatially across the horizontal sections of the domain. The heights for the analysis of swirls in the surface layers of the convection zone (z = −200 km), in the photosphere (z = 100 km), and in the low chromosphere (z = 700 km) are indicated by blue dashed lines.

The simulation started from a relaxed, purely hydrodynamical model, to which a unipolar, vertical magnetic field of 50 G was added. The lateral boundary conditions are periodic for both the plasma and the magnetic field, while at the top and bottom of the box the magnetic field is forced to be vertically oriented. More details on the simulation setup can be found in Calvo (2018, Sect. 2) and in Battaglia et al. (2021).

This choice of initial magnetic field configuration is roughly representative of a predominantly unipolar magnetic network patch of a quiet Sun region. The configuration and top boundary condition of the magnetic field favor the production of vertically oriented vortex tubes in the chromosphere, as was demonstrated by Battaglia et al. (2021, Appendix A). A stronger initial field would yield a more homogeneous structure of vortices, while vanishing magnetic fields would lead to less and rather isotropically distributed vortices.

For this study, we analyzed 30 time instances of the CO5BOLD simulation with a cadence of 4 min, which cover a total of 2 h in physical time. This cadence period is two-thirds of the mean granular lifetime (Hirzberger et al. 1999).

2.2. Identification algorithm

In this paper, we employed the EVC method presented in Paper I. It can be considered an extension of the curvature center method proposed by Sadarjoen & Post (1999), where the velocity field and its derivatives are used instead of streamlines. In more detail, the method consists in accurately estimating the center of rotation of every rotating fluid particle (grid cell) from the instantaneous horizontal velocity field alone. Fluid particles belonging to the same vortex share a common axis of rotation (Lugt 1979), and therefore their estimated centers of rotation, dubbed EVCs, should cluster around the true core of the vortical structure. Consequently, vortices are identified through clusters of EVCs.

To accurately compute the EVC of any given grid cell that presents some degree of curvature in the velocity field, one has to estimate the radius of curvature and the radial direction of the local flow. For this purpose, we employ the Rortex criterion, R, proposed by Tian et al. (2018) and Wang et al. (2019). The Rortex criterion is a mathematical criterion, like the vorticity, and it is defined as

$$\begin{array}{c}\hfill R=\mathit{\omega }·{\mathit{u}}_{\mathrm{r}}-\sqrt{{(\mathit{\omega }·{\mathit{u}}_{\mathrm{r}})}^2-{\lambda }^2},\end{array}$$(1)

where ω is the vorticity vector, u_r is the normalized, real eigenvector of the velocity gradient tensor, and λ is the swirling strength criterion. For more details on these quantities, we refer the reader to Paper I. However, whereas the vorticity and other mathematical criteria are affected by the presence of shear flows, the Rortex criterion measures the rigid-body rotational part of the flow alone. Therefore, it is the optimal quantity to extract physical information on the curvature of the flow from the velocity field and allows unprecedented accuracy in the estimation of the center of rotation (Paper I).

Given a map containing all the computed EVCs, clusters indicating the presence of vortices can, in principle, be identified by eye. Nevertheless, in Paper I we proposed a modified version of the clustering by fast search and find of density peaks (CFSFDP) algorithm (Rodriguez & Laio 2014) to automatize the identification process. Moreover, a cleaning procedure is proposed to remove misidentifications caused by noise or by coherently nonspiraling curvatures in the flow.

The EVC method and the associated automated algorithm are implemented in an open-source Python package called SWIRL. For more details on the method, the clustering algorithm, and the test cases, we refer the reader to Paper I.

3. Results and discussion

In this section, we test the applicability of the SWIRL algorithm in automatically identifying swirls in photospheric and chromospheric horizontal sections of the simulation introduced in Sect. 2.1. We then present the results of a statistical study performed over the full set of data cubes that addresses the properties of small-scale swirls and their relation with the surface magnetic field of the Sun in numerical simulations.

The SWIRL algorithm requires careful tuning of several parameters based on the specific characteristics of the flows being analyzed. Detailed descriptions of these parameters can be found in Paper I and in the GitHub repository for the SWIRL code. The values used in this study are listed in Table 1.

We find the identification process to be particularly sensitive to the number of “stencils” used, as well as the “noise” and “kink” parameters. Increasing the number of stencils increases the robustness of the identification process to small-scale turbulence. However, using too many stencils (typically more than ∼10) can lead to lower computational performance without significantly improving the results.

The noise and kink parameters are responsible for cleaning up false detections, and their adjustment depends on the level of noise and turbulence in the flow. Higher values can lead to false detections, while excessively low values can cause true vortices to be missed. Empirically, the parameter values that have shown good performance in CO5BOLD simulations of the solar atmosphere range from about 0.5 to 1.5. However, we encourage users of the SWIRL algorithm to experiment with different values.

The clustering parameters can also be adjusted based on decision graphs (see Fig. 6 in Paper I). However, the values given in Table 1 should generally lead to satisfactory results for most applications.

3.1. Validation of the SWIRL algorithm on CO5BOLD simulations

3.1.1. Photosphere

We started from a photospheric, two-dimensional, horizontal subsection of a time instance of the simulation data, which is shown in Fig. 2. The chosen height is z = 100 km as shown in Fig. 1. We notice the granular pattern of the flow with integranular lanes harboring magnetic flux concentrations (top panel). The magnetic field is predominantly of positive polarity because of the initial condition of the simulation.

Fig. 2. Two-dimensional horizontal subsection of the simulated photosphere. The section measures 4.0 × 4.0 Mm2 and is taken at z = 100 km. The simulation time instance corresponds to t = 5774 s. The horizontal velocity field of the subsection is depicted using a vector plot. The length of each arrow corresponds to the magnitude of the horizontal flow, and a reference scale is included in the bottom-right corner. Top: Vertical magnetic field Bz at z = 100 km. Middle: Rortex criterion R. Bottom: G-EVC map. Contours where R ≠ 0 are shown in gray in the middle panel.

The middle panel of Fig. 2 shows the Rortex criterion, R, computed from the horizontal velocity field. Positive values of R (green) indicate counterclockwise curvature in the flow, while clockwise curvatures are characterized by negative R (purple). Small-scale patches where R ≠ 0 appear to roughly track the photospheric magnetic flux concentrations, but form a rather chaotic pattern of different rotational strengths and orientations².

It would be difficult (if not impossible) to discern coherent vortical structures using the Rortex map alone. A priori, we do not know if a two-dimensional region of R ≠ 0 is part of a coherent vortical structure or simply stems from the turbulent nature of the flow. Indeed, the Rortex criterion is a local criterion, defined on a small stencil of only very few grid cells. However, to distinguish turbulent, local rotations from actual vortical flows, additional information about the large-scale properties of the flow is needed. For example, if a single fluid parcel is deflected, it may exhibit local rotation and therefore be identified by a mathematical criterion such as Rortex. However, the presence of a vortex flow requires several fluid parcels to rotate coherently about a common axis. It is the coherent behavior of many parcels that indicates the presence of a vortical structure.

This conundrum can be partially solved by considering the G-EVC map, which is shown in the bottom panel of Fig. 2. The G-EVC map is obtained by counting the number of EVCs in every grid cell³. Clockwise and counterclockwise EVCs count as −1 and +1, respectively, and their sum determines the grid cardinality, s, in each grid cell. In principle, a cluster of EVCs indicates the location of a vortex core, and therefore high absolute values of the grid cardinality, |s|, can be used to infer the presence of a vortex. For example, by inspecting the bottom panel of Fig. 2, we expect a counterclockwise vortex to be found around (x,  y) = (3.1 Mm,  3.25 Mm) and a clockwise one close to (x,  y) = (0.7 Mm,  2.4 Mm). On the other hand, we can presume that the Rortex patches around (x,  y) = (0.2 Mm,  3.7 Mm) visible in the middle panel of Fig. 2 do not represent a swirl, because the grid cardinality is relatively low in that region.

The SWIRL algorithm automatically finds clusters of G-EVCs, and thus detects candidate vortex centers. The vortices identified on the 4.0 × 4.0 Mm² photospheric subsection of the CO5BOLD simulation are shown in Fig. 3. For simplicity, we represent vortices with colored disks centered on the estimated vortex core. However, it is important to note that the SWIRL algorithm returns a collection of grid cells that form the vortex for each identification. As a result, the true shape of the identified vortices, while generally exhibiting a roundish appearance, tends to be more irregular than what is presented here. The radius of the disk corresponds to the effective radius of the vortex, r_eff, which is computed as

$$\begin{array}{c}\hfill r_{\mathrm{eff}}=\sqrt{\frac{N_{\mathrm{c}}}{\pi }}\mathrm{\Delta }x,\end{array}$$(2)

Fig. 3. Vortices identified by the SWIRL algorithm in the two-dimensional, horizontal velocity field of Fig. 2. The location and effective size of the identified vortices are indicated by colored disks. Clockwise vortices are represented by purple disks, while counterclockwise ones are shown in green. The vertical magnetic field Bz is color coded and saturates at ±1000 G. The gray squares denote the 0.6 × 0.6 Mm2 regions shown in Fig. 4.

where N_c is the number of EVCs belonging to the cluster and Δ is the grid spacing. Here, the effective radius is defined through the effective area occupied by the grid cells that form that vortex. The color of the disks indicates the rotation direction: green for counterclock wise vortices and purple for clockwise ones. In total, 21 vortices have been identified by the code with an average effective radius of ∼50 km.

Most of the detected vortices lie within or nearby strong magnetic flux concentrations. This is expected, because photospheric swirling motions are known to be tightly coupled to small-scale surface magnetic fields (see, e.g., Moll et al. 2012; Battaglia et al. 2021). There are also a few exceptions: for example, the two clockwise vortices around (x,  y) = (4.0 Mm,  3.7 Mm). These apparently nonmagnetic events could be related to the footpoints of vortex arches in high-plasma-β regions or to nonmagnetic bright points, such as those reported in numerical simulations by Muthsam et al. (2010), Moll et al. (2011), Battaglia et al. (2021), and Calvo et al. (2016).

Figure 4 shows zoom-in plots of six different 0.6 × 0.6 Mm² regions of the photospheric section shown in Fig. 3. The horizontal velocity field, which is represented by instantaneous streamlines, is particularly turbulent in magnetic flux concentrations, resulting in multiple spiraling configurations within the same magnetic structure. In general, the identified vortices correlate well with the spiraling instantaneous streamlines.

Fig. 4. Zoom-in plots of the photospheric regions outlined in Fig. 3. The location and effective size of the identified vortices are indicated by colored disks. The vertical magnetic field Bz is color coded with the same scale as in Fig. 3, while the horizontal velocity field is represented by instantaneous streamlines.

In panel A, no vortices have been identified despite the large negative value of the grid cardinality, s, in that same region (see bottom panel of Fig. 2). That cluster of EVCs is caused by the semi-circular clockwise configuration of the flow visible at coordinates (x,  y) = (0.3 Mm,  0.3 Mm) in the center of panel A of Fig. 4. In this case, the SWIRL algorithm identifies the cluster of (G-)EVCs as a possible candidate vortex, but correctly discards it during the cleaning procedure because the flow is not fully spiraling. For details on the cleaning procedure, we refer the reader to Paper I.

The only two misidentifications are found in panel C at coordinates (x,  y)  =  (0.1 Mm,  0.25 Mm) and (x,  y)  =  (0.25 Mm, 0.45 Mm). The SWIRL algorithm identified two clockwise vortices in these locations. However, the instantaneous velocity streamlines do not indicate the presence of spiraling flows. Generally, the code proved to be reliable in identifying vortical motions at the photospheric level. Out of 20 identified swirls in the snapshot of Fig. 2, only two were misidentified, giving an estimated accuracy of ∼90%.

3.1.2. Chromosphere

To further assess the reliability of the SWIRL algorithm when applied to realistic numerical simulations of the solar atmosphere, we repeated the identification analysis on a chromospheric section of the simulation box. The chosen subsection covers the same horizontal domain at the same time instance as taken in Sect. 3.1.1, but at z = 700 km above the average surface of optical depth τ₅₀₀ = 1, which corresponds to the bottom of the chromosphere (see Fig. 1). The identified vortices are shown in Fig. 5, and are also shown in more detail in the zoom-in plots of Fig. 6.

Fig. 5. Vortices identified by the SWIRL algorithm in the same horizontal domain as that of Fig. 3 but at the base of the chromosphere (z = 700 km). The location and effective size of the identified vortices are indicated by colored disks. Clockwise vortices are shown in purple, while counter-clockwise ones are shown in green. The vertical magnetic field Bz is color coded and saturated at ±200 G. The gray squares denote the 0.6 × 0.6 Mm2 regions shown in Fig. 6.

Fig. 6. Zoom-in plots of the chromospheric regions outlined in Fig. 5. The location and effective size of the identified vortices are indicated by colored disks. The vertical magnetic field Bz is color coded, while the horizontal velocity field is represented by instantaneous streamlines.

At first sight, we notice that the swirls identified in the chromospheric layer appear to be more numerous and larger than the photospheric ones. Indeed, in Fig. 5 there are 74 vortices, with the largest one measuring 266 km in diameter. Multiple swirls are found in the magnetic region around (x,  y) = (1.0 Mm,  2.0 Mm), which stems from the strong photospheric magnetic flux concentration visible at the same coordinates in Fig. 3. Panel C of Fig. 6 shows a 1.0 × 1.0 Mm² close-up view of that region with streamlines derived from the horizontal velocity field and multiple spiraling patterns of different orientation can be seen. Battaglia et al. (2021) found that multiple swirls typically coexist in strong and complex magnetic flux concentrations in numerical simulations, dubbing this type of formation “superposition of swirls”.

Overall, the SWIRL algorithm identified most of the swirls in the chromospheric section of the simulation, as we can infer from the horizontal velocity field streamlines shown in Fig. 6. The estimated effective radii also correlate well with the visual size of the spiraling streamlines. There are nonetheless a few exceptions. For example, a small-scale clockwise vortex at (x, y)∼(0.6 Mm,  0.35 Mm) of panel D appears to have been missed, while a misidentification probably occurred around (x, y)∼(0.65 Mm,  0.85 Mm) of the same panel. Moreover, the radius of the relatively large counterclockwise vortical system shown in the left of panel B is presumably underestimated. An analysis of the radial profile of the tangential velocity (as done by, e.g., Silva et al. 2020) would be necessary to draw robust conclusions, but based on visual inspection of closed instantaneous streamlines, we can estimate that the size of the vortex, as estimated visually, could be as much as four times larger than that computed by SWIRL.

3.1.3. Three-dimensional structures

To investigate the three-dimensionality of the vortical structures self-consistently emerging in the simulation, we applied the SWIRL algorithm to the full 9.6 × 9.6 Mm² horizontal domain at all heights between z = −300 km (surface layers of the convection zone) and z = 1000 km (middle chromosphere). For this analysis, we used the same parameters (Table 1) at all heights.

As the automated identification is carried out on two-dimensional horizontal slices, only vertically extending vortices will be identified by our approach. Horizontal small-scale swirls have also been observed in the solar atmosphere (see, e.g., Steiner et al. 2010; Fischer et al. 2020), but they probably do not impact the upward transport of energy and mass as they do not reach the upper atmospheric layers.

To construct three-dimensional swirling structures, we search for vertical alignments between two-dimensional swirls identified at different heights in the simulation box. For this purpose, we consider two swirls with the same orientation to be part of the same vortical structure if the distance between their centers is smaller than a certain threshold. For this study, we chose the threshold to be 40 km in the horizontal direction over a vertical distance of 20 km, which corresponds to four grid cells horizontally and two grid cells vertically. In this way, a missed identification in one plane between two adjacent planes with corresponding identification does not preclude the identification of the full three-dimensional structure.

Moreover, as the horizontal threshold is smaller than the swirl average radius (see Sect. 3.2), the risk of two swirls being improperly connected is minimal. Using larger thresholds would increase the risk of erroneously connecting two separate swirls. Using excessively small thresholds carries the danger of missing a three-dimensional structure when the SWIRL algorithm misses the detection of a vortex in a single plane.

We started with the two-dimensional swirls identified in the horizontal plane located at z = 700 km. We then looked for horizontally aligned swirls in the plane 20 km below and above it. Whenever such an alignment was found, we reiterated the process starting from the previously connected two-dimensional swirl. In this way, we can construct swirls of coherent vertical extension that represent the three-dimensional extension of the two-dimensional swirls identified by the SWIRL algorithm on the different horizontal planes.

The three-dimensional vortices identified as above are shown in Fig. 7 for the time instance t = 5774 s of the simulation. We note that only the vortices reaching the height of z = 700 km are displayed in this figure, because this was the starting point for the three-dimensional stacking process. Vortical structures that are restricted to the surface layers of the convection zone or photosphere are omitted, as are purely chromospheric vortices that do not extend to the photosphere. Therefore, Fig. 7 shows only three-dimensional swirls that connect the photosphere to the chromosphere.

Fig. 7. Three-dimensional vortical structures identified for the time instance t = 5774 s of the CO5BOLD simulation. The displayed structures are obtained by stacking two-dimensional vortices in different height levels that are sufficiently well aligned with each other in the vertical direction. Only those structures that are rooted in the photosphere and reach the plane at z = 700 km are displayed. The three-dimensional vortices are colored according to the Rortex value R averaged over their surface at each height z. The surface of optical depth τ500 = 1 is shown with the vertical magnetic field Bz color coded on it. The black box outlines the 4.0 × 4.0 Mm2 horizontal domain used for Figs. 3 and 5, while zoomed-in renderings of the two red boxes labeled A and B are shown in Figs. 8 and 9, respectively.

The majority of the vertically extending swirls stem from photospheric magnetic flux concentrations, as we can see from the vertical magnetic field B_z color coded on the τ₅₀₀ = 1 surface of Fig. 7. Moreover, multiple swirls coexist in strong and complex magnetic foot points, which is in agreement with the results obtained in Figs. 3 and 5 from the two-dimensional sections.

Figure 8 shows an example of a superposition of swirls in more detail. The three-dimensional domain, which encloses the large magnetic flux concentration located at (x, y)∼(1.0 Mm,  2.0 Mm) in Figs. 3 and 5, is outlined by the red box labeled A in Fig. 7. Multiple vortices are identified in this patch and can be visualized by the instantaneous velocity field streamlines shown in the right panel of Fig. 8. The magnetic field lines shown in the left panel are mostly vertically oriented. This is typical in strong magnetic flux concentrations with plasma-β ≪ 1, where plasma-β is the ratio between the gas pressure, p_g, and the magnetic pressure, p_m = B²/8π. We recall that swirling motions and essentially untwisted, vertically oriented magnetic fields are not mutually exclusive. Such a configuration can be thought of as a quasi-rigidly rotating, stiff magnetic flux concentration.

Fig. 8. Three-dimensional rendering of a superposition of swirls stemming from a relatively large and complex small-scale photospheric magnetic flux concentration. Left: Identified three-dimensional swirls colored according to the mean Rortex value R, as in Fig. 7. Thick tubes represent magnetic field lines with the intensity of the magnetic field color coded on them. The corrugated surface near z = 0 km represents the τ500 = 1 surface. Right: Instantaneous streamlines of the velocity field belonging to the vortical structures. The strength of the velocity field is color coded on the streamlines.

Also visible in Fig. 7 are isolated swirls stemming from relatively small and weak magnetic footpoints. An example of such an event is shown in Fig. 9, which is outlined by the red box labeled B in Fig. 7. The magnetic field is weaker in this case; hence the plasma-β is closer to unity in that region, in particular in the photospheric layers. Under these conditions, the flow dominates the magnetic field and the frozen-in magnetic field lines are dragged by the rotating plasma. Indeed, we observe slightly twisted magnetic field lines in the proximity of the three-dimensional swirl. The orientation of the twist in the magnetic field lines (counterclockwise) is contrary to the rotation of the flow (clockwise) when thinking of upwardly directed positive polarity. Such a configuration is compatible with the propagation of an Alfvénic pulse (see Liu et al. 2019a; Battaglia et al. 2021) and the event shown in Fig. 9 is structurally similar to the one analyzed by Battaglia et al. (2021), which proved to be Alfvénic in nature.

Fig. 9. Three-dimensional rendering of an isolated swirling structure stemming from a relatively small photospheric magnetic flux concentration. Left: Identified three-dimensional swirls colored according to the mean Rortex value R, as in Fig. 7. Thick tubes represent magnetic field lines with the intensity of the magnetic field color coded on them. The corrugated surface near z = 0 km represents the τ500 = 1 surface. Right: Instantaneous streamlines of the velocity field belonging to the vortical structures. The strength of the velocity field is color coded on the streamlines.

3.2. Statistics

In this section, we investigate the properties of small-scale swirls in the simulated solar atmosphere from a statistical point of view. For this purpose, we ran the SWIRL algorithm on the horizontal planes of 30 time instances of the CO5BOLD simulation, covering a total physical time interval of 2 h. We used the full 9.6 × 9.6 Mm² horizontal extent of the simulation domain between z = −300 km and z = 1000 km. Our statistical analysis assesses the properties of swirls in the surface layers of the convection zone, the photosphere, and the low chromosphere of the numerical simulation. We note that all the swirls identified in the 30 time instances have been taken into account for the analysis presented in this section, regardless of whether or not they are part of a three-dimensional structure.

3.2.1. Vertical profile of swirl properties

The distributions of the number density per unit area, n_2D, the effective radius, r_eff, and the effective rotational period, P_eff, of the identified swirls as a function of height, z, are shown in Fig. 10. The data at each height z are generally not symmetrically distributed and are best fitted by a generalized extreme value distribution. Therefore, we show the 68.2% and 95.4% percentile areas around the median, labeled 1σ and 2σ, respectively.

Fig. 10. Statistical distributions as a function of height z of the number density of swirls per unit area, n2D (top), the effective radius, reff (middle), and the effective period of rotation, Peff (bottom). The median and the 1σ and 2σ deviations of the distributions are shown at each height z. The average optical surface τ500 = 1 (z = 0 km) is marked by a dotted line, while the heights of the surface layers of the convection zone, photosphere, and low chromosphere used in the analysis are indicated by dashed blue lines.

The number density of swirls, n_2D, is computed as the ratio between the average number of identified swirls at each height, N_s, and the area of a horizontal plane through the simulation domain, A_box = 9.6 × 9.6 Mm². The top panel of Fig. 10 shows that the number density decreases from n_2D ∼ 4 Mm⁻² in the surface layers of the convection zone (z = −300 km), reaching a minimum value of n_2D ∼ 1 Mm⁻² at around z = 200 km. Turbulent convection is a natural source of vortices, which explains the large abundance of identified swirls below the average optical surface τ₅₀₀ = 1. Into the upper photosphere and chromosphere, the number density increases again up to n_2D ∼ 5 Mm⁻² around z = 1000 km. The ratio between the statistical numbers of chromospheric and photospheric swirls corroborates the visual impression we note when comparing Figs. 3 and 5. This scenario is also in agreement with three-dimensional renderings of the swirling strength criterion⁴ shown by Moll et al. (2012) and Battaglia et al. (2021) in numerical simulations. However, the origin of this difference is still not well understood.

The middle panel of Fig. 10 shows the distribution of the effective radius of the swirls, r_eff, computed via Eq. (2), as a function of height z. The median profile, as well as the 1σ and 2σ percentiles, are roughly flat in the surface layers of the convection zone and in the chromosphere, while a slight rise characterizes the low photosphere. The median value is r_eff ∼ 50 km throughout the surface layers of the convection zone and r_eff ∼ 60 km in the upper photosphere. The distribution is skewed towards larger values, but 97% of the radii of the identified swirls measure less than 100 km and 150 km in the subsurface region and in the low chromosphere, respectively. Swirls are therefore statistically larger in the upper layers of the simulated domain. This growth can be explained by the expansion of the plasma ascending into the photosphere caused by the steep decrease in mass density (Nordlund et al. 1997).

The distribution of the effective rotational period, P_eff, of the identified swirls is shown in the bottom panel of Fig. 10. The effective rotational period of a swirl is computed as

$$\begin{array}{c}\hfill P_{\mathrm{eff}}=\frac{4\pi }{{|⟨R⟩}_{\mathrm{swirl}}|},\end{array}$$(3)

where ⟨R⟩_swirl is the average Rortex criterion computed over the swirl area. The median of the distribution reaches its peak of P_eff ∼ 140 s around z = 400 km with a marked skewness towards larger values, that is, towards slower swirls. This result is compatible with the growth of the typical swirl radius seen in the middle panel, which also reaches its maximum value at around the same height. Indeed, the growth in size caused by the expansion of the photospheric plasma causes the swirls to rotate slower because of the conservation of angular momentum. Moreover, the structure of the distribution is also in agreement with the vertical profiles of average vorticity and swirling strength presented by Moll et al. (2011), Canivete Cuissa & Steiner (2020), and Battaglia et al. (2021).

3.2.2. Swirls and magnetic fields

Next, we investigate the relation between the vertical magnetic field, B_z, and the properties of the identified swirls. Figure 11 shows the bivariate distribution of the average Rortex criterion ⟨R⟩_swirl and the average vertical magnetic field ⟨B_z⟩_swirl of the identified swirls. The averages are taken over the area of the swirls in the three horizontal sections corresponding to the surface layers of the convection zone, the photosphere, and the chromosphere. The different panels correspond to the swirls identified at the bottom of the chromosphere (z = 700 km), in the photosphere (z = 100 km), and in the surface layers of the convection zone (z = −200 km), as shown in Fig. 1. The effective radius of the swirl is color coded.

Fig. 11. Bivariate distribution of rotational and magnetic characteristics of vortices at z = 700 km (chromosphere, top), z = 100 km (photosphere, middle), and z = −200 km (surface layers of the convection zone, bottom). Every identified vortex in these layers is represented by a scatter point according to the Rortex criterion ⟨R⟩swirl and the vertical magnetic field ⟨Bz⟩swirl averaged over their area. The effective radius reff of the vortex is color coded.

The distributions are symmetric with respect to the sign of the Rortex criterion in all three panels. Therefore, there is no preferred orientation for small-scale swirls in the simulated solar atmosphere. This result is in agreement with the observations reported by, for example, Giagkiozis et al. (2018) and Liu et al. (2019a).

In the surface layers of the convection zone, swirls are almost homogeneously distributed with respect to the vertical magnetic field for |⟨B_z⟩_swirl|≲10² G. The properties of these swirls generated by turbulence are expected to be independent of the magnetic field in weakly magnetized regions, because in such circumstances the magnetic field does not affect the dynamics of the plasma. However, an over-density of swirls can be found for |⟨B_z⟩_swirl|≳10² G, especially in the positive-polarity end, meaning that swirls tend to be particularly associated with hecto-Gauss magnetic fields. Magnetic flux concentrations can impact the convective dynamics below the solar surface and couple it to the photosphere (see, e.g., Battaglia et al. 2021, Fig. 2), so that swirls in highly magnetized subsurface regions are part of three-dimensional photospheric vortical structures. We do not find any particular pattern regarding the effective radius of the swirls. The majority of the identified swirls measure between ∼30 km and ∼100 km, which is in accordance with the results of Fig. 10.

A clear asymmetry towards positive-polarity magnetic fields is noticeable in the middle and top panels of Fig. 11, which correspond to photospheric and chromospheric layers, respectively. We already encountered this asymmetry in the polarity of the vertical magnetic field in Figs. 3 and 5, and its origin can be traced back to the initial conditions of the numerical simulation. Indeed, even after relaxation, the polarity of the initial magnetic field persists in most of the photospheric magnetic flux concentrations and within the magnetic canopy of the low chromosphere. The fact that most photospheric and chromospheric vortices are found in regions with positive-polarity magnetic field is therefore a consequence of our choice of the initial condition to mimic a magnetic network patch of preferred polarity; we expect this asymmetry to be lifted in numerical simulations with no preferred initial magnetic configuration. In the surface layers of the convection zone, the initial imbalance is leveled out by the action of a subsurface small-scale turbulent dynamo (see, e.g., Rempel 2014), which ultimately generates the negative-polarity magnetic flux concentrations hosting the less frequent swirls on the left side of the top and middle panels of Fig. 11.

On average, chromospheric swirls are larger than photospheric ones. We discuss this difference in Sect. 3.2.1. However, we notice that the largest swirls in the photosphere and in the chromosphere are found in correspondence with strong magnetic fields. In Fig. 7, we see that most of the coherent three-dimensional vortical structures in the simulated atmosphere are anchored in photospheric magnetic flux concentrations. Therefore, magnetically dominated regions appear to provide preferable conditions for the creation and preservation of large and coherent vortical structures extending throughout the photosphere and into the chromosphere.

The distribution in the top panel of Fig. 11 forms a “butterfly” pattern, which hints at the existence of a relation between the Rortex criterion and the vertical magnetic field in chromospheric swirls. The stronger the magnetic field hosting the vortex, the faster it rotates. Moreover, the relation also depends on the radius of the swirl, as the growth of ⟨R⟩_swirl as a function of ⟨B_z⟩_swirl is reduced for larger swirls. In the following, we propose a simple analytical model to explain this relation.

Figure 12 is analogous to Fig. 11 but shows the average plasma-β over the swirl area, ⟨β⟩_swirl, instead of the average vertical magnetic field, ⟨B_z⟩_swirl. Similar to Fig. 11, the distribution is symmetric with respect to the sign of the average Rortex criterion, ⟨R⟩_swirl, in all panels, showing no preferred direction of rotation for the identified swirls. In the surface layers of the convection zone, the vast majority of swirls are found in β > 1 conditions, which means that the gas dominates over the magnetic pressure. These swirls are most certainly induced by the turbulent dynamics of the convection zone.

Fig. 12. Bivariate distribution of the Rortex criterion, ⟨R⟩swirl, and plasma-β averaged over the area of the identified vortices at z = 700 km (chromosphere, top), at z = 100 km (photosphere, middle), and at z = −200 km (surface layers of the convection zone, bottom). The effective radius reff of the vortex is color coded.

In the photosphere, we observe the emergence of two different populations of swirls. The first group, characterized by plasma-β > 1, has the same convective origin as those featured in the bottom panel. A second collection of swirls is instead found in β ≲ 1 conditions, where the magnetic field dictates the dynamics of the plasma. The swirls belonging to the second group are embedded in strong magnetic flux concentrations and can represent the footpoints of the coherent three-dimensional structures observable in Fig. 7.

Finally, in the top panel of Fig. 12, we notice a butterfly pattern similar to that seen in Fig. 11. The fastest and largest swirls are characterized by low β values, which correspond to the ones with high ⟨B_z⟩_swirl in Fig. 11. The dynamics of these swirls are dominated by the magnetic field and the model proposed in Sect. 3.2.3 qualitatively applies.

There are also a large number of chromospheric swirls for which the local β is larger than one. These swirls populate weakly magnetized areas of the chromosphere. In these regions, purely hydrodynamical mechanisms, such as baroclynic forces or shocks, could be at the origin of these chromospheric swirls, which appear to be locally produced.

3.2.3. A simple model of magnetic chromospheric swirls

Let us consider a chromospheric swirl coupled to a strong magnetic flux tube, such as those identified above. For simplicity, we assume their shape to be a cylinder. We further assume the system to be in stationary magnetohydrodynamic radial equilibrium and to rotate as a rigid body. The MHD momentum equation in the radial coordinate, r, can then be written as

$$\begin{array}{c}\hfill {\partial }_{r}p-(\mathit{B}·\mathbf{\nabla })B_r=\rho \frac{v_{\varphi }^2}{r},\end{array}$$(4)

where p = p_g + p_m is the total pressure, that is the sum of the gas pressure, p_g, and the magnetic pressure, p_m = B²/8π; B = (B_r, B_ϕ, B_z) is the magnetic field in cylindrical coordinates; ρ is the plasma density; and v_ϕ is the plasma rotational velocity. In the chromosphere, a strong magnetic flux tube is characterized by plasma-β ≲ 1, and therefore it is safe to further assume that the total pressure p is dominated by the magnetic field component, that is, p ≃ p_m. Moreover, the rigid-body rotational velocity of the plasma is related to the angular velocity by v_ϕ = Ωr. Taking this into account, Eq. (4) becomes

$$\begin{array}{c}\hfill \frac{1}{8\pi }{\partial }_r\left({\mathit{B}}^2\right)-(\mathit{B}·\mathbf{\nabla })B_r=\rho {\mathrm{\Omega }}^{2}r.\end{array}$$(5)

We model a chromospheric section of the magnetic flux tube with purely vertical magnetic field B_z = B(r) and density ρ. Therefore, the magnetic field inside the cylinder is B = B(r)e_z. In this scenario, Eq. (5) can be further simplified into

$$\begin{array}{c}\hfill \frac{1}{8\pi }{\partial }_r\left(B^2\right)=\rho {\mathrm{\Omega }}^{2}r,\end{array}$$

and integrated over the radius of the swirl r, leading to

$$\begin{array}{c}\hfill \frac{1}{4\pi \rho }(B{(r)}^2-B{(0)}^2)={\mathrm{\Omega }}^2{r}^2.\end{array}$$(6)

A physical solution to the equation above exists only if B(r)² > B(0)², that is, if the rotation of the plasma is supported by a negative magnetic pressure gradient toward the vortex core. Indeed, the cores of vortices within large magnetized regions are often found to be associated with reduced magnetic pressure. We provide an example of such an event in Appendix A.

We recognize the local Alfvén speed, $v_{\mathrm{A}}(r)=B(r)/\sqrt{4\pi \rho }$, and the rotational velocity of the swirl, v_ϕ = Ωr, on the left- and right-hand sides of Eq. (6), respectively,

$$\begin{array}{c}\hfill v_{\mathrm{A}}{(r)}^2-v_{\mathrm{A}}{(0)}^2=v_{\varphi }^2,\end{array}$$(7)

where v_A(0) is the Alfvén speed computed in the vortex center. The above equation states that the Alfvén speed in the magnetic flux tube is the upper limit to the swirl rotational velocity. If we assume the magnetic field in the vortex core to be weak enough compared to the bulk of the flux tube, then the swirl rotates at approximately the local Alfvén speed, v_ϕ ≈ v_A.

For r = 0, Eq. (7) predicts v_ϕ = 0, which is consistent with the assumption of rigid-body rotation. Moreover, in Appendix A, we show that the structure of the chromospheric swirl shown in Fig. 9 qualitatively agrees with the model presented above.

For a consistency test of Eq. (7), we estimate the expected rotational velocity, $v_{\varphi }^{\exp }$, of a chromospheric swirl and compare it to the Alfvén speed in chromospheric swirls of low plasma-β conditions. We base our estimation on typical values of the effective radius, r_eff, and the Rortex criterion, R, for chromospheric swirls in plasma-β ≪ 1. From Fig. 12, we infer that such values are r_eff ∼ 120 km and R ∼ 0.2 Hz.

Using the formula

$$\begin{array}{c}\hfill v_{\varphi }=\mathrm{\Omega }r\approx \frac{R}{2}{r}_{\mathrm{eff}},\end{array}$$(8)

we calculate an expected rotational velocity for swirls at z = 700  in low plasma-β conditions of $v_{\varphi }^{\exp }\sim 12\mathrm{km}{\mathrm{s}}^{-1}$.

Figure 13 shows the probability distribution of the average Alfvén speeds computed over the effective area of swirls identified in the low chromosphere (600 km < z < 900 km). There, we differentiate between swirls in magnetic conditions (plasma-β ≤ 0.5), mixed conditions (0.5 < plasma-β ≤ 2.0), and nonmagnetic conditions (plasma-β > 2.0), for a succinct labeling. The expected rotational velocity of chromospheric swirls inferred from Fig. 12 and Eq. (8) and the distribution of Alfvén speeds averaged over the effective area of swirls in low plasma-β conditions are consistent with the analytical model and Eq. (7). For comparison, the distribution of the average sound speeds computed over the effective area of the swirls is outlined by its median together with the 5 − 95 percentile range. The estimated rotational speed and the rotational speed derived from the model are clearly above the sound speed.

Fig. 13. Probability distributions of the average Alfvén speed, vA, computed over the effective area of swirls identified between z = 600 km and z = 900 km. The total kernel density estimate is shown in black. The distributions are divided according to the average plasma-β conditions computed over the effective area of the swirls. The nonmagnetic category corresponds to plasma-β > 2.0, the mixed category corresponds to 0.5 < plasma-β ≤ 2.0, and the magnetic category corresponds to plasma-β ≤ 0.5. The red area and vertical dashed line correspond to the 5–95 percentile range and the median of the average sound speed computed over the effective area of the swirls, respectively.

Figure 14 shows a scatter plot of all the identified vortices at z = 700 km as a function of the Alfvén speed ⟨v_A⟩_swirl and the estimated rotational velocity ${⟨v_{\varphi }⟩}_{\mathrm{swirl}}=\frac{1}{2}{r}_{\mathrm{eff}}{⟨R⟩}_{\mathrm{swirl}}$, both averaged over the swirl area. A large dispersion characterizes the distribution, which is expected given the rough approximations made in deriving Eq. (7), but a linear trend is perceivable. The solid black line shows the relation v_ϕ = v_A. For β ≲ 1, we see that v_ϕ ≲ v_A for almost all swirls, confirming that v_A is an upper limit for v_ϕ. However, this limit applies less well for weak field swirls with β ≳ 1.

Fig. 14. Bivariate distribution of the local average Alfvén speed, ⟨vA⟩swirl, and the average estimated rotational velocity of the swirl, ${⟨v_{\varphi }⟩}_{\mathrm{swirl}}=\frac{1}{2}{r}_{\mathrm{eff}}{⟨R⟩}_{\mathrm{swirl}}$ at the bottom of the chromosphere (z = 700 km). The averages are computed over the swirl area. The average plasma-β over the swirl area is color coded. A power-law fit of the type y = axb is shown in dashed black. The fitted parameters are a = 0.16,  b = 1.13. Density contours of the scattered points are shown in gray.

We fitted a power-law function of the type y = ax^b to the data. The resulting curve is represented by the black dashed line in the log–log plot of Fig. 14. The fitted exponent is b = 1.16, which is quite close to the modeled linear exponent b = 1. Another measure of the linear correlation between ⟨v_A⟩_swirl and r_eff⟨R⟩_swirl for the identified chromospheric swirls can be obtained in the form of the Pearson’s correlation coefficient r_P. For the dataset shown in Fig. 14, we obtain r_P = 0.45, which demonstrates a discrete degree of linear correlation between these two quantities.

3.2.4. Torsional Alfvénic waves

We also investigated the correlation between swirls in the simulated solar atmosphere and perturbations in the magnetic field lines. Battaglia et al. (2021) reported that a toroidal perturbation in the predominantly vertically directed magnetic field can be found in upwardly propagating pulses of swirling plasma. The same authors introduced the magnetic swirling strength, λ^B, which is a measure for the toroidal components, or twists, in magnetic flux tubes. The simultaneous presence of a twist in the magnetic field lines and a vortical motion in the plasma may hint at the presence of torsional Alfvénic waves propagating cojointly with the rotating magnetic flux concentration.

Figure 15 shows the bivariate distribution of the swirls identified at the bottom of the chromosphere (z = 700 km), in the photosphere (z = 100 km), and in the surface layers of the convection zone (z = −200 km) as a function of the Rortex, ⟨R⟩_swirl, and magnetic swirling strength, ⟨λ^B⟩_swirl, averaged over the swirl area. The average strength of the vertical magnetic field over the swirl area, ⟨B_z⟩_swirl, is color coded.

Fig. 15. Bivariate distribution of the rotational characteristics of vortices at z = 700 km (chromosphere, top), z = 100 km (photosphere, middle), and z = −200 km (surface layers of the convection zone, bottom). Every identified vortex in these layers is represented by a scatter point according to the Rortex criterion ⟨R⟩swirl and the magnetic swirling strength criterion ⟨λB⟩swirl averaged over the swirl area. The vertical magnetic field ⟨Bz⟩swirl averaged over the swirl area is color coded.

If we do not consider the polarity of the vertical magnetic field, the swirls identified in the surface layers of the convection zone (bottom panel) are distributed almost symmetrically with respect to ⟨R⟩_swirl and ⟨λ^B⟩_swirl. Once more, we explain this symmetry with the isotropical turbulence that dominates the surface layers of the convection zone. However, the magnetic field orientation reveals a pattern that is even more prominent in the photosphere (middle panel): most of the swirls embedded in positive-polarity magnetic fluxes are concentrated in the top-left and bottom-right quadrants of the bivariate distribution, while the ones associated with negative-polarity magnetic fields are found in the top-right and bottom-left quadrants.

The excess of red points (⟨B_z⟩_swirl > 0 G) in the photospheric and chromospheric distributions is due to the initial conditions of the present simulations. Moreover, the pattern is less pronounced in the surface layers of the convection zone because of the swirls that are randomly generated by turbulence and are not part of coherent photospherical structures.

The pattern revealed by the middle and bottom panels of Fig. 15 can be explained if we consider the swirls to be Alfvénic in nature, as proposed by Liu et al. (2019a) and Battaglia et al. (2021). For a vertical magnetic field B = B_ze_z and an incompressible plasma in magneto-hydrostatic equilibrium in the ideal MHD approximation, a torsional Alfvén wave is characterized by velocity and magnetic field perturbations, v and b, that obey (see, e.g., Priest 2014, Chap. 4)

$$\begin{array}{c}\hfill \mathit{v}=-\frac{\omega \mathit{b}}{\mathit{k}·\mathit{B}},\end{array}$$(9)

where ω is the angular frequency of the plane wave and k is the wave-vector indicating the propagation direction. For a vertically propagating torsional Alfvén wave, that is, k = ke_z with k > 0, Eq. (9) can be simplified as

$$\begin{array}{c}\hfill \mathit{v}=-\frac{v_{\mathrm{A}}}{B_z}\mathit{b},\end{array}$$(10)

where v_A > 0 is the local Alfvén speed, while v and b are perturbations in the horizontal plane.

From Eq. (10) we conclude that the perturbations v and b are parallel or anti-parallel depending on the polarity of the vertical magnetic field, that is, on the sign of B_z. If we use the Rortex and the magnetic swirling strength criteria as proxies to quantify such perturbations, then we can write

$$\begin{array}{c}\hfill \mathrm{sign}\left(R{\lambda }^{\mathrm{B}}\right)=-\mathrm{sign}\left(B_z\right),\end{array}$$(11)

and the distributions in the surface layers of the convection zone and photosphere of Fig. 15 appear to statistically follow this relation. The clockwise vortex associated with the counterclockwise twist of the positive-polarity magnetic field lines shown in Fig. 9 is a practical illustration of Eq. (11).

The Alfvénic pattern encountered in the lower and middle layers of the simulation box seems to disappear in the chromosphere (top panel of Fig. 15). A high degree of symmetry is restored in the distribution, although the polarity of the magnetic field is predominantly positive. If the chromospheric swirls were associated with upwardly propagating Alfvén waves, we would expect Eq. (11) to be respected and the scatter points to populate mainly the top-left and bottom-right quadrants of the plot.

We find two explanations for the systematic violation of Eq. (11) in the simulated chromosphere. First, in Sect. 3.2.1, we show that chromospheric swirls are more abundant than photospheric ones. Therefore, a large fraction of swirls must be generated locally in the simulated chromosphere and, as they are not linked to a photospheric coherent structure, they do not share the same properties. Second, the boundary conditions of the simulation force the magnetic field to be strictly vertical at the top boundary, which may cause upwardly directed Alfvénic waves to be reflected. In that case, the wave-vector becomes k = −ke_z, Eq. (10) picks up a minus sign on the right-hand side, and therefore Alfvénic waves propagating downward are characterized by parallel perturbations v and b when embedded in positive-polarity magnetic fields. Consequently, swirls associated with downwardly directed torsional Alfvén pulses would populate the top-right and bottom-left quadrants of the top panel of Fig. 15. We note that the narrower distribution of data points over ⟨λ^B⟩_swirl in the chromosphere compared to the photosphere and surface layers of the convection zone is due to the fact that the magnetic swirling strength is proportional to the magnetic field strength, which in turn is much smaller in the chromosphere; it does not indicate a smaller twist angle.

To characterize the abundance of swirls exhibiting imprints of upwardly propagating Alfvénic waves, we computed the fraction of them, f_Alfévnic, for which Eq. (11) is respected. Figure 16 shows the results obtained by considering all the swirls identified in any time instances of the simulation (green curve) and those only forming three-dimensional swirling structures (blue curve). In the latter case, only the structures that reach both the surface layers (z = 0 km) and the low chromosphere (z = 700 km) are taken into account.

Fig. 16. Fraction of swirls obeying Eq. (11), fAlfévnic, as a function of height z. The green profile refers to all identified swirls, while the blue curve takes into account only those swirls that form coherent structures connecting the surface layers (z = 0 km) to the low chromosphere (z = 700 km). Shaded areas represent statistical standard deviations. The average optical surface τ500 = 1 (z = 0 km) is marked by a dotted line, while the heights of the surface layers of the convection zone, photosphere, and low chromosphere used in the analysis are indicated by dashed blue lines.

In the photosphere, approximately 80% of all the swirls show perturbations in the plasma and in the magnetic field that are compatible with torsional Alfvénic waves, which is in accordance with the pattern observed in Fig. 15. However, this fraction decreases as we move upward in the simulation box and falls below 50% in the chromosphere. Regarding swirls that belong to coherent three-dimensional structures, we notice a higher fraction at all heights, reaching ∼90% in the photosphere and ∼80% at z = 700 km.

Therefore, Fig. 16 suggests that a significant fraction of the identified coherent three-dimensional swirls present characteristics compatible with torsional Alfvénic pulses propagating upward in the simulated solar atmosphere. On the other hand, vortical structures that do not couple the photosphere to the chromosphere appear to show fewer imprints of these waves, and are therefore probably of a different nature and likely have a different origin.

4. Summary and conclusions

In this paper, we employed the recently developed SWIRL algorithm to investigate small-scale swirls in radiative MHD numerical simulations of the solar atmosphere. The methodology at the core of this algorithm considers both the local and global properties of the velocity field in the detection process. Therefore, the SWIRL algorithm is specifically tailored to identifying coherent vortical structures, whereas conventional methods, such as the vorticity or the swirling strength, can only recognize local curvatures in the flow. The identification process is automatized through the implementation of a state-of-the-art clustering algorithm. This approach requires minimal interaction with the user, which reduces the risk of human bias in the identification process, and ensures a high level of precision and consistency. In Paper I, we validated the robustness of the SWIRL algorithm against noise and turbulence.

In a first stage of the present paper, we tested the reliability of the code in identifying swirls that emerge self-consistently from the simulated photospheric and chromospheric flows. The interplay between magnetic fields, turbulence, convective flows, and shocks significantly increases the complexity of the flow with respect to the tests carried out in Paper I. In addition, fine-tuning of the algorithm parameters is necessary, especially for the number of stencils, the “noise” parameter, and the “kink” parameter. We provide the list of parameter values used in this study in Table 1; we consider these to be suitable default values for applying the SWIRL algorithm to numerical simulations of the solar atmosphere.

The algorithm detected photospheric and chromospheric swirls with high accuracy and precision based on the instantaneous streamlines of the horizontal component of the velocity field. Occasional misidentifications can occur, as shown in Figs. 4 and 6. Moreover, the identification method implemented in the SWIRL algorithm is not Galilean invariant and therefore swirls that are advected at speeds comparable to their rotational velocity could be missed. This shortcoming should not affect photospheric swirls, which are predominantly rooted in intergranular lanes and tightly coupled to magnetic flux concentrations, but it could be relevant to swirls in a dynamical environment such as the chromosphere. Further investigation of this aspect is required in order to improve the performance of the algorithm and to reduce such inaccuracies. Together with the tests carried out in Paper I, we conclude that the present SWIRL algorithm is a reliable tool for the identification of swirls in the solar atmosphere and astrophysical flows in general.

The SWIRL algorithm is currently limited to the identification of vortices in two-dimensional planes. Therefore, in order to investigate the presence of coherent three-dimensional vortical structures extending vertically in the simulation domain, we ran the algorithm on all the horizontal sections of a particular time instance of the CO5BOLD simulation. Successively, we stacked up the identified vortices that were approximately vertically aligned in order to reconstruct the three-dimensional structures. Figure 7 shows an example result of this procedure and demonstrates that the vast majority of small-scale swirls that reach the chromosphere stem from photospheric magnetic flux concentrations. Depending on the intensity and complexity of these magnetic regions, we can observe isolated vortex tubes or multiple swirls that coexist and interact within a magnetic element (see also Battaglia et al. 2021).

The procedure adopted in this paper to find three-dimensional swirls is relatively basic and vortical structures in the simulation domain may have been missed. Moreover, this method can outline vertically extending structures only, while horizontally directed vortex tubes and arches have been shown to populate the solar atmosphere as well. An identification code that could handle the three-dimensional flow of a simulation domain (or subdomain) would, in principle, be required to properly characterize and study these structures. As we argue in Paper I, the SWIRL methodology and algorithm can, in principle, be extended to three dimensions. However, the computational costs that such an upgrade would entail complicate matters considerably, especially regarding the automated clustering task. The identification process on a 960  ×  960 plane of the CO5BOLD simulations typically takes around ∼2 min on a single CPU. Of this time, approximately 90% of the computational time is dedicated to the clustering step. As the size of the data set increases, the percentage of time spent on clustering is expected to rise due to the inherent computational complexity of the clustering algorithm. Consequently, the overall time required for the complete identification process will also increase.

In the second part of this paper, we present a statistical analysis of the properties of small-scale swirls in numerical simulations of the solar atmosphere and near-surface convection zone. Our study indicates that, statistically, around one small-scale swirl can be found in each Mm² of the photosphere, while in the low chromosphere the number density of swirls grows to approximately 4 Mm⁻². Because of these different abundances, approximately three out of four chromospheric swirls must be generated locally in the chromosphere, but the physical mechanism responsible for this generation is still unknown.

Shelyag et al. (2012) and Canivete Cuissa & Steiner (2020) presented an analysis of the generation of vortical motions based on the evolution equation of the vorticity and of the swirling strength, respectively. Both papers concluded that the origin of swirling motions in the chromosphere must be traced to the action of magnetic fields. However, in Sect. 3.2.2 we show that a considerable number of chromospheric swirls are found in high plasma-β regions, which indicates that hydrodynamical forces or shocks may be responsible for part of the small-scale swirls in the upper solar atmosphere.

If we extrapolate the obtained number densities to the whole Sun, our results hint at the steady presence of ∼6 × 10⁶ and ∼2 × 10⁷ swirls in the photosphere and the chromosphere, respectively. These numbers greatly exceed previous estimations based on simulations and observations reported in the literature (see, e.g., Wedemeyer-Böhm et al. 2012; Kato & Wedemeyer 2017; Giagkiozis et al. 2018; Liu et al. 2019a,b; Dakanalis et al. 2022). On the other hand, our analysis reveals that the average size of the swirls in the simulated atmosphere settles down to around 50 − 60 km in radius, although larger vortices can be systematically found in the chromosphere. In summary, swirls may be more numerous and smaller than previously thought.

For comparison, we mention in the following a few studies addressing the number densities and typical radii of small-scale swirls in the solar atmosphere. A more comprehensive list of these values can be found in Tziotziou et al. (2023). Wedemeyer-Böhm et al. (2012) counted on average ∼2.0 × 10⁻³ Mm⁻² (3.8 arcmin⁻²) long-lived chromospheric swirls with a typical radius of 1.4 × 10³ km (2.0 arcsec) from observations obtained with the CRisp Imaging SpectroPolarimeter (CRISP) instrument of the Swedish 1m Solar Telescope (SST). Automated surveys have been carried out by Giagkiozis et al. (2018) and Liu et al. (2019b) on photospheric observations obtained with CRISP/SST and with the Solar Optical Telescope (SOT) on board the Hinode satellite, respectively. The first study identified on average 2.7 × 10⁻² Mm⁻² swirls with a mean radius of 290 km, while the second one found number densities that are closer to our result, namely 2.4 × 10⁻¹ Mm⁻², but with an average radius of 280 km. Using a new automated identification method based on the morphological characterization of Hα spectral lines (Dakanalis et al. 2021), Dakanalis et al. (2022) found a number density of chromospheric swirls of 8 × 10⁻² Mm⁻² and an average radius of 1.3 × 10³ km from CRISP/SST observations. From numerical simulations, Kato & Wedemeyer (2017) also detected a relatively high number of chromospheric swirls with an average number density of 8.6 × 10⁻¹ Mm⁻². However, in this case, the average radius of the identified swirls was 338 km.

Nevertheless, we would not recommend a blunt comparison between the results presented in this paper and previous results found in the literature. First, the identifications performed on observational data heavily rely on the methods used to estimate the horizontal velocity fields. For example, LCT techniques should be used with caution as they present several limitations, especially in estimating granular and subgranular flows (Verma et al. 2013; Tremblay et al. 2018). To our knowledge, the only study that is not affected by this issue is the one presented by Dakanalis et al. (2022), because these authors detected swirls directly from chromospheric filtergrams. One possible solution to consider is a deep learning approach, as proposed by Asensio Ramos et al. (2017). However, it is important to be aware that the simulations on which the models are trained may introduce bias into the results if their vortical flows are not consistent with the real ones.

Second, the properties of vortical motions appear to be heavily dependent on the spatial resolution available, as shown by Yadav et al. (2020) in numerical simulations. Other details regarding the simulations, such as the initial and boundary conditions or the strength of the magnetic field, can also deeply affect the characteristics of vortical motions (see, e.g., Appendix A of Battaglia et al. 2021 or Canivete Cuissa et al. 2022 for simulations with different initial magnetic fields). A comprehensive investigation of the influence of the numerical setup on the characteristics of vortices is an essential step towards a deeper understanding of their formation and evolution in numerical simulations of the solar atmosphere.

Finally, different datasets and different automated algorithms have been used for the identification of swirls in the solar atmosphere. A comparative study between the available algorithms would be necessary to assess strengths and weaknesses of the different detection methods.

Given the clear correlation between magnetic flux concentrations and vortical motions, we investigated how the properties of the small-scale swirls vary as a function of the vertical component of the magnetic field. We find indications of a relation between the vertical magnetic field, the angular velocity, and the size of chromospheric swirls. We explain this relation with a simple model of a homogeneously dense magnetic flux tube in a low-plasma-β environment with a magnetic pressure gradient that supports its rotation.

This model assumes stationary radial equilibrium and rigid-body rotation. We acknowledge that these assumptions are only very basic and do not accurately capture the complex nature of chromospheric swirls. For example, swirls do not rigidly rotate (see, e.g., Silva et al. 2020) and flows in the highly dynamical chromosphere are not stationary. Nevertheless, they allow a straightforward analytical analysis and interpretation of the swirl properties and of our statistical results. The model suggests that chromospheric swirls can rotate at maximum speeds that approach the local Alfvén speed, and the data gathered from the simulation support this conclusion. To our knowledge, this result represents a new property of chromospheric swirls that has not yet been investigated, and that could have profound implications for the total energy transport associated with small-scale vortical motions in the solar atmosphere.

Park et al. (2016) observed a chromospheric swirl measuring 0.5 − 1.0 Mm and rotating at an average speed of 13 km s⁻¹ with CRISP/SST. Although there is no available information regarding the magnetic field for this particular event, the observed rotational velocity is in the range of typical Alfvén speeds in chromospheric conditions. Other observational studies, such as those of Wedemeyer-Böhm & Rouppe van der Voort (2009), Morton et al. (2013), Liu et al. (2019a), and Murabito et al. (2020), suggest slower average rotational velocities for chromospheric swirls, rarely exceeding 2 km s⁻¹. We expect future high-resolution observations with the Daniel K. Inouye Solar Telescope (DKIST) to shed light on the real rotational velocities of swirls in the solar atmosphere.

Finally, we carried out a statistical analysis in order to investigate the possible Alfvénic nature of photospheric and chromospheric swirls. We find a clear relation between the orientation of the identified swirls, the orientation of the toroidal magnetic perturbations in the swirling area, and the polarity of the vertical magnetic field emerging from the data, in particular in the photosphere. In 80% of the identified photospheric swirls, this relation is compatible with the propagation of torsional Alfvén waves according to Eq. (10). In the chromosphere, the correlation between swirls and Alfvénic waves seems to vanish, probably because of the local generation of vortical motions, which do not have a photospheric counterpart. However, when considering only those swirls that extend from the photosphere to the chromosphere, we find 90% and 80% of them to show imprints of Alfvénic waves in the photosphere and chromosphere, respectively. Together with the rotational speeds reported for chromospheric swirls, our study strongly suggests a strong connection between coherent vortical structures in the solar atmosphere and Alfvénic waves.

In conclusion, this paper demonstrates the reliability and the capability of the SWIRL algorithm in identifying vortical motions in magnetized, turbulent, and highly dynamical astrophysical flows such as those characterizing the solar atmosphere. We believe that, combined with state-of-the-art methods for the estimation of horizontal velocities and high-resolution observational campaigns, the SWIRL algorithm can provide reliable information from which rigorous conclusions can be drawn as to the statistical properties and nature of swirls in the solar atmosphere.

Acknowledgments

The authors acknowledge support by Swiss National Science Foundation under grant ID 200020_182094 and the University of Zürich under grant UZH Candoc 2022 ID 7104. This work has profited from discussions with the team of K. Tziotiou and E. Scullion (conveners) “The Nature and Physics of Vortex Flows in Solar Plasma” and with the team of P. Keys (convener) “WaLSA: Waves in the Lower Solar Atmosphere at High Resolution” (https://walsa.team/), both at the International Space Science Institute (ISSI).

References

Appendix A: Modeling chromospheric swirls

We show that the model proposed in Sect. 3.2.3 for chromospheric swirls is qualitatively consistent with the swirling structures observable in the simulation. We chose to analyze the swirl shown in the three-dimensional rendering of Fig. 9, because it represents a nearly ideal case of an isolated swirl. We recall that our model consists of a magnetic flux tube that swirls as a rigid body and is characterized by plasma-β ≲ 1. We find that the magnetic flux tube (or vortex) center must correspond to a minimum in magnetic pressure and that the rotational velocity is comparable to the local Alfvén speed, v_A.

The top panel of Fig. A.1 shows the plasma-β in a 1.0 × 1.0 Mm² section centered on the swirl at the height of 700 km, that is, at the bottom of the chromosphere. The swirling plasma is depicted through the vector plot. The size of the chromospheric swirl is probably underestimated by the SWIRL algorithm. We notice that the plasma-β is consistently lower than 1 in the vortical region, with minimal values reaching plasma-β ∼ 0.05 in the outer layers of the swirl. Therefore, the dynamics of the swirl under consideration are mostly dominated by the magnetic field, in agreement with our assumption.

Fig. A.1. Physical conditions in the chromospheric horizontal section (z = 700 km) of the isolated vortical structure visualized in Fig. 9. The plasma-β (top), magnetic pressure, pm (middle), the local Alfvén speed, vA (bottom) are shown. The horizontal velocity field is represented by a vector field, while the purple circle denotes the position and effective size of the chromospheric swirl.

The magnetic pressure, p_m, is shown in the middle panel of Fig. A.1. The swirling plasma is highly magnetized, while regions further away show a very low magnetic pressure. Nonetheless, the vortex core is found in a local minimum of magnetic pressure, which means that a negative pressure gradient —which opposes the centrifugal force— is directed towards the center of the swirl. Finally, we show in the bottom panel of the same figure that the rotational speed of the chromospheric swirl (∼1 km s⁻¹) is of approximately the same magnitude as the local Alfvén speed (∼1 − 4 km s⁻¹). Moreover, we also find a local minimum of the Alfvén speed in proximity to the vortex core, which corroborates Eq. (7) and our findings regarding Fig. 14.

In conclusion, the simple model of a chromospheric swirl embedded in a strong magnetic flux tube we put forward in Sect. 3.2.2 seems to qualitatively capture the main physical properties of an isolated vortical structure that self-consistently emerged from the CO5BOLD simulation.

All Tables

All Figures

Fig. 1. Average stratifications of density, ρ, temperature, T, and rms of the vertical component of the velocity field, vz,  rms. The profiles represent averages, both temporally across the 30 time instances of the CO5BOLD simulation, and spatially across the horizontal sections of the domain. The heights for the analysis of swirls in the surface layers of the convection zone (z = −200 km), in the photosphere (z = 100 km), and in the low chromosphere (z = 700 km) are indicated by blue dashed lines.

In the text

Fig. 2. Two-dimensional horizontal subsection of the simulated photosphere. The section measures 4.0 × 4.0 Mm2 and is taken at z = 100 km. The simulation time instance corresponds to t = 5774 s. The horizontal velocity field of the subsection is depicted using a vector plot. The length of each arrow corresponds to the magnitude of the horizontal flow, and a reference scale is included in the bottom-right corner. Top: Vertical magnetic field Bz at z = 100 km. Middle: Rortex criterion R. Bottom: G-EVC map. Contours where R ≠ 0 are shown in gray in the middle panel.

In the text

Fig. 3. Vortices identified by the SWIRL algorithm in the two-dimensional, horizontal velocity field of Fig. 2. The location and effective size of the identified vortices are indicated by colored disks. Clockwise vortices are represented by purple disks, while counterclockwise ones are shown in green. The vertical magnetic field Bz is color coded and saturates at ±1000 G. The gray squares denote the 0.6 × 0.6 Mm2 regions shown in Fig. 4.

In the text

	Fig. 4. Zoom-in plots of the photospheric regions outlined in Fig. 3. The location and effective size of the identified vortices are indicated by colored disks. The vertical magnetic field Bz is color coded with the same scale as in Fig. 3, while the horizontal velocity field is represented by instantaneous streamlines.
In the text

Fig. 5. Vortices identified by the SWIRL algorithm in the same horizontal domain as that of Fig. 3 but at the base of the chromosphere (z = 700 km). The location and effective size of the identified vortices are indicated by colored disks. Clockwise vortices are shown in purple, while counter-clockwise ones are shown in green. The vertical magnetic field Bz is color coded and saturated at ±200 G. The gray squares denote the 0.6 × 0.6 Mm2 regions shown in Fig. 6.

In the text

	Fig. 6. Zoom-in plots of the chromospheric regions outlined in Fig. 5. The location and effective size of the identified vortices are indicated by colored disks. The vertical magnetic field Bz is color coded, while the horizontal velocity field is represented by instantaneous streamlines.
In the text

Fig. 7. Three-dimensional vortical structures identified for the time instance t = 5774 s of the CO5BOLD simulation. The displayed structures are obtained by stacking two-dimensional vortices in different height levels that are sufficiently well aligned with each other in the vertical direction. Only those structures that are rooted in the photosphere and reach the plane at z = 700 km are displayed. The three-dimensional vortices are colored according to the Rortex value R averaged over their surface at each height z. The surface of optical depth τ500 = 1 is shown with the vertical magnetic field Bz color coded on it. The black box outlines the 4.0 × 4.0 Mm2 horizontal domain used for Figs. 3 and 5, while zoomed-in renderings of the two red boxes labeled A and B are shown in Figs. 8 and 9, respectively.

In the text

Fig. 8. Three-dimensional rendering of a superposition of swirls stemming from a relatively large and complex small-scale photospheric magnetic flux concentration. Left: Identified three-dimensional swirls colored according to the mean Rortex value R, as in Fig. 7. Thick tubes represent magnetic field lines with the intensity of the magnetic field color coded on them. The corrugated surface near z = 0 km represents the τ500 = 1 surface. Right: Instantaneous streamlines of the velocity field belonging to the vortical structures. The strength of the velocity field is color coded on the streamlines.

In the text

Fig. 9. Three-dimensional rendering of an isolated swirling structure stemming from a relatively small photospheric magnetic flux concentration. Left: Identified three-dimensional swirls colored according to the mean Rortex value R, as in Fig. 7. Thick tubes represent magnetic field lines with the intensity of the magnetic field color coded on them. The corrugated surface near z = 0 km represents the τ500 = 1 surface. Right: Instantaneous streamlines of the velocity field belonging to the vortical structures. The strength of the velocity field is color coded on the streamlines.

In the text

Fig. 10. Statistical distributions as a function of height z of the number density of swirls per unit area, n2D (top), the effective radius, reff (middle), and the effective period of rotation, Peff (bottom). The median and the 1σ and 2σ deviations of the distributions are shown at each height z. The average optical surface τ500 = 1 (z = 0 km) is marked by a dotted line, while the heights of the surface layers of the convection zone, photosphere, and low chromosphere used in the analysis are indicated by dashed blue lines.

In the text

Fig. 11. Bivariate distribution of rotational and magnetic characteristics of vortices at z = 700 km (chromosphere, top), z = 100 km (photosphere, middle), and z = −200 km (surface layers of the convection zone, bottom). Every identified vortex in these layers is represented by a scatter point according to the Rortex criterion ⟨R⟩swirl and the vertical magnetic field ⟨Bz⟩swirl averaged over their area. The effective radius reff of the vortex is color coded.

In the text

	Fig. 12. Bivariate distribution of the Rortex criterion, ⟨R⟩swirl, and plasma-β averaged over the area of the identified vortices at z = 700 km (chromosphere, top), at z = 100 km (photosphere, middle), and at z = −200 km (surface layers of the convection zone, bottom). The effective radius reff of the vortex is color coded.
In the text

Fig. 13. Probability distributions of the average Alfvén speed, vA, computed over the effective area of swirls identified between z = 600 km and z = 900 km. The total kernel density estimate is shown in black. The distributions are divided according to the average plasma-β conditions computed over the effective area of the swirls. The nonmagnetic category corresponds to plasma-β > 2.0, the mixed category corresponds to 0.5 < plasma-β ≤ 2.0, and the magnetic category corresponds to plasma-β ≤ 0.5. The red area and vertical dashed line correspond to the 5–95 percentile range and the median of the average sound speed computed over the effective area of the swirls, respectively.

In the text

Fig. 14. Bivariate distribution of the local average Alfvén speed, ⟨vA⟩swirl, and the average estimated rotational velocity of the swirl, ${⟨v_{\varphi }⟩}_{\mathrm{swirl}}=\frac{1}{2}{r}_{\mathrm{eff}}{⟨R⟩}_{\mathrm{swirl}}$ at the bottom of the chromosphere (z = 700 km). The averages are computed over the swirl area. The average plasma-β over the swirl area is color coded. A power-law fit of the type y = axb is shown in dashed black. The fitted parameters are a = 0.16,  b = 1.13. Density contours of the scattered points are shown in gray.

In the text

Fig. 15. Bivariate distribution of the rotational characteristics of vortices at z = 700 km (chromosphere, top), z = 100 km (photosphere, middle), and z = −200 km (surface layers of the convection zone, bottom). Every identified vortex in these layers is represented by a scatter point according to the Rortex criterion ⟨R⟩swirl and the magnetic swirling strength criterion ⟨λB⟩swirl averaged over the swirl area. The vertical magnetic field ⟨Bz⟩swirl averaged over the swirl area is color coded.

In the text

Fig. 16. Fraction of swirls obeying Eq. (11), fAlfévnic, as a function of height z. The green profile refers to all identified swirls, while the blue curve takes into account only those swirls that form coherent structures connecting the surface layers (z = 0 km) to the low chromosphere (z = 700 km). Shaded areas represent statistical standard deviations. The average optical surface τ500 = 1 (z = 0 km) is marked by a dotted line, while the heights of the surface layers of the convection zone, photosphere, and low chromosphere used in the analysis are indicated by dashed blue lines.

In the text

	Fig. A.1. Physical conditions in the chromospheric horizontal section (z = 700 km) of the isolated vortical structure visualized in Fig. 9. The plasma-β (top), magnetic pressure, pm (middle), the local Alfvén speed, vA (bottom) are shown. The horizontal velocity field is represented by a vector field, while the purple circle denotes the position and effective size of the chromospheric swirl.
In the text