11institutetext: Max-Planck-Institut für Sonnensystemforschung, Justus-von-
Liebig Weg 3, 37077 Göttingen

# Small-scale dynamo in cool stars

III. Changes in the photospheres of F3V to M0V stars

Tanayveer Bhatia, 11<EMAIL_ADDRESS>Robert Cameron 11 Hardi Peter 11 Sami
Solanki 11

(Received / Accepted )

###### Abstract

Context. Some of the quiet solar magnetic flux could be attributed to a small-
scale dynamo (SSD) operating in the convection zone. An SSD operating in cool
main-sequence stars is expected to affect the atmospheric structure, in
particular the convection, and should have observational signatures.

Aims. We aim to investigate the distribution of SSD magnetic fields as well as
their effect on bolometric intensity characteristics, vertical velocity and
spatial distribution of kinetic energy (KE) and magnetic energy (ME) in the
lower photosphere of different spectral types.

Methods. We analyse the SSD and purely hydrodynamic simulations of the near
surface layers of F3V, G2V, K0V and M0V stars. We compare the time-averaged
distributions and power spectra in SSD setups relative to the hydrodynamic
setup. Properties of the individual magnetic fields are also considered.

Results. PDFs of field strength at the $\tau=1$ surface are quite similar for
all cases. The M0V star displays the strongest fields, but relative to the gas
pressure, the fields on the F3V star reach the largest values. All stars
display an excess of horizontal field relative to vertical field in the middle
photosphere, with this excess becoming increasingly prominent towards later
spectral types. These fields result in a decrease in upflow velocities,
slightly smaller granules as well as the formation of bright points in
intergranular lanes. The spatial distribution of KE and ME is also similar for
all cases, implying a simple pressure scale height proportionality of
important scales.

Conclusions. SSD fields have rather similar effects on the photospheres of
cool main-sequence stars, namely, a significant reduction in convective
velocities as well as a slight reduction in granule size, and concentration of
field to kG levels in intergranular lanes associated with the formation of
bright points. The distribution of field strengths and energies is also rather
similar.

###### Key Words.:

Stars: atmospheres – Stars: magnetic fields – Stars: late-type – Convection –
Dynamo

##  1 Introduction

Magnetism in cool stars is ubiquitous. In addition, a significant number of
cool stars show a solar-like activity cycle (Wilson, 1978). The magnetic
fields associated with these cycles are expected to arise from a large-scale
dynamo operating in the interiors of cool stars (Brandenburg & Subramanian,
2005; Charbonneau, 2014). However, there is also an additional, cycle-
independent component of stellar fields, the quiet-star magnetism. From
detailed observations of the quiet Sun, (see, e.g., reviews by Solanki (1993);
de Wijn et al. (2009); Sánchez Almeida & Martínez González (2011); Bellot
Rubio & Orozco Suárez (2019)) as well as state-of-the-art simulations (Vögler
& Schüssler, 2007; Rempel, 2014), this component was realised to be
substantial and could, in part, be explained by invoking a small-scale dynamo
(SSD) mechanism which would amplify magnetic fields via turbulent motions of
the plasma. In fact, recent global SSD simulations (Hotta & Kusano, 2021;
Hotta et al., 2022) showed that the field generated can be significantly
super-equipartition at small scales in the deep convection zone, being strong
enough to affect the meridional circulation and the differential rotation
profile. For stars other than the sun, the influence of SSD fields on quiet
star phenomena like granulation, pressure oscillations and basal chromospheric
activity remain yet to be studied.

Quiet star magnetic fields also have the potential to affect interpretation of
radial velocity (RV) observations of stellar spectra. RV measurements allow
detection of exoplanets by accounting for Doppler shifts in stellar spectral
lines due to a gravitationally-induced ”wobble” caused by a planet’s orbital
motion. However RVs can be affected by stellar magnetism. Starspots can affect
RVs in multiple ways, e.g. via the Evershed effect (Solanki, 2003; Rimmele &
Marino, 2006), while faculae reduce the granular blueshift (Brandt & Solanki,
1990). Granulation, whose presence is not just visible in short time-scale
”noise”, but also as a net blueshift in photospheric spectral lines of most
solar-like stars (Dravins, 1987), could also potentially affect RVs. For a
comprehensive list of factors potentially affecting high precision RV
measurements, see Table A-4 in Crass et al. (2021). Shporer & Brown (2011)
demonstrated the impact convective blueshift can have on RV measurements
during transits at m/s accuracy level via a simple model. With RV measurements
reaching sub-m/s precision, allowing detection of Earth-like rocky exoplanets,
it becomes imperative to understand the sources of stellar ”noise” properly,
including the contribution from magnetic fields.

Stellar light curves also show variation on the order of granulation
timescales. The amplitude of variations over such short timescales are well-
correlated with the stellar surface gravity (Bastien et al., 2013, 2016),
allowing an independent estimation of the latter (as opposed to asteroseismic
measurements). In addition, Shapiro et al. (2017) showed that, for the Sun,
the total solar irradiance (TSI) could be reliably reconstructed just from the
consideration of granulation noise from simulations and solar magnetograms in
a forward model. This is encouraging for modelling stellar variability at
shorter timescales.

In our recent work (Bhatia et al. (2022), hereafter Paper I), we showed that
SSD magnetic fields significantly reduce the convective velocities and can be
strong enough in earlier spectral types (F-star and earlier) to even affect
the stratification and scale heights near the surface (which could influence
scales of granulation). In addition, it was also shown that changing the
metallicity also leads to a difference in SSD-associated field strengths and
properties of momentum transport near the surface (Witzke et al., 2022).
Hence, it becomes imperative to understand the effect of SSD fields on the
photospheres of different stars.

In this paper, we describe the distribution of photospheric quiet-star
magnetic fields as expected to arise from an SSD mechanism operating in the
near-surface convection of late-type dwarfs. We also look at the effects of
this magnetic field on the velocities, the bolometric intensity as well as the
energy distribution in the photospheres of these stars.

##  2 Methods

We use the models described in Paper I, namely, the purely hydrodynamic (HD)
models and the models with SSD fields. The setup, number of snapshots, time
range etc. are similar to those in Paper 1. Briefly, we consider four sets of
a local box-in-a-star simulation of F3V, G2V, K0V and M0V stars, with each set
consisting of a time series of an HD run and an SSD run. The boxes have 512
grid points in both the horizontal directions and 500 grid points in the
vertical direction. The resolution (and the physical size) is such that all
simulations cover a similar number of pressure scale heights
($n_{H_{p}}=\log(p_{\mathrm{gas}}/p_{\mathrm{gas}(\langle\tau\rangle=1)})$)
and horizontally scaled to maintain the aspect ratio. The G2V star is used as
a reference, with $n_{H_{p}}\sim 7.5$ below the surface and $\sim 6-8$ above
the surface (corresponding to 4 Mm below, 1 Mm above). The horizontal extent
is 9 Mm $\times$ 9 Mm 111The K-star SSD and HD setups were rerun for this
paper with a taller atmosphere so as to have a similar $n_{H_{p}}$ above the
surface for all cases. The new time series was used to redo all analyses
mentioned in Paper I and the resulting plots were practically the same.. This
corresponds to a resolution of approximately 17.6(10) km in the horizontal
(vertical) direction. For the F, K and M-stars, the corresponding resolutions
are approximately 45(26), 8.2(4.6) and 4(2.3) km, respectively. The stellar
type is determined by specifying the entropy of inflows (effectively setting
the $T_{\rm eff}$) and by setting a constant gravitational acceleration in the
$z-$direction. We note that the upper boundary condition for magnetic field is
set to be vertical. We simply made a choice between vertical and potential
upper boundary222A test run with potential upper boundary for the G2V star
resulted in more horizontal fields above the surface, and a small increase in
the magnitude of magnetic field (possibly a consequence of more low-lying
loops available for recirculation and amplification). At and below the
$\tau=1$ surface, there was practically no change in the thermodynamic
structure.. We refer the reader to Paper 1 for further details about the
setup.

Before we start describing the results, we provide some information on the
meaning of symbols and conventions we follow. All averages over time are
denoted by an overline $\overline{q}$. All averages over space of 2D data are
denoted by angular brackets $\langle q\rangle$. The standard deviation for
bolometric intensity in a single snapshot is calculated as
$\sigma_{I}=\sqrt{\langle(I-\langle I\rangle)^{2}\rangle}$. The calculation of
spatial power spectra is covered in Appendix A. All plots have error bars
corresponding to standard error $\epsilon=\sigma/\sqrt{N}$, where $N$ is the
number of snapshots and $\sigma$ is the standard deviation over time averages.
The colour coding is the same as that in Paper I, with blue for the F3V, black
for the G2V, green for the K0V, and red for the M0V star. Dashed lines with
lighter corresponding colours refer to the hydrodynamic case, unless stated
otherwise.

##  3 Results

###  3.1 Distribution of the magnetic fields

Figure 1: Magnetic field characteristics. Top: Horizontally averaged magnitude
of magnetic field. Bottom: Horizontally averaged inclination of magnetic
field, defined as $\sqrt{\langle B_{x}^{2}+B_{y}^{2}\rangle/\langle
B_{z}^{2}\rangle}$. The horizontal axis is the number of pressure scale
heights below the surface (positive is below, negative is above).

Figure 2: PDF of the strength of magnetic field $B$ for the SSD F- (blue), G- (black), K- (green) and M-star (red). Top: PDF of $B$ calculated for the geometric surface $z_{\langle\tau\rangle=1}$ corresponding to the height at which $\langle\tau\rangle=1$. Bottom: Same plot as the top, but for the iso-$\tau=1$ surface. The vertical lines correspond to the pressure equipartition field $B_{\rm eqp,gas}=\sqrt{8\pi(p_{\mathrm{gas}})}$ (dotted) and $B_{\rm eqp,tot}=\sqrt{8\pi(p_{\mathrm{gas}}+\rho v_{h}^{2}})$ (dashed), for each case (see Sect. 4.1 for details on the description of $B_{\rm eqp,tot}$). The labels contain the mean value for each case in gauss. Table 1: From left to right: Characteristic of magnetic field near the surface for each stellar type: average field strength ($\overline{B}$), average vertical field strength ($\overline{|B_{z}|}$), average horizontal field strength ($\overline{|B_{h}|}$) and the percent fraction of area occupied by kilogauss fields. ($\overline{A_{B>1\mathrm{kG}}/A_{\mathrm{tot}}}$). All values are averages over time along with 1$\sigma$ standard deviation. Simulation |  $\overline{B}$ [G] |  $\overline{|B_{z}|}$ [G] |  $\overline{B_{h}}$ [G] |  $\overline{A_{B>1\mathrm{kG}}/A_{\mathrm{tot}}}$ [%]
---|---|---|---|---
| $z_{\langle\tau\rangle=1}$ | $\tau=1$ | $z_{\langle\tau\rangle=1}$ | $\tau=1$ | $z_{\langle\tau\rangle=1}$ | $\tau=1$ | $z_{\langle\tau\rangle=1}$ | $\tau=1$
F3V | 130 $\pm$ 10 | 187 $\pm$ 14 | 66 $\pm$ 6 | 94 $\pm$ 8 | 98 $\pm$ 7 | 143 $\pm$ 11 | 0.07 $\pm$ 0.06 | 0.80 $\pm$ 0.31
G2V | 113 $\pm$ 10 | 126 $\pm$ 11 | 59 $\pm$ 5 | 66 $\pm$ 6 | 85 $\pm$ 8 | 95 $\pm$ 9 | 0.46 $\pm$ 0.14 | 0.47 $\pm$ 0.14
K0V | 117 $\pm$ 9 | 118 $\pm$ 9 | 60 $\pm$ 4 | 60 $\pm$ 4 | 89 $\pm$ 7 | 90 $\pm$ 7 | 0.53 $\pm$ 0.12 | 0.50 $\pm$ 0.11
M0V | 106 $\pm$ 6 | 106 $\pm$ 6 | 55 $\pm$ 3 | 55 $\pm$ 3 | 81 $\pm$ 5 | 81 $\pm$ 4 | 0.43 $\pm$ 0.17 | 0.42 $\pm$ 0.16

Fig. 1 shows the horizontally-averaged magnitude (top panel) and horizontally-
averaged inclination (bottom panel) of the magnetic field near the surface.
All cases show similar value of magnetic field strength, except for the
F-star, which shows a somewhat higher value at and above the surface (which is
marked by the dotted vertical line). The field inclination shows the field
becoming more horizontal for all cases as one goes higher up in the
atmosphere. Between
$-2>\log(p_{\mathrm{gas}}/p_{\mathrm{gas}(\langle\tau\rangle=1)})>-4$, all
cases show a peak in $B_{h}/B_{z}$, which probably corresponds to low lying
magnetic field loops. Higher up, $B_{h}/B_{z}$ tends towards zero, in
accordance with the upper boundary condition of a vertical field.

Table 1 lists the magnetic field characteristics for all cases for
$\langle\tau\rangle=1$ horizontal slice as well as for $\tau=1$ iso-
surface333The $\tau=1$ iso-surface refers to the surface where $\tau=1$ for
each vertical column in the 3D cube. The data points are calculated by
interpolating logarithmically against the corresponding $\tau$ column to where
$\tau=1$.. Here, $\tau$ refers to the optical depth corresponding to a
reference wavelength of 500 nm. The $\tau=1$ iso-surface, hence, provides an
observational point of view for understanding the results. We also consider
the horizontal slice because the thermodynamic stratification is expected to
be quite uniform in the horizontal direction, allowing a better understanding
of the physics of the magnetic field distribution. The columns show that for
the horizontal slice, the average field strength is quite similar (100 to 130
G) for all cases, but the value increases significantly for the F-star (almost
190 G) and G-star (almost 130 G) if the $\tau=1$ iso-surface is considered.
For the other cases, the change is almost within the standard deviation and
decreases with $T_{\rm eff}$. In addition, the area fraction of kG fields for
the F-star is around 0.8% for the $\tau=1$ surface, whereas for the
$\langle\tau\rangle=1$ horizontal slice, it is almost zero. We note here that
this result is a consequence of the field strength corresponding to pressure
equipartition (including the $\rho v_{h}^{2}$ term) at this height being just
about one kilogauss (kG); this result is discussed in detail in Sect. 4.1. For
the other cases, the area fraction is roughly 0.5% both ways.

Fig. 2 shows the probability density function (PDF) of the magnitude of the
magnetic field for the horizontal slice at $\langle\tau\rangle=1$ (top panel)
and for the $\tau=1$ iso-surface (bottom panel). All stars show a rather
similar distribution of photospheric magnetic fields $B$, with most of the
field between 20 to 100 G (slightly higher for the F-star) and a more rapid
drop-off (”knee”) in the kG regime. For the horizontal slice, the field
strength for the kG-knee shows an inverse trend with $T_{\mathrm{eff}}$. These
kG fields form mostly in the downflow lanes (see appendix B for separate PDFs
of magnetic field in upflows and downflows). The vertical dotted and dashed
lines mark the field strength corresponding to equipartition with $p_{\rm
gas}$ (that is, $B_{\rm eqp,gas}$) and $p_{\rm gas}+\rho v_{h}^{2}$ (that is,
$B_{\rm eqp,tot}$), respectively, for all the cases. We note that there seems
to be a trend in the location of the kG-knee and the equipartition field
strength. We also note that, if compared against just gas pressure, the F-star
seems to have super-equipartition fields. This tendency towards having super-
equipartition fields relative to $p_{\rm gas}$ decreases towards later
spectral types, with the M-star fields being decidedly sub-equipartition. We
discuss this relation further in Sect. 4.1. Briefly, we expect the strongest
fields to be situated in intergranular lanes, with the strength of the field
being such that it balances against the external pressure.

When one considers the $\tau=1$ iso-surface (bottom panel), the distribution
of fields is roughly similar and does not show the trend for kG fields from
the $\langle\tau\rangle=1$ slice. In addition, the field strengths for the
F-star are generally higher. For the fields concentrated by flux expulsion,
i.e. the fields in equipartition with the flows (around a few 100 G), this is
due to a depression in the $\tau=1$ iso-surface in the intergranular lanes. As
opacity in photospheres of cool stars depends steeply on temperature, the
$\tau=1$ surface generally dips below $z_{\langle\tau\rangle=1}$ in the cooler
intergranular lanes. For the kG fields, however, the dominant mechanism is
evacuation of plasma in the flux tubes leading $\tau=1$ forming deeper below
and radiation escaping from lower levels. This depression in optical depth is
somewhat similar to the Wilson depression (WD) (Wilson & Maskelyne, 1774) of
the optical surface in sunspots. Because of increasing pressure due to
stratification and conservation of flux, the field strength is stronger deeper
down. The magnitude of this depression scales with $T_{\mathrm{eff}}$ as well
as pressure scale height (Beeck et al., 2015), and is the strongest for the
F-star and weakest for the M-star. These results are consistent with
simulations of weak stellar magnetism by Salhab et al. (2018), where the
authors also reported similar values of kG field concentrations for all cases
as well as the scaling of Wilson depression with pressure scale height. The
factors influencing the kG field distribution are discussed later in the
context of convective collapse (Spruit, 1979) as well as pressure balance in
individual flux concentrations in Sect. 4.1.

###  3.2 Bolometric intensity

In Paper 1, we considered the changes in thermodynamic structure due to SSD-
generated magnetic fields. We showed that these changes resulted in a decrease
in density scale height $H_{\rho}$, as well as in convective velocities, near
the stellar surface. Here, we consider how these changes affect the intensity
structure.

Figure 3: Snapshots of the bolometric intensity $I$ (in $10^{10}$ erg/cm2/s)
for the SSD case (first row) and the corresponding vertical magnetic field
$B_{z}$ (in kG) (second row) and the horizontal magnetic field $B_{h}$ (in kG)
(third row) at the iso-$\tau=1$ surface for spectral types (from left to
right) F, G, K and M, respectively. Snapshots of $I$ for the HD case for
comparison (last row). The colour bars are identical for the respective SSD
and the HD case. Green circles indicate magnetic bright points.

Figure 4: Spatial power spectra $P_{k}$ of bolometric intensity $I_{\mathrm{bol}}$ for the SSD (solid) and HD (dashed) cases plotted against spatial frequency $1/x$ (normalized by the horizontal box size $x_{0}$ for each star). Top: $P_{k}$ for all cases normalized by the SSD F-star total power ($\sum_{k}P_{k}$). Bottom: Relative change in power at different scales between the SSD and HD cases. The gray shaded region refer to the approximate scales corresponding to range in granule sizes, for which the centre of gravity is calculated in Table 3. Table 2: Average values for various quantities related to $v_{z}$ and $I_{\mathrm{bol}}$. For each quantity, left column is SSD, right is HD. Here, contrast is defined as $\sigma_{I}/\langle I_{\rm bol}\rangle$. All units are cgs. Simulation | upflow frac. [%] | $\overline{v_{z,\mathrm{rms}}}\,[10^{5}]$ | $\overline{\langle I_{\mathrm{bol}}\rangle}\,[10^{10}]$ | $\overline{\sigma_{I}}\,[10^{10}]$ | contrast [%]
---|---|---|---|---|---
| SSD | HD | SSD | HD | SSD | HD | SSD | HD | SSD | HD
F3V | 57.33 | 57.14 | 5.167 | 5.568 | 4.653 | 4.662 | 1.002 | 1.026 | 21.54 | 22.01
G2V | 56.64 | 56.83 | 2.518 | 2.561 | 2.543 | 2.535 | 0.440 | 0.444 | 17.30 | 17.50
K0V | 55.38 | 55.69 | 1.316 | 1.335 | 1.043 | 1.038 | 0.120 | 0.121 | 11.51 | 11.65
M0V | 60.31 | 60.05 | 0.632 | 0.660 | 0.425 | 0.425 | 0.011 | 0.011 | 2.63 | 2.63

The presence of SSD magnetic fields affects the bolometric intensity
$I_{\mathrm{bol}}$ in multiple ways. The evacuation of plasma due to
concentrated magnetic fields in intergranular lanes leads to formation of
bright points. This is generally attributed to the ”hot-wall” effect (Spruit,
1976), where the low density plasma in the intergranular lanes is heated up by
the surrounding hot, dense upflows, which causes the former to appear bright.
This effect is easily seen in the modelled F, G, K and (to a lesser degree) M
stars (see Fig. 3).

In addition, there are changes in the spatial distribution of
$I_{\mathrm{bol}}$ due to the SSD. The top panel of Fig. 4 shows the magnitude
of the spatial power spectra $P_{k}$ of $I_{\mathrm{bol}}$ for all the HD and
SSD cases. See appendix A for details on how the power spectrum is calculated.
We note that the spatial frequency has been scaled by the box size, which
essentially corresponds to a pressure scale height scaling (see Sect. 2.2 of
Paper 1 for details). This ensures that all plots have the same range on the
$x$-axis. The bottom panel shows the relative change in the power between SSD
and HD cases $(P_{k,\mathrm{SSD}}/P_{k,\mathrm{HD}})-1$, with positive values
corresponding to an increase in power in the SSD case.

The usual interpretation of $P_{k}$, as calculated here, is in terms of level
of contrast at different spatial scales (Nordlund et al., 1997). We take the
peaks of these spectra to be an indication of the granulation scales, as the
largest contrast is expected to exist at the scale of granulation (bright
granule centres vs. dark intergranular lanes). Due to the variation of granule
sizes over a range of scales, we calculate the centre of gravity (CG) over a
range of $k$ (marked by the gray shaded region) to estimate an average spatial
wave number $k_{\mathrm{CG}}$ corresponding to an average granule size. The
peaks of the power spectra lie between $x_{0}/x\in(4,9)$. With that in mind,
we restrict the limits for the CG calculation between $x_{0}/x\in(2,20)$. For
reference, this corresponds to a spatial frequency between 0.22 Mm-1 (or 4.5
Mm in length) and 2.2 Mm-1 (or 0.45 Mm in length) for the G-star, where the
typical granule size is $\sim 1.5-2$ Mm, corresponding to a spatial frequency
of $\sim 0.67-0.5$ Mm-1. The results are presented in Table 3. We interpret
the positive change in $k_{\mathrm{CG}}$ for all cases as an indication for a
slight decrease in average granule size for SSD cases (as $x\sim 1/k$),
relative to the HD cases. The results indicate a positive $\delta k$ (smaller
scales) for granulation in the presence of an SSD compared to the HD case.
Modifying the limits ($k_{1}$ from 2 to 3 and $k_{2}$ from 10 to 20) affects
$\delta k$ but not its sign. The value of $\delta k/k_{HD}$ ranges from +1.1%
to +3.8% for all the stellar types considered (M, K, G and F) for $k_{1}=2$ or
3 and $k_{2}=10$, 15 or 20.

Table 3: Center-of-gravity (CG) for average granulation scale in Fig. 4, calculated as $(\int_{k_{1}}^{k_{2}}P_{k}kdk)/(\int_{k_{1}}^{k_{2}}P_{k}dk)$, where $k_{1},k_{2}$ are the left and right bounds of the shaded region. Simulation | $k_{\mathrm{CG,SSD}}[x_{0}/x]$ | $k_{\mathrm{CG,HD}}[x_{0}/x]$ |  $\delta k/k_{\mathrm{HD}}$ (%)
---|---|---|---
F3V | 10.5 | 10.08 | +4.09
G2V | 8.65 | 8.41 | +2.77
K0V | 7.57 | 7.30 | +3.68
M0V | 9.36 | 9.15 | +2.24

Figure 5: PDF of the bolometric intensity for the SSD (solid) and HD (dashed)
with vertical axis in semilog.

In Fig. 5, we plot the PDF of the bolometric intensity for all the SSD as well
as HD cases. The F, G and K-star show a clearly bimodal distribution, with the
bright peak corresponding to mean granular intensities and downflows
corresponding to mean intergranular lanes intensities. This is consistent with
distributions obtained from a variety of other simulations of stellar
photospheres (Beeck et al., 2012; Magic et al., 2013; Beeck et al., 2013a;
Salhab et al., 2018). The M-star also has a somewhat bimodal distribution, but
it is not as prominent as the other cases, since the contrast between granules
and intergranular lanes is very low. All SSD cases show an extended bright
tail (right side of the PDFs) which corresponds to the formation of magnetic
bright points in intergranular lanes. In addition, the G and K-stars show
slight excess intensities in the bright flank (above the bright tail,
corresponding to bright granules) for the SSD cases, whereas F-star shows a
slight decrease. The F and G-stars also show a steeper fall-off at the dark
flank (left side of the PDFs), which contributes a reduction in contrast for
the SSD cases, as noted in Table 2.

The effect of SSD fields on intensity at sub-granular scales is more varied
between spectral types. Fig. 4, bottom panel, shows a prominent increase in
power for the K-star at the smallest spatial scales, corresponding to the
high-contrast bright points present in the intergranular lanes. This is
consistent with the bright tail in Fig. 5. A similar interpretation also holds
for the G-star. Visually, we also see bright points in the F-star case (see
first column in Fig. 3). However, these bright points do not lead to an
increase in power. This is because magnetic fields also restrict convective
velocities, acting as an effective ”viscosity”, which makes the flow more
laminar and leads to granules with a smoother appearance, thus decreasing the
amount of small-scale sub-structure and contrast. We caution against
interpreting the results in the $k>100$ range, since here the effects of
numerical diffusion are non-negligible.

###  3.3 Vertical velocity

In Paper 1, we showed that there is a general reduction in vertical velocity
$v_{z}$ for SSD cases near the photosphere. Here we examine in more detail how
the distribution of $v_{z}$ changes relative to the HD cases when accounting
for magnetic field generation.

Figure 6: Vertical velocity $v_{z}$ at the iso-$\tau=1$ surface normalized by
the r.m.s vertical velocity in the HD case $(v_{z_{r}ms})_{\mathrm{HD}}$ (as
noted in Table 2). Top: PDF of $v_{z}/(v_{z_{r}ms})_{\mathrm{HD}}$ for the SSD
(solid) as well as HD (dashed) cases. Bottom: The difference between the PDF
for SSD and HD cases $v_{z,\mathrm{SSD}}-v_{z,\mathrm{HD}}$. Positive
(negative) values on horizontal axis correspond to upflows (downflows).

The top panel of Fig. 6 shows the PDF of $v_{z}$ at the $\tau=1$ iso-surface,
normalized by $(v_{z,\mathrm{rms}})_{\mathrm{HD}}$. This normalization allows
us to compare the shapes of the PDF between the different stars and to examine
the changes between the SSD and HD cases (Fig. 6, bottom panel). First of all,
we note that all cases show a similar PDF, with a sharp high peak for upflows
and a broad low peak for downflows. There are a couple of exceptions to the
general trend: the downflow peak for the M-star is lower than for the others.
This might partly be a consequence of a somewhat higher upflow fraction for
the M-star ($\sim$60%) compared to the other stars ($\sim$55-57%). Another
difference is the upflow peak for the F-star, which is offset to relatively
smaller velocities compared to the other stellar types. This is consistent
with a thicker tail for the high upflow velocities $v_{z}$, and probably
reflects the larger spread in $v_{z}$ for the F-star. Notwithstanding these
small differences, the distribution of velocities is remarkably similar for
all the four stars.

In the presence of SSD magnetic fields, we see that there is a decrease in the
mean upflow velocities. This is possibly a consequence of reduced kinetic
energy in upflows due to presence of SSD fields, with the fields acting
effectively as enhanced viscosity. The mean downflow velocities remain
relatively unchanged. This may be related to the fact that magnetic fields in
downflows are close to vertical, allowing downflowing plasma to remain
relatively unhindered, whereas the magnetic field above granules is largely
horizontal, which impedes upflowing plasma (Schüssler & Vögler, 2008; Rempel,
2014).

###  3.4 Spatial distribution of energy

Figure 7: Spatial power spectrum of magnetic and kinetic energies at the
iso-$\tau=1$ surface. Top: Power spectrum of magnetic energy for all the SSD
cases. Middle: Power spectrum of the kinetic energy for all the SSD (solid,
dark) and HD (dashed, light) cases. Bottom: Percent change in the kinetic
energy power spectrum for SSD cases, relative to HD cases. The top and middle
plots are normalized by the total kinetic energy for the F-star (solid blue).

The spatial power spectrum plot for the magnetic energy (ME) in Fig. 7 (top
panel) shows a fairly similar distribution for the G, K and M-star, whereas
for the F-star, the spectrum is slightly steeper at the larger scales and has
higher power than the spectra for other stars at all wavenumbers. The power
spectra for kinetic energy (KE) (middle panel) are also very similar for all
the stars at smaller wavenumbers and roughly similar for the larger
wavenumbers. In fact, the relative changes in the KE power spectra between the
SSD and the HD cases (bottom panel) are remarkably similar for all the models,
with a decrease in energy at the largest scales (smallest wavenumbers) and the
smallest scales (largest wave numbers). On the other hand, there is no
significant change in the power at scales roughly corresponding to granule
sizes (see Sect. 3.2 for details on granulation scales). Since the dimensions
of all the stars are scaled to have a similar number of vertical pressure
scale heights (and the horizontal size is scaled accordingly to maintain
aspect ratio), the similarity in all the power spectra point to a simple
pressure scale height scaling of the relevant dynamics.

##  4 Discussion

###  4.1 kG magnetic fields

Figure 8: Horizontal force balance in a vertical slice of a magnetic field
concentration for the F-star. Left: Plot of bolometric intensity with the
selected magnetic element highlighted by the green rectangle. Red contours
represent areas with $B^{2}/8\pi>p_{\mathrm{gas}}$. Right: Map of $B_{z}$ in
the $x$-plane with horizontal axis in Mm and vertical axis in number of
pressure scale heights above (negative) and below (positive) the
$\langle\tau\rangle=1$ height. The dark red, red and yellow contours show the
iso-$\log_{10}\tau=$1,0,-1 surfaces, respectively. The line plot on top shows
the diagonal terms of the total stress tensor (refer to Eq. 1 and Sec. 4.1 for
details). Animation available. Figure 9: Same as Fig. 8, but for K-star.
Animation available.

As mentioned in the introduction, the distribution of small-scale magnetic
field strengths is well-studied for the solar case. The weak (sub-kG) field
distribution is explained in terms of an equipartition between KE and ME: weak
turbulent magnetic fields get carried upward and outwards from granule centres
and get collected in intergranular lanes up to kinetic energy equipartition in
a process called flux expulsion (Weiss, 1966). However, the field strength
corresponding to KE equipartition is substantially sub-kG. To explain the
presence of kG fields, the convective collapse mechanism (Parker, 1978; Spruit
& Zweibel, 1979; Spruit, 1979) has usually been invoked. A qualitative picture
of the mechanism is as follows: magnetic field gets accumulated in downflow
lanes due to flux expulsion from granule centres. This field can intensify up
to equipartition with kinetic energy. At this point, a nascent flux tube
forms. While plasma in the tube keeps flowing downwards, the flow of plasma
into the lanes is restricted by the Lorentz force, leading to the tube getting
evacuated. Finally, to maintain horizontal force balance, the tube gets
compressed, causing the magnetic field to amplify, potentially up to pressure
equipartition ($p_{\mathrm{mag}}=B^{2}/8\pi\approx p_{\mathrm{gas}}$). The
efficiency of this mechanism can be thought of in terms of the minimum plasma
$\beta=p_{\mathrm{gas}}/p_{\mathrm{mag}}$ (or the strongest fields) for which
this instability can work in an idealized flux tube. Rajaguru et al. (2002)
showed that this minimum $\beta$ shows an increasing trend with
$T_{\mathrm{eff}}$ and decreasing $g$, implying that convective collapse is
more efficient for hotter stars.

However, this mechanism can only explain field strength up to pressure
equipartition, whereas our simulations reveal presence of super-equipartition
fields for the F and G-stars. The reason for this is that convective collapse
is an idealized model that assumes flux concentrations exist as stable thin
vertical flux tubes and only accounts for the importance of magnetic pressure
and gravity. Realistic simulations of flux concentrations have shown that the
assumptions are not quite satisfied, especially near the surface where
$\beta\sim 1$ (Yelles Chaouche et al., 2009).

Figure 10: Average over time of the line plots described in Fig. 8 (top) and
Fig. 9 (bottom), respectively. Error bars are 1$\sigma$ standard error.

Some idealized convection studies have indeed shown to result in magnetic
field concentrations well above $B_{\mathrm{eqp}}$, for example, in Bushby et
al. (2008), where they cite dynamic pressure as a major factor leading to
super equipartition fields. This becomes especially important for hotter stars
like our F-star, where velocities in the photosphere are, on average, sonic,
and the turbulent pressure is non-negligible compared to gas pressure. To
understand which factors are important for the formation and evolution of flux
concentrations we consider the momentum balance equation in the following
form,

$\displaystyle\frac{d(\rho v_{i})}{dt}$
$\displaystyle=-\nabla_{j}\left(p_{\mathrm{gas}}\delta_{ij}+\rho
v_{i}v_{j}+\frac{1}{4\pi}\left(\frac{B^{2}}{2}\delta_{ij}-B_{i}B_{j}\right)\right)+\rho
g_{i}$ (1)

Here, the indices $i,j$ correspond to the $x,y,z$ directions and
$g_{i}=(0,0,-g_{\mathrm{surf}})$. The terms on the right hand side within the
parentheses are the gas pressure and, when averaged over time, the Reynolds
and Maxwell stresses, respectively. The Maxwell stresses themselves can be
decomposed into an isotropic magnetic pressure term and a general magnetic
tension term, respectively.

We consider individual field concentrations for the F-star (Fig. 8) and the
K-star (Fig. 9). We note that Eq. 1 must always be satisfied at any given
point of time. For horizontal force balance, however, we only consider the
diagonal terms ($i=j$) of the total stress tensor (first bracket on the right
hand side) for simplicity. For a cut along the $y$-axis, these terms (gas
pressure $p_{\rm gas}$, magnetic pressure $B^{2}/8\pi$, $\rho v_{y}^{2}$ and
$B_{y}^{2}/4\pi$, plus the sum of all the components) are plotted for the
highlighted field concentration in the left panel for both Fig. 8 and 9 at a
given instant.

If there was perfect horizontal pressure balance, the thick blue line in Fig.
8 and 9 would be flat. That is clearly not the case, however, even for the
K-star. This means we are possibly missing contributions from the time-
dependent $d(\rho v_{i})/dt$ term, as well as the off-diagonal components
($i\neq j$ in Eq. 1) of the stress tensor. To account for the former term, at
least, we average over the lifetime of the magnetic field concentrations for
both the F- and K-star. The corresponding averages are plotted in Fig. 10.
Luckily, the horizontal balance is reasonably well-maintained (that is, the
thick blue line is relatively flat) with just the diagonal components for the
F-star (top) as well as the K-star (bottom). The difference between the two
cases is the strength of the magnetic pressure relative to the gas pressure:
in the F-star, $B^{2}/8\pi>p_{\rm gas}$ in the centre of the field
concentration, whereas for the K-star, $B^{2}/8\pi\sim p_{\rm gas}/3$. The
reason for this difference, in this case, is the extra contribution from the
$\rho v_{y}^{2}$ term, which is roughly $p_{\rm gas}/4$ outside of the field
concentration. We note that the value of $B_{\rm eqp,tot}$ is around 1 kG for
the F-star for the $z_{\langle\tau\rangle=1}$ slice, which explains why the kG
fraction in this case is essentially zero.

The degree of evacuation of a flux tube is dependent on the level of
superadiabaticity (Spruit & Zweibel, 1979), with hotter stars having higher
superadiabaticity near the surface (Rajaguru et al., 2002; Beeck et al.,
2013a). This explains the trend of kG fields in G, K and M-stars relative to
pressure-equipartition field strength. However, based on the analysis above,
for the F-star, an additional contribution from the diagonal component of the
Reynolds stresses (here, $\rho v_{y}^{2}$) must be accounted for at the very
least. When this is done, the equipartition field strength rises to $\sim 1$
kG and the fraction of field strength above this drops to zero essentially, as
noted in section 3.1. Another thing to note here is the depth of the
iso-$\tau=1$ surface in terms of $n_{H_{p}}$: For the F-star, the dip in the
surface in the intergranular lane can be $>1.5n_{H_{p}}$, whereas for the
K-star, this dip is much less than $1.5n_{H_{p}}$. This explains why the PDFs
for the iso-$\tau=1$ surface look rather similar for all cases: the horizontal
force balance is maintained relatively deeper down for hotter stars (we ”see”
stronger fields), whereas for cooler stars, the depth difference is not so
significant.

###  4.2 Granulation and intensity distribution

All models show slight changes in the apparent granulation with the inclusion
of SSD fields. In Paper I, we showed that the inclusion of SSD fields results
in a reduction in the ratio of horizontal to vertical velocities
$v_{h,\mathrm{rms}}/v_{z,\mathrm{rms}}$ as well as in the density scale height
$H_{\rho}$. Above the surface, $H_{\rho}$ reduces for the SSD F-star, but
stays largely the same for the SSD G, K and M-star, as the decrease in
turbulent pressure is roughly compensated by an increase in magnetic pressure.
However, $v_{h,\mathrm{rms}}/v_{z,\mathrm{rms}}$ decreases for all cases.
Based on simple momentum conservation arguments, the diameters of convection
cells at any depth is given by $D\approx 4(v_{h}/v_{z})H_{\rho}$ (Nordlund et
al., 2009), so one would expect the granule size to decrease accordingly. To
check that, we base the average granule size on the peak of the spatial power
spectrum $P_{k}$ of $I_{\mathrm{bol}}$. Previous simulations have shown a
tight correlation between granule diameters derived from this relation and
$P_{k}$ (Magic et al., 2013). A decrease in granule size would imply a shift
to smaller spatial scales (higher spatial frequency) of the peak. This is, in
fact, the case for all the stars, as shown in Table 3. This is also supported
by observational indications of the relation between magnetic field and
granule size for the Sun. Various studies of variation in granule size within
an active region (Title et al., 1992; Narayan & Scharmer, 2010) and over the
solar activity cycle (Ballot et al., 2021), show a general inverse correlation
between the granule size and the magnetic field strength.

Figure 11: Mean values of normalized intensity $\hat{I}=(I-\overline{\langle
I\rangle_{\mathrm{HD}}})/\overline{\sigma_{I_{\mathrm{HD}}}}$ for each bin of
$B_{z}$ (top) and $B_{h}$ (bottom) at the iso-$\tau=1$ surface, plotted over
the respective joint histograms

At the smaller scales, the features are more varied between spectral types. As
mentioned before, all cases exhibit magnetic bright points in intergranular
lanes. We can more clearly see the relation between magnetic field
concentrations and intensity by considering the joint histogram of $B_{z}$ and
normalized intensity $\hat{I}=(I-\overline{\langle
I\rangle_{\mathrm{HD}}})/\overline{\sigma_{I_{\mathrm{HD}}}}$ and plotting the
mean values of $\hat{I}$ values for each $B_{z}$ bin (Fig. 11, top panel).
Here we see that F, G and K-star have quite similar mean $\hat{I}$ values,
with a correlation between bright features and strong $B_{z}$ values,
qualitatively similar to observational findings for the Sun (Kahil et al.,
2017). The trend is also similar for the M-star, but the correlation is not as
strong. The fact that the plots for F, G and K-star basically overlap is a
consequence of the normalization of intensity by $\sigma_{I}$. In addition,
the horizontal fields $B_{h}$ show no such correlation with intensity (Fig.
11, bottom panel), firmly connecting bright points to vertical field
concentrations, confirming the result found for the Sun by, e.g., Riethmüller
& Solanki (2017).

The presence of magnetic bright points is expected to enhance power in the
$I_{\mathrm{bol}}$ spatial power spectra at smaller scales, but this is
clearly the case only for G and K stars. This can be understood in terms of
where the $\tau=1$ surface is formed with respect to where the energy transfer
shifts from convective to radiative (Nordlund & Dravins, 1990) for different
stellar types. As discussed in Sect. 4.2 of Paper I, for the F and G-star, the
$\tau=1$ layer forms below where most of overturning of plasma takes place,
leading to naked granules whereas for the K and M-star, it forms above,
leading to hidden granules (see also Sect. 3.2 of Beeck et al. (2013b) for a
more comprehensive description of what constitutes a naked vs. hidden
granule). Especially for the F-star, the turbulent structures within granules
are seen very clearly in the HD case. However, with SSD fields, this turbulent
appearance smooths out significantly as the magnetic field ends up acting like
an effective viscosity hindering the flow. This affects not only the intensity
distribution but also the overall radiative flux: For the F-star, there is a
slight decrease in the bolometric intensity (see column 4 of Table 2), i.e.
the decrease in intensity due to the magnetic field’s effect on convection
dominates over the enhancement due to bright points. This effect is not so
strong for the G and K star, since the increase in contrast due to bright
points dominates. In addition, there is a slight increase in bolometric
intensity. For the M-star, there is practically no change as bright points are
relatively infrequent and their contrast is relatively low as well.

###  4.3 Energy distribution and vertical velocities

Our simulations show remarkably similar spatial distributions of not just KE
and ME, but also the change in KE between SSD and HD simulations. Since all
boxes are scaled by number of pressure scale heights, this implies a simple
pressure scale height proportionality between important scales in the
considered stellar atmospheres. The presence of SSD fields results in
reduction of KE at sub-granular scales as well as the scale of the whole box.
Energy for the magnetic fields is extracted from the KE reservoir at small
(sub-granular) scales, which leads to fields with strength near kinetic energy
equipartition, and this ME cascades to the largest scales, resulting in a net
reduction of average KE. Global simulations of SSD also show a similar
reduction in KE at the largest scale (Hotta et al., 2022), pointing to
possibly significant global consequences in stellar convection with SSD
fields.

An observational signature associated with this reduction in convective
velocities is change in the convective blueshift of photospheric lines
(Dravins et al., 1981). Convective blueshift is one of the few measurable
quantities that encodes information about stellar granulation. A reduction in
$v_{z}$ would imply shifts in line profiles used for estimating convective
blueshift. Since these simulations are gray, we refrain from calculating line
bisectors etc., as these are likely to be inaccurate. We plan to carry on this
analysis in a subsequent paper in this series.

The reduction of KE at smallest scales of around $40\%$ (see Fig. 7 bottom
panel near $k{\sim}10^{2}$) reflects the near-equipartition division of energy
between KE and ME at the scales where field amplification takes place. The SSD
mechanism and the corresponding scale-dependent transfer of energy between the
kinetic and magnetic energy reservoirs has been extensively studied and was
shown to be quite universal (Moll et al., 2011).

Figure 12: Difference between the SSD and HD cases for the mean normalized
intensity $\hat{I}$ as a function of normalized vertical velocity
$v_{z}/(v_{z,rms})_{\mathrm{HD}}$ at the iso-$\tau=1$ surface. Positive values
on the vertical axis mean a higher mean value for the intensity from the SSD
case at the corresponding velocity bin. Positive values on the horizontal axis
represent upflows.

Another way to study the effect of SSD fields on intensity is to consider the
joint PDFs of $\hat{I}$ and $v_{z}/(v_{z,rms})_{\mathrm{HD}}$ and examine the
difference between the mean $\hat{I}$ for SSD and HD cases. In brief, we
compute the mean value of intensity in every velocity bin for the SSD and HD
case and plot the difference in Fig. 12. We see that, for the F, G and K-star
SSD cases, $\hat{I}$ is enhanced in downflows, as would be expected from
bright points forming in downflows. Interestingly, there is a slight increase
in $\hat{I}$ for upflows as well for all cases (except the F-star), which
could mean that granules are in general a little brighter for SSD cases.

##  5 Conclusion

The presence of SSD fields in our simulations affects the photosphere in a
rather similar manner for the cool-star spectral types considered here: fields
are amplified due to turbulent plasma motions at sub-granular spatial scales.
These fields then get collected in intergranular lanes, where they get
concentrated to kG levels while (roughly) maintaining horizontal force
balance. For the F-star, the horizontal force balance is satisfied only after
inclusion of Reynolds stresses, which leads to magnetic pressure being higher
than gas pressure in the strongest flux concentrations. Magnetic bright points
are also occasionally visible in the downflow lanes, with a clear correlation
of $B_{z}$ with $\hat{I}$, as well as an increase in the bright tail of
intensity distribution. There is also an overall slight decrease in granule
size. Because of the SSD fields, the upflow velocities also decrease, again
with a similar signature in PDFs of $v_{z}/(v_{z,rms})_{\mathrm{HD}}$ for all
cases. This decrease in upflow velocities signals a possible reduction in
expected convective blueshift. In summary, an SSD acting in stellar
photospheres is expected to have an effect on spectral line shifts, limb
darkening as well as stellar variability at short time scales. We plan to
investigate these possibilities in subsequent papers in this series.

###### Acknowledgements.

We thank the anonymous referee for their comments and careful consideration,
which significantly helped improve the presentation of this paper. TB is
grateful for access to the supercomputer Cobra at Max Planck Computing and
Data Facility (MPCDF), on which all the simulations were carried out. This
project has received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (grant
agreement No. 695075).

## References

  * Ballot et al. (2021) Ballot, J., Roudier, T., Malherbe, J. M., & Frank, Z. 2021, arXiv e-prints, arXiv:2106.03556
  * Bastien et al. (2013) Bastien, F. A., Stassun, K. G., Basri, G., & Pepper, J. 2013, Nature, 500, 427
  * Bastien et al. (2016) Bastien, F. A., Stassun, K. G., Basri, G., & Pepper, J. 2016, ApJ, 818, 43
  * Beeck et al. (2013a) Beeck, B., Cameron, R. H., Reiners, A., & Schüssler, M. 2013a, A&A, 558, A48
  * Beeck et al. (2013b) Beeck, B., Cameron, R. H., Reiners, A., & Schüssler, M. 2013b, A&A, 558, A49
  * Beeck et al. (2012) Beeck, B., Collet, R., Steffen, M., et al. 2012, A&A, 539, A121
  * Beeck et al. (2015) Beeck, B., Schüssler, M., Cameron, R. H., & Reiners, A. 2015, A&A, 581, A42
  * Bellot Rubio & Orozco Suárez (2019) Bellot Rubio, L. & Orozco Suárez, D. 2019, Living Reviews in Solar Physics, 16, 1
  * Bhatia et al. (2022) Bhatia, T. S., Cameron, R. H., Solanki, S. K., et al. 2022, A&A, 663, A166
  * Brandenburg & Subramanian (2005) Brandenburg, A. & Subramanian, K. 2005, Phys. Rep, 417, 1
  * Brandt & Solanki (1990) Brandt, P. N. & Solanki, S. K. 1990, A&A, 231, 221
  * Bushby et al. (2008) Bushby, P. J., Houghton, S. M., Proctor, M. R. E., & Weiss, N. O. 2008, MNRAS, 387, 698
  * Charbonneau (2014) Charbonneau, P. 2014, ARA&A, 52, 251
  * Crass et al. (2021) Crass, J., Gaudi, B. S., Leifer, S., et al. 2021, arXiv e-prints, arXiv:2107.14291
  * de Wijn et al. (2009) de Wijn, A. G., Stenflo, J. O., Solanki, S. K., & Tsuneta, S. 2009, Space Sci. Rev., 144, 275
  * Dravins (1987) Dravins, D. 1987, A&A, 172, 200
  * Dravins et al. (1981) Dravins, D., Lindegren, L., & Nordlund, A. 1981, A&A, 96, 345
  * Hotta & Kusano (2021) Hotta, H. & Kusano, K. 2021, Nature Astronomy, 5, 1100
  * Hotta et al. (2022) Hotta, H., Kusano, K., & Shimada, R. 2022, ApJ, 933, 199
  * Kahil et al. (2017) Kahil, F., Riethmüller, T. L., & Solanki, S. K. 2017, ApJS, 229, 12
  * Magic et al. (2013) Magic, Z., Collet, R., Asplund, M., et al. 2013, A&A, 557, A26
  * Moll et al. (2011) Moll, R., Pietarila Graham, J., Pratt, J., et al. 2011, ApJ, 736, 36
  * Narayan & Scharmer (2010) Narayan, G. & Scharmer, G. B. 2010, A&A, 524, A3
  * Nordlund & Dravins (1990) Nordlund, A. & Dravins, D. 1990, A&A, 228, 155
  * Nordlund et al. (1997) Nordlund, A., Spruit, H. C., Ludwig, H. G., & Trampedach, R. 1997, A&A, 328, 229
  * Nordlund et al. (2009) Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Living Reviews in Solar Physics, 6, 2
  * Parker (1978) Parker, E. N. 1978, ApJ, 221, 368
  * Rajaguru et al. (2002) Rajaguru, S. P., Kurucz, R. L., & Hasan, S. S. 2002, ApJ, 565, L101
  * Rempel (2014) Rempel, M. 2014, ApJ, 789, 132
  * Riethmüller & Solanki (2017) Riethmüller, T. L. & Solanki, S. K. 2017, A&A, 598, A123
  * Rimmele & Marino (2006) Rimmele, T. & Marino, J. 2006, ApJ, 646, 593
  * Salhab et al. (2018) Salhab, R. G., Steiner, O., Berdyugina, S. V., et al. 2018, A&A, 614, A78
  * Sánchez Almeida & Martínez González (2011) Sánchez Almeida, J. & Martínez González, M. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 437, Solar Polarization 6, ed. J. R. Kuhn, D. M. Harrington, H. Lin, S. V. Berdyugina, J. Trujillo-Bueno, S. L. Keil, & T. Rimmele, 451
  * Schüssler & Vögler (2008) Schüssler, M. & Vögler, A. 2008, A&A, 481, L5
  * Shapiro et al. (2017) Shapiro, A. I., Solanki, S. K., Krivova, N. A., et al. 2017, Nature Astronomy, 1, 612
  * Shporer & Brown (2011) Shporer, A. & Brown, T. 2011, ApJ, 733, 30
  * Solanki (1993) Solanki, S. K. 1993, Space Sci. Rev., 63, 1
  * Solanki (2003) Solanki, S. K. 2003, A&A Rev., 11, 153
  * Spruit (1976) Spruit, H. C. 1976, Sol. Phys., 50, 269
  * Spruit (1979) Spruit, H. C. 1979, Sol. Phys., 61, 363
  * Spruit & Zweibel (1979) Spruit, H. C. & Zweibel, E. G. 1979, Sol. Phys., 62, 15
  * Title et al. (1992) Title, A. M., Topka, K. P., Tarbell, T. D., et al. 1992, ApJ, 393, 782
  * Vögler & Schüssler (2007) Vögler, A. & Schüssler, M. 2007, A&A, 465, L43
  * Weiss (1966) Weiss, N. O. 1966, Proceedings of the Royal Society of London Series A, 293, 310
  * Wilson & Maskelyne (1774) Wilson, A. & Maskelyne, N. 1774, Philosophical Transactions of the Royal Society of London Series I, 64, 1
  * Wilson (1978) Wilson, O. C. 1978, ApJ, 226, 379
  * Witzke et al. (2022) Witzke, V., Duehnen, H. B., Shapiro, A. I., et al. 2022, arXiv e-prints, arXiv:2211.02722
  * Yelles Chaouche et al. (2009) Yelles Chaouche, L., Solanki, S. K., & Schüssler, M. 2009, A&A, 504, 595

##  Appendix A Spatial Power spectrum

To compute the spatial power spectrum, we use the following procedure:

  1. 1.

For a given 2D quantity $q$, take the 2D FFT (with the zero-frequency mode
shifted to the center, e.g., with numpy.fft.fftshift in python using numpy)
giving $\tilde{q}$.

  2. 2.

Multiply $\tilde{q}$ with its complex conjugate $\tilde{q}^{*}$ to get
$P_{k_{x},k_{y}}=\tilde{q}\tilde{q}^{*}$.

  3. 3.

For each radial wave number $k=\sqrt{k_{x}^{2}+k_{y}^{2}}$, construct a
1-pixel wide mask.

  4. 4.

Take the mean of all the data in each mask and multiply by the radius of the
mask to get power $P_{k}$.

For bolometric intensity, we use $q=\sqrt{I}$. For kinetic energy,
$q=\sqrt{\rho/2}v$. For magnetic energy, $q=B/\sqrt{8\pi}$.

##  Appendix B Additional plots

Figure 13: PDF of the magnitude of magnetic field in upflows, for the
geometric surface $z_{\langle\tau\rangle=1}$ (top) and for the $\tau=1$ iso-
surface (bottom).

Figure 14: PDF of the magnitude of magnetic field in downflows, for the
geometric surface $z_{\langle\tau\rangle=1}$ (top) and for the $\tau=1$ iso-
surface (bottom).


THE DEFLECTION OF LIGHT INDUCED BY THE SUN’S GRAVITATIONAL FIELD AND MEASURED
WITH GEODETIC VLBI

O.A. TITOV1, A.A. GIRDIUK2
1 Geoscience Australia
PO Box 378, Canberra, ACT 2601, Australia
e-mail<EMAIL_ADDRESS>
2 Institute of Applied Astronomy, Russian Academy of Sciences
Kutuzov Quay 10, 191187, Saint-Petersburg, Russia
e-mail<EMAIL_ADDRESS>

ABSTRACT.

The Sun’s gravitational field deflects the apparent positions of close objects
in accordance with the formulae of general relativity. Optical astrometry is
used to test the prediction, but only with the stars close to the Sun and only
during total Solar eclipses. Geodetic Very Long Baseline Interferometry (VLBI)
is capable of measuring the deflection of the light from distant radio sources
anytime and across the whole sky. We show that the effect of light deflection
is equivalent to the gravitational delay calculated during the reduction of
VLBI data. All reference radio sources display an annual circular motion with
the magnitude proportional to their ecliptic latitude. In particular, radio
sources near the ecliptic pole draw an annual circle with magnitude of 4 mas.
This effect could be easily measured with the current precision of the
geodetic VLBI data.

1\. INTRODUCTION

Very Long Baseline Interferometry (VLBI) is capable of measuring precise group
delays - the difference in arrival times of radio waves at two radio
telescopes (Schuh & Behrend, 2012) from distant extragalction radio sources
(quasars). Accurate positions of these radio sources are obtained with an
accuracy of 40 microarcsec (Ma et al., 2009).

Gravitational time delay caused by the Solar gravitational field is calculated
during the reduction of geodetic VLBI data (Shapiro, 1964, 1967). The
conventional formula for calculating gravitational delay is formulated in
terms of the positions of the radio telescopes within the barycentric
reference frame of the Solar System (Kopeikin, 1990; Eubanks et al., 1991;
Soffel et al., 1991; Klioner, 1991), rather than the baseline length between
the radio telescopes.

We propose an alternate gravitational delay formula using a Taylor series
expansion. We show that the conventional formula can be split into a sum of
several terms, and the major term links the gravitational delay and the well-
known formula for the light deflection angle. The light deflection angle can
be considered equivalent for all baselines and estimated for each radio source
at times of interest. We develop a new approach to probe the formula for the
light deflection angle at an arbitrary elongation from the Sun. Finally,
estimates of the light deflection for several reference radio sources based on
VLBI observations in 1991-2001 are presented.

2\. GRAVITATIONAL DELAY VS LIGHT DEFLECTION

The light deflection angle $\alpha$ at an arbitrary elongation $\theta$ from
the Sun, is given by

$\alpha=\frac{(\gamma+1)GM}{c^{2}r}\frac{\sin\theta}{1-\cos\theta}$ (1)

where $G$ is the gravitational constant, $M$ is the mass of a gravitational
body, $c$ is the speed of light, $r$ is the distance from the Earth to the
Sun, $\gamma$ is the parameter of the Parametrised Post-Newtonian formalism
(PPN), $\theta$ is the elongation angle (i.e. the angular distance between the
Sun and the gravitational body, as measured by the Earth observer) (Shapiro,
1967; Ward, 1970). The conventional gravitational delay is calculated as
follows

$\tau_{grav}=\frac{(\gamma+1)GM}{c^{3}}\ln\frac{|\vec{r_{1}}|+(\vec{s}\cdot\vec{r_{1}})}{|\vec{r_{2}}|+(\vec{s}\cdot\vec{r_{2}})},$
(2)

The VLBI total delay model also comprises a term due to the transformation
from the barycentric to the geocentric reference frames.

$\tau_{coord}=\frac{(\gamma+1)GM}{c^{2}r}\frac{(\vec{b}\cdot\vec{s})}{c}$ (3)

It can be shown that the three formulae are approximately linked as (Titov &
Girdiuk, 2015)

$\tau_{GR}\approx\tau_{grav}+\tau_{coord}=\alpha\frac{b}{c}\sin\varphi\cos A$
(4)

where angles $\varphi,\theta,A$ are linked by the standard spherical triangle
formula

$\begin{array}[]{c}\cos\psi=-\cos\varphi\cos\theta-\sin\varphi\sin\theta\cos
A,\end{array}$ (5)

Converting the three dot products with angles $\varphi,\;\psi,\;\theta$ (Fig
1), where

$(\vec{b}\cdot\vec{s})=|\vec{b}|\cos\varphi,\;\;\;(\vec{b}\cdot\vec{r_{2}})=|\vec{b}||\vec{r_{2}}|\cos\psi,\;\;\;(\vec{r_{2}}\cdot\vec{s})=-|\vec{r_{2}}|\cos\theta$
(6)

Figure 1: Angles $\varphi$, $\psi$, $\theta$ and $A$, originated by position
of gravitational mass (B), quasar (Q) and baseline vector ($\vec{b}$). If the
Sun plays the role of the gravitational mass, then, the point B on Fig 1 is
also the position of the Solar system barycentre.

More accurately, the total contribution of the post-Newtonian (PN) general
relativity effect to the total VLBI delay is given by (we recall now
$\gamma=1$ for the case of general relativity)

$\begin{array}[]{ll}\tau_{GR}\approx&\frac{2GMb}{r_{2}c^{3}}\frac{\sin\varphi\sin\theta\cos
A}{1-\cos\theta}+\frac{GM}{c^{3}}\frac{b^{2}(1-\cos^{2}\varphi\cos^{2}\theta)}{r_{2}^{2}(1-\cos\theta)}-\frac{GM}{c^{3}}\frac{b^{2}\sin^{2}\varphi\sin^{2}\theta\cos^{2}A}{r_{2}^{2}(1-\cos\theta)^{2}}\end{array}$
(7)

The two additional terms appear because the elongation angles $\theta_{1}$ and
$\theta_{2}$ between the direction to the observed radio source and the centre
of the gravitational body, as measure from each radio telescope, are not
equivalent. The parallactic effect for a baseline of 6,000 km is about 8
arcsec, and the deflection angles $\alpha_{1}$ and $\alpha_{2}$ for the
”first” and the ”second” radio telescope differing. At a large elongation the
additional terms are negligible, so in the small angle approximation this
effect should be considered.

$\begin{array}[]{ll}\tau_{GR}&=\frac{4GM}{c^{2}R_{2}}\frac{b}{c}\sin\varphi(\cos
A-\frac{b}{2R_{2}}\sin\varphi\cos
2A)=\alpha_{2}\cdot\frac{b}{c}\sin\varphi(\cos
A-\frac{b}{2R_{2}}\sin\varphi\cos 2A)\end{array}$ (8)

where $\alpha_{2}$ is the classical deflection angle for light propagated
through a gravitation field (Einstein, 1916) for the second station
$\alpha_{2}=\frac{4GM}{c^{2}R_{2}}$ and $R_{2}=r_{2}\cos\theta_{2}$. In an
alternative form

$\tau_{GR}=(\alpha_{2}-\frac{GM}{c^{2}R_{2}}\frac{b}{R_{2}}\frac{\sin\varphi\cos
2A}{\cos A})\cdot\frac{b}{c}\sin\varphi\cos
A=(\alpha_{2}+\alpha^{\prime}_{2})\cdot\frac{b}{c}\sin\varphi\cos A$ (9)

where $\alpha^{\prime}_{2}$ is the additional deflection angle
$\alpha^{\prime}_{2}=\frac{GM}{c^{2}R_{2}}\frac{b}{R_{2}}\frac{\sin\varphi\cos
2A}{\cos A}$

3\. OBSERVATIONAL DATA

We processed the available VLBI data between 1991 and 2001 using the OCCAM
software (Titov et al., 2004). A fraction of well-established reference radio
sources were assigned as astrometrically unstable, i.e. their coordinates were
not fixed. The observational group delays were approximated by the theoretical
values, and the O–C (observed minus calculated) differences were used to
obtain a set of estimates of the standard daily parameters (Earth orientation,
station positions, wet troposphere delays and gradients, etc.) and corrections
to the selected radio source coordinates (right ascension and declination)
regardless of its elongation from the Sun. The consensus gravitational delay
(3) and the coordinate term (4) were not applied for the calibration of the
group delay for the radio sources. The ionosphere fluctuations were calibrated
in the conventional way, and, as it will be shown later, this makes possible
VLBI observations of radio sources in the range of 1∘.5 to 3∘ from the Sun.

Figure 2: The variations of the coordinates of the radio sources 1606+106
(left) and 0229+131 (right) in 1991-2001. A smaller minimum angle $\theta$ for
0229+131 results in larger size of the circle.

We used the Least Squares collocation method to estimate the wet troposphere
fluctuations with the OCCAM software. The mutual correlation between
observables are introduced in the data adjustment process. The difference
between the VLBI estimates of the wet troposphere delays and independent
radiometer data is typically within 3-6 mm, or 10-20 ps (Titov & Stanford,
2013). Thus, the impact of the wet troposphere delay on the astrometric light
deflection angle estimated near the Sun is negligible.

Figure 3: Variations of the daily estimates of the PPN parameter $\gamma$ for
radio source 0229+131 in 1991-2001 at the different scale of angle $\theta$.

Fig 2 shows the variations of the coordinates of two selected radio source
over 10 years (1991-2001). Each point on the plot corresponds to a daily
position estimated from a single 24-hour session. About 10 annual circles are
drawn by each radio source on the plots. The catalogue (non-deflected)
positions are at the reference origin, and the total deflection angle is equal
to distance to the point from the reference origin. The ecliptic latitude of
the four sources are different, therefore, elongation angles vary over
different ranges. As a result, the deflection angle does not change
essentially for the former radio sources (the reference origin is in the
middle of the ring) with a magnitude about 4 mas. In contrast, for the latter
case, the deflection angle varies from zero to 0′′.2 (for the minimum
elongation of 3∘), and the reference origin lies near the edge of the ring.

The variations of the deflection angle are easily converted with formula (1)
to the variations of the PPN parameter $\gamma$. Fig 3 shows the variations of
the deflection angle as a function of $\theta$ for the reference radio source
0229+131 at different scales. This radio source is almost on the ecliptic
plane (ecliptic latitude $\epsilon\approx$ -1∘.5), and the elongation angle is
limited by the technical capability of observing near the Sun. On 29 April,
1997 (session code 97APR29XE) the mean elongation angle was 1∘.55, and the
measured deflection reached a maximum among all selected radio sources,
$\alpha$ = 0′′.3100 +/- 0′′.0009. However, this estimate was obtained from
only 7 observations. The best accuracy was achieved on 1 May, 1996 (session
code 96MAY01XO) with 122 observations. The deflection angle estimate is
$\alpha$ = 0′′.1588 +/- 0′′.0002, and the relative accuracy is 0.0012, in
spite of a larger mean elongation angle for this date.

4\. CONCLUSION

The effects of general relativity are explicitly contained in both components
of the total VLBI delay model – the gravitational and geometric delays. While
the former component uses the individual barycentre positions of the radio
telescopes to calculate the effect, the latter component is expressed in terms
of the baseline between the radio telescopes. Coupling between both parts has
not been investigated until now.

The Shapiro effect as measured by radars, and the deflection of light measured
with traditional astronomical instruments are considered two independent tests
of general relativity. In this paper we show that the total group delay model
joins the two tests within one observational technique - geodetic VLBI. The
gravitational delay that traditionally originated from the Shapiro effect is
linked to the light deflection angle. Therefore, the two approaches, VLBI
delay and angular, are absolutely equivalent.

For almost all realistic situations this angle does not depend on the baseline
length, thus, a standard geodetic VLBI interferometer acts as a traditional
astronomical instrument. In addition, the coordinate term explicitly presented
in the conventional geometric delay model ceases to exist because it is
compensated by the same effect in the gravitational delay with the opposite
sign. Thus, the proposed alternative version of the general relativity
contribution to the total VLBI group delay model is free from coordinate
effects. Therefore, the two approaches, time delay and angular deflection, are
absolutely equivalent.

The final equation of the general relativity contribution also comprises two
smaller terms which are significant at very small angular separation between
the deflecting body and distant radio source. These terms may be considered as
an increment in the light deflection angle due to the additional time delay
for propagation of the light from station 1 to station 2. This effect become
significant at $\frac{b}{R}>0.1$.

acknowledgements

This paper is published with the permission of the CEO, Geoscience Australia.

5\. REFERENCES

Einstein, A., 1916, Annalen der Physik, 354, 769.

Eubanks, T.M., Carter, M.S., Josties, F.J., Matsakis, D.N., McCarthy, D.D.,
1991, IAU Colloq. 127: Reference Systems, 256.

Klioner, S.A., 1991, In: Proc. of AGU Chapman Conference on Geodetic VLBI:
Monitoring Global Change, Washington DC, 188.

Kopeikin, S.M., 1990, Sov. Astron., 34, 5.

Ma, C., Arias, E.F., Bianco, G., et al., 2009, IERS Technical Note, 35.

Schuh, H., Behrend, D., 2012, J. Geodyn., 61, 68.

Shapiro, I.I., 1964, Phys. Rev. Lett., 13, 789.

Shapiro, I.I., 1967, Science, 157, 806.

Soffel, M.H., Wu, X., Xu, C., Mueller, J., 1991, AJ, 101, 2306.

Titov, O., Girdiuk, A. 2015, A&A, 574, AA128

Titov, O., Stanford, L. 2013, 21st Meeting of EVGA, Reports of the Finnish
Geodetic Institute, p. 151-154

Titov, O., Tesmer, V., Boehm, J. 2004, IVS for Geodesy and Astrometry 2004 GM
Proceedings, 267

Ward, W.R., 1970, ApJ, 162, 345.


[1em]1em1em

# Generative Pre-Trained Transformer for Cardiac Abnormality Detection

Pierre Louis Gaudilliere1 Halla Sigurthorsdottir1 Clémentine Aguet1
Jérôme Van Zaen1 Mathieu Lemay1 Ricard Delgado-Gonzalo1

1 Swiss Center for Electronics and Microtechnology (CSEM SA)  Neuchâtel
Switzerland

###### Abstract

ECG heartbeat classification plays a vital role in diagnosis of cardiac
arrhythmia. The goal of the Physionet/CinC 2021 challenge was to accurately
classify clinical diagnosis based on 12, 6, 4, 3 or 2-lead ECG recordings in
order to aid doctors in the diagnoses of different heart conditions.
Transformers have had great success in the field of natural language
processing in the past years. Our team, CinCSEM, proposes to draw the parallel
between text and periodic time series signals by viewing the repeated period
as “words” and the whole signal as a sequence of such words. In this way, the
attention mechanisms of the transformers can be applied to periodic time
series signals. In our implementation, we follow the Transformer Encoder
architecture, which combines several encoder layers followed by a dense layer
with linear or sigmoid activation for generative pre-training or
classification, respectively. The use case presented here is multi-label
classification of heartbeat abnormalities of ECG recordings shared by the
challenge. Our best entry, not exceeding the challenge’s hardware limitations,
achieved a score of 0.12, 0.07, 0.10, 0.10 and 0.07 on 12-lead, 6-lead,
4-lead, 3-lead and 2-lead test set respectively. Unfortunately, our team was
unable to be ranked because of a missing pre-print.

##  1 Introduction

In the recent years transformers [1] have shown to perform strongly in the
field of natural language processing. By solving the problem of long-term
dependencies and effectively paying attention to the right words for the
context, very impressive transformer-based language generation and translation
models have been created including Google’s BERT [2] and OpenAI’s GPT-3 [3].

A sentence, or a piece of text, is a series of words that represents a finite
amount of repeated patterns with a semantic meaning. Similarly, certain time
series signals such as electrocardiogram (ECG signals) can be thought of as
series of a finite amount of repeated patterns (heartbeats). Here we exploit
this similarity by modeling single heartbeats as words and full ECG signals as
sentences.

The contribution of the present work is two-fold:

  * •

Drawing the conceptual parallel between periodic time series signals and text

  * •

Using that to create a simple, yet powerful, representation of ECG signals,
representing each heartbeat as a “word embedding” and an ECG signal as a
sequence of such “words”

The paper is organized as follows. Section 2 exposes the background. In
Section 3, we describe the pipeline from raw ECG signals to multi-label
classification, as well as the model architecture. Then, in Section 4, we
expose the obtained results. Finally, in Section 5, we give a critical review
of the implementation and propose further lines of research.

##  2 Previous Work

In the recent years, transformers have been used on time-series based tasks
[4, 5]. Regarding ECG classification, Yan et al. [6] got good results on the
MIT-BIH database, i.e. a per heartbeat classification task. Furthermore, [7]
used a convolutional neural network (CNN)+transformer architecture alongside
hand-crafted features to win the CinC2020 challenge [8], where each recording
(duration between 5 s and 30 min) was weakly labelled (a multi-label
classification problem). However, none of these works drew the parallel
between their time-series signals and text. Furthermore, none of the ECG-based
works used a heartbeat as a word embedding, nor used their transformer on a
sequence of periodic “words”.

Figure 1: Overview of the complete system: we split the raw ECG signal into a
sequence of heartbeats that is fed into a transformer. The last fully
connected layer takes the output of the transformer to produce a prediction of
the next heartbeat or to output probabilities for multi-label classification
of cardiac arrhythmia.

##  3 Methods

In this work, we present a deep transformer network designed for multi-label
classification of $d_{classes}=28$ cardiac arrhythmia including atrial
fibrillation, atrial flutter or premature atrial contraction. Our network uses
sequences of heartbeats embedded as “words”. The latter are fed into a
Transformer architecture [1] that relies entirely on a parallelizable self-
attention mechanism. We first pre-trained the model so that it learns about
features of the ECG waveforms and to make the model more generalizable by
adding external data
sources111https://physionet.org/content/edb/1.0.0/222https://physionet.org/content/mitdb/1.0.0/.
During pre-training, the model receives as input a sequence of single
heartbeats and it outputs a guess of the next heartbeat. Therefore, we fed the
output of the transformer into a fully connected layer with $d_{model}=1000$
output neurons and linear activation. Finally, for the classification task
(supervised learning), we used transfer learning and replaced the last layer
with a fully connected layer with 28 output probabilities and sigmoid
activation. Figure 1 gives an overview of the complete system.

###  3.1 Dataset

The database combines eight worldwide datasets, for a total of 88k ECG
recordings shared for training [9]. we left aside recordings from St
Petersburg since they are 30 min long and more susceptible to artifacts such
as movement. Moreover, we restricted the total number of heartbeats in a
sequence to be $maximum\\_position\\_encoding=50$. We truncated longer
sequences, whereas we zero padded shorter sequences at the end for the missing
heartbeats. We separated the remaining recordings into five folds, stratified
by labels and datasets. We assigned four of the folds to training while we
preserved the fifth for validation. To cope for possible inconsistencies in
labeling between datasets, we trained on different combination of datasets.
All the challenge datasets were included for validation. Some labels (SA,
TInv) are underrepresented which make it harder for the model to learn a good
representation. We combined the labels that were considered equivalent by the
challenge. Finally, we did not assign a label to recordings that have only
non-scored labels.

###  3.2 Signal Pre-Processing

To be able to properly view the signal as a series of repeated heartbeats, we
need to perform some pre-processing steps. First, we scaled the ECG leads by
their respective ADC gain. Then, we applied an IIR (infinite impulse response)
high-pass filter with a cutoff frequency of 0.5 Hz. We then tested four R-peak
detectors333https://github.com/berndporr/py-ecg-detectors (Christov, Hamilton,
Two-Average and Pan Tompkins detector) in order to separate the raw ECG
signals into heartbeats. Finally, to uniformize the sampling frequency between
the datasets, we re-sampled each recording to 500 Hz, and we updated the
timing of the R-peaks accordingly.

###  3.3 Word Embeddings

The first step in a transformer model is to create word embeddings. These
embeddings represent the word in an embedding space. Here, we tested the most
simplistic way of representing a heartbeat as a “word embedding”. We
considered 1/3 of the R-R interval before the R-peak and 2/3 of the R-R
interval after the R-peak to be a part of the heartbeat. We then aligned the
heartbeats at the R-peaks. To create embeddings of the same length, we defined
a maximum normal length of a heartbeat to be $d_{model}=1000$ samples, or 2 s.
We zero padded any heartbeat below this length and we truncated any heartbeat
above this length. Therefore, one can infer the variability of the R-R
interval by looking at the number of added zeros: the shorter the interval the
more we padded with zeros. Finally, we computed the root mean square (RMS) of
the active ECG leads to obtain a single ECG signal (Figure 2).

Figure 2: We aligned the heartbeats at the R-peaks and we zero padded on the
left and right for a total length of 1000 samples. Then, we computed the RMS
of the active leads.

###  3.4 Transformer Model

As transformers rely uniquely on self-attention, we use positional encoding
[1] in order to infer the relative timing of the heartbeats within a sequence.
We added together the embedding representations $(x_{0},...,x_{n})$, where
$x_{i}\in R^{d_{model}}$, with the positional encodings $(p_{0},...,p_{n})$
and the result is given as input to the transformer.

Our transformer uses a stack of $N=5$ encoders. The input of each encoder is a
masked multi-head attention layer followed by a feed forward network which
combines two dense layers with $\mathit{dff}$ and $d_{model}$ output neurons
respectively. For regularization, we applied a dropout at the output of each
layer before it is added to the input of the layer and normalized.

The masked multi-head attention layer at the beginning of the transformer is
composed of $h=8$ heads. Each embedding representation $x_{i}$ is passed to a
dense linear layer of $d_{model}$ units to create a query, key and value
vector ($q,k,$ and $v$ respectively). These vectors are split equally between
the heads and stacked into matrices $\mathbb{Q},\mathbb{K}$ and $\mathbb{V}$,
whose depth is then equal to $d_{q,k,v}=d_{model}//h=125$. Finally, self-
attention is computed using equation (1).

$\displaystyle
Attention(\mathbb{Q},\mathbb{K},\mathbb{V})=softmax\left(\frac{\mathbb{QK}^{T}}{\sqrt{d_{k}}}\right)\mathbb{V}$
(1)

The output of the transformer is fed into a last fully connected layer either
with $d_{model}$ output neurons and linear activation, or with $d_{class}$
output probabilities and sigmoid activation.

###  3.5 Training

We deemed a class positive if the output probability of the last layer for
that class is greater than a fixed $threshold=0.5$. When pre-training the
model, the loss function is the mean squared error whereas during supervised
training, the loss function is the standard binary cross entropy loss averaged
across the classes. During pre-training, for a sequence of $n$ heartbeats, we
took the first $i$ ($1\leq i\leq n-1$) heartbeats and the model had to
generate a prediction of the next $i+1$ heartbeat. We used the Adam optimizer
[1] ($\beta_{1}=0.9,\beta_{2}=0.98$ and $\epsilon=10^{-9}$) whose learning
rate varies following equation (2), with $warmup\\_steps=4000$.

$\displaystyle
lr=\frac{1}{\sqrt{d_{model}}}\cdot\min\left(\frac{1}{\sqrt{step\\_num}},\frac{step\\_num}{warmup\\_steps^{1.5}}\right)$
(2)

The complete model is composed of $90,536,240$ trainable parameters,
initialized using Xavier uniform initialization. we trained the models over 50
epochs on subsets of the data made publicly available for the 2021
Physionet/CinC challenge. We ran the different models using
Tensorflow:2.5.0-gpu on three RTX 2080 Ti Turbo GPUs. We selected the
hyperparameters described in Table 1 from the references.

Table 1:  List of selected hyperparameters to train the model

Hyperparameters Value Sampling Frequency Hz 500 Batch size 128 Number of
classes, $d_{class}$ 28 Threshold 0.5 Dropout 0.1 Number of encoders, $N$ 5
Maximum position encoding 50 Embedding size, $d_{model}$ 1000 Number of heads,
$h$ 8 Feed forward layer, $\mathit{dff}$ 2048 Depth of query, key, and value
vectors, $d_{q,k,v}$ 125

Table 2:  CinC 2021 validation and final test scores and running times

Leads Challenge Datasets Running Time (minutes) Validation set CPSC test G12EC
test Undisclosed test UMich test Test set Train Validation Test All-lead 0.10
0.15 0.10 0.08 0.10 0.10 677 29 115 12-lead 0.13 0.20 0.12 0.09 0.13 0.12 677
29 115 6-lead 0.07 0.04 0.08 0.09 0.05 0.07 677 26 112 4-lead 0.12 0.18 0.11
0.08 0.11 0.10 677 26 120 3-lead 0.11 0.18 0.10 0.07 0.10 0.10 677 28 111
2-lead 0.07 0.07 0.08 0.07 0.07 0.07 677 26 112

##  4 Results

We pre-trained the model on external data sources as well as St Petersburg.
However, because of their length (30 min long), the data was often corrupted
and lead to poor prediction. After 20 epochs the loss converged to
$2.99\mathrm{E}{-5}$. We then decided to pre-train the model on Georgia and
PTB_XL and the loss converged to $1.90\mathrm{E}{-5}$ after 5 epochs. However,
neither linear probe nor transfer learning of the pre-trained model improved
the score as opposed to supervised learning alone. Therefore, pre-training was
not included in the final model.

We first trained the model on Georgia and PTB_XL. Then, adding Chapman and
Ningbo data sources positively impacted the performance. However, when
training on the full database, our model exceeded the hardware limitation.

The best model utilized the Two-Average R-peak detector and was trained on
Georgia, PTB_XL, Chapman and Ningbo. Unfortunately, we did not manage to
obtain a successful entry with that model. Our best entry, the model trained
on Geogia and PTB_XL with the Christov detector, achieved a score of 0.12,
0.07, 0.10, 0.10 and 0.07 on 12-lead, 6-lead, 4-lead, 3-lead and 2-lead test
set respectively (Table 2).

##  5 Discussion and Conclusions

Although the use of transformers and the interpretation of ECG signal as
sentences and heartbeats as words seem promising, our model did not manage to
obtain a satisfying score. Neither linear probe nor transfer learning of the
pre-trained model improved the score as opposed to supervised training alone.
Our approach is highly dependent on the R-peak detectors used. Open-source
R-peak detectors are not perfect and wrong R-peak detection could be a major
source of error. Finally, although the RMS provides good generalization, it is
interesting to explore other word embedding methods such as the CNN developed
by Sigurthorsdottir et al. 2020 [10].

## References

  * [1] Vaswani A, Shazeer N, Parmar N, Uszkoreit J, Jones L, Gomez A, Kaiser L, Polosukhin I.  Attention is all you need Dec. 2017;.
  * [2] Devlin J, Chang MW, Lee K, Toutanova K.  Bert: Pre-training of deep bidirectional transformers for language understanding May. 2019;.
  * [3] Brown T, Mann B, Ryder N, Subbiah M, Kaplan J, Dhariwal P, Neelakantan A, Shyam P, Sastry G, Askell A, Agarwal S, Herbert-Voss A, Krueger G, Henighan T, Child R, Ramesh A, Ziegler D, Wu J, Winter C, Hesse C, Chen M, Sigler E, Litwin M, Gray S, Chess B, Clark J, Berner C, McCandlish S, Radford A, Sutskever I, Amodei D.  Language models are few-shot learners Jul 2020;.
  * [4] Song H, Rajan D, Thiagarajan J, Spanias A.  Attend and diagnose: Clinical time series analysis using attention models Nov 2017;.
  * [5] Wu N, Green B, Ben X, O’Banion S.  Deep transformer models for time series forecasting: The influenza prevalence case Jan 2020;.
  * [6] Yan G, Liang S, Zhang Y, Liu F.  Fusing transformer model with temporal features for ecg heartbeat classification.  In 2019 IEEE International Conference on Bioinformatics and Biomedicine (BIBM). Jan 2019; 898–905.
  * [7] Natarajan A, Chang Y, Mariani S, Rahman A, Boverman G, Vij S, Rubin J.  A wide & deep transformer neural network for 12-lead ecg classification.  In Proceedings of the 2020 Computing in Cardiology CinC Conference 2020; .
  * [8] Perez Alday EA, Gu A, Shah AJ, Robichaux C, Wong AI, Liu C, Liu F, Rad AB, Elola A, Seyedi S, Li Q, Sharma A, Clifford GD, Reyna MA.  Classification of 12-lead ecgs: the Physionet/Computing in Cardiology Challenge 2020.  In Physiol Meas 2020. Nov 11 2020; Doi: 10.1088/1361-6579/abc960.
  * [9] Reyna M, Sadr N, Perez Alday EA, Gu A, Shah AJ, Robichaux C, Rab AB, Elola A, Seyedi S, Ansari S, Ghanbari H, Li Q, Sharma A, Clifford GD.  Will two do? Varying dimensions in electrocardiography: the Physionet/Computing in Cardiology Challenge 2021.  Computing in Cardiology 2021;48:1–4.
  * [10] Sigurthorsdottir H, Van Zaen J, Delgado-Gonzalo R, Lemay M.  ECG classification with a convolutional recurrent neural network.  In Proc. CinC 2020. Sep 2021; Rimini, Italy.

Pierre Louis Gaudilliere
Rue Jaquet-Droz 1, Neuchâtel, NE, Switzerland
<EMAIL_ADDRESS>