11institutetext:  Istituto ricerche solari Aldo e Cele Daccò (IRSOL), Faculty
of Informatics, Università della Svizzera italiana, CH-6605 Locarno,
Switzerland 22institutetext: Euler Institute, Università della Svizzera
italiana (USI), CH-6900 Lugano, Switzerland 33institutetext: Kiepenheuer-
Institut für Sonnenphysik (KIS), Schöneckstrasse 6, 79104 Freiburg i.Br.,
Germany

# Small-scale vortical motions in cool stellar atmospheres

J. R. Canivete Cuissa 11 F. Riva 1122 O. Steiner 1133

(Received xxx; accepted yyy)

###### Abstract

Context. Small-scale vortices in the solar atmosphere have received
considerable attention in recent years. These events are considered potential
conduits for the exchange of energy and mass between the solar atmospheric
layers from the convective surface to the corona. Similar events may occur in
the atmospheres of other stars and play a role in energy transfer within their
atmospheres.

Aims. Our aim is to study the presence and properties of small-scale swirls in
numerical simulations of the atmospheres of cool main-sequence stars. Our
particular focus is on understanding the variations in these properties for
different stellar types and their sensitivity to the surface magnetic field.
Furthermore, we aim to investigate the role of these events in the energy
transport within the simulated atmospheres.

Methods. We analyze three-dimensional, radiative-magnetohydrodynamic, box-in-
a-star, numerical simulations of four main-sequence stars of spectral types
K8V, K2V, G2V, and F5V. These simulations include a surface small-scale dynamo
responsible for amplifying an initially weak magnetic field. Thus, we can
study models characterized by very weak, or, magnetic fields in near
equipartition. To identify small-scale vortices in horizontal layers of the
simulations, we employ the automated algorithm SWIRL.

Results. Small-scale swirls are abundant in the simulated atmospheres of all
the investigated cool stars. The characteristics of these events appear to be
influenced by the main properties of the stellar models and by the strength of
the surface magnetic field. In addition, we identify signatures of torsional
Alfvénic pulses associated with these swirls, which are responsible for a
significant vertical Poynting flux in the simulated stellar photospheres.
Notably, this flux is particularly significant in the K8V model, suggesting a
possible link to the enhanced basal Ca ii H and K fluxes observed in the range
of $B-V$ color index $1.1\leq B-V\leq 1.4$. Finally, we present a simple
analytical model, along with an accompanying scaling relation, to explain a
peculiar result of the statistical analysis that the rotational period of
surface vortices increases with the effective temperature of the stellar
model.

Conclusions. Our study shows that small-scale vortical motions are not unique
to the solar atmosphere and may be relevant in explaining the chromospheric
activity properties of main sequence cool dwarf stars.

###### Key Words.:

Stars: atmospheres – Stars: magnetic field – Magnetohydrodynamics (MHD)

##  1 Introduction

Given the ubiquitous presence of small-scale swirls in the solar atmosphere
(see Tziotziou et al., 2023, for a review), it is reasonable to expect similar
vortical flows in other stellar atmospheres. However, the unique properties of
individual stars, such as the gravitational acceleration at the stellar
surface, the effective temperature, and the mass density, may influence the
properties of these small-scale features.

In addition, magnetic fields play a crucial role in the generation and
dynamics of swirls in the solar atmosphere as was demonstrated, e.g., with
simulations by Canivete Cuissa & Steiner (2020) or Battaglia et al. (2021).
This suggests that the surface magnetic field of other stars are likewise
important for stellar swirls (Wedemeyer et al., 2013). Specifically, magnetic
fields appear to be tightly related to the transport of vortical motion and
Poynting flux from the photosphere to the outer solar atmosphere within so
called magnetic tornadoes, which possibly contribute to the heating of the
chromosphere or corona (Wedemeyer-Böhm et al., 2012).

The precise mechanism driving the upward transport of energy within magnetic
tornadoes remains a topic of debate, with torsional Alfvénic pulses111The term
Alfvénic (instead of Alfvén) is used to remind that these pulses or waves are
not purely but mainly of Alfvénic nature. They may also exhibit a compressive
component. identified as the most probable candidate (Liu et al., 2019;
Battaglia et al., 2021). In this scenario, small-scale swirls are the plasma
counterparts to twist in the field lines of the vertically directed
photospheric magnetic flux concentrations. These perturbations, originating in
the surface layers of the atmosphere, can propagate upward along the flux tube
into the chromosphere as pulses of torsional Alfvén waves. Moreover, numerical
simulations indicate that about $80\,\%$ of the swirls in the photosphere are
consistent with this scenario (Canivete Cuissa & Steiner, 2024).

Of particular interest is the possible contribution of magnetic tornadoes to
the chromospheric basal flux as a function of spectral type (Boro Saikia et
al., 2018). Observations of main-sequence stellar populations reveal that
stars with $1.1\leq B-V\leq 1.4$ exhibit a significantly enhanced
chromospheric basal flux, the cause of which is not yet fully understood.
Therefore, we aim to investigate the presence and properties of small-scale
vortices in the simulated atmospheres of different main sequence stellar
models, and their role in the transport of energy and mass through the
different layers.

For this purpose, three-dimensional radiative-magnetohydrodynamic (R-MHD) box-
in-a-star numerical models have been analyzed, specifically four main-sequence
stellar models of spectral types K8V, K2V, G2V, and F5V. Originally, these
simulations have been performed to study the surface small-scale dynamo (SSD)
in cool main-sequence stars by Riva et al. (2024). Thus, with these models at
hand, it became possible to investigate the properties of vortices in the
presence of a magnetic field that self-consistently emerges through SSD
action. We carry out a statistical analysis of the fundamental properties of
the vortices during both the kinematic and saturation phases of the SSD in
different heights within the simulated atmosphere, and across different
stellar models.

The paper is organized as follows: Sect. 2 describes the simulations and the
vortex identification algorithm used in this study. In Sect. 3, we present
results on the qualitative and quantitative properties of the vortices in the
four stellar models. In addition, we introduce a simple analytical model for
photospheric vortices that replicates one of the results of our statistical
analysis. Finally, Sect. 4 provides a comprehensive summary and the
conclusions.

##  2 Methods

###  2.1 Numerical simulations

The three-dimensional R-MHD simulations that serve as the basis of the present
investigation have been carried out with the CO5BOLD code (Freytag et al.,
2012) in the so called box-in-a-star setup. This corresponds to simulating a
small partial volume near the optical surface of a star, encompassing the
layers where energy transport changes from predominantly convective to purely
radiative. A constant and vertical gravitational acceleration is applied. The
simulations use the long-characteristic radiative transfer module of CO5BOLD
with Rosseland mean opacities and the diffusion approximation in the regime of
high optical depths ($\log\tau_{\rm R}\gtrsim 2$). Details on the numerical
setup can be found in Riva et al. (2024).

Table 1: Parameters describing the four simulated stellar models. See text in
Sect. 2.1 for more details.

Spectral Type | K8V | K2V | G2V | F5V
---|---|---|---|---
$T_{\rm eff}\,[{\rm K}]$ | 4005 | 5000 | 5766 | 6506
$\log{(g)}$ | 4.66 | 4.59 | 4.44 | 4.24
$H_{\rm p}(z=0)\,[{\rm km}]$ | 60 | 100 | 182 | 348
$L_{x},L_{y}\,[{\rm km}]$ | 2200 | 3500 | 6000 | 15000
$L_{z}\,[{\rm km}]$ | 968 | 1540 | 3072 | 7800
Below $z=0$ [km] | 748 | 1190 | 2412 | 6150
Above $z=0$ [km] | 220 | 350 | 660 | 1650
$\Delta x,\,\Delta y,\,\Delta z\,[{\rm km}]$ | 2.93 | 4.67 | 8.00 | 20.00
$L_{\rm gran}\,[{\rm km}]$ | 354 | 609 | 990 | 2490
$N_{\rm kin}$ | 27 | 20 | 18 | 20
$N_{\rm sat}$ | 20 | 32 | 36 | 34
222$T_{\rm eff}$, $L_{\rm gran}$, and $H_{P}$ are evaluated as time averages
during the kinematic phase of the simulations. For more details, the reader
can refer to Riva et al. (2024).

The stellar models are primarily distinguished by their effective temperature,
$T_{\rm eff}$, gravitational acceleration, $g$, and the chemical composition
(metallicity), which here is taken to be the solar one. The main parameters
$g$ and the average $T_{\rm eff}$ of the models are provided in the top set of
rows in Tab. 2 together with the average pressure scale height at the surface
of the model, $H_{\rm p}(z=0)$. Here, surface denotes the space and time
averaged optical depth surface $\tau_{\rm R}=1$, which is also taken as the
origin of the geometrical $z$-scale, $z=0$. The surface gravity $g$ was chosen
to be that of a dwarf star placed on the main sequence.

The middle set of rows of Tab. 2 specifies the horizontal and vertical sizes
of the simulated domains, $L_{x}$, $L_{y}$, and $L_{z}$, along with the grid
spacing, indicated by $\Delta x$, $\Delta y$, and $\Delta z$. The vertical
extension, $L_{z}$, is chosen such that the bottom boundary is at a depth
where the entropy of the downflowing plasma approaches the essentially
constant entropy of the upflowing plasma (see, Salhab et al., 2018, Fig. 11),
whereas the top boundary is located well above the entropy minimum in the
photosphere. The height extents of the computational domains, below and above
the average Rosseland optical depth unity, are also displayed in this section
of the table. Also given is the average size of a granule, $L_{\rm gran}$,
wich is computed as the integral length scale of the vertical velocity
spectrum according to Eq. (5) of Riva et al. (2024). The horizontal extent is
specifically chosen to encompass approximately six granular scales within the
computational domain. Vertically, the simulation domains extend from the near
surface layers of the convective zone to beyond the entropy minimum of the
photosphere, corresponding to approximately the top layers of the photosphere.

Pre-existing and relaxed purely hydrodynamic simulations were seeded with a
vertical magnetic field of $B_{z}=1\,{\rm mG}$. These simulations were then
further evolved and underwent a kinematic phase in which the seed magnetic
field was exponentially amplified by the action of a sub-surface SSD until it
reached saturation. The bottom set of rows indicates the number of analyzed
snapshots during the kinematic and saturated phases, $N_{\rm kin}$ and $N_{\rm
sat}$, respectively. The selection of snapshots for the kinematic and
saturated phases is based on the mean ratio between kinetic and magnetic
energy densities in the convection zone, $\langle E_{\rm kin}/E_{\rm
mag}\rangle_{\rm conv}$. We consider a time instance to be in the kinematic
phase when $\langle E_{\rm kin}/E_{\rm mag}\rangle_{\rm conv}<2.5\cdot
10^{-3}$, while saturation333For the four simulations, the maximum values of
$\langle E_{\rm kin}/E_{\rm mag}\rangle_{\rm conv}$ are around 0.10. For
further details, the reader can refer to Riva et al. (2024). is defined when
$\langle E_{\rm kin}/E_{\rm mag}\rangle_{\rm conv}>2.5\cdot 10^{-2}$.

###  2.2 Identification algorithm

We employ the method proposed by Canivete Cuissa & Steiner (2022, 2024) to
identify swirls in numerical simulations. This method accurately estimates the
coordinates of the center of rotation, that is, the estimated vortex centers
(EVC), for each fluid particle (grid cell) that exhibits some degree of
curvature in the instantaneous velocity field. Since vortices can be viewed as
groups of fluid particles coherently rotating about a common axis of rotation
(Lugt, 1979), clusters of EVCs form in the neighborhood of the core of a
vortex structure. Therefore, the method uses clusters of EVCs to identify
vortices. In this respect, the EVC method can be seen as a more accurate
version of the curvature center method proposed by Sadarjoen & Post (1999),
since it uses the velocity field and its derivatives to extract more detailed
information about the curvature of the flow.

The accuracy of the method hinges on an accurate estimate of the radius of
curvature and radial direction of the local flow. This requires the
measurement of the rigid-body rotational component of the flow. Traditional
methods have used mathematical quantities such as vorticity or swirling
strength (Zhou et al., 1999) for this purpose, but Canivete Cuissa & Steiner
(2022) showed that the optimal quantity is the Rortex criterion $R$ (Tian et
al., 2018; Wang et al., 2019). This advanced mathematical criterion measures
alone the rigid body rotational component of a rotating flow and can be
computed as,

$R=\boldsymbol{\omega}\cdot\boldsymbol{u}_{\rm
r}-\sqrt{\left(\boldsymbol{\omega}\cdot\boldsymbol{u}_{\rm
r}\right)^{2}-\lambda^{2}}\,,$ (1)

where $\boldsymbol{\omega}$ is the vorticity vector, $\boldsymbol{u}_{\rm r}$
is the normalized, real eigenvector of the Jacobian of the velocity vector,
and $\lambda$ is the swirling strength criterion, which formally is twice the
imaginary part of the complex conjugated eigenvalues of the Jacobian of the
velocity vector. For more details on these quantities, the reader can refer to
Canivete Cuissa & Steiner (2020, 2022).

Clusters of EVCs can in principle be identified by eye. However, Canivete
Cuissa & Steiner (2022) adapted the algorithm named clustering by fast search
and find of density peaks (CFSFDP, Rodriguez & Laio, 2014) to automate the
process and reduce human bias in the identification. In addition, the proposed
algorithm includes a cleaning procedure to remove misidentifications caused by
noise or coherent but non-spiraling curved flows.

The algorithm outputs a list of identified vortices, each characterized by its
center coordinates, the region it occupies in a horizontal section (in grid-
cells), and its rotation direction (clockwise or counterclockwise). To
quantify the effective radius of the swirl, $r_{\rm eff}$, the following
formula is employed,

$r_{\rm eff}=\sqrt{\frac{N_{\rm c}}{\pi}}\Delta x\,,$ (2)

where $N_{\rm c}$ represents the number of grid cells occupied by the swirl
and $\Delta x$ is the grid spacing. Additionally, meaningful properties can be
computed by combining this information with the values of physical variables
on the swirl surface. For instance, the effective period of rotation, $P_{\rm
eff}$, can be estimated using,

$P_{\rm eff}=\frac{4\pi}{\langle R_{z}\rangle_{\rm swirl}}\,,$ (3)

where $\langle\,\cdot\,\rangle_{\rm swirl}$ denotes the spatial average over
the swirl surface, and $R_{z}$ represents the vertical component of the Rortex
vector, $\boldsymbol{R}$. Another noteworthy property is the effective plasma
rotational velocity, $v_{\phi,\,{\rm eff}}$, determined by,

$v_{\phi,\,{\rm eff}}=\frac{2\pi r_{\rm eff}}{P_{\rm eff}}\,.$ (4)

A Python implementation of the resulting algorithm is open source on
GitHub444https://github.com/jcanivete/swirl and named SWirl Identification by
Rotation-centers Localization (SWIRL, Canivete Cuissa, 2022). For more details
on the method, the clustering algorithm, and the test cases, the reader can
refer to Canivete Cuissa & Steiner (2022).

##  3 Results and discussion

In this section, we present and discuss the results of our analysis of the
properties of the swirls in the simulated stellar atmospheres. First, we
qualitatively validate the vortex identification algorithm, SWIRL, by
comparing its output with emerging bolometric intensity maps and horizontal
velocity fields. We then perform a comprehensive statistical analysis of the
swirls properties, exploring their tight relationship with magnetic fields and
study their energetics. Finally, we present a simple model of a photospheric
vortex that supports a surprising result revealed by the statistical analysis.

###  3.1 Qualitative validation of SWIRL

Figure 1: Emerging bolometric intensity, $I_{\rm bol}$, for the K8V, K2V, G2V,
and F5V models during the kinematic phase. Clockwise and counterclockwise
vortices identified by the SWIRL algorithm are indicated by pink and green
disks, respectively. The blue squares in the K8V and F5V panels mark the
boundaries of the zoomed plots shown in Fig. 2.

Figure 2: Identified swirls in the areas marked in blue in Fig. 1 for the K8V
and F5V models. Clockwise and counterclockwise vortices identified by the
SWIRL algorithm are indicated by pink and green disks, respectively. The
horizontal component of the velocity field is represented by blue
instantaneous streamlines.

Figure 1 shows the emerging bolometric intensity, $I_{\rm bol}$, at a specific
time during the kinematic phase for the four stellar models. At this time, the
magnetic field remains relatively weak and exerts only a feeble influence on
the flow dynamics. The granulation, characterized by cells of bright (warm)
upflowing plasma surrounded by darker (cooler) downflow in the intergranular
space, is clearly visible in all models. In particular, the granules in the
cooler stars (K8V and K2V) exhibit a smoother surface compared to the more
vigorous dynamics visible in the G2V and F5V models.

From the horizontal velocity field at the surface ($z=0\,{\rm km}$) and at the
corresponding specific times, the SWIRL algorithm identifies the vortices
shown in the four panels of Fig. 1. Each identified vortex is represented by a
circle of size $r_{\rm eff}$, denoting the effective radius of the vortex, and
color-coded according to the direction of rotation (clockwise or
counterclockwise). In Fig. 2, the identified vortices, accompanied by
instantaneous streamlines of the horizontal velocity field, are shown in
detail in two zoomed regions for the K8V and the F5V model.

The majority of the vortices are located within intergranular lanes,
consistent with the results reported for simulations (Canivete Cuissa &
Steiner, 2024) and observations (Bonet et al., 2008) of the solar photosphere
. This behavior is expected due to the conservation of angular momentum in
intergranular downflows (Nordlund, 1985; Wedemeyer & Steiner, 2014). An
examination of the instantaneous streamlines in the zoomed plots of Fig. 2
shows that the SWIRL algorithm consistently identifies vortical patterns in
the velocity field. This is true regardless of whether the flow appears more
laminar-like, as in the K8V model, or more vigorous and highly turbulent, as
in the F5V model.

In Fig. 1 it can be seen that certain swirls, though not all, exhibit
increased radiative intensity compared to their surroundings. An illustration
of such a case is in the left panel of the Fig. 2 located at
$(x,y)=(1.10,0.25)\,{\rm Mm}$. For comparison, in Appendix A we provide plots
of the plasma mass density $\rho$ at $z=0\,{\rm km}$ during the same time
instances as shown in Figs. 1 and 2. We notice that areas of enhanced
bolometric intensity in the integranular space that are associated with
swirling motions exhibit significantly reduced mass densities, as is
particularly well visible in Fig. 10. Partial evacuation and enhanced
radiative intensity in connection with swirling motions are the characteristic
ingredients of so called non-magnetic bright points. Non-magnetic bright
points occurring in magnetic field free simulations of the solar atmosphere
have been investigated by Calvo et al. (2016). Here, they frequently occur
under weak magnetic conditions in the kinematic phase of all stellar models.

In the saturation phase of the SSD, the basic characteristics and distribution
of the swirls in the bolometric intensity maps, and the reliability of their
detection with the SWIRL algorithm carries over from the kinematic phase.
Therefore, we have chosen not to show further similar plots here. The primary
difference lies in the increased prevalence of bright features within
intergranular regions, attributed to magnetic flux concentrations (see, Riva
et al., 2024, Fig. 2), that appear to correlate with the presence of vortices.
This aspect is quantitatively investigated in Sect. 3.2.1.

###  3.2 Statistical Analysis

Table 2: Location of the selected horizontal planes for the statistical
analysis.

Plane | $\ln(p_{0}/p)$ | $z_{K8V}\,[{\rm km}]$ | $z_{K2V}\,[{\rm km}]$ | $z_{G2V}\,[{\rm km}]$ | $z_{F5V}\,[{\rm km}]$
---|---|---|---|---|---
Convection zone | -0.5 | -31 | -57 | -110 | -206
Surface | 0.0 | 0 | -1 | -2 | -6
Middle photosphere | 2.0 | 101 | 153 | 270 | 475
High photosphere | 4.0 | 183 | 284 | 486 | 835
555The values of $z$ for each model correspond to the plane that best
approximates the mean value of $\ln(p_{0}/p)$ relative to the mean $\tau_{\rm
R}=1$ surface.

In this subsection, we examine the characteristics of the vortices identified
by the SWIRL algorithm in both the kinematic and saturated phases. The number
of snapshots analyzed in each regime for each stellar model is detailed in the
bottom set of rows of Tab. 2.

Our analysis refers to four different heights in the simulated atmospheres,
roughly corresponding to the surface layer of the convection zone, the stellar
surface, the middle photosphere, and the high photosphere. The specific
location is determined by the quantity $\ln(p_{0}/p)$, where $p_{0}$ is the
horizontally averaged pressure at the mean optical depth $\tau_{\rm R}=1$
($z\approx 0\,{\rm km}$), and $p$ is the horizontally averaged pressure at the
targeted plane. Table 5 gives the assumed values of $\ln(p_{0}/p)$ along with
the corresponding heights $z$ in the simulation domains of the four models.

####  3.2.1 Surface distributions

Figure 3: Bivariate distribution of the bolometric intensity $I_{\rm bol}$ and
the absolute vertical magnetic field $|B_{z}|$ in the center of the swirls
identified on the surface layer at $z=0\,{\rm km}$ of the four stellar models.
The distributions during the kinematic and saturation phases are shown in blue
and orange, respectively. The vertical black dotted lines indicate the average
bolometric intensity in the intergranular (lower value) and granular (higher
value) regions. The red shaded area indicates the identification range of
magnetic bright feature (MBF) as defined by Riva et al. (2024).

To further investigate the distribution of vortices with respect to the
granular configuration, Fig. 3 shows the bivariate distribution of the
bolometric intensity $I_{\rm bol}$ and the absolute vertical magnetic field
$|B_{z}|$ for all the swirls identified at $z=0\,{\rm km}$. The values of
$I_{\rm bol}$ and $|B_{z}|$ are taken at the center of each swirl. Moreover,
we distinguish swirls identified during the kinematic phase (blue) from swirls
in the saturated phase (orange) of the sub-surface SSD.

The vertical dashed lines in Fig. 3 denote the average bolometric intensity in
the intergranular and granular regions for the four stellar models. These
regions are defined based on the vertical velocity of the plasma, $v_{z}$, at
the surface ($z=0\,{\rm km}$). Specifically, grid cells within the
intergranular and granular regions satisfy the conditions
$v_{z}<-\sigma_{v,z}$ and $v_{z}>\sigma_{v,z}$, respectively, where
$\sigma_{v,z}$ is the standard deviation of the distribution of $v_{z}$ at
$z=0\,{\rm km}$.

During the kinematic phase (blue), the majority of the swirls in all models
are located in regions of low bolometric intensity, corresponding to
intergranular lanes. This finding confirms the visual impression conveyed by
Fig. 1. The high end tails of the four distributions extend into the range of
average granular bolometric intensity, indicating the presence of swirls
associated with non-magnetic bright points (see left panel of Fig. 2) or
within granules. As for the distribution in the vertical magnetic field, it
spans several orders of magnitude ($10^{-4}\,{\rm G}\lesssim|B_{z}|\lesssim
10^{2}\,{\rm G}$) as the magnetic field is amplified by the sub-surface SSD,
starting from an initial value of $1\,{\rm mG}$. Despite the large range, the
distribution appears symmetrical in all four models. We conclude that the
location of the swirls is essentially independent of the magnetic fields
during this phase.

In the saturated phase (red), the distributions in $I_{\rm bol}$ remain
centered around the intergranular bolometric intensities for the four models.
At the same time, the distributions of $|B_{z}|$ become narrower and shift to
higher values ($1\,{\rm G}\lesssim|B_{z}|\lesssim 10^{4}\,{\rm G}$) due to
stronger magnetic fields in the saturation phase.

A notable feature emerges in all panels: the bivariate distributions exhibit a
distinct presence of swirls in regions of high $I_{\rm bol}$ and near
$\max(|B_{z}|)$, particularly evident in the K2V and G2V models. This
observation strongly suggests the presence of swirls associated with magnetic
bright features (MBFs), a phenomenon previously documented in both solar
observations (see, e.g., Bonet et al., 2008) and simulations (see, e.g., Moll
et al., 2011). Notably, the red boxes in the top-right corners of the four
panels indicate the identification range for MBFs, as defined in Riva et al.
(2024, Appendix B). This indicates that, for a sufficiently strong surface
magnetic field, swirls are found in correspondence with MBFs across different
stellar models.

####  3.2.2 Physical properties

Figure 4: Properties of the swirls during the kinematic and saturation phases
of the SSD and at different heights in the atmosphere of the four stellar
models. The average number of swirls per granule in the simulation domain
(top), the average ratio between radius and pressure scale height, $r_{\rm
eff}/H_{\rm p}$ (middle-top), rotational period, $P_{\rm eff}$ (middle-
bottom), and plasma rotational velocity $v_{\phi,\,{\rm eff}}$ (bottom), are
shown. The heights of the cross sections labelled as convection zone, surface,
middle photosphere, and higher photosphere are given in Tab. 5 for the four
stellar models.

Next, we investigate the main physical characteristics of swirls across the
four stellar models, the four heights in the simulation domain, and the two
phases of the sub-surface SSD. The physical properties we are interested in
are: the average number of swirls per granule, the average ratio between the
effective radius and the pressure scale height at $z=0\,{\rm km}$, $r_{\rm
eff}/H_{\rm p}$, the average effective period of rotation $P_{\rm eff}$, and
the average rotational velocity $v_{\phi,\,{\rm eff}}$. The approximate size
of the granules and pressure scale heights at $z=0\,{\rm km}$ for the four
models are given in Tab. 2, while for computing $r_{\rm eff}$, $P_{\rm eff}$,
and $v_{\phi,\,{\rm eff}}$ we employ Eqs. (2), (3), and (4), respectively.

The results of the statistical analysis are shown in Fig. 4. Focusing on the
kinematic phase, the top panel of the figure highlights that the convection
zone has the highest number of vortices per granule in all four models. This
number is computed as the total number of vortices in the corresponding
horizontal plane divided by the area of this plane, $L_{x}\cdot L_{y}$, times
the approximate area of a granule $L_{\rm gran}^{2}$. As we ascend through the
atmosphere, the number of vortices per granule gradually decreases. These
results are consistent with simulations of the solar atmosphere (see Moll et
al., 2012; Canivete Cuissa & Steiner, 2024). The prevalence of highly
turbulent convective flows in the sub-surface layers of cool stars naturally
leads to the formation of vortex motions that explain the observed trend.
Interestingly, within the convection zone, the G2V model stands out as having
the highest number of vortices among the different models.

Conversely, the F5V model has the highest number of vortices in the high
photosphere. This can be attributed to the energetic nature of the flows
observed in the F5V model, which leads us to expect the presence of shocks
already in the photosphere. These shocks could serve as a source of vortex
motions in the upper atmospheres, especially in the virtual absence of
magnetic fields (Moll et al., 2012).

During the saturation phase, there are significant changes in the number of
vortices in the convection zone and in the high photosphere. In the convection
zone, the presence of strong magnetic fields leads to a significant
suppression of vortex formation. This suppression results from the
“stiffening” effect of magnetic fields on the plasma as it approaches
equipartition – a phenomenon known as magnetic quenching (see, e.g.,
Brandenburg & Subramanian, 2005).

In contrast, the number of vortices in the high photosphere increases
significantly, by a factor of about 3 to 4, except for the F5V model where the
increase is less strong. This observation underscores the central role of
magnetic fields in the generation and dynamics of vortex features in the upper
atmosphere. The results are consistent with simulations of the solar
atmosphere performed by Moll et al. (2012) and Battaglia et al. (2021).
Moreover, the similar number of vortices in the high photosphere across models
indicates that magnetic effects outweigh shocks in the formation of vortices
within the magnetized upper photosphere. These results are in agreement with a
similar conclusion presented in Canivete Cuissa & Steiner (2020).

The second panel of the Fig. 4 shows the mean ratio between the effective
radius of the swirl, $r_{\rm eff}$, and the pressure scale height at
$z=0\,{\rm km}$, $H_{\rm p}$. This ratio allows a more meaningful comparison
of results between different models, taking into account the different scales
inherent in each simulation. In stars characterized by larger granules, such
as the F5V model, one would expect the swirls to be larger due to the broader
scales of the flows. Therefore, by normalizing the radius by the surface
pressure scale height – an essential length scale in the models – we obtain a
normalized radius that allows a direct and fair comparison between the
different models.

The results show minimal variation in the mean normalized radius of vortices
across spectral type, model layers, and magnetic field strengths, with values
around $r_{\rm eff}/H_{\rm p}\sim 0.25$. Notably, the only discernible trend
indicates a subtle increase in the photosphere. This result is consistent with
the statistical analysis of the simulated solar atmosphere presented in
Canivete Cuissa & Steiner (2024). The result suggests an intrinsic
relationship between the average size of the swirls and the typical length
scale of the simulated atmosphere. Given that the pressure scale height is
related to the characteristic size of granular flows (see Tab. 2), this
observed relationship seems reasonable.

The bottom two panels of Fig. 4 show the average rotational period, $P_{\rm
eff}$, and the rotational velocity, $v_{\phi,\,{\rm eff}}$, of the swirls. An
interesting observation is the correlation between the rotational velocity and
the effective temperature of the stellar model at all heights. This
correlation is expected, since higher effective temperatures generally
correspond to higher plasma velocities. However, in contrast to this, we also
observe the counterintuitive result that the rotational period also increases.
Obviously, the increasing length scale with increasing effective temperature
compensates for the higher plasma velocity. Predicting which effect dominates
is not trivial. We further explore this intriguing finding with a simple
analytical model in Sect. 3.3.

In the saturation phase of the SSD, a noticeable trend emerges: swirls in the
sub-surface and surface layers of the simulated atmospheres show slower
rotation than in the kinematic phase. This phenomenon can be attributed to the
conversion of turbulent kinetic energy into magnetic energy by the SSD,
resulting in a deceleration of the rotation of vortices generated by turbulent
convective plasma. However, an interesting contrast occurs in the high
photosphere, where vortices tend to rotate faster in the saturation phase of
the SSD.

####  3.2.3 Alfvénic properties

Figure 5: Bivariate distribution of the rotational characteristics of vortices
in the middle photosphere during the saturated phase of the SSD for the four
stellar models. Every identified vortex is represented by a scatter point
according to the Rortex criterion $R_{z}$ and the magnetic Rortex criterion
$R_{z}^{\rm B}$ averaged over its area. The vertical magnetic field $B_{z}$
averaged over the swirl area is color-coded.

Figure 6: Same as Fig. 5 but for the high photosphere.

Sections 3.2.1 and 3.2.2 underscore the fundamental role played by magnetic
fields in shaping the properties of the swirls in the four models of stellar
photospheres. Therefore, the goal of this subsection is to investigate this
relationship in more detail. Of particular interest is the possible
association of swirls in magnetized regions with torsional Alfvénic pulses
(see, e.g., Wedemeyer-Böhm et al., 2012; Liu et al., 2019; Battaglia et al.,
2021). These pulses provide a mechanism for energy transport from the surface
atmospheric layer to the chromosphere and above.

Canivete Cuissa & Steiner (2024) conducted a thorough analysis of the swirl
properties identified in simulations of the solar atmosphere. Their
investigation focused on the imprints of torsional Alfvénic waves in
atmospheric swirls, showing that about $80\,\%$ of photospheric swirls are
consistent with the scenario of a torsional Alfvénic pulse. Next, we reproduce
part of this analysis for the four stellar models.

We begin by considering an incompressible plasma in magneto-hydrostatic
equilibrium within the ideal MHD approximation. Then, for a torsional Alfvén
wave propagating in the vertical direction along a vertically oriented
magnetic field $\boldsymbol{B}=B_{z}\boldsymbol{e_{z}}$, the dispersion
relation of the wave results in,

$\boldsymbol{v}^{\prime}=-\frac{v_{\rm A}}{B_{z}}\boldsymbol{b}^{\prime}\,,$
(5)

where $\boldsymbol{v}^{\prime}$ and $\boldsymbol{b}^{\prime}$ denote plasma
velocity and magnetic field perturbations in the horizontal plane,
respectively, while $v_{\rm A}=|\boldsymbol{B}|/\sqrt{4\pi\rho}\geq 0$ is the
local Alfvén speed. A detailed derivation can be found in Priest (2014, Sect.
4.3.1).

From Eq. (5), we deduce that the perturbations $\boldsymbol{v}^{\prime}$ and
$\boldsymbol{b}^{\prime}$ are either parallel or anti-parallel, depending on
the polarity of the vertical magnetic field, determined by the sign of
$B_{z}$. Furthermore, assuming that the velocity perturbation manifests itself
as a vortex, that is a rotational motion of the plasma in the horizontal
plane, the magnetic field perturbation corresponds to the horizontal component
of a twist in the predominantly vertically oriented magnetic field lines.

To quantify the horizontal velocity perturbations in the plasma,
$\boldsymbol{v}^{\prime}$, we use the vertical component of the Rortex vector,
$R_{z}$. Similarly, we define the magnetic Rortex vector, $\boldsymbol{R}^{\rm
B}$, which measures the equivalent to the rigid-body rotational component for
the torsional component of the magnetic field. The magnetic Rortex vector is
defined as the Rortex vector, Eq. (1), but with the magnetic field
$\boldsymbol{B}$ as the input field, resulting in,

$\boldsymbol{R}^{\rm B}=\boldsymbol{\tilde{J}}\cdot\boldsymbol{u}_{\rm r}^{\rm
B}-\sqrt{\Big{(}\boldsymbol{\tilde{J}}\cdot\boldsymbol{u}_{\rm r}^{\rm
B}\Big{)}^{2}-\Big{(}\lambda^{\rm B}\Big{)}^{2}}\,,$ (6)

where $\boldsymbol{\tilde{J}}=\boldsymbol{\nabla}\times\boldsymbol{B}$ is
proportional to the current density, $\boldsymbol{u}_{\rm r}^{\rm B}$ is the
normalized, real eigenvector of the Jacobian of the magnetic field, and
$\lambda^{\rm B}$ is the magnetic swirling strength (Battaglia et al., 2021).

Using the vertical components of the Rortex and magnetic Rortex vectors,
$R_{z}$ and $R_{z}^{\rm B}$, as proxies to quantify the perturbations
associated with upward propagating torsional Alfvénic pulses, we can express
Eq. (5) as,

${\rm sign}{\Big{(}R_{z}R_{z}^{\rm B}\Big{)}}=-{\rm
sign}{\Big{(}B_{z}\Big{)}}\,.$ (7)

This equation can be used to detect imprints of torsional Alfvénic pulses in
the identified swirls of the stellar models. As we are interested in events
that potentially transport energy from the photosphere into the upper
atmosphere, our analysis targets the lower and upper photospheric layers of
the simulations during the saturation phase of the SSD.

Figure 5 shows the bivariate distribution of the average $R_{z}$ and
$R_{z}^{\rm B}$ for swirls identified in the middle photosphere of the four
stellar models during the saturated phase of the SSD. The averages are
computed over the area of each swirl, and the color indicates the polarity and
strength of the vertical component of the magnetic field in which the swirl is
embedded.

The pattern revealed by the four plots indicates that the vast majority of
swirls obey Eq. (7) in all models, especially in regions where the vertical
magnetic field is sufficiently strong. Considering only swirls with an average
$|B_{z}|>200\,{\rm G}$ over their surface, we find compatibility with Eq. (7)
in percentages of $83.1\,\%,89.3\,\%,92.1\,\%$, and $90.2\,\%$ for the K8V,
K2V, G2V, and F5V models, respectively. Notably, the results for the solar-
like model are in agreement with those of Canivete Cuissa & Steiner (2024),
who performed a similar analysis for a solar atmosphere simulation with an
initial vertical magnetic field of $50\,{\rm G}$.

The same analysis is performed for the high photospheric layer of the four
stellar models, and the results are shown in Fig. 6. The pattern indicating
compliance with Eq. (7) persists in all four panels of the figure, but a
higher level of noise is evident, especially in the K2V model. In this case,
the compatibility with Eq. (7) decreases to $69.1\,\%,58.0\,\%,66.3\,\%$, and
$80.2\,\%$ for the K8V, K2V, G2V, and F5V models, respectively, with a
threshold of $|B_{z}|>50\,{\rm G}$. This trend mirrors the results obtained by
Canivete Cuissa & Steiner (2024) in simulations of the solar atmosphere. The
cause may be a higher proportion of vortices locally generated by shocks that
are not associated with torsional Alfvén pulses, or due to the upper boundary
conditions of the simulations.

####  3.2.4 Energetics

Figure 7: Kernel density estimation (KDE) of the distribution of the vertical component of the Poynting flux vector associated with Alfvénic swirls, $S_{z}^{\rm swirl}$ (green). For each entry into the distribution, the corresponding quantity is the mean over the area of those swirls that are compatible with torsional Alfvénic waves in the middle photosphere of the four stellar models. The KDEs of the horizontal and vertical terms of $S_{z}^{\rm swirl}$, $S_{z}^{\rm h}$ and $S_{z}^{\rm v}$, are shown in orange and blue, respectively. The red vertical dashed lines represent the mean vertical component of the Poynting flux vector computed over the entire domain at the middle photosphere, $S_{z}^{\rm tot}$, while the green vertical dashed lines represent the mean values of $S_{z}^{\rm swirl}$. Table 3: Mean vertical Poynting flux in the middle photosphere computed over the surface of Alfvénic swirls, $S_{z}^{\rm swirl}$, and over the entire simulation domain, $S_{z}^{\rm tot}$. See text in Sect. 3.2.4 for more details. Poynting flux | K8V | K2V | G2V | F5V
---|---|---|---|---
$S_{z}^{\rm swirl}$ [${\rm kW}\,{\rm m}^{-2}$] | 229.6 | 169.1 | 213.3 | 296.2
$S_{z}^{\rm tot}$ [${\rm kW}\,{\rm m}^{-2}$] | 11.5 | 18.4 | 26.9 | 36.5
Table 4: Estimated chromospheric flux and basal flux proxy attributed to
Alfvénic swirls for the four stellar models. The last row provides the
observed lower values of $\log{R^{\prime}_{\rm HK}}$ for stars with similar
effective temperatures. See text in Sect. 3.2.4 for more details.

| K8V | K2V | G2V | F5V
---|---|---|---|---
$F_{\rm swirl}$ [${\rm erg}\,{\rm cm}^{-2}\,{\rm s}^{-1}$] | $3.5\cdot 10^{5}$ | $2.2\cdot 10^{5}$ | $4.7\cdot 10^{5}$ | $5.5\cdot 10^{5}$
$F_{\rm bol}$ [${\rm erg}\,{\rm cm}^{-2}\,{\rm s}^{-1}$] | $1.5\cdot 10^{10}$ | $3.5\cdot 10^{10}$ | $6.3\cdot 10^{10}$ | $1.0\cdot 10^{11}$
$\log{(F_{\rm swirl}/F_{\rm bol})}$ | -4.6 | -5.2 | -5.1 | -5.3
$\log{R^{\prime}_{\rm HK}}$ | -4.8 | -5.2 | -5.2 | -5.0
666The $\log{R^{\prime}_{\rm HK}}$ values are estimated from Boro Saikia et
al. (2018, Fig. 3).

As we observe a strong correlation between swirls and twist of their magnetic
field, indicating compatibility with torsional Alfvénic waves, we investigate
whether these events have the potential to transport energy upward in the
simulated models. To do this, we use the Poynting flux vector, defined as
$\boldsymbol{S}=\boldsymbol{B}\times(\boldsymbol{v}\times\boldsymbol{B})/4\pi$.
Since our focus is on the vertical flux, we consider only the vertical
component of the Poynting flux vector, $S_{z}$. Following Shelyag et al.
(2013), $S_{z}$ can be decomposed as

$S_{z}=\underbrace{\frac{1}{4\pi}v_{z}\Big{(}B_{x}^{2}+B_{y}^{2}\Big{)}}_{S_{z}^{\rm
v}}-\underbrace{\frac{1}{4\pi}B_{z}\Big{(}v_{x}B_{x}+v_{y}B_{y}\Big{)}}_{S_{z}^{\rm
h}},$ (8)

where $S_{z}^{\rm v}$ is a term associated with vertical motions of the
plasma, while $S_{z}^{\rm h}$ is associated with horizontal motions of the
plasma.

In the following, we restrict our analysis to swirls in the middle photosphere
that obey Eq. (7) and are embedded in a sufficiently strong vertical magnetic
field ($|B_{z}|>200\,{\rm G}$ averaged over the swirl area). We will refer to
this selection as Alfvénic swirls.

Figure 7 shows the distributions of the $S_{z}^{\rm v}$ and $S_{z}^{\rm h}$
terms averaged over the area of the Alfvénic swirls in the middle photosphere
of the four stellar models. To better visualize the distributions, a kernel
density estimation (KDE) is provided. We find that for all four models, the
term associated with the vertical motions of the plasma, $S^{\rm v}_{z}$,
(blue curve) peaks at negative values and has a negative skew, while the other
term, $S^{\rm h}_{z}$, (orange curve) is characterized by a positive peak and
skew. This behavior is consistent with results from simulations of the solar
atmosphere by Battaglia et al. (2021).

The sum of the two terms giving the distributions of the vertical component of
the Poynting flux vector associated with Alfvénic swirls in the middle
photosphere of the four models, $S_{z}^{\rm swirl}$, is also shown in Fig. 7
(green curve). These distributions are more centered than the other two, but
retain a slight tendency toward positive values, both in terms of peak and
skewness. As a result, the means of the four distributions have positive
values, as indicated by the green vertical dashed lines in the Fig. 7 and
reported in the first row of Tab. 3.

For comparison, in the second row of Tab. 3 and in Fig. 7 (vertical red dashed
line), we present the average vertical component of the Poynting flux vector
computed over the entire domain in the middle photosphere, $S_{z}^{\rm tot}$.
The mean vertical Poynting flux associated with Alfvénic swirls is
approximately an order of magnitude larger than the mean vertical Poynting
flux in the simulation box at the same height. This implies that Alfvénic
swirls are associated with a substantial upward transport of energy relative
to the surrounding environment. In the following considerations, we assume
that only the net vertical Poynting flux of Alfvénic swirls reaches the
chromosphere, the rest being reflected or dissipated already in the
photosphere.

An important observation is that the mean vertical Poynting flux over the
computational domain increases with the effective temperature of the model,
while the mean vertical Poynting flux associated with Alfvénic swirls is at a
minimum in the K2V stellar model and increases again for model K8V. This
result is particularly interesting in the light of studies of the ratio of the
chromospheric Ca ii H and K flux to the bolometric flux, $\log{R^{\prime}_{\rm
HK}}$, for main-sequence stellar populations (see, e.g., Boro Saikia et al.,
2018, Fig. 3). These studies show that the lower bound of chromospheric
emission, known as the basal flux, is minimal and constant for stars within
the range $0.5\leq B-V\leq 1.1$ ($6000\,{\rm K}\gtrsim T_{\rm eff}\gtrsim
4500\,{\rm K}$) and increases linearly for cooler stars down to $B-V\sim 1.4$
($T_{\rm eff}\sim 4000\,{\rm K}$). The reason for the high basal fluxes in
such cool stars remains unclear.

The contribution of Alfvénic swirls to the vertical Poynting flux over the
total horizontal domain is,

$F_{\rm swirl}=\dfrac{1}{A^{\rm tot}}\sum_{\rm i}{A^{\rm
swirl}_{i}S_{z,\,i}^{\rm swirl}}\,,$ (9)

where we sum over each Alfvénic swirl contribution to the vertical Poynting
flux: $A^{\rm swirl}_{i}=\pi(r_{\rm eff}^{\,i})^{2}$ is the effective area of
the swirl, and $S_{z,\,i}^{\rm swirl}$ is the vertical Poynting flux averaged
over the swirl area. We also divide by the total horizontal area of the
simulation domain, $A^{\rm tot}=L_{x}^{2}$.

For the bolometric flux, we use Stefan-Boltzmann’s law,

$F_{\rm bol}=\sigma T_{\rm eff}^{4}\,,$ (10)

where $\sigma=5.7\cdot 10^{-5}\,{\rm erg}\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm
K}^{-4}$ is the Stefan–Boltzmann constant and $T_{\rm eff}$ is the effective
temperature of the star.

Assuming that a fraction $r$ of $F_{\rm swirl}$ reaches and dissipates in the
low chromosphere where the dominant radiative losses take place in the
emission lines of Ca ii H and K, we can compute its contribution to the
chromospheric Ca ii H and K flux $F_{\rm HK}^{\rm swirl}=rF_{\rm swirl}$. By
taking the logarithm of the ratio between the chromospheric H and K flux due
to the Alfvénic swirls and the bolometric flux, $\log{R^{\prime,\,{\rm
swirl}}_{\rm HK}}=\log{(F^{\rm swirl}_{\rm HK}/F_{\rm bol})}$, we get an
estimate of the basal flux due to Alfvénic swirls (see Linsky et al., 1979).

In Table 6, we present the computed values for $F_{\rm swirl}$, $F_{\rm bol}$,
and the basal flux proxy with $r=1$ across the four stellar models. The final
row includes approximate lower limits of $\log{R^{\prime}_{\rm HK}}$ for stars
characterized by $B-V$ color indices of $1.37$ (K8V), $0.90$ (K2V), $0.65$
(G2V), and $0.47$ (F5V). Notably, the basal flux proxy for the K8V model
surpasses those of the other stellar models. Furthermore, an approximately
$80\%$ fraction of the flux attributed to Alfvénic swirls would suffice to
account for the observed basal flux in stars with $B-V\sim 1.37$. While
acknowledging the hypothetical nature of this result, it introduces a novel
and intriguing explanation for the significant basal fluxes observed in stars
with $1.1\leq B-V\leq 1.4$. Concerning the remaining stellar models, the
contribution of Alfvénic swirls may feasibly account for the basal flux in
stars resembling the $K2V$ and $G2V$ types. However, such a contribution
proves insufficient for the $F5V$ model.

###  3.3 Period - Temperature scaling relation

In Sect. 3.2.2 we found that the rotational period of surface vortices
increases with the effective temperature of stellar models. At first glance,
this result may seem counterintuitive. In order to better understand this
phenomenon, we construct a simplified analytical model of a stellar surface
vortex. For simplicity, we assume a purely hydrostatic model, but a
generalization can be obtained by considering the total pressure $p_{\rm
tot}=p_{\rm gas}+p_{\rm mag}$ instead of just the gas pressure.

Let us consider a vertically oriented, axially symmetric vortex in rotational
hydrostatic equilibrium in the surface layers of a star. We assume a constant
angular velocity $\Omega$ and a uniform plasma density $\rho$. To describe the
dynamics of this system, we use Euler’s momentum equation,

$\rho\frac{D\boldsymbol{v}}{Dt}=-\boldsymbol{\nabla}p+\boldsymbol{F}\,,$ (11)

where $D\boldsymbol{v}/Dt$ represents the material derivative applied to the
velocity field $\boldsymbol{v}$, $p$ is the gas pressure, and $\boldsymbol{F}$
accounts for external forces.

Equilibrium along the radial and vertical directions, $\hat{r}$ and $\hat{z}$,
implies that,

$\displaystyle\frac{\partial p}{\partial r}$
$\displaystyle=\rho\Omega^{2}r\,,$ $\displaystyle\frac{\partial p}{\partial
z}$ $\displaystyle=-\rho g\,.$ (12)

Furthermore, due to axial symmetry, the gas pressure $p$ depends only on the
radius $r$ and height $z$. Consequently, we can formulate the total derivative
using the following differential expressions,

${\rm d}p=\frac{\partial p}{\partial r}\,{\rm d}r+\frac{\partial p}{\partial
z}\,{\rm d}z=\rho\Omega^{2}r\,{\rm d}r-\rho g\,{\rm d}z\,,$ (13)

where we have substituted the partial derivatives with the right-hand sides of
Eq. (12). After integration, we obtain a typical pressure equation for a
rotating fluid in hydrostatic equilibrium,

$p-p_{0}=\frac{1}{2}\rho\Omega^{2}r^{2}-\rho g(z-z_{0})\,,$ (14)

where $p$ and $p_{0}$ are the gas pressure at coordinates $(r,z)$ and
$(0,z_{0})$, respectively. In particular, Eq. (14) indicates that pressure
iso-surfaces in vortex flows exhibit a parabolic shape, with the vortex core
characterized by lower pressure compared to the periphery. This behavior is
illustrated in Fig. 8.

Figure 8: Two-dimensional representation of an axially symmetric flow rotating
at angular velocity $\Omega$. The surface of constant pressure $p=p_{0}$ is
shown, where $p_{0}$ is the pressure at $(r,z)=(0,z_{0})$.

If we consider a pressure isosurface ($p=p_{0}$), we can derive the following
scaling relation,

$H_{\rm p}\sim\frac{\Omega^{2}r_{\rm eff}^{2}}{g}\,,$ (15)

where the pressure scale height, $H_{\rm p}$, is used as the typical vertical
length scale, replacing $z-z_{0}$, and the radius of the vortex, $r_{\rm
eff}$, as the typical radius, $r$. Rearranging Eq. (15) and introducing the
definition $v_{\phi,\,{\rm eff}}=\Omega r_{\rm eff}$, we obtain,

$v_{\phi,\,{\rm eff}}^{2}\sim H_{\rm p}g\sim T_{\rm eff}\,,$ (16)

where we exploit the property of a stratified atmosphere where the pressure
scale height $H_{\rm p}$ is proportional to $T_{\rm eff}/g$. This relationship
allows the model to accurately capture the observed trend shown in the bottom
panel of Fig. 4. Specifically, the model predicts that the rotational velocity
of the vortices increases with the effective temperature of the star.

This scaling relation is consistent with the results of Calvo et al. (2016),
who found a similar relationship by considering the pressure, density, and
temperature contrasts between the vortex core and the surrounding plasma.

Considering the empirical observation shown in the second panel of Fig. 4,
which indicates that $r_{\rm eff}/H_{\rm p}\sim{\rm const}$, we can refine Eq.
(15) as follows,

$\frac{1}{\Omega^{2}}\sim\frac{H_{\rm p}}{g}\,.$ (17)

Finally, using the relations $P_{\rm eff}=2\pi/\Omega$ and $H_{\rm p}\sim
T_{\rm eff}/g$, we can deduce that,

$P_{\rm eff}\sim\sqrt{\frac{T_{\rm eff}}{g^{2}}}\,.$ (18)

This equation establishes a correlation between the period of rotation of the
vortex and key properties of the stellar model, namely $g$ and $T_{\rm eff}$.

To derive a scaling relation for the rotational period that depends only on
the effective temperature, we can use the Stefan-Boltzmann law and rely on
empirical scaling relations established for cool, main-sequence stars
(Demircan & Kahraman, 1991),

$\displaystyle L$ $\displaystyle\sim R^{2}T_{\rm eff}^{4}\,,$ $\displaystyle
L$ $\displaystyle\sim M^{3.7}\,,$ $\displaystyle R$ $\displaystyle\sim
M^{0.9}\,,$ (19)

where $L$, $M$, and $R$ represent the luminosity, mass, and radius of the
star, respectively. When combined with $g\sim M/R^{2}$ at the star’s surface,
these relations lead to the following empirical scaling relation,

$P_{\rm eff}\sim T_{\rm eff}^{2.2}\,.$ (20)

This equation confirms that the rotational period for a simplified vortex
model increases with the effective temperature of the star, which supports the
statistical results obtained from the simulations.

##  4 Summary and conclusions

In this paper, we carried out an analysis on three-dimensional, radiative-MHD,
numerical simulations of the atmospheres of main sequence dwarfs of the
spectral types K8V, K2V, G2V, and F5V. The focus was on the properties of
vortex motions within different layers of the simulation domain. In
particular, we investigated how these properties are affected by the strength
of the magnetic field, which is amplified by the action of a sub-surface SSD.

As expected, vortex motions are found to be ubiquitous in numerical
simulations of stellar atmospheres also other than the solar one. However, we
found that different properties of the vortices depend on the characteristics
of the stellar model and on the strength of the surface magnetic field. The
presence of a strong magnetic field led to a reduction in the number and
rotational speed of vortices in the surface layers of the convection zone, due
to magnetic quenching. However, in the photosphere, an increased number of
vortices showing higher rotational velocities were found in all four stellar
models with respect to the kinematic phase, highlighting the fundamental role
of the magnetic field in the formation and dynamics of small-scale swirls in
stellar atmospheres.

The size of the vortices was observed to be larger in hotter stars. However,
when normalized to the average pressure scale height at the surface, the
average radius of the swirls became nearly constant across the different
stellar models and layers, being about one fourth of the pressure scale
height. This result suggests a direct correlation between the size of the
vortices and the scale of the flows, particularly the granular and
intergranular flows. However, to firmly establish this correlation, a
resolution study is still needed to assess how this result is being affected
by the numerical spatial resolution of the simulations (see, e.g., Yadav et
al., 2020).

We investigate the relation between stellar swirls and perturbations in the
surface magnetic fields generated by saturated SSD action. Our investigation
reveals a robust correlation between swirls and twists in the magnetic field,
suggesting that over $80\,\%$ of the swirls in the middle photosphere are
consistent with torsional Alfvénic waves. These events, especially when
associated with a sufficiently strong vertical magnetic field, contribute to
an average, upward energy flux that significantly exceeds the average Poynting
flux over the entire simulation box. This result highlights that swirls
associated with Alfvénic pulses, termed Alfvénic swirls, have the potential to
play a significant role in the transport of energy to the upper layers of
stellar atmospheres.

Moreover, our investigation reveals that the average vertical Poynting flux
associated with Alfvénic swirls is higher in the K8V model compared to the K2V
model. This observation holds particular significance in the context of
studies on chromospheric activity in main sequence stars. Our proxy for the
basal flux originating from Alfvénic swirls in diverse stellar models aligns
with recent observational data. Specifically, if approximately $80\%$ of the
vertical Poynting flux linked to Alfvénic swirls in the K8V model reaches the
low chromosphere and dissipates, it could potentially explain the heightened
basal flux observed in stars with $B-V\sim 1.4$.

If this vertical Poynting flux is associated with torsional Alfvénic waves
propagating along magnetic flux tubes connecting the photosphere and the low
chromosphere, the fraction reaching the low chromosphere might be close to
unity due to an efficient transport mechanism. Conversely, the majority of the
vertical Poynting flux in the photosphere is a result of random motions, with
most of it being either reflected or dissipating still in the photosphere
(see, Battaglia et al., 2021, Fig. 11). Therefore, we propose that only a
minute fraction of the Poynting flux generated in the photosphere manages to
reach the upper layers. Within this small proportion, Alfvénic swirls likely
play a predominant role. This hypothesis introduces a novel framework for
interpreting the heightened basal fluxes observed in stars within the range
$1.1\leq B-V\leq 1.4$, offering valuable insights into understanding
chromospheric activity in cool stars.

Furthermore, we discovered that the rotational period of swirls increases with
the effective temperature of the stellar model. This finding is not trivial,
considering that hotter stars typically exhibit faster swirl rotations due to
higher flow velocities. Based on a simple model of a vortex in rotational
hydrostatic equilibrium and incorporating empirical mass-luminosity and mass-
radius scaling relations for cool, main-sequence stars, we derived a scaling
relation connecting the rotational period with the effective temperature of
the stellar model, $P_{\rm eff}\propto T_{\rm eff}^{2.2}$. This relation
supports the results obtained from the statistical analysis.

In conclusion, our study demonstrates the existence of swirls in numerical
models of stellar atmospheres, with their specific properties depending on the
properties of the stellar models. In addition, we propose that torsional
Alfvénic pulses associated with small-scale swirling motions may be an
important contributor to the basal flux in Ca ii H and K (S-index) of cool
stars, in particular with respect to the enhanced basal flux in the range
$1.1\leq B-V\leq 1.4$. Although this hypothesis involves some crude
assumptions, it effectively reproduces the trend of the observed fluxes and is
therefore worth considering. To make further progress, simulations of a wider
parameter space and a finer model grid is needed for a more detailed
comparison with measurements of the S-index. Also for these simulations, the
chromospheric layers and corresponding physics should be taken into account,
which would make it possible to compute the Ca ii H and K fluxes directly from
the models and to compute a synthetic S-index and $R^{\prime}_{\rm HK}$
conversion. Including and excluding the areas of Alfvénic swirls would then
provide more insights into the role of these kind of swirls for the
chromospheric basal flux of cool stars. By continuing to study these
phenomena, we can gain a more complete understanding of the intricate
processes that take place in stellar atmospheres and their impact on stellar
systems.

###### Acknowledgements.

The authors acknowledge support by the 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” (www.walsa.com), both at the
International Space Science Institute (ISSI). Part of the numerical
simulations were carried out on Piz Daint at CSCS under project IDs s1059,
s1172, and u14, and the rest was carried out on the HPC ICS cluster at USI.

## References

  * Battaglia et al. (2021) Battaglia, A. F., Canivete Cuissa, J. R., Calvo, F., Bossart, A. A., & Steiner, O. 2021, A&A, 649, A121
  * Bonet et al. (2008) Bonet, J. A., Márquez, I., Sánchez Almeida, J., Cabello, I., & Domingo, V. 2008, ApJ, 687, L131
  * Boro Saikia et al. (2018) Boro Saikia, S., Marvin, C. J., Jeffers, S. V., et al. 2018, A&A, 616, A108
  * Brandenburg & Subramanian (2005) Brandenburg, A. & Subramanian, K. 2005, Phys. Rep, 417, 1
  * Calvo et al. (2016) Calvo, F., Steiner, O., & Freytag, B. 2016, A&A, 596, A43
  * Canivete Cuissa (2022) Canivete Cuissa, J. R. 2022, SWIRL, https://zenodo.org/doi/10.5281/zenodo.10016646
  * Canivete Cuissa & Steiner (2020) Canivete Cuissa, J. R. & Steiner, O. 2020, A&A, 639, A118
  * Canivete Cuissa & Steiner (2022) Canivete Cuissa, J. R. & Steiner, O. 2022, A&A, 668, A118
  * Canivete Cuissa & Steiner (2024) Canivete Cuissa, J. R. & Steiner, O. 2024, A&A, in press
  * Demircan & Kahraman (1991) Demircan, O. & Kahraman, G. 1991, Ap&SS, 181, 313
  * Freytag et al. (2012) Freytag, B., Steffen, M., Ludwig, H. G., et al. 2012, Journal of Computational Physics, 231, 919
  * Linsky et al. (1979) Linsky, J. L., Worden, S. P., McClintock, W., & Robertson, R. M. 1979, ApJS, 41, 47
  * Liu et al. (2019) Liu, J., Nelson, C. J., Snow, B., Wang, Y., & Erdélyi, R. 2019, Nature Communications, 10, 3504
  * Lugt (1979) Lugt, H. J. 1979, The Dilemma of Defining a Vortex, ed. U. Müller, K. G. Roesner, & B. Schmidt (Berlin, Heidelberg: Springer Berlin Heidelberg), 309–321
  * Moll et al. (2011) Moll, R., Cameron, R. H., & Schüssler, M. 2011, A&A, 533, A126
  * Moll et al. (2012) Moll, R., Cameron, R. H., & Schüssler, M. 2012, A&A, 541, A68
  * Nordlund (1985) Nordlund, A. 1985, Sol. Phys., 100, 209
  * Priest (2014) Priest, E. 2014, Magnetohydrodynamics of the Sun
  * Riva et al. (2024) Riva, F., Steiner, O., & Freytag, B. 2024, A&A, in press
  * Rodriguez & Laio (2014) Rodriguez, A. & Laio, A. 2014, Science, 344, 1492
  * Sadarjoen & Post (1999) Sadarjoen, I. A. & Post, F. H. 1999, in Data Visualization ’99, ed. E. Gröller, H. Löffelmann, & W. Ribarsky (Vienna: Springer Vienna), 53–62
  * Salhab et al. (2018) Salhab, R. G., Steiner, O., Berdyugina, S. V., et al. 2018, A&A, 614, A78
  * Shelyag et al. (2013) Shelyag, S., Cally, P. S., Reid, A., & Mathioudakis, M. 2013, ApJ, 776, L4
  * Tian et al. (2018) Tian, S., Gao, Y., Dong, X., & Liu, C. 2018, Journal of Fluid Mechanics, 849, 312
  * Tziotziou et al. (2023) Tziotziou, K., Scullion, E., Shelyag, S., et al. 2023, Space Sci. Rev., 219, 1
  * Wang et al. (2019) Wang, Y.-q., Gao, Y.-s., Liu, J.-m., & Liu, C. 2019, Journal of Hydrodynamics, 31, 464
  * Wedemeyer et al. (2013) Wedemeyer, S., Ludwig, H. G., & Steiner, O. 2013, Astronomische Nachrichten, 334, 137
  * Wedemeyer & Steiner (2014) Wedemeyer, S. & Steiner, O. 2014, PASJ, 66, S10
  * Wedemeyer-Böhm et al. (2012) Wedemeyer-Böhm, S., Scullion, E., Steiner, O., et al. 2012, Nature, 486, 505
  * Yadav et al. (2020) Yadav, N., Cameron, R. H., & Solanki, S. K. 2020, ApJ, 894, L17
  * Zhou et al. (1999) Zhou, J., Adrian, R. J., Balachandar, S., & Kendall, T. M. 1999, Journal of Fluid Mechanics, 387, 353

##  Appendix A Surface distribution of swirls: density

Figure 9: Density $\rho$ at the average $\tau_{\rm R}=1$ surface for the K8V,
K2V, G2V, and F5V models during the kinematic phase. Pink and green disks
correspond to the clockwise and counter-clockwise vortices identified by the
SWIRL algorithm. The blue squares in the K8V and F5V panels represent the
boundaries of the zoom-in plots shown in Fig. 10.

Figure 10: Identified swirls in the areas depicted in blue in Fig. 9 for the
K8V and F5V models. Pink disks represent clockwise swirls, green counter-
clockwise swirls. The horizontal component of the velocity field is
represented by blue instantaneous streamlines.

Figure 9 shows the density $\rho$ at the average $\tau_{\rm R}=1$ surface for
the K8V, K2V, G2V, and F5V models at the same time instances as in Fig. 1. The
figure shows a granulation pattern consistent with that of Fig. 1,
characterized by cells of bright and warm upflowing plasma surrounded by
denser, cooler intergranular downdrafts. Figure 10 zooms in on two regions of
the K8V and F5V models, accompanied by instantaneous streamlines of the
horizontal velocity field. Similar to Figs. 1 and 2, the vortices identified
by the SWIRL algorithm are represented by disks which size and color is
defined by their effective radius and direction of rotation, respectively.

These figures confirm the visual impression conveyed by Fig. 1, indicating
that the majority of vortices are located within or in close proximity to
intergranular flows. It is also apparent that certain bright features observed
at bolometric intensity, commonly referred to as non-magnetic bright spots,
are associated with low density regions and swirling motions. A prominent
example is shown in the left panel of Figs. 2 and 10 at
$(x,y)=(1.10,0.25)\,{\rm Mm}$. These events, previously observed in numerical
simulations of the solar atmosphere by Calvo et al. (2016), are here
identified in simulations of other stellar types as well.