# Evidence that cuprate superconductors form an array of nanoscopic Josephson
junctions

Hércules S. Santana and E. V. L. de Mello  Instituto de Física, Universidade
Federal Fluminense, 24210-346 Niterói, RJ, Brazil<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Recent measurements of charge instabilities in overdoped compounds rekindled
the proposal that cuprates become superconductors by long-range order through
Josephson coupling between nanoscopic charge domains. We use the theory of
phase-ordering dynamics to show that incommensurate charge density waves
(CDWs) are formed in the CuO planes by a series of free-energy wells separated
by steep barriers. Charge oscillations in these domains give rise to a net
hole-hole attraction proportional to the height of these barriers.
Concomitantly, the self-consistent calculations yield localized
superconducting amplitudes in the CDWdomains characterizing a granular
superconductor. We show that a transition by long-range phase order promoted
by Josephson coupling elucidates many well-known features of cuprates like the
high magnetic penetration depth anisotropy and the origin of the pseudogap,
among others. Furthermore, the average Josephson energy reproduces closely the
planar superfluid density temperature dependence of La-based films and the
superconducting giant proximity effects of cuprates, a 20-year-old open
problem.

##  I Introduction

One of the main challenges of condensed-matter physics is a complete theory
for high critical temperature superconductors (HTSs). Such a theory has been
hindered by many issues like the absence of clear Fermi surfaces, the
pseudogap, and the prominence of various forms of collective
fluctuationsKeimer  _et al._ (2015). To characterize whether these distinct
orders compete with or strengthen each other, many experimen tal techniques
have been refined to detect even overlapping fluctuations. After imaging a
granular structure with high spatial resolution scanning tunneling microscopy
(STM) in underdoped Bi2Sr2CaCu2O8+d (Bi2212), Lang et alLang  _et al._ (2002)
proposed that superconducting (SC) long-range order could be achieved by
Josephson coupling between nanoscopic domains. This idea gained more
recognition after the measurements of CDW in YBa2Cu3O6+x (YBCO) single
crystals but was still not considered as a general theory of cuprates for two
reasons: first, the absence of incommensurate charge ordering (CO) or CDW data
in the overdoped region and, second, the lack of a theoretical model that
could justify and describe the physical formation of the small Josephson
junctions between the CO domains.

In this paper, we use the theory of phase-ordering dynamics to show that the
CDW or CO may be formed by a two-dimensional array of free energy potential
wells with similar properties of granular superconductors. In recent years
there was a great improvement in the precision of the CO wavelength
$\lambda_{\rm CO}$ measurements by STM, x-ray, and Resonant X-Ray Scattering
(REXS)Comin and Damascelli (2016). The very fine variation of $\lambda_{\rm
CO}$ with the doping $p$ (or hole per CuO unit cell) revealed in these
experiments can be reproduced theoretically by a phase separation formalism
based on the time-dependent nonlinear Cahn-Hilliard (CH) differential
equationCahn and Hilliard (1958). In this approach, the charge modulations
maybe tuned up to reproduce the measured $\lambda_{\rm CO}(p)$ on 100% volume
fraction of the simulationsde Mello  _et al._ (2009); de Mello (2012); de
Mello and Kasal (2012); de Mello and Sonier (2014, 2017).

The phase separation free energy reproduces the CDW structure and acts as a
crystal field that promotes the Cooper pair formation, leading to a direct
connection between the SC interaction and the charge modulations. Some
critique the CDW-mediated superconductivity, arguing that CDW is limited to
the underdoped region. However recent measurementsWu  _et al._ (2017); Chen
_et al._ (2019); Fei  _et al._ (2019) and new x-ray diffraction demonstrated
CDW correlations in overdoped La2-xSrxCuO4 (LSCO) up to compounds of at least
$x\equiv p=0.21$ and possibly up to $p=0.25$Miao  _et al._ (2021). Here we
perform CDW simulations and develop a SC theory of cuprates based on Josephson
coupling between local SC order parameters in these domains and their long-
range phase order (LRO). The Josephson coupling is closely related with the
superfluid densitySpivak and Kivelson (1991); de Mello (2021) $\rho_{\rm sf}$
and reproduces its measured temperature variation $\rho_{\rm sf}(T)$ of
several overdoped LSCO films with great accuracy.

The CDW-LRO approach is also appropriate to describe the giant proximity
effects (GPE) experiments in YBCO S-I-S wiresDecca  _et al._ (2000) and in
LSCO Josephson type trilayers S-N’-S junctionsBozovic  _et al._ (2004);
Morenzoni  _et al._ (2011) where I is an insulator with nonzero doping, S is a
superconductor, and N’ is a superconductor layer in the normal phase.
According to conventional theory, the critical current should diminish
exponentially with the size $d$ of a barrier made of non-SC materialsBozovic
_et al._ (2004). For traditional low temperature Josephson junctions the
thickness of the barrier $d$ is generally comparable with the barrier
coherence lengthBozovic  _et al._ (2004) $\xi_{\rm N}$. However, despite the
very small SC coherence lengths $\xi_{\rm SC}$ of cuprates, GPE was measured
with several large barriersDecca  _et al._ (2000); Bozovic  _et al._ (2004)
with $d\gg\xi_{\rm N}$. Thus, technically the YBCO non-SC spacer with $d\sim
100$ nm and La-based trilayers with $d\sim 100-1000$ nm are too large to carry
a critical current, in clear contradiction with the experimentsDecca  _et al._
(2000); Bozovic  _et al._ (2004); Morenzoni  _et al._ (2011). On the other
hand, in our approach, S, N’, and I all have CDW and differ only by the
presence or absence of LRO, which is very sensitive to external perturbations
like an applied current or magnetic field.

##  II The CDW simulations

The starting point is the definition of the time-dependent phase separation
(PS) order parameter associated with the local electronic density, $u({\bf
r},t)=[p({\bf r},t)-p]/p$, where $p$ is the average hole density and $p({\bf
r},t)$ is the charge density at a position ${\bf r}$ in the CuO plane and at a
time of simulation $t$. The CH equation is based on the electronic phase
separation Ginzburg-Landau (GL) free energy expansion in terms of the
conserved charge order parameter $u$ de Mello and da Silveira Filho (2005); de
Mello  _et al._ (2009); de Mello (2012); de Mello and Kasal (2012); de Mello
and Sonier (2014, 2017).

$f(u)={{\frac{1}{2}\varepsilon|\nabla u|^{2}+V_{\rm GL}(u,T)}}\;,$ (1)

where $\varepsilon$ is the parameter that controls the charge modulations
scale and ${V_{\rm GL}}(u,T)=-\alpha[T_{\rm PS}-T]u^{2}/2+B^{2}u^{4}/4+...$ is
a temperature-dependent double-well potential that characterizes the rise of
charge oscillations below the onset of phase separation temperature $T_{\rm
PS}$. We do not know $T_{\rm PS}$, but there are indications that it is close
to the pseudogap temperature $T^{*}$. In the simulations, when $T\leq T^{*}$,
the values of $\alpha$ and $B$ are always one. This free energy in terms of
the phase separation order parameter is much simpler than the Ginzburg-Landau-
Wilson free energy in terms of SC and pair density wave fields (PDW)Fradkin
_et al._ (2015), but it suitably reproduces the details of the CO structure of
distinct compounds and their localization energy ${V_{\rm GL}}$.

Figure 1:  (a) The $V_{\rm GL}(t)$ evolution with time. Initially, $V_{\rm
GL}$ is flat, which corresponds to a system above the onset phase separation
temperature $T_{\rm PS}$. The phase separation potential wells or valleys
increase with time, and this behavior is correlated with the decrease of the
temperature below $T_{\rm PS}$. As the temperature goes below $T^{*}$, we
assume the stable large time behavior $V_{\rm GL}(t\rightarrow\infty)$ to be
that of low temperature $V_{\rm GL}(T\leq T^{*})$. (b) The Cooper pair
formation in the CDW free energy valleys. At the top, we represent some planar
Cu atoms (blue-filled circles) attracted to hole-poor domains represented by
black lines as an illustration. Hole fluctuations in these domains produces
atomic fluctuations that affect also the other holes, promoting an atomic
mediated interaction (represented by the springs as an illustration). At low
temperature ($T\leq T^{*}$) the Cooper pairs may be formed in the CDW valleys
(the encircled pair of black dots), and at $T\leq T_{c}$ they superflow and
become uniform on the CuO plane (in agreement with the CO x-ray scattering
decreasing signal below $T_{c}$Wise  _et al._ (2008); Chang  _et al._ (2012)).

The CH equation can be written derived by a continuity equation of the local
free energy current density ${\bf J}=M{\bf\nabla}(\delta f/\delta u)$,Bray
(1994)

$\displaystyle\frac{\partial u}{\partial t}$ $\displaystyle=$
$\displaystyle-{\bf\nabla.J}$ (2) $\displaystyle=$
$\displaystyle-M\nabla^{2}[\varepsilon^{2}\nabla^{2}u-\alpha^{2}(T)u+B^{2}u^{3}],$

where $M$ is the mobility or the charge transport coefficient that sets both
the phase separation time scale and the contrast between the values of $u$ for
the two phases.

The equation is solved by a stable and fast finite difference scheme with free
boundary conditionsde Mello and da Silveira Filho (2005), yielding the phase
separation order parameter $u({\bf r},t=n\delta t)$, function of position
${\bf r}$ and $n$ simulation time step $\delta t$. The limiting cases are
$u({\bf r}_{i},t)\approx 0$, corresponding to homogeneous systems above or
near $T_{\rm PS}$ or small charge variations like the observed CDW, and
$u({\bf r}_{i},t\rightarrow\infty)=\pm\>1$, corresponding to the extreme case
(at low temperatures) of complete phase separation. The local charge density
is derived from $p({\bf r},t)=p\times(u({\bf r},t)+1)$, and the latter case
(strong phase separation) applies to static stripesTranquada  _et al._ (1995);
Thampy  _et al._ (2017), while the former (weak phase separation) to very
small $\Delta p\approx 10^{-2-3}$ variations around $p$, like that measured in
YBa2Cu3O6+δ (Y123)Kharkov and Sushkov (2016).

Figure 1(a) shows that below the $T_{\rm PS}$ temperature the order parameter
$u({\bf r},t)$ evolves in time and $V_{\rm GL}({\bf r},t)$ valleys become
deeper, which is expected to occur when the temperature decreases, favoring
the mesoscopic phase separation. As the temperature goes below $T^{*}$, we
assume the large time behavior $V_{\rm GL}({\bf r},t\rightarrow\infty)$ shown
in Fig. 1(a) to become the low temperature $V_{\rm GL}({\bf r},T\rightarrow
0)$ that generates the CDW [Fig. 1(b)]. In this large time regime, $V_{\rm
GL}({\bf r})$ depends on the temperature through the usual form of the first
GL coefficient defined after Eq. 1, that is, $(1-T/T^{*})^{2}$. This
temperature dependence is the only relevant dependence of $V_{\rm GL}({\bf
r},T)$ in our calculations.

Figure 2:  Top and three dimensional perspective of different properties at a
later time or low temperature: (a) A checkerboard CDW map for a LSCO system
with average hole density $p=0.12$ on a $100\times 100$ unit cell simulations
and out of plane view. (b) The corresponding GL free energy potential $V_{\rm
GL}({\bf r}_{i})$. The inset shows the three-dimensional perspective, which
reveals the array of deep wells. (c) The BdG SC amplitude $\Delta_{d}({\bf
r}_{i})$ map in a $28\times 28$ unit cell portion of the density map of (a).
The three-dimensional perspective plot demonstrates that the amplitudes
$\Delta_{d}({\bf r}_{i})$ are completely localized inside the CDW domains. (d)
Typical histogram of the amplitudes shows the $\Delta_{d}({\bf r}_{i})$
variation with the local $p_{i}$ concentration.

A typical later time and low-temperature $V_{\rm GL}(u({\bf r}))$ used in the
calculations is shown in Fig. 2(b), which leads to the LSCO checkerboard
structure. We also show the three-dimensional side view (in the inset) with
its wells or valleys in form of “ice cream cones”. Notice that the free energy
is defined over the CuO plane, and the third dimension is its values in each
site, which demonstrates that the minima occur at the center of the charge
domains and the maxima at the borders. The detailed structure and strength of
this potential is shown in the Supplemental MaterialSup .

For YBCO, the CDW domains are formed in puddles or patches with different
modulation in either a or b-direction in the CuO planeComin  _et al._ (2015);
Kim  _et al._ (2021). The LRO superconducting calculations with $V_{\rm
GL}(u({\bf r}))$ and alternating stripe-like puddles follow in the same way of
the checkerboard of Fig. 2 since the method requires only the formation of
such finite-charge domainsde Mello and Sonier (2017). We emphasize that the
calculations described and performed here are pure two dimensional in order to
model the observed charge arrangements on the CuO plane and does not take into
account possible out-of-plane influences, like the quenched disorder which
arises from oxygen interstitials discovered in the in HgBa2CuO4+y systemCampi
_et al._ (2015).

##  III Superconducting Calculations

The new form of analysis here is the three-dimensional view of the planar CDW
maps shown in Fig. 2(a), the GL free energy potential $V_{\rm GL}$ in Fig.
2(b), and the SC pair amplitude $\Delta_{d}({\bf r}_{i})$ in a $28\times 28$
site portion of the $100\times 100$ site density map from Fig 2(a) in Fig.
2(c). These illustrative plots are for an LSCO sample with an average hole
density $p=0.12$. These plots show clearly that alternating rich and poor
charge regions develop in the same kind of valleys in the form of “ice cream
cones”, as it is demonstrated by the out-of-plane view of the insets. The SC
pair amplitudes inside the wells in Fig. 2(c) have the same properties of
isolated SC grains, which motivates the proposal of an array of Josephson
junctions.

While Fig. 1(a) shows the evolution of $V_{\rm GL}(x_{i})$ that leads to the
CDW, Fig. 1(b) shows the low-temperature $V_{\rm GL}(x_{i})$ that leads to the
formation of the SC interaction in the CuO planes. At the top we represent
some planar atoms (blue circles) attracted (repelled) by hole-poor (hole-rich)
CDW domains represented schematically by black lines. High-energy x-ray
diffractionChang  _et al._ (2012) revealed that CDW modulations displace the
Cu and O atoms whose oscillations around their equilibrium positions are
sensed by the holes and may give rise to a net hole-hole attractionde Mello
and Sonier (2017); de Mello (2021) illustrated by the springs.

We argued before that this induced hole-hole SC interaction is proportional to
the depth of the $V_{\rm GL}$ wells averaging over the whole system, that is,
${\left<V_{\rm GL}(p)\right>}$de Mello and Sonier (2017); de Mello (2021). The
mean-field self-consistent Bogoliubov-de Gennes (BdG) calculations with this
attractive pair interaction over the whole system yield the local pair
amplitude map $\Delta_{d}({\bf r}_{i})$. The results also depend on the local
hole density $p(i)$ and, consequently, have the same CDW modulations
($\lambda_{\rm CO}$), leading to what is known as pair density wave Keimer
_et al._ (2015).

Since the $\Delta_{d}({\bf r}_{i})$ fit perfectly inside the wells, each CDW
domain has an independent local SC order parameter phase $\phi_{i}$, exactly
as a granular superconductor. Therefore, each CDW domain may behave like
mesoscopic SC grains interconnected by Josephson junctions with energies
$E_{\rm J}({\bf r}_{lm})$ between $l$ and $m$ domains like those shown with
alternate red and blue colors in Fig. 2(c). In the bottom of Fig. 1(b) we show
that, according to previous calculationsde Mello and Sonier (2017); de Mello
(2021), the Cooper pairs are present up to $T\leq T^{*}$ and are represented
by two encircled black dots. At $T\leq T_{c}$, LRO sets in, and the Cooper
pairs superflow through the system, establishing a uniform charge density on
the CuO plane. This leads to a decrease of the CDW signal observed by high-
energy x-ray diffraction below $T_{c}$Wise  _et al._ (2008); Chang  _et al._
(2012) and was widely interpreted by due to the competition between CDW and
superconductivity. We emphasize that our calculations establish the opposite:
The CDW hosts the charge domains and the local Cooper pairs in alternating
domains, but they spread through the system when LRO sets in at $T\leq T_{c}$.
This is shown schematically at the bottom of Fig. 1(b).

Following the above arguments, $T_{\rm c}$ is determined by the competition
between thermal disorder and the average planar Josephson energy
${\left<E_{\rm J}(p,T)\right>}$ that depends also on the spatial averaged
$d$-wave pair amplitude

${\left<\Delta_{d}(p,T)\right>}=\sum_{i}^{N}\Delta_{d}({\bf r}_{i},p,T)/N,$
(3)

where the sum, like in the case of ${\left<V_{\rm GL}(p)\right>}$, is over the
N unit cells of a CuO plane. As explained previouslyde Mello and Sonier (2014,
2017), the $d$-wave relation for the average Josephson coupling energy is
proportional to the $s$-wave expressionAmbegaokar and Baratoff (1963):

${\left<E_{\rm
J}(p,T)\right>}=\frac{\pi\hbar{\left<\Delta_{d}(p,T)\right>}}{4e^{2}R_{\rm
n}(p)}{\rm tanh}\bigl{[}\frac{\left<\Delta_{d}(p,T)\right>}{2k_{\rm
B}T}\bigr{]}.$ (4)

Here $R_{\rm n}(p)$ is the average planar tunneling resistance between the
grains that is assumed to be proportional to the in-plane normal state
resistance just above $T_{\rm c}$. In the array of Josephson junction model,
the current is composed of Cooper pairs tunneling between the CDW domains and
by the normal carriers or quasi-particle plane currentBruder  _et al._ (1995).
For a $d$-wave HTS near $T_{\rm c}$ the supercurrent is dominantBruder  _et
al._ (1995), which justifies the use of the experimental $R_{\rm n}(p)$ in the
above equation. Therefore, LRO is attained when the average ${\left<E_{\rm
J}(p,T)\right>}$ is strong enough to overcome thermal phase disorder, or
${\left<E_{\rm J}(p,T_{\rm c})\right>}=k_{\rm B}T_{\rm c}$, and that is how we
derive $T_{\rm c}$de Mello and Sonier (2014, 2017); de Mello (2020a, b).

These in-plane calculations are the basic pillars to the three-dimensional LRO
in the whole system, which we may infer from transport measurements. For low
$p$, just above $T_{\rm c}$, the $z$-direction resistivity $\rho_{c}$ is
$\approx 10^{3}-10^{6}$ larger than the $a$ or $b$-axis resistivityOno and
Ando (2003); Komiya  _et al._ (2002). Despite this large difference, it is
surprising that both $\rho_{c}(T)$, and $\rho_{ab}(T)$ fall to zero at the
same temperature ($T_{\rm c}$) and this puzzling behavior can be understood in
terms of the Josephson coupling of Eq. 4. The smaller planar resistances yield
larger $E_{\rm J}$ that promote first LRO in the planes, but each plane $i$
would have its own SC phase $\theta_{i}$ if it was not for the weaker inter-
plane $E_{\rm J}$ coupling. Thus the $z$-direction Josephson coupling connects
the planes, leading to only a single-phase $\theta$ in the whole system, and
both $z$ and $ab$ resistivity drop off together at $T_{\rm c}$, which is a
plausible explanation for this long-known non-conventional result. Again using
Eq. 4, a smaller low temperature superfluid density is expected along the
$z$-direction than along the plane, which is confirmed by the measured large
anisotropy of $ab$ and $c$-axis magnetic penetration depthPanagopoulos  _et
al._ (2000, 1999). Thus, although $z$-direction coupling is fundamental, the
SC properties like the local SC amplitudes $\Delta_{d}({\bf r}_{i},p,T)$ and
$T_{\rm c}(p)$ develop and depend entirely on the CuO planes.

##  IV the Superfluid Density

For low-temperature superconductors, the temperature $T_{\theta}$ at which LRO
disappears is very large compared with $T_{\rm c}$, but they are estimated to
be comparable for HTSEmery and Kivelson (1995). Our basic point is that
$T_{\rm c}=T_{\theta}$, and at higher temperatures, the Cooper pairs in the
nanoscopic grains have all different $\theta_{i}$ and, consequently are in an
incoherent state. This approach leads us to infer a close connection between
Josephson coupling and superfluid phase stiffness.

The first thing to notice is that $J_{\rm c}(T)$ is proportional to the two-
dimensional (2D) superfluid density $n_{\rm s}(T)$ or the local Josephson
currentSpivak and Kivelson (1991) that is also proportional to the phase
stiffness $\rho_{\rm sc}$Božović  _et al._ (2016). Along these lines, we made
${\left<E_{\rm J}(p,0)\right>}$ equal to the 2D zero-temperature superfluid
phase stiffnessde Mello (2020a, 2021) $\rho_{\rm sc}(p,0)$ and reproduced the
measured scale relation between the zero temperature superfluid density and
$T_{\rm c}(p)$. To underdoped compounds, this relation is known as Uemura’s
lawUemura  _et al._ (1989), and a similar relation for La-based overdoped
films was obtained more recentlyBožović  _et al._ (2016).

Along these lines, we show in Fig. 3 that $({\left<E_{\rm
J}(p,T)\right>}/k_{\rm B}-T)$ reproduces the measuredBožović  _et al._ (2016)
2D phase stiffness $\rho_{\rm sc}(p,T)$ temperature dependence. The plots show
very good agreement with $p=0.16$ and 0.19 films according to Ref. [Božović
_et al._ , 2016]. For $p=0.22$ we could extract with certainty only the
experimental values at $T=0$ and $T_{\rm c}$, but the theory seems to
reproduce closely also the intermediate data. In the plots we used
${\left<E_{\rm J}(p,T)\right>}$ from Eq. 4 with
${\left<\Delta_{d}(p,T)\right>}$ and $R_{\rm n}(p)$ derived previously in
Refs.[de Mello, 2021,de Mello, 2020a].

Figure 3:  Comparison between the measured phase stiffness $\rho_{\rm
sc}(p,T)$ for La-based films with $p=0.16,0.19$, and 0.22 from Ref.[Božović
_et al._ , 2016]) and $({\left<E_{\rm J}(p,T)\right>}/k_{\rm B}-T)$ (dashed
lines), with an almost perfect agreement. We show also the Josephson coupling
${\left<E_{\rm J}(0.16,T)\right>}$ that sets the SC phase order scale
$T_{\theta}$ derived in Ref. [de Mello, 2021] and the temperature $T$ to make
clear how to derive $T_{\rm c}(p=0.16)$.

Another indication that ${\left<E_{\rm J}(p,0)\right>}/k_{\rm B}-T$ is equal
to $\rho_{\rm sc}(p,T)$ comes from the proximity effect experiment on YBCO
wiresDecca  _et al._ (2000). They measured the critical current $I_{\rm c}$
through a small insulator region of size $d$ inserted in the SC wire. In their
Fig. 2(b), they showed that the product $I_{\rm c}R_{\rm n}$ at low-
temperature is constant for several values of $d$ from 40 to 110 nm, which is
much larger than the estimated coherence length $\xi_{\rm SC}$ of 9 nmDecca
_et al._ (2000). We interpret this result recalling that the low-temperature
Josephson current isKetterson and Song (1999)

$I_{\rm J}=\frac{E_{\rm J}(0)}{\hbar/(2e)}=\frac{\pi\Delta_{d}(0)}{2eR_{\rm
n}}$ (5)

Following the above discussion, taking $I_{\rm J}=I_{\rm
c}\propto{\left<E_{\rm J}(p,0)\right>}$ and according to Eq. 4, the product
$I_{\rm c}R_{\rm n}$ is independent of the insulator spacer $d$, in perfect
agreement with the measurements with the S-I-S YBCO junctions. This result, in
the framework of our theory, suggests that the CDW is also present in the
underdoped insulator part of the YBCO wire, in agreement with the connection
between CDW and the antinodal PG mentioned in the previous sectionWise  _et
al._ (2009); Comin  _et al._ (2014).

##  V Results and Discussion

In general, the critical current density is the product of the superfluid
velocity and density, that is, $J_{\rm c}=v_{s}\rho_{\rm sc}$. We have shown
above that ${\left<E_{\rm J}(p,T)\right>}/k_{\rm B}-T$ is equal to $\rho_{\rm
sc}(T)$ and that the low-temperature critical current $I_{\rm
c}(0)\propto{\left<E_{\rm J}(p,0)\right>}$ on the YBCO wiresDecca  _et al._
(2000). Based on these two results, we assume that the average superfluid
velocity is $v_{s}\propto{\left<E_{\rm J}(p,T)\right>}/{\left<E_{\rm
J}(p,0)\right>}$ and that,

$J_{\rm c}(T)\propto\left({\left<E_{\rm J}(p,T)\right>}/k_{\rm
B}-T\right)\frac{{\left<E_{\rm J}(p,T)\right>}}{{\left<E_{\rm
J}(p,0)\right>}}\;.$ (6)

This expression will be used to compare with the GPE measured critical
currentsDecca  _et al._ (2000); Bozovic  _et al._ (2004).

The first system to deal here is the trilayers along the $z$ direction
composed of La1.85Sr0.15CuO4 films with $T_{\rm c}\approx 45$ K (S) and the
underdoped superconductor La2CuO4+δ with $T^{\prime}_{\rm c}\approx 25$ K (N’)
sandwiched according to Ref.[Bozovic  _et al._ , 2004]. Initially we study S
and N’ separately, with the average amplitudes
${\left<\Delta_{d}(p=0.12,T)\right>}$ and
${\left<\Delta_{d}(p=0.15,T)\right>}$ de Mello (2020a, 2021) and the factor of
three in their normal resistances $R_{\bf n}$ just above their $T_{\rm c}$
according to Ref. [Bozovic  _et al._ , 2002]. For $T\leq 25$ K there is LRO in
both S and N’. Therefore, in principle, $J_{\rm c}$ would be limited by the
lower $T^{\prime}_{\rm c}\approx 25$ K, but they measured a persisting $J_{\rm
c}$ almost to $T=35$ K independently of the number of underdoped N’
layersBozovic  _et al._ (2004).

Figure 4:  S and N’ layers plots after they are connected in the S-N’-S
structure both with the average hole density $p^{\prime}=0.135$.
${\left<E_{\rm J}(p^{\prime},T)\right>}$ for S (green) and N’ (purple). We
also plot $[{\left<E_{\rm J}(p^{\prime},T)\right>}/k_{\rm B}-T]=\rho_{\rm
sc}(p^{\prime},T)$ for S (red) and N’ (yellow). The smaller $\rho_{\rm
sc}(p^{\prime},T)$ of N’ is the maximum critical current through the S-N’-S
system [the superfluid density through S (red line) is larger] and is assumed
to be proportional to the theoretical critical current $J_{\rm c}$ (dashed
line).

The main reason to have supercurrents above $T^{\prime}_{\rm c}=25$ K is the
presence of SC amplitudes $\Delta_{d}$ in the CDW charge domains without LRO,
what characterizes the PG region of the N’ superconductor. Another important
factor is the rearrangement of the holes between the layers according to
subsequent experiments: using similar La-based materials but combining undoped
insulators with overdoped metallic La2-xSrxCuO4 with $x=0.38$Smadici  _et al._
(2009) and $0.44$Suter  _et al._ (2012). The authors of those studies verified
that the conducting holes pile up at the interface and redistribute themselves
from the hole-rich to the insulator layers. Analysis with REXSSmadici  _et
al._ (2009) and $\mu$-SRSuter  _et al._ (2012) revealed that this
reorganization of the holes occurs near the interface and they are not
followed by the Sr dopant atoms. The rearrangement of the hole densities add
on or subtract carriers and affect directly the local pair amplitudes
$\Delta_{d}$, which depend on the local densities $p(i)$ as discussed above.
The redistribution and uniformization of the charges will enhance
${\left<\Delta_{d}\right>}$ in N’ with a concomitant decrease in S.

Figure 5:  The solid red line is the S-N’-S experimental results from Fig. 2
of Ref. [Bozovic  _et al._ , 2004], and the dashed line is our calculations
that are proportional to $J_{\rm c}$. We made both curves equal near $T=5$ K
because our estimates are more precise at low temperatures. The inset shows
schematically the coupled systems used in the calculations. In the trilayer
experimentsBozovic  _et al._ (2004); Morenzoni  _et al._ (2011) the layers are
on top of each other and connected by a small Josephson coupling, but in the
text, we explain that the calculations are similar.

This enhancement or weakening was observed in Meissner effects studies by low
energy $\mu$-SRMorenzoni  _et al._ (2011). They showed that the local magnetic
field ${\left<B_{x}\right>}$ along the $z$ direction in a similar
heterostructure S-N’-S is excluded like in a single uniform superconductor
even at temperatures more than 3 times larger than $T^{\prime}_{\rm c}$. This
reinforces our point that S and N’ at $T^{\prime}_{\rm c}<T<T_{\rm c}$ have
both the same structure differing only by the presence or absence of LRO in
the charge domains.

Now we will apply these ideas to the GPE in the La-based trilayer
systemBožović  _et al._ (2016) . The effect of the checkerboard charge
inhomogeneities in the local SC properties were studied in detailde Mello
(2020a) in a system similar to N’ with $p=0.12$, and the simulations were
shown above in Fig. 2. The $\Delta_{d}({\bf r}_{i})$ amplitude histogram is
shown in Fig. 2(d), and the calculations show that CDW higher (lower) local
hole densities have larger (smaller) local SC amplitudesde Mello and Sonier
(2017). Therefore, when a supercurrent flows, the S-N’-S system acquires a
mean hole density close to their average; that is, the mean hole density in N’
and S become $p^{\prime}\approx 0.135$. This rearrangement increases
${\left<\Delta_{d}(N^{\prime},T)\right>}$ by a factor of approximately 1.5 and
decreases ${\left<\Delta_{d}(S,T)\right>}$ by a factor of approximately 0.6.
Taking these changes into account, we calculate the new Josephson couplings in
N’ and S keeping their original $R_{\rm n}$ unchanged. The results of these
calculations are shown in Fig. 4, where we plot these renormalized
${\left<E_{\rm J}(p^{\prime},T)\right>}$, that is with $p^{\prime}=0.135$
(green for S and purple for N’) and $\rho_{\rm sc}(p^{\prime},T)$ for the N’
(yellow) and S (red) coupled layers. Since the critical current must be
constant through the layers and N’ holds the smaller values, we take exactly
this (yellow) lower current as the theoretical critical current through the
entire system (dashed line).

Since $J_{\rm c}=v_{s}\rho_{\rm sc}$, $v_{s}(p,T)\propto{\left<E_{\rm
J}(p,T)\right>}/{\left<E_{\rm J}(p,0)\right>}$, and using the dashed
$\rho_{\rm sc}(p^{\prime},T)$ curve on Fig. 4, we obtain our estimation to the
critical current measurementsBozovic  _et al._ (2004) shown in Fig. 5. To
perform this comparison, we made our results equal to the lowest measured
temperature $T\approx 5$ K, where our estimation to $v_{s}(T)$ is more
precise. Our results follow near the experimental $J_{\rm c}(T)$ but do not
have the round features displayed by the red curve that was taken from Ref.
[Bozovic  _et al._ , 2004]. This is mainly because of the poor estimation of
$v_{s}(T)$ at finite temperatures. Another reason is that the S-N’-S system is
formed with stacked layers along the $z$ direction and our calculations are
with layers in the CuO plane as shown in the inset of Fig. 5. But as discussed
before, the calculations follow along the same lines in the sense that LRO is
attained first in CuO planes and, afterward, by interlayer Josephson coupling,
the whole system becomes superconducting. Another possible correction is that
the $z$-direction current may enhance influences of out-of-plane dopants or
oxygen interstitialsCampi  _et al._ (2015).

We argued that the formation of a granular superconductivity in the CO domains
is the most basic property of cuprate superconductors. We showed that phase-
ordering kinetics may describe electronic phase separation transitions by the
GL potential $V_{\rm GL}$ with incommensurate free energy wells that originate
the observed CDW. Charge fluctuations inside the domains may induce hole-hole
SC interaction proportional to the depth of the $V_{\rm GL}$ wells averaging
over the whole system, that is, ${\left<V_{\rm GL}(p)\right>}$de Mello and
Sonier (2017); de Mello (2021). The average Josephson coupling ${\left<E_{\rm
J}(p,T)\right>}$ between the local SC order parameter yields the LRO or phase-
ordering temperature that is made equal to $T_{\rm c}$. Furthermore,
${\left<E_{\rm J}(p,T)\right>}-k_{\rm B}T$ provides perfect agreement with the
measured superfluid density temperature dependence of La-based overdoped
filmsBožović  _et al._ (2016). The granular superconductivity is also the key
mechanism behind the GPE because both S, N’, and I (with nonzero doping) used
in the experiments have all localized SC order parameters and differ only by
the absence or presence of LRO.

We should emphasize that the Josephson coupling of Eq. 4 is inversely
proportional to the resistivity, which leads to much smaller $E_{\rm J}$
inter-plane coupling, and its identification with the local superfluid density
is confirmed by the anisotropic magnetic penetration length
measurementsPanagopoulos  _et al._ (1999, 2000). We finish pointing out that
the presence of CDW in electron-doped compoundsda Silva Neto  _et al._ (2018)
is a further indication that the method applies to all cuprates, and which we
will explore further in the future.

##  VI Acknowledgements

We are grateful to David Möckli for critical reading the manuscript and
acknowledge partial support by the Brazilian agencies CNPq and FAPERJ.

## References

  * Keimer  _et al._ (2015) B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179–186 (2015).
  * Lang  _et al._ (2002) K. M. Lang, V. Madhavan, J. E. Hoffman, E. W. Hudson, H. Eisaki, S. Uchida, and J. C. Davis, Nature 415, 412 (2002).
  * Comin and Damascelli (2016) R. Comin and A. Damascelli, Ann. Rev. of Cond. Mat. Phys. 7, 369 (2016).
  * Cahn and Hilliard (1958) J. W. Cahn and J. E. Hilliard, J. Chem. Phys. 28, 258 (1958).
  * de Mello  _et al._ (2009) E. V. L. de Mello, R. B. Kasal, and C. A. C. Passos, J. Phys.: Condens. Matter 21, 235701 (2009).
  * de Mello (2012) E. V. L. de Mello, Europhys. Lett. 99, 37003 (2012).
  * de Mello and Kasal (2012) E. V. L. de Mello and R. B. Kasal, Physica C: Superconductivity 472, 60 (2012).
  * de Mello and Sonier (2014) E. V. L. de Mello and J. E. Sonier, J. Phys.: Condens. Matter 26, 492201 (2014).
  * de Mello and Sonier (2017) E. V. L. de Mello and J. E. Sonier, Phys. Rev. B 95, 184520 (2017).
  * Wu  _et al._ (2017) J. Wu, A. T. Bollinger, X. He, and I. Božović, Nature 547, 432 (2017).
  * Chen  _et al._ (2019) S.-D. Chen, M. Hashimoto, Y. He, D. Song, K.-J. Xu, J.-F. He, T. P. Devereaux, H. Eisaki, D.-H. Lu, J. Zaanen, and Z.-X. Shen, Science 366, 1099 (2019).
  * Fei  _et al._ (2019) Y. Fei, Y. Zheng, K. Bu, W. Zhang, Y. Ding, X. Zhou, and Y. Yin, Science China Physics, Mechanics & Astronomy 63, 227411 (2019).
  * Miao  _et al._ (2021) H. Miao, G. Fabbris, R. J. Koch, D. G. Mazzone, C. S. Nelson, R. Acevedo-Esteves, G. D. Gu, Y. Li, T. Yilimaz, K. Kaznatcheev, E. Vescovo, M. Oda, T. Kurosawa, N. Momono, T. Assefa, I.K. Robinson, E.S. Bozin, J.M. Tranquada, P.D. Johnson, M.P.M. Dean, npj Quantum Materials 6, 31 (2021).
  * Spivak and Kivelson (1991) B. I. Spivak and S. A. Kivelson, Phys. Rev. B 43, 3740 (1991).
  * de Mello (2021) E. V. L. de Mello, J. of Phys.: Cond. Matter 33, 145503 (2021).
  * Decca  _et al._ (2000) R. S. Decca, H. D. Drew, E. Osquiguil, B. Maiorov, and J. Guimpel, Phys. Rev. Lett. 85, 3708 (2000).
  * Bozovic  _et al._ (2004) I. Bozovic, G. Logvenov, M. A. J. Verhoeven, P. Caputo, E. Goldobin, and M. R. Beasley, Phys. Rev. Lett. 93, 157002 (2004).
  * Morenzoni  _et al._ (2011) E. Morenzoni, B. M. Wojek, A. Suter, T. Prokscha, G. Logvenov, and I. Bozovic, Nature Communications 2, 272 (2011).
  * de Mello and da Silveira Filho (2005) E. de Mello and O. T. da Silveira Filho, Physica A 347, 429 (2005).
  * Fradkin  _et al._ (2015) E. Fradkin, S. A. Kivelson, and J. M. Tranquada, Rev. Mod. Phys. 87, 457 (2015).
  * Wise  _et al._ (2008) W. D. Wise, M. C. Boyer, K. Chatterjee, T. Kondo, T. Takeuchi, H. Ikuta, Y. Wang, and E. W. Hudson, Nature Physics 4, 696 (2008).
  * Chang  _et al._ (2012) J. Chang, E. lackburn, . T. Holmes, N. B. Christensen, J. Larsen, J. Mesot, R. Liang, D. A. Bonn, W. N. Hardy, A. Watenphul, M. V. Zimmermann, E. M. Forgan, and S. M. Hayden, Nature Physics 8, 871 (2012).
  * Bray (1994) A. Bray, Adv. Phys. 43, 357 (1994).
  * Tranquada  _et al._ (1995) J. M. Tranquada, B. J. Sternlieb, J. D. Axe, Y. Nakamura, and S. Uchida, Nature 375, 561 (1995).
  * Thampy  _et al._ (2017) V. Thampy, X. M. Chen, Y. Cao, C. Mazzoli, A. M. Barbour, W. Hu, H. Miao, G. Fabbris, R. D. Zhong, G. D. Gu, J. M. Tranquada, I. K. Robinson, S. B. Wilkins, and M. P. M. Dean, Phys. Rev. B 95, 241111 (2017).
  * Kharkov and Sushkov (2016) Y. A. Kharkov and O. P. Sushkov, Scientific Reports 6, 34551 (2016).
  * (27) “See supplemental material at [url will be inserted by publisher]for potential simulations and discussions.” .
  * Comin  _et al._ (2015) R. Comin, R. Sutarto, E. H. da Silva Neto, L. Chauviere, R. Liang, W. N. Hardy, D. A. Bonn, F. He, G. A. Sawatzky, and A. Damascelli, Science (New York, N.Y.) 347, 1335 (2015).
  * Kim  _et al._ (2021) H.-H. Kim _et al._ , Phys. Rev. Lett. 126, 037002 (2021).
  * Campi  _et al._ (2015) G. Campi, A. Bianconi, N. Poccia, G. Bianconi, L. Barba, G. Arrighetti, D. Innocenti, J. Karpinski, N. D. Zhigadlo, S. M. Kazakov, M. Burghammer, M. v. Zimmermann, M. Sprung, and Ricci, Nature 525, 359 (2015).
  * Ambegaokar and Baratoff (1963) V. Ambegaokar and A. Baratoff, Phys. Rev. Lett. 10, 486 (1963).
  * Bruder  _et al._ (1995) C. Bruder, A. van Otterlo, and G. T. Zimanyi, Phys. Rev. B 51, 12904 (1995).
  * de Mello (2020a) E. V. de Mello, J. Phys.: Condens. Matter 32, 40LT02 (2020a).
  * de Mello (2020b) E. de Mello, J. Phys.: Condens. Matter 32, 38LT01 (2020b).
  * Ono and Ando (2003) S. Ono and Y. Ando, Phys. Rev. B 67, 104512 (2003).
  * Komiya  _et al._ (2002) S. Komiya, Y. Ando, X. F. Sun, and A. N. Lavrov, Phys. Rev. B 65, 214535 (2002).
  * Panagopoulos  _et al._ (2000) C. Panagopoulos, J. R. Cooper, T. Xiang, Y. S. Wang, and C. W. Chu, Phys. Rev. B 61, R3808 (2000).
  * Panagopoulos  _et al._ (1999) C. Panagopoulos, J. L. Tallon, and T. Xiang, Phys. Rev. B 59, R6635 (1999).
  * Emery and Kivelson (1995) V. J. Emery and S. A. Kivelson, Nature 374, 434 (1995).
  * Božović  _et al._ (2016) I. Božović, X. He, J. Wu, and A. T. Bollinger, Nature 536, 309 (2016).
  * Uemura  _et al._ (1989) Y. J. Uemura, _et al._ Phys. Rev. Lett. 62, 2317 (1989).
  * Ketterson and Song (1999) J. B. Ketterson and S. Song,  _Superconductivity_ (Cambridge University Press, London, 1999).
  * Wise  _et al._ (2009) W. D. Wise, K. Chatterjee, M. C. Boyer, T. Kondo, T. Takeuchi, H. Ikuta, Z. Xu, J. Wen, G. D. Gu, Y. Wang, and E. W. Hudson, Nature Physics 5, 213 (2009).
  * Comin  _et al._ (2014) R. Comin, A. Frano, M. M. Yee, Y. Yoshida, H. Eisaki, E. Schierle, E. Weschke, R. Sutarto, F. He, A. Soumyanarayanan, He Yang, M.Le Tacon, I.S Elfimov, E. Hoffman Jennifer, G.A Sawatzky, B. Keimer, A. Damascelli, Science (New York, N.Y.) 343, 390 (2014).
  * Bozovic  _et al._ (2002) I. Bozovic, G. Logvenov, I. Belca, B. Narimbetov, and I. Sveklo, Phys. Rev. Lett. 89, 107001 (2002).
  * Smadici  _et al._ (2009) S. Smadici, J. C. T. Lee, S. Wang, P. Abbamonte, G. Logvenov, A. Gozar, C. D. Cavellin, and I. Bozovic, Phys. Rev. Lett. 102, 107004 (2009).
  * Suter  _et al._ (2012) A. Suter, E. Morenzoni, T. Prokscha, H. Luetkens, B. Wojek, G. Logvenov, A. Gozar, and I. Božovič, Physics Procedia 30, 271 (2012).
  * da Silva Neto  _et al._ (2018) E. H. da Silva Neto, M. Minola, B. Yu, W. Tabis, M. Bluschke, D. Unruh, H. Suzuki, Y. Li, G. Yu, D. Betto,K. Kummer, F. Yakhou, N.B Brookes, M. Le Tacon, M. Greven, B. Keimer, A. Damascelli, Phys. Rev. B 98, 161114 (2018).


11institutetext: Max-Planck-Institut für Sonnensystemforschung, Justus-von-
Liebig Weg 3, 37077 Göttingen 22institutetext: School of Space Research, Kyung
Hee University, Yongin, Gyeonggi-Do, 446-701, Korea

# Small-scale Dynamo in Cool Stars

I. Changes in stratification and near-surface convection for main-sequence
spectral types

Tanayveer S. Bhatia, 11<EMAIL_ADDRESS>Robert H. Cameron 11 Sami K.
Solanki, 1122 Hardi Peter 11 Damien Przybylski 11 Veronika Witzke 11 Alexander
Shapiro 11

(Submitted: 22 Mar 2022, Proposed for acceptance: 31 May 2022)

###### Abstract

Context. Some of the small-scale solar magnetic flux can be attributed to a
small-scale dynamo (SSD) operating in the near-surface convection. The SSD
fields have consequences for solar granular convection, basal flux, as well as
chromospheric heating. A similar SSD mechanism is expected to be active in the
near-surface convection of other cool main-sequence stars, but this has never
been investigated.

Aims. We aim to investigate changes in stratification and convection due to
inclusion of SSD fields for F3V, G2V, K0V and M0V spectral types in the near-
surface convection.

Methods. 3D magnetohydrodynamic (MHD) models of the four stellar boxes,
covering the subsurface convection zone up to the lower photosphere in a small
cartesian box, are studied using the MURaM radiative-MHD simulation code. The
SSD runs are compared against reference hydrodynamic runs.

Results. An SSD is found to efficiently produce magnetic field with energies
ranging between 5% to 80% of the plasma kinetic energy at different depths.
This ratio tends to be larger for larger $T_{\mathrm{eff}}$. The relative
change in density and gas pressure stratification for the deeper convective
layers due to SSD magnetic fields is negligible, except for the F-star. For
the F-star, there is a substantial reduction in convective velocities due to
Lorentz force feedback from magnetic fields, which, in turn, reduces the
turbulent pressure.

Conclusions. SSD in near-surface convection for cool main-sequence stars
introduces small but significant changes in thermodynamic stratification
(especially for the F-star) due to reduction in convective velocities.

###### Key Words.:

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

##  1 Introduction

The interpretation of data from stellar observations requires a stellar model
to compare against. Traditionally, these models have been 1D global models
(Carbon & Gingerich, 1969) that use formulations of mixing-length theory (MLT)
(Böhm-Vitense, 1958). Later models accounted for line-blanketing effects and
used an opacity distribution function (ODF) approach to calculate opacities
(Strom & Kurucz, 1966). Among others, the MARCS code (Gustafsson et al., 1975,
2008), the ATLAS code (Kurucz, 1979; Castelli & Kurucz, 2003) and the PHOENIX
code (Allard & Hauschildt, 1995; Hauschildt et al., 1999) have enabled the
calculation of synthetic stellar spectra with a detailed accounting for the
relevant physics. These models enabled, for example, the accurate
determination of abundances, stellar evolution tracks, and even for
constraining the chemical evolution of galaxies (Edvardsson et al., 1993).

However, convection is a 3D process. A phenomenological 1D approach is
insufficient to characterize properties of granulation and plasma turbulence.
Radiative aspects of granulation, in particular, are not fully captured by an
MLT-like approach. 1D models also require free parameters, like the mixing
length parameter $\alpha_{MLT}$ and turbulent velocities for line broadening
calculations. Hence, 3D stellar atmosphere models become important for a
realistic interpretation of stellar characteristics from observations. The
first 3D simulations of solar convection were pioneered by Nordlund (1982);
Stein & Nordlund (1998); Nordlund & Stein (1990). These simulations are
realistic in the sense that they are directly comparable to solar
observations: they reproduce granulation pattern and associated spectral line
widths and asymmetries quite well (Asplund et al., 2000).

Early stellar models (Nordlund & Dravins, 1990) showed the limitations of MLT-
based models in accurately reproducing the near-surface temperature gradient,
which affects radiative properties in the lower photospheres. Subsequently,
various models have been constructed over a grid of effective temperature,
surface gravity and metallicity. Currently, the most comprehensive 3D grids
(Magic et al., 2013; Ludwig et al., 2009) cover a wide range of stellar type
on the HR-diagram, but they are purely hydrodynamic.

However, stellar convection is not a purely hydrodynamic process. Most cool
stars are expected to have magnetic fields. Hence, a complete description of
their photospheres should also take into account the effects of such fields.
The best studied star in this context is the Sun. There is a rich variety of
solar magnetic field-related phenomena ranging from sunspots and active
regions to network fields, forming plages and faculae all the way down to
small-scale mixed polarity turbulent magnetic field filling the rest of the
solar surface. State-of-the-art solar simulations reproduce all of these
features, from sunspots (Rempel et al., 2009) and plages (Vögler et al., 2005;
Yadav et al., 2021) to the quiet-sun magnetism (Rempel, 2014).

The quiet-sun small-scale field, that is, the field associated with regions of
the sun not showing any activity, is partly attributable to a small-scale
dynamo (SSD) operating in the convection zone (Vögler & Schüssler, 2007;
Rempel, 2014; Pietarila Graham et al., 2010). This field is believed to have a
significant magnitude of around $\sim$ 130 G, based on Hanle depolarization
(Trujillo Bueno et al., 2004). Additional evidence that a fair fraction of the
small-scale field is a result of an SSD comes from the fact that internetwork
magnetic flux does not follow the solar cycle (Buehler et al., 2013; Lites et
al., 2014). In realistic radiative-MHD simulations, the effect of quiet-sun
magnetic fields (self consistently generated via an SSD mechanism) has
previously been shown to be important, e.g. to reproduce the correct solar
intensity contrast (Danilovic et al., 2010) and account for inferred
photospheric magnetic field strength based on Hanle-effect diagnostics
(Shchukina & Trujillo Bueno, 2011). In addition, one often obtains a rough
equipartition between kinetic and magnetic energy in SSD simulations (Hotta et
al., 2015; Haugen et al., 2004; Schekochihin et al., 2004), which implies a
substantial reduction in plasma velocities since most of the energy in the
magnetic fields is obtained from the plasma motions. The importance of
magnetic fields generated from an SSD in other stellar types and its effect on
the intensity characteristics, however, remains yet to be explored. Hence, it
becomes imperative to investigate the effect of quiet-star small-scale
magnetic fields on these quantities.

This paper is a part of a project aimed at constructing a grid of magneto-
convective stellar atmospheres, ranging in temperature ($3500<T<7000$ K),
gravity ($4.3<\log_{10}g<4.8$ in cgs units) and metallicity (in this paper,
only solar metallicities are considered). We investigate four stellar cases:
F3V, G2V, K0V and M0V, and study the relative change in stratification,
convection and intensity from purely hydrodynamic setups.

In §2, we outline the simulation code and the setup. Then, we present the
results of the simulations in §3, followed by interpretation of the results in
§4. Lastly, we summarize the results and corresponding discussion in §5.

##  2 Methods

###  2.1 Simulation Code

The code we use throughout this work is MURaM (Vögler et al., 2005; Rempel,
2014, 2017). It is a 3D radiative-MHD code that solves the conservative MHD
equations for compressible, partially ionized plasma. It uses a multi-group
radiative transfer scheme with short characteristics (Nordlund, 1982). The
equations for mass ($\rho$ \- density), momentum ($\@vec{v}$ \- velocity, $p$
\- pressure) and energy ($\epsilon_{h}$ \- enthalpy density) conservation are
solved, along with the induction equation ($\@vec{B}$ \- magnetic field)

$\displaystyle\partial_{t}\rho=$ $\displaystyle-\nabla\cdot(\rho\@vec{v})$ (1)
$\displaystyle\partial_{t}(\rho\@vec{v})=$
$\displaystyle-\nabla\cdot(\rho\@vec{v}\@vec{v})-\nabla
p+\rho\@vec{g}+\@vec{F}_{\rm{SR}}+\@vec{F}_{\rm{L}}$ (2)
$\displaystyle\begin{split}\partial_{t}(\epsilon_{h}+\rho
v^{2}/2)=&-\nabla\cdot(\@vec{v}(\epsilon_{h}+p+\rho v^{2}/2))+\\\
&\@vec{v}\cdot(\@vec{g}+\@vec{F}_{\rm{L}}+\@vec{F}_{\rm{SR}})+Q_{\rm{rad}}+Q_{\rm{res}}\end{split}$
(3) $\displaystyle\partial_{t}\@vec{B}=$
$\displaystyle\nabla\times(\@vec{v}\times\@vec{B})$ (4)

Here $\@vec{F}$ refers to forces and $Q$ refers to source terms. The subscript
$\rm{SR}$ refers to semi-relativistic ”Boris correction”-related terms (Boris,
1970; Gombosi et al., 2002) (negligible for our setups) and $\rm{L}$ refers to
Lorentz force. The two $Q$ terms in the energy equation account for the
radiative heating/cooling and resistive heating (since the hydrodynamic energy
is conserved, instead of the total energy). For details, we refer the reader
to Rempel (2014, 2017). In this work, the grey approximation is used for
solving the radiative transfer equations, where the frequency dependence of
the opacity is replaced by an average value. This is an acceptable
approximation for this work since we are mainly interested in the structure
below and just above the surface (Vögler et al., 2004). Finally, the FreeEOS
equation of state (Irwin, 2012) with solar abundances (Asplund et al., 2009)
is used to close the set of equations.

The effective temperature ($T_{\mathrm{eff}}$) (related to the radiative
output), the surface gravitational acceleration ($g$) (related to the
hydrostatic balance) and the metallicity ($Z$) uniquely specify the spectral
class of a star. The MURaM code uses a constant gravitational acceleration
$g$, the gas pressure at bottom boundary ($p_{\mathrm{bot}}$), and the inflow
entropy at the bottom boundary ($s_{\mathrm{bot}}$) as free parameters. The
$p_{\mathrm{bot}}$ and $s_{\mathrm{bot}}$ parameters determine the height of
the $\tau=1$ surface and the $T_{\mathrm{eff}}$.

###  2.2 Setup and Parameters

We consider four stellar cases: F3V, G2V, K0V, and M0V. These cases were
chosen to cover a broad range of $T_{\mathrm{eff}}$ for stellar types with
convective envelopes. All boxes have the same number of gridpoints (512
$\times$ 512 in the horizontal direction and 500 in the vertical direction).
The scaling for the horizontal and the vertical geometric extent is done such
that the number of granules in each box is roughly the same and the number of
pressure scale heights below the photosphere is also similar ($\sim 7.5$). For
the G-star, this corresponds to 4 Mm below the surface and a horizontal extent
of 9 Mm $\times$ 9 Mm.

Table 1: Parameters for the simulation setup

Type |  $z_{\downarrow}$**$*$These quantities are determined after running the simulations. Note that the change in $T_{\mathrm{eff}}$ will not influence the total radiative output on long timescales ($>10^{5}$ yr) corresponding to the Kelvin-Helmholtz timescale. See Spruit (1982) for details.$(z_{\mathrm{tot}})$ | $x_{0},y_{0}$ | $\Delta x,y$ | $\Delta z$ | $\log_{10}g$ |  $(T_{\mathrm{eff}})_{HD}$**$*$These quantities are determined after running the simulations. Note that the change in $T_{\mathrm{eff}}$ will not influence the total radiative output on long timescales ($>10^{5}$ yr) corresponding to the Kelvin-Helmholtz timescale. See Spruit (1982) for details. |  $(T_{\mathrm{eff}})_{\mathrm{SSD}}$**$*$These quantities are determined after running the simulations. Note that the change in $T_{\mathrm{eff}}$ will not influence the total radiative output on long timescales ($>10^{5}$ yr) corresponding to the Kelvin-Helmholtz timescale. See Spruit (1982) for details. |  $\langle||B||\rangle_{\tau=1}$**$*$These quantities are determined after running the simulations. Note that the change in $T_{\mathrm{eff}}$ will not influence the total radiative output on long timescales ($>10^{5}$ yr) corresponding to the Kelvin-Helmholtz timescale. See Spruit (1982) for details. |  $\langle|B_{z}|\rangle_{\tau=1}$**$*$These quantities are determined after running the simulations. Note that the change in $T_{\mathrm{eff}}$ will not influence the total radiative output on long timescales ($>10^{5}$ yr) corresponding to the Kelvin-Helmholtz timescale. See Spruit (1982) for details. |  $\langle||B||\rangle_{\langle\tau\rangle=1}$**$*$These quantities are determined after running the simulations. Note that the change in $T_{\mathrm{eff}}$ will not influence the total radiative output on long timescales ($>10^{5}$ yr) corresponding to the Kelvin-Helmholtz timescale. See Spruit (1982) for details.
---|---|---|---|---|---|---|---|---|---|---
| (Mm) | (Mm) | (km) | (km) | (cm/s2) | (K) | (K) | (G) | (G) | (G)
F3V | 11.11 (13.00) | 23 | 45 | 26 | 4.301 | 6817$\pm$7 | 6807$\pm$8 | 188$\pm$15 | 93$\pm$8 | 132$\pm$10
G2V | 4.09 (5.00) | 9 | 17.5 | 10 | 4.438 | 5834$\pm$9 | 5840$\pm$9 | 127$\pm$12 | 66$\pm$6 | 113$\pm$10
K0V | 2.05 (2.31) | 4.62 | 8.2 | 4.62 | 4.609 | 4668$\pm$5 | 4671$\pm$5 | 103$\pm$6 | 50$\pm$3 | 103$\pm$6
M0V | 0.90 (1.14) | 2.043 | 3.99 | 2.28 | 4.826 | 3825$\pm$1 | 3827$\pm$2 | 106$\pm$6 | 54$\pm$4 | 107$\pm$6

111

The boundaries are periodic in the horizontal $x,y$, direction. The top
boundary ($z_{\mathrm{top}}$) is open to outflows and closed to inflows, with
vertical magnetic fields. The bottom boundary ($z_{\mathrm{bot}}$) is
symmetric222This refers to the way the derivative is handled across the ghost
cells. Symmetric boundary implies the same value in the ghost cell ($q_{g}$)
next to the boundary domain cell ($q_{d}$), such that the derivative across
the boundary is zero ($q_{g}=q_{d}$), and anti-symmetric implies a value with
the opposite sign ($q_{g}=-q_{d}$) for mass flux ($\rho\@vec{v}$), entropy
downflows and magnetic fields. Note that the magnetic field boundary condition
also allows horizontal field to be advected across the bottom boundary. This
’mimics’ the presence of magnetic field deeper in the convection zone, as
considered previously based on equipartition arguments (Rempel, 2014; Hotta et
al., 2015). The magnetic field BC may not necessarily preserve the
$\nabla\cdot\@vec{B}=0$ constraint. However, the hyperbolic divergence
cleaning approach (Dedner et al., 2002) employed in MURaM takes care of the
$\nabla\cdot\@vec{B}$ errors reasonably well
($(\nabla\cdot\@vec{B})_{\mathrm{rms}}/(||B||/\Delta z)\sim O(10^{-3})$, see
Fig. 7).

For each star, we performed simulations with magnetic fields (SSD), and purely
hydrodynamic (HD) simulations. Purely hydrodynamic simulations were
initialized with density and internal energy profiles generated using the 1D
stellar code MESA (Paxton et al., 2019) for the F-, K- and M-star and using
the standard solar model from Christensen-Dalsgaard et al. (1996) for the
G-star. These were then run for several hours in stellar time till convection
became relaxed and there were no transients visible in velocity and pressure
vertical slices. Then the simulation box was seeded with net zero-flux
magnetic field with a negligibly small field strength ($10^{-5}$ G) and run
till photospheric magnetic field strength reached saturation. The results
presented in the subsequent sections are averaged over a few hours of stellar
time (after saturation), and over a number of snapshots, for all the eight
cases (see Table 2 for further details for each setup).

Table 1 describes the detailed setup for all the simulations: for all four
stellar types, it gives the height of the $\tau=1$ surface above the bottom
boundary (and the total vertical extent), the horizontal extent, the
horizontal resolution, the vertical resolution, the log constant surface
gravity, the effective temperature of the SSD and HD cases, the average
magnetic field magnitude at the $\tau=1$ iso-surface333This is the surface
corresponding to where $\tau=1$ in each vertical column of the simulation
cube, and this is usually somewhat corrugated because downflows are cooler
than upflows and opacity is extremely sensitive to temperature in the relevant
ranges. $\langle||B||\rangle_{\tau=1}$, the average unsigned vertical field at
the $\tau=1$ iso-surface $\langle|B_{z}|\rangle_{\tau=1}$, and the average
magnetic field magnitude at the $\langle\tau\rangle=1$ horizontal slice
$\langle||B||\rangle_{\langle\tau\rangle=1}$ of the SSD cases. The effective
temperature is calculated by averaging the angle-averaged bolometric
luminosity over time.

##  3 Results

Table 2: Color-coding for the plots, along with number of snapshots considered
and total time in hours for each stellar box

Star |  | SSD |  |  | HD |
---|---|---|---|---|---|---
| Color | $N$ |  $t$ (h) | Color | $N$ |  $t$ (h)
F3V | Blue | 264 | 12.4 | Light Blue | 124 | 14.8
G2V | Black | 178 | 7.7 | Gray | 117 | 15.0
K0V | Green | 125 | 10.6 | Lime | 193 | 16.8
M0V | Red | 263 | 5.4 | Light Red | 197 | 10.1
444All HD plots are dashed lines.

Figure 1: Emergent intensity and surface vertical velocities in different
stellar types for models with and without magnetic field. From top to bottom:
Snapshot of the bolometric intensity and $v_{z}$ at $\tau=1$ for the HD case
(row 1 and 2), bolometric intensity and $v_{z}$ at $\tau=1$ for the SSD case
(row 3 and 4) and the corresponding vertical magnetic field at $\tau=1$ (row
5, from left to right) for spectral types F, G, K and M, respectively. The
green circles indicate the bright points and corresponding magnetic field
concentrations.

All the magnetic simulations exhibit dynamo action and develop magnetic fields
with energy within an order of magnitude of the kinetic energy through most of
the simulation box. The change in the partition of energy influences the
internal structure as well as convective velocities. Snapshots of the
bolometric intensity and vertical velocity at the $\tau=1$ surface for the SSD
and the HD setups (along with the vertical magnetic field for the SSD cases)
are shown in Fig. 1. The SSD cases show distinct intergranular bright points
which correspond with strong magnetic field concentrations. In the following
subsections, we examine the horizontally averaged structure of the magnetic
field and its effects on the stratification as well as convection for these
stellar types. The analysis of the magnetic fields in the lower photosphere
and their effect on the bolometric intensity and vertical velocity will be
covered in the next paper in this series.

All 1D plots are averages over a number of snapshots spanning a few hours of
stellar time (see Table 2 for exact numbers). The error bars are standard
error (standard deviation normalized by the square root of the number of
snapshots $\sigma/\sqrt{N}$) of the average 1D structure, the assumption being
that over this time span, the snapshots are statistically independent.

All quantities are plotted as a function of number of pressure scale heights
relative to the height where $\langle\tau\rangle=1$,
$n_{H}=\log(p_{\mathrm{gas}}/p_{\mathrm{gas}(\langle\tau\rangle=1)})$ (with
this definition, positive values correspond to the interior). Since the non-
magnetic bottom boundary conditions are identical ($p_{\mathrm{bot}}$ and
$s_{\mathrm{bot}}$ are the same) for HD and SSD runs, the deviations are
calculated geometrically and plotted against the corresponding HD pressure
scale axis.

This extent in terms of pressure scale heights ranges from 7.5 (bottom) to -5
(near top) for all the simulations. Since the focus of this paper is on the
near surface convection zone, we exclude from our analysis the portion of the
box corresponding to $n_{H}<-1$. We also exclude the region corresponding to
$n_{H}>6$, to possible numerical bottom boundary effects.

###  3.1 Magnetic field structure

Figure 2: Top: Horizontally averaged magnetic field magnitude. Bottom: The
ratio of the horizontal r.m.s field strength to the vertical r.m.s field
strength. The horizontal axis is the number of pressure scale heights
$\log_{10}(p_{\mathrm{gas}}/p_{\mathrm{gas}(\tau=1)})$, calculated for the HD
cases, below the surface (dotted vertical black line). The shaded regions
correspond to 1-$\sigma$ standard error ($e=\sigma/\sqrt{N}$, $N$ is the
number of snapshots) computed over time averaging of snapshots.

As mentioned in the introduction to this section, the dynamo action results in
a significant amount of magnetic field, with the overall magnitude roughly
similar for all the cases, and with a somewhat decreasing trend with
$T_{\mathrm{eff}}$ near the surface (Fig. 2, top panel and last column of
Table 1). The relation of the magnetic energy to kinetic energy is discussed
in §3.3.

The bottom panel of Fig. 2 shows the ratio $B_{h}/B_{z}$, which gives an
indication of the 3D structure of the magnetic fields. For a fully isotropic
distribution of magnetic field, one would expect $B_{x}^{2}\approx
B_{y}^{2}\approx B_{z}^{2}$. This implies that the ratio of the horizontal
r.m.s. component of the magnetic field fluctuations $B_{\mathrm{h,rms}}$ to
the vertical component $B_{\mathrm{z,rms}}$ should be $\sim\sqrt{2}$. For G-,
K- and M-stars, this ratio is slightly less than $\sqrt{2}$ in the middle of
the box ($2<n_{H}<5$), indicating near-isotropy, whereas for the F-star, it is
significantly lower ($\sim$ 1).

Near the $\tau=1$ surface, $B_{\mathrm{h,rms}}/B_{\mathrm{z,rms}}$ for all the
stars is lower because of intensification of vertical magnetic fields in the
intergranular lanes (e.g. Spruit (1979)).

###  3.2 Changes in Stratification

Figure 3: Deviations in $\rho$, $T$ and $p_{\mathrm{gas}}$ for F (blue), G
(black), K (green) and M-star (red) cases. The vertical axis gives the
geometric deviations as a percentage relative to the hydrodynamic case. The
horizontal axis is the number of pressure scale heights
$\log_{10}(p_{\mathrm{gas}}/p_{\mathrm{gas}(\tau=1)})$, calculated for the HD
cases, below the surface (dotted vertical black line). The shaded regions
correspond to 1-$\sigma$ standard error ($e=\sigma/\sqrt{N}$, $N$ is the
number of snapshots) of the mean solid curve.

Fig. 3 shows the plots of deviation from the mean HD stratification. The
deviations from the HD simulations are presented as the relative percent
change in the horizontally averaged 1D structure. For any quantity of interest
$q$ (e.g, density, temperature etc.), these deviations are calculated as
$(q_{\mathrm{SSD}}-q_{\mathrm{HD}})/q_{\mathrm{HD}}$. This means that a
positive value for the deviation corresponds to a higher value for the SSD
case relative to the HD case.

All simulations show slight ($\leq$ 2%) changes in thermodynamic
stratification relative to the corresponding HD simulations, with the
magnitude of deviations below the surface roughly increasing from coolest to
hottest stellar type. Below 3 pressure scale heights, the G, K and M
simulations show negligible $(\leq 0.1\%)$ deviations in thermodynamic
quantities.

Closer to the surface, the F-star simulation shows up to 1.5% reduction in
density and up to 1% reduction in gas pressure (Fig. 3, top and bottom panel,
blue line). This trend is opposite to that seen for other cases, all of which
show a slight ($<1\%$) enhancement in density and pressure. These results are
analyzed in §4.1. For context, a 1.5% deviation in temperature for an F-star
would correspond to a temperature change of $\sim$100 K.

###  3.3 Distribution of energies

Figure 4: a) Ratio of magnetic to kinetic energy, b) percent change in kinetic
energy, c) ratio of kinetic to internal energy, and d) ratio of magnetic to
internal energy (bottom) for F-, G-, K- and M-stars.

In Fig. 4, the panel a) shows the ratio of magnetic to kinetic energy. For all
cases, this ratio is within an order of magnitude throughout the box. This
ratio for the F-star is significantly higher compared to the other stars. In
addition, the ratio for M-star has a minima deeper down in the box than for
other stars. As before, the general trend shows a decrease in this ratio with
decreasing $T_{\mathrm{eff}}$.

Fig. 4b shows the relative change in kinetic energy between the SSD and the HD
setups. All stars show a marked decrease in kinetic energy in the SSD case. As
with the top panel, F-star shows a significantly stronger reduction in kinetic
energy compared to the other cases. The M-star as well shows a stronger
reduction in KE near the surface compared to the G and K-stars.

Fig. 4c shows the ratio of kinetic to internal energy. For the F-star, this
ratio is within an order of magnitude near the surface, whereas for the other
stars, it is less than $1\%$.

Fig. 4d shows the ratio of the magnetic to the internal energy. There is a
clear trend in this ratio with stellar type, with it being lowest for M-star
($\sim 5\times 10^{-5}$ near bottom, to $\sim 10^{-4}$ near surface) and
highest for F-star ($\sim 10^{-3}$ near bottom, to $\sim 10^{-2}$ near
surface).

###  3.4 Changes in Velocities

Figure 5: Relative decrease in convective velocities $v_{z}^{2}$ (top),
relative decrease in ratio of horizontal to vertical r.m.s velocities
$v_{\mathrm{h,rms}}/v_{\mathrm{z,rms}}$ (middle), and the actual
$v_{\mathrm{h,rms}}/v_{\mathrm{z,rms}}$ ratio (bottom).

All SSD cases exhibit a decrease in vertical velocities $v_{z}^{2}$ as well as
the ratio of horizontal to vertical r.m.s. velocities
$v_{\mathrm{h,rms}}/v_{\mathrm{z,rms}}$ (Fig. 5), relative to the
corresponding HD cases. The reduction in $v_{z}^{2}$ follows similar trend for
all the four stars, with a decrease of 4-8% near the surface, going up to 20%
near the bottom boundary. For G-, K- and M-star, the decrease in
$v_{\mathrm{h,rms}}/v_{\mathrm{z,rms}}$ is similar (5% near the surface, going
up to 12% near the bottom) but is more pronounced for the F-star case (10%
near the surface, going up to 20% near the bottom boundary). This implies a
change in the horizontal extent of sub-surface granulation, which follows
Nordlund et al. (2009), who showed using simple mass conservation that the
horizontal extent of granules is proportional to $H(v_{h}/v_{z})$, where $H$
is the local density scale height.

The ratio $v_{\mathrm{h,rms}}/v_{\mathrm{z,rms}}$ gives an idea of the 3D
velocity structure. Close to the $\tau=1$ surface, the ratio increases
suddenly. This increase corresponds to where convective flows turn over, as
the atmosphere becomes convectively stable and the flows above that point are
mainly due to convective overshoot. The exact pressure scale depth where this
turning over takes place depends of the effective temperature: For the F- and
G-star, this turnover takes place within half a pressure scale height of the
$\tau=1$ surface whereas, for the M-star, it takes place well below the
$\tau=1$ surface (around $n_{H}=3$). This is also probably why the minima of
the ME/KE ratio in Fig. 4 for M-star is significantly below the surface
compared to the G and K-star, as density (and, consequently, KE) is lower
above the height where most of the overturning takes place.

The changes in $v_{z}^{2}$ as well as the
$v_{\mathrm{h,rms}}/v_{\mathrm{z,rms}}$ with depth seem to follow a regular
trend, and is consistent with the results obtained previously in whole-
convection zone simulations with a small-scale dynamo set up (Hotta et al.,
2015).

##  4 Discussion

###  4.1 Turbulent pressure

Figure 6: Hydrostatic force balance with gas pressure $p_{\mathrm{gas}}$ (top)
and total pressure $p_{\mathrm{tot}}$ (middle) gradient term and $\rho g$
term, normalized by $\langle\rho g\rangle$ as $\langle dp/dz+\rho
g\rangle/\langle\rho g\rangle$, and ratio of turbulent pressure to gas
pressure (bottom). All quantities are plotted for SSD (solid) and HD (dashed)
cases. Note that the extent of $y$-axis is larger by more than an order of
magnitude for force balance with $p_{\mathrm{gas}}$ compared to force balance
with $p_{\mathrm{tot}}$

We first discuss the deviations of the pressure and density in the SSD models
from a purely HD model. The pressure and density deviations (Fig. 3) are very
small ($\sim 0.1\%$) for the G-, K- and M-star in the convection zone (and not
significantly larger near the surface for the K- and M-star). The deviations
are more prominent for the F-star ($\sim 1\%$). This discrepancy can be
understood by considering the contribution of turbulent pressure in the
overall hydrostatic balance.

Turbulent pressure becomes important for hydrostatic balance when kinetic
energy is within an order of magnitude or so of the internal energy. This is
indeed the case for our F-star models. This implies that turbulent pressure
(which is largely due to plasma motions) can be a significant fraction of the
gas pressure $p_{\mathrm{gas}}$. Quantitatively, this can be seen from a crude
MLT calculation of the Mach number (see appendix A for derivation)

$\displaystyle M$ $\displaystyle\approx 0.138\left(\frac{T}{\text{10${}^{3}$
K}}\right)^{5/6}\left(\frac{\rho}{\text{10${}^{-7}$
g/cm${}^{3}$}}\right)^{-1/3}.$ (5)

Using the above equation, we calculate the F-star photospheric Mach number to
be about 0.75 (see Table 3), slightly lower than simulation value in Beeck et
al. (2013) who found it to be 0.9. For the other stars, the velocity is
decidedly subsonic.

Table 3: Mach numbers computed using MLT

Simulation |  $T_{\mathrm{surf}}$ (103 K) |  $\rho_{\mathrm{surf}}$ (10-7 g/cm3) | $M$
---|---|---|---
F3V | 6.19 | 0.59 | 0.75
G2V | 6.17 | 2.39 | 0.47
K0V | 4.94 | 6.67 | 0.28
M0V | 3.96 | 21.3 | 0.16
555$\rho_{\mathrm{surf}}$ and $T_{\mathrm{surf}}$ are obtained from HD
simulation data

Now, the hydrostatic balance is expressed in terms of the balance between
force due to pressure gradient $p^{\prime}=dp/dz$ and gravity $\rho g$, acting
in the vertical (viz. radial) direction. These two terms should be
approximately equal in magnitude. If just $p_{\mathrm{gas}}$ is considered,
this balance does not hold very well, with deviations increasing strongly with
$T_{\mathrm{eff}}$ (Fig. 6, top panel).

The turbulent pressure consists of terms from the total stress tensor. From
Reynold stresses, the $\rho v_{i}v_{j}$ term and from Maxwell stresses, the
$((B^{2}/2)\delta_{ij}-B_{i}B_{j})/4\pi$ term is obtained. As mentioned in the
introduction, the presence of SSD magnetic fields implies a reduction in
kinetic energy, where energy is redistributed between the plasma motions and
magnetic fields via Lorentz force feedback. In Fig. 4, the ratio of magnetic
to kinetic energy (top panel) corresponds remarkably well with the reduction
in kinetic energy relative to the HD case (middle panel), implying most of the
energy in the magnetic fields is drawn from the kinetic energy reservoir. The
magnetic energy is within an order of magnitude of the kinetic energy for
subsurface plasma. This is consistent with the results on equipartition of
energy in several SSD simulations (Hotta et al., 2015; Haugen et al., 2004;
Schekochihin et al., 2004).

This order of magnitude equipartition results in a reduction of plasma
velocities, and, consequently, the magnitude of the Reynolds stress term. The
contribution from the Maxwell stress, on the other hand, can be either
negative or positive, depending on whether the effects of magnetic tension
dominate over that of magnetic pressure in the vertical direction. More
quantitatively (see appendix B for a derivation), the total turbulent pressure
can be expressed as

$\displaystyle p_{\mathrm{turb}}$ $\displaystyle=\rho
v_{z}^{2}+\frac{B_{h}^{2}-B_{z}^{2}}{8\pi}.$ (6)

With this included in the calculation of the pressure gradient term,
hydrostatic balance is satisfied, as can be seen in middle panel of Fig. 6
(note the difference in the extent of $y$-axis for the top and middle panels).
In the bottom panel of the same figure, the reduction in
$p_{\mathrm{turb}}/p_{\mathrm{gas}}$ is most prominent for the F-star. Since
the changes in $v_{z}^{2}$ are similar (Fig. 5, top panel) for all stars,
there must be a reduction in density for the F-star to compensate for the
significant change in $p_{turb}$.

Based on the expression for $p_{\mathrm{turb}}$, we introduce an effective
turbulent velocity $v=\sqrt{p_{\mathrm{turb}}/\rho}$

With this, it becomes possible to relate changes in density stratification to
changes in total pressure gradient. A straightforward consideration of the
density scale height $H_{\rho}$ (see appendix B) yields

$\displaystyle H_{\rho}$ $\displaystyle=RT/(\mu g)+v^{2}/g$ (7)

This tells us that where $B_{h}^{2}<B_{z}^{2}$ and $(\rho
v_{z}^{2})_{\mathrm{SSD}}<(\rho v_{z}^{2})_{\mathrm{HD}}$, the scale height of
the SSD model is smaller than for the HD model,
$(H_{\rho})_{\mathrm{SSD}}<(H_{\rho})_{\mathrm{HD}}$. This is valid if
assuming that the change in $T$ and $\mu$ is relatively small. From the
proportionality between density scale height and $v_{z}^{2}$ (Eq. 7), the
decrease in vertical velocity is associated with a decrease in the local
density scale height. This can be inferred from Fig. 6, bottom panel, where
the ratio $p_{\mathrm{turb}}/p_{\mathrm{gas}}$ ($\propto v^{2}/(RT/\mu)$) is
noticeably lower for F-star near the surface.

These conditions are satisfied relatively well for the F-star case. A decrease
in $H_{\rho}$ implies a steepening of the density stratification for the SSD
case relative to the HD case. This is exactly the case in the top panel of
Fig. 3: the reduction in density goes from a $\sim$0% to $\sim 1.5$% near the
surface. Below the surface, it is reasonably good to assume a perfect gas
equation of state, with $p_{\mathrm{gas}}\propto\rho T$. Assuming changes in
$T$ are small compared to $\rho$ and $p_{\mathrm{gas}}$, this, in turn,
implies that the changes in $p_{\mathrm{gas}}$ follows a similar trend as
changes in $\rho$. Near the surface, this picture breaks down since the mode
of energy transfer changes from convective to radiative and this simplified
analysis is no longer valid.

The effect of turbulent pressure on stratification and convection has
previously been considered for stellar model envelopes from an MLT perspective
(Henyey et al., 1965) as well as for 3D HD simulations (Ludwig & Kučinskas,
2012; Jørgensen & Weiss, 2019), but, to the best of our knowledge, the effect
of magnetic fields from a SSD have not been considered before.

###  4.2 Changes in velocity structure

As mentioned in §3.4, the ratio $v_{\mathrm{h,rms}}/v_{\mathrm{z,rms}}$ (Fig.
5, bottom panel), gives an idea of where convection overturns relative to
pressure scale height for a given effective temperature. The trend in
overturning for different spectral types fits the usual picture of ”hidden”
granulation below the $\tau=1$ surface for cooler, denser stars like M and K
vs. the ”naked” granulation above the $\tau=1$ surface of hotter, more
rarefied stars like G and F, as discussed in Nordlund & Dravins (1990).

This ratio is close to 1 for the HD G-, K- and M-stars (dashed lines). For the
HD F-star, however, it is significantly higher (dashed light blue line),
indicating a higher degree of isotropy in the velocity structure666Fully
isotropic flow requires $v_{x}^{2}=v_{y}^{2}=v_{z}^{2}$. Hence, this would
imply $v_{\mathrm{h,rms}}/v_{\mathrm{z,rms}}=\sqrt{2}$. This implies
relatively higher horizontal velocities for the HD F-star. The nearly-
isotropic velocity profile for the HD F-star can be attributed to the kinetic
energy being a non-negligible fraction of the internal energy. We speculate
that these stronger horizontal velocities contribute to the stronger magnetic
fields in the vertical direction for the F-star (as inferred from the lower
$B_{\mathrm{h,rms}}/B_{\mathrm{z,rms}}$ for the F-star in Fig. 2, bottom
panel).

In the case of SSD (solid lines), all stars follow a similar trend for
$v_{\mathrm{h,rms}}/v_{\mathrm{z,rms}}$, with the ratio being $\leq 1$ up to
the point where convection overturns. The change in velocity field structure
in the presence of SSD fields is a hard problem that depends on the SSD
saturation mechanism. Note that the $B_{\mathrm{h,rms}}/B_{\mathrm{z,rms}}$
ratio is lower for F-star compared to other cases. Perhaps the relatively
stronger fields in the vertical direction restrict horizontal flows more for
F-star and lead to a larger change in $v_{\mathrm{h,rms}}/v_{\mathrm{z,rms}}$.
An understanding of this behavior requires a more detailed analysis which is
beyond the scope of this paper.

##  5 Conclusions

We have investigated the magnetic field self-consistently generated by an SSD
acting in the near-surface layers of main-sequence stars of spectral types
F3V, G2V, K0V and M0V. The SSD mechanism operates in all cases to amplify
magnetic fields from a seed field of negligible strength and zero net flux.
The magnetic fields from the SSD have an energy density that is a non-
negligible fraction of the kinetic energy density. These fields act back on
the plasma to reduce the convective velocities, which in turn reduces the
turbulent pressure. This becomes substantial for the F-star as it is hot
enough to have kinetic and internal energy within an order of magnitude near
the surface, which gives magnetic fields stronger than those in G-, K- and
M-stars, especially in the vertical direction. The equation for hydrostatic
balance for total pressure and the reduction of convective velocities implies
a reduction in the density scale height itself. This is significant enough for
the F-star to result in reduced density and gas pressure throughout the box.
This effect tends to get smaller towards later spectral types.

This paper only covers the near-surface convection zone. Of particular
interest to observational studies would be the magnetic field structure in the
lower photosphere and changes in the intensity characteristics. In addition,
the change in scale height and the changes in $v_{h}/v_{z}$ imply changes in
granulation scale. All these points are investigated in a subsequent paper.

###### Acknowledgements.

TB ran the simulations, performed the analysis and wrote most of the text. DP
helped with running the simulations. All authors contributed to the final
version of the text. 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. We also acknowledge Alan Irwin for providing the open source
FreeEOS equation of state. 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 and 715947).

## References

  * Allard & Hauschildt (1995) Allard, F. & Hauschildt, P. H. 1995, ApJ, 445, 433
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, Annual Review of Astron and Astrophys, 47, 481
  * Asplund et al. (2000) Asplund, M., Nordlund, Å., Trampedach, R., Allende Prieto, C., & Stein, R. F. 2000, A&A, 359, 729
  * Beeck et al. (2013) Beeck, B., Cameron, R. H., Reiners, A., & Schüssler, M. 2013, A&A, 558, A48
  * Böhm-Vitense (1958) Böhm-Vitense, E. 1958, ZfA, 46, 108
  * Boris (1970) Boris, J. P. 1970, NRL Memorandum Report 2167
  * Buehler et al. (2013) Buehler, D., Lagg, A., & Solanki, S. K. 2013, A&A, 555, A33
  * Carbon & Gingerich (1969) Carbon, D. F. & Gingerich, O. 1969, in Theory and Observation of Normal Stellar Atmospheres, ed. O. Gingerich, 377
  * Castelli & Kurucz (2003) Castelli, F. & Kurucz, R. L. 2003, in Modelling of Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray, Vol. 210, A20
  * Christensen-Dalsgaard et al. (1996) Christensen-Dalsgaard, J., Dappen, W., Ajukov, S. V., et al. 1996, Science, 272, 1286
  * Danilovic et al. (2010) Danilovic, S., Schüssler, M., & Solanki, S. K. 2010, A&A, 513, A1
  * Dedner et al. (2002) Dedner, A., Kemm, F., Kröner, D., et al. 2002, Journal of Computational Physics, 175, 645
  * Edvardsson et al. (1993) Edvardsson, B., Andersen, J., Gustafsson, B., et al. 1993, A&A, 500, 391
  * Gombosi et al. (2002) Gombosi, T. I., Tóth, G., De Zeeuw, D. L., et al. 2002, Journal of Computational Physics, 177, 176
  * Gustafsson et al. (1975) Gustafsson, B., Bell, R. A., Eriksson, K., & Nordlund, A. 1975, A&A, 500, 67
  * Gustafsson et al. (2008) Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951
  * Haugen et al. (2004) Haugen, N. E., Brandenburg, A., & Dobler, W. 2004, Phys. Rev. E, 70, 016308
  * Hauschildt et al. (1999) Hauschildt, P. H., Allard, F., & Baron, E. 1999, ApJ, 512, 377
  * Henyey et al. (1965) Henyey, L., Vardya, M. S., & Bodenheimer, P. 1965, ApJ, 142, 841
  * Hotta et al. (2015) Hotta, H., Rempel, M., & Yokoyama, T. 2015, ApJ, 803, 42
  * Irwin (2012) Irwin, A. W. 2012, FreeEOS: Equation of State for stellar interiors calculations
  * Jørgensen & Weiss (2019) Jørgensen, A. C. S. & Weiss, A. 2019, MNRAS, 488, 3463
  * Kurucz (1979) Kurucz, R. L. 1979, ApJS, 40, 1
  * Lites et al. (2014) Lites, B. W., Centeno, R., & McIntosh, S. W. 2014, PASJ, 66, S4
  * Ludwig et al. (2009) Ludwig, H. G., Caffau, E., Steffen, M., et al. 2009, Mem. SA.It., 80, 711
  * Ludwig & Kučinskas (2012) Ludwig, H. G. & Kučinskas, A. 2012, A&A, 547, A118
  * Magic et al. (2013) Magic, Z., Collet, R., Asplund, M., et al. 2013, A&A, 557, A26
  * Nordlund (1982) Nordlund, A. 1982, A&A, 107, 1
  * Nordlund & Dravins (1990) Nordlund, A. & Dravins, D. 1990, A&A, 228, 155
  * Nordlund & Stein (1990) Nordlund, Å. & Stein, R. F. 1990, Computer Physics Communications, 59, 119
  * Nordlund et al. (2009) Nordlund, Å., Stein, R. F., & Asplund, M. 2009, Living Reviews in Solar Physics, 6, 2
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10
  * Pietarila Graham et al. (2010) Pietarila Graham, J., Cameron, R., & Schüssler, M. 2010, ApJ, 714, 1606
  * Rempel (2014) Rempel, M. 2014, ApJ, 789, 132
  * Rempel (2017) Rempel, M. 2017, ApJ, 834, 10
  * Rempel et al. (2009) Rempel, M., Schüssler, M., Cameron, R. H., & Knölker, M. 2009, Science, 325, 171
  * Schekochihin et al. (2004) Schekochihin, A. A., Cowley, S. C., Taylor, S. F., Maron, J. L., & McWilliams, J. C. 2004, ApJ, 612, 276
  * Shchukina & Trujillo Bueno (2011) Shchukina, N. & Trujillo Bueno, J. 2011, ApJ, 731, L21
  * Spruit (1979) Spruit, H. C. 1979, Sol. Phys., 61, 363
  * Spruit (1982) Spruit, H. C. 1982, A&A, 108, 348
  * Stein & Nordlund (1998) Stein, R. F. & Nordlund, Å. 1998, ApJ, 499, 914
  * Stix (2004) Stix, M. 2004, The sun : an introduction
  * Strom & Kurucz (1966) Strom, S. E. & Kurucz, R. 1966, AJ, 71, 181
  * Trujillo Bueno et al. (2004) Trujillo Bueno, J., Shchukina, N., & Asensio Ramos, A. 2004, Nature, 430, 326
  * Vögler et al. (2004) Vögler, A., Bruls, J. H. M. J., & Schüssler, M. 2004, A&A, 421, 741
  * Vögler & Schüssler (2007) Vögler, A. & Schüssler, M. 2007, A&A, 465, L43
  * Vögler et al. (2005) Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335
  * Yadav et al. (2021) Yadav, N., Cameron, R. H., & Solanki, S. K. 2021, A&A, 645, A3

##  Appendix A Mach number derivation

We derive the expression for Mach number used in Eq. (5) from §4.1. Consider
the hydrostatic pressure balance equation, the expression for pressure-scale
height ($H_{p}$), the expression for sound speed ($c_{s}$) and the ideal gas
equation of state:

$\displaystyle p^{\prime}$ $\displaystyle=-g\rho$ (8) $\displaystyle H_{p}$
$\displaystyle=-p/p^{\prime}$ (9) $\displaystyle c_{s}^{2}$
$\displaystyle=\gamma p/\rho$ (10) $\displaystyle p$
$\displaystyle=\rho(R/\mu)T$ (11)

Here, $H_{p}$ is the pressure scale height, $p^{\prime}=dp/dz$, $\mu$ is the
mean molar mass, and $1/\mu=(1+E)(X+Y/4+Z/2)$, where $E$ is ionization
fraction and $X,Y,Z$ are the H, He and metal abundances. In the Böhm-Vitense
MLT (Böhm-Vitense, 1958), convective flux and velocities are expressed as,

$\displaystyle F_{conv}$ $\displaystyle={\alpha c_{p}}\rho
vT(\nabla-\nabla_{a})/2$ (12) $\displaystyle v^{2}$
$\displaystyle={\alpha^{2}\delta}gH_{p}(\nabla-\nabla_{a})/8$ (13)

Here, $\alpha$ is the mixing length parameter (usually taken to be between 1.5
to 2), $\delta=1-(\partial\ln\mu/\partial\ln T)_{p}$ (abundance gradient with
temperature for an ideal gas), $c_{p}=(\partial U/\partial
T)_{p}+p\delta/(\rho T)$ (heat capacity at constant pressure) and
$\nabla=d\ln{T}/d\ln{p}$, with the subscript $a$ referring to ‘adiabatic’. See
chapter 6 of Stix (2004) for details.

Near the surface, as most energy is carried by radiation,
$F_{conv}\approx\sigma T^{4}$. We take the ratio of Eq. 12 and Eq. 13 to
eliminate $\nabla-\nabla_{a}$ and use $\sigma T^{4}$ instead of $F_{conv}$.
Taking Mach number $M=v/c_{s}$ and $c_{s}=\sqrt{\gamma RT/\mu}$ and
eliminating $H_{p}$ and $p$ using Eqs. 8,9,11,

$\displaystyle M$
$\displaystyle=\left(\frac{\alpha\delta\sigma\mu^{1/2}}{4c_{p}\gamma^{3/2}R^{1/2}}\right)^{1/3}T^{5/6}\rho^{-1/3}$
(14)

Using relevant values for the terms in the parentheses (all in cgs), we get
Eq. (5) in terms of $T$ and $\rho$ as above

$\displaystyle\alpha$ $\displaystyle=1.8$ (average literature value)
$\displaystyle\delta$ $\displaystyle=1$ (assume $\mu$ is constant)
$\displaystyle\sigma$ $\displaystyle=\text{5.67e-5
erg/cm${}^{2}$/s/K${}^{4}$}$ (Stefan-Boltzmann constant) $\displaystyle\gamma$
$\displaystyle=5/3$ (monoatomic adiabatic index) $\displaystyle\mu$
$\displaystyle=1.2$ (solar surface abundance) $\displaystyle R$
$\displaystyle=\text{8.314e7 erg/K/mol}$ (universal gas constant)
$\displaystyle c_{p}$ $\displaystyle=\gamma R/(\mu(\gamma-1))$ (isobaric
specific heat capacity)

Note that this value of $\mu$ assumes $E\sim 0$ near $\tau=1$. For an M-star,
this may not necessarily be the case. However, even if $E\to 1$, $\mu$ changes
by a factor of 2 and, accordingly, $M$ cahnges by a factor of $2^{1/2}\approx
1.4$ and the qualitative result still holds.

##  Appendix B Turbulent pressure and pressure scale height derivation

In §4.1, we use an expression for $p_{turb}$ in Eq. (6). This can be derived
as follows: consider the total (Reynolds and Maxwell) stress tensor
$\sigma_{ij}$ for ideal MHD

$\displaystyle\sigma_{ij}$
$\displaystyle=\rho(v_{i}v_{j})+\left(p+\frac{B^{2}}{8\pi}\right)\delta_{ij}-\frac{B_{i}B_{j}}{4\pi}$
(15)

Here $i,j$ represent the $x,y$ and $z$ directions.

The pressure balance along the vertical ($z$) direction involves the term
$\langle\nabla_{z}\cdot\sigma_{iz}\rangle_{z}$ which is equal to $\partial
p_{\mathrm{tot}}/\partial z$. We can now compare these two terms to get an
expression for total pressure ($p_{tot}$). Assuming all the off-diagonal terms
($i\neq j$) are negligible (which is akin to saying there is no cross-
correlation between the vertical and the horizontal velocity and magnetic
field components), we obtain the following:

$\displaystyle\langle p_{\mathrm{tot}}\rangle_{z}$
$\displaystyle=\biggl{<}{\rho
v_{z}^{2}+\left(p+\frac{B_{x}^{2}+B_{y}^{2}+B_{z}^{2}}{8\pi}\right)-\frac{\left(B_{z}^{2}+\cancelto{0}{B_{z}B_{x}}+\cancelto{0}{B_{z}B_{y}}\right)}{4\pi}}\biggr{>}_{z}$
(16) $\displaystyle\langle p_{\mathrm{tot}}\rangle_{z}$
$\displaystyle=\langle\rho v_{z}^{2}\rangle+\langle p\rangle_{z}+\frac{\langle
B_{h}^{2}-B_{z}^{2}\rangle_{z}}{8\pi}$ (17)

The next step is to consider how using $p_{tot}$ instead of $p$ affects the
density scale height $H_{\rho}$. In Eq. (8), we substitute expression for
$p_{tot}$ instead of $p$ and use Eq. (11) to eliminate $p$. With this, we
obtain,

$\displaystyle\frac{d}{dz}\left(\rho(RT/\mu+v^{2})\right)$
$\displaystyle=-\rho g$ (18)

Note that we take $v^{2}$ to be $v_{z}^{2}+(B_{h}^{2}-B_{z}^{2})/(8\pi\rho)$
The above equation can be rearranged to get

$\displaystyle\frac{d\rho}{\rho}$
$\displaystyle=-\frac{gdz}{RT/\mu+v^{2}}-\frac{d(RT/\mu+v^{2})}{RT/\mu+v^{2}}$
(19)

From this, one gets the expression for density scale height $H_{\rho}$ as

$\displaystyle H_{\rho}$ $\displaystyle=RT/(\mu g)+v^{2}/g$ (20)

The second term above includes the contribution from the turbulent pressure.
Since this term is almost always smaller for the SSD case relative to the HD
case (mainly because the convective velocities are lower for SSD), the
corresponding pressure scale height for the SSD cases is also almost always
smaller.

##  Appendix C Diagnostics

For any MHD simulation, the computation and evolution of the magnetic field
must be divergence free. To ensure this, MURaM uses a hyperbolic divB cleaning
algorithm (Dedner et al., 2002). Here we show the horizontally averaged
divergence of magnetic field across the box for all 4 SSD cases. Since the
units for $\nabla\cdot B$ are field/length, a proper comparison requires a
normalization. We do so with $||B||/\Delta z$. As mentioned in §2.2, at best,
the error is $O(10^{-3})$.

Figure 7: div$\cdot$B error for the four SSD cases F, G, K and M. The quantity
plotted is horizontally averaged
$(\nabla\cdot\@vec{B})_{\mathrm{rms}}/(||B||/\Delta z)$.