# Overscreening–free electron–phonon interaction in realistic materials

Fulvio Paleari  CNR-ISM, Division of Ultrafast Processes in Materials
(FLASHit), Area della Ricerca di Roma 1, Via Salaria Km 29.3, I-00016
Monterotondo, Scalo, Italy  Andrea Marini  CNR-ISM, Division of Ultrafast
Processes in Materials (FLASHit), Area della Ricerca di Roma 1, Via Salaria Km
29.3, I-00016 Monterotondo, Scalo, Italy

(August 27, 2024)

###### Abstract

State–of–the–art model Hamiltonians, like Fröhlich’s, or density functional
theory approaches to electron–phonon interaction suffer from a severe
overscreening error. This is due to the incorrect treatment of the screening
of the ionic potential. We derive a correct formulation of the electron–phonon
interaction and demonstrate its validity by numerically implementing the new
scheme in a paradigmatic material: MgB2, a system whose double-gap,
low-$T_{c}$ superconductivity depends on exceptionally high phonon linewidths.
We find that the present, correct treatment enhances the linewidths by $57\%$
with respect to what has been previously reported for the anomalous $E_{2g}$
mode. We also discover that the $A_{2u}$ mode is also anomalous (its strong
coupling being completely quenched by the overscreened expression). Our
results deeply question methods based on state–of–the–art approaches and
impact a wide range of fields such as thermal conductivity, phononic
instabilities and non–equilibrium lattice dynamics.

Introduction. The research on the physics induced by electron-phonon
interaction is one of the most prolific topics in materials science and solid
state physicsW. Schäfer (2002); H. Haug (2008); Giustino (2017); van Leeuwen
(2004). The electron-phonon interaction is at the basis of the description of
such disparate phenomena as BCS superconductivityCampi  _et al._ (2021); Choi
_et al._ (2002a); Profeta  _et al._ (2012); Schrodi  _et al._ (2020); Zhang
_et al._ (2021); Nosarzewski  _et al._ (2021), lattice thermal conductivity
and transportLiao  _et al._ (2015); Wang  _et al._ (2016); Tong  _et al._
(2019); Zhou  _et al._ (2020); Gold-Parker  _et al._ (2018), structural
stability and phase transitionsPisana  _et al._ (2007); Roy  _et al._ (2021);
Nguyen  _et al._ (2020); Caruso  _et al._ (2017), ultrafast lattice
dynamicsTong and Bernardi (2020); Caruso (2021), spectral shape of Raman,
x-ray, neutron and related spectroscopiesNovko  _et al._ (2020); Reichardt and
Wirtz (2020); Shukla  _et al._ (2003); Ruf  _et al._ (2001), and more. State-
of-the-art methods are based either on first principles, with density
functional perturbation theory (DFPT) schemesBaroni  _et al._ (1987, 2001);
Gonze (1995a, b) or model electron–phonon approaches like the Fröhlich
HamiltonianFröhlich (1954); Giustino (2017); W. Schäfer (2002). Yet, both
approaches, and more generally any approach introducing an effective electron-
phonon Hamiltonian, suffer from a systematic overscreening error. This error
is due to the different way electrons and phonons are affected by the
electron-phonon interaction, and it has been noted by several authorsKeating
(1968); Giustino (2017); Calandra  _et al._ (2010). Although some attempts to
overcome it have been madeCaruso  _et al._ (2017), the overscreening of the
phonon quasiparticle corrections has mostly been either ignored or consciously
neglected in the literature, including very recent calculationsCampi  _et al._
(2021).

In practice, the overscreening error follows from the assumption of a starting
electron–phonon Hamiltonian written asMarini and Pavlyukh (2018)

$\hat{H}=\hat{H}_{ph}+\hat{W}_{e-ph}+\hat{H}_{e},$ (1)

where $\hat{H}_{ph}$ is the noninteracting phonon system taken in the harmonic
approximation with eigenvalues $\Omega_{\nu}$,
$\hat{W}_{e-ph}=\sum_{ij\nu}g^{SCR}_{\nu ij}\hat{\rho}_{ij}\hat{Q}_{\nu}$ is
the electron–phonon interaction, and $\hat{H}_{e}$ the electronic Hamiltonian.
Here, $g^{SCR}_{\nu ij}$ is the screened electron–phonon matrix element
between electronic states $i$ and $j$ and phonon mode $\nu$, $\hat{\rho}_{ij}$
is the electron density operator and $\hat{Q}_{\nu}$ the phonon displacement.

Eq.(1) is used in countless fields in physics. Its key property is that the
renormalization of $\Omega_{\nu}$ and the screening of $g_{\nu ij}$ appear
simultaneously due to electronic correlations. The most elementary and known
treatment is based on the Hartree approximation, which corresponds to the
well–known Random–Phase Approximation Mattuck (1992); Giustino (2017).
Regardless of the approximation employed, it is generally not possible to
obtain meaningful phonon frequencies without including screening. The problem
is that physically, screening arises dynamically from the response of the
electron system to the perturbation induced by lattice vibrations, while in
Eq. (1) it is embedded in the approximated Hamiltonian from the start. This
procedure causes any scattering event between phonons and electrons to be
preemptively screened: when the dynamical response of the phonon system is
computed at various orders in perturbation theory, the screening is then
counted multiple times in an uncontrollable way. This implies that there is no
formal way to avoid the overscreening error when Eq.(1) is used.

In this manuscript, we provide an overscreening error–free formulation of the
electron-phonon interaction that avoids the use of Eq.(1). By taking the
paradigmatic case of MgB2 as a testbed, we demonstrate that the overscreening
error can lead to an underestimation of the phonon linewidths as large as
$57\%$ (the well-known $E_{2g}$ mode). At the same time we also discover
instances where the ovescreening error leads to a large reduction of the
linewidths by more than one order of magnitude (the $A_{2u}$ mode). The
present overscreening error–free approach also improves the agreement with
experimental results. These results provide new theoretical insights on
electron-phonon physics, while also being especially relevant to all the
fields where electron-phonon interaction is customarily modelled with the
overscreened Hamiltonian, Eq.(1).

In the following, we will first explain our approach, then we will describe
the results on MgB2 and finally discuss the wider implications of our
findings.

Figure 1: Phonon linewidths in MgB2 computed at the A point in the Brillouin
zone. Only the relevant phonon modes are shown. Orange: overscreened (OE)
linewidths. Teal: linewidths without overscreening (no OE). The full and empty
bars represent the contributions to the linewidths stemming from intraband and
interband processes, respectively.

Screened and bare electron–phonon interaction. The phonon-related quantities
appearing in Eq.(1) depend on the treatment of the electron–electron
interaction, $\hat{V}_{e-e}$, embodied in $\hat{H}_{e}$. On the one hand, if
$\hat{V}_{e-e}$ were neglected, then the phonon quantities would be bare,
i.e., unscreened: we would have $\Omega_{\nu}={\Omega}^{BARE}_{\nu}$ and
${g}^{BARE}_{\nu ij}$, with

${g}^{BARE}_{\nu ij}=\int\mathrm{d}^{3}\mathbf{r}\
\phi_{j}^{*}(\mathbf{r})\partial_{\nu}\bigr{|}_{eq}V_{e-ion}(\mathbf{r})\phi_{i}(\mathbf{r})$
(2)

In Eq.(2) the $\phi_{i}(\mathbf{r})$ are eigenfunctions of $\hat{H}_{e}$ and
$V_{e-ion}$ is the electron-ion interaction. On the other hand, if the effect
of $\hat{V}_{e-e}$ is included by using either a fully many–body approach
Giustino (2017) or any mean-field (mf) approximation, like DFPT, then the
phonon frequencies, $\Omega^{SCR}_{\lambda q}$, and the electron–phonon matrix
element, $g^{SCR}_{\nu ij}$, become screened. The matrix elements
$g^{SCR}_{\nu ij}$ are given by:

$g^{SCR}_{\nu ij}=\int\mathrm{d}^{3}\mathbf{r}\mathbf{r}^{\prime}\
\epsilon^{-1}(\mathbf{r},\mathbf{r}^{\prime})\phi_{j}^{*}(\mathbf{r}^{\prime})\partial_{\nu}\bigr{|}_{eq}V_{e-ion}(\mathbf{r}^{\prime})\phi_{i}(\mathbf{r}^{\prime}),$
(3)

where $\epsilon$ is the corresponding static electronic dielectric function
responsible for the screening. Since unscreened phonon frequencies are
unphysical Mattuck (1992), any approach that aims at being predictive and
accurate must include screening. If, however, the approach is based on Eq.(1),
then the overscreening error cannot be avoided, and thus it is always present
from the start when dealing with electron-phonon interactionsTong  _et al._
(2019); Klemens and Williams (1986); Marsiglio (2020); Migdal (1958);
Eliashberg (1960); Lüders  _et al._ (2005); Marques  _et al._ (2005); Allen
(1987); Reichardt and Wirtz (2019); Peter Y. Yu (2010); Burkel (2002); Mahan
(1990). In the language of many-body perturbation theory, the overscreening
error takes the form of a double counting of the terms originating from the
perturbative expansion of $\epsilon^{-1}$.

Expression for the phonon linewidths. The phonon linewidths are given by twice
the imaginary part of the retarded phonon-electron self-energy evaluated at
the phonon energy itselfMahan (1990)111In this work we use the well-known
Allen approximationAllen (1972); Giustino (2017) since it has proven to work
well for low-temperature MgB2Calandra and Mauri (2005). This allows us both to
ease the computational burden by accelerating numerical convergence, and to
properly compare with the results of Ref. Shukla  _et al._ (2003), where the
Allen formula was also used.. The resulting expressions (with superscript $OE$
denoting the overscreened case) are:

$\gamma_{\lambda q}^{OE}=\frac{4\pi\Omega_{\lambda
q}}{N_{k}}\sum_{nmk}|g^{SCR}_{\lambda
q,mnk}|^{2}\delta(\varepsilon_{nk})\delta(\Omega_{\lambda
q}-\Delta^{q}_{nmk})$ (4)

and

$\gamma_{\lambda q}=\frac{4\pi\Omega_{\lambda q}}{N_{k}}\sum_{nmk}\
\mathcal{G}^{\lambda q}_{mnk}\ \
\delta(\varepsilon_{nk})\delta(\Omega_{\lambda q}-\Delta^{q}_{nmk}).$ (5)

These expressions describe energy-conserving scatterings between phonons and
electron-hole pairs. In particular, $\Omega_{\lambda q}$ is the phonon energy
of branch $\lambda$ and momentum $q$,
$\Delta_{nmk}^{q}=\varepsilon_{nk}-\varepsilon_{mk-q}$ is the energy of an
electronic transition from band $m$ with momentum $k-q$ to band $n$ of
momentum $k$, and $N_{k}$ is the total number of $k$-points included in the
discrete $k$-integral. Eq.(4) is affected by the overscreening error while
Eq.(5) is not. The difference in the two expressions is the intensity of the
coupling strength, which in the overscreened case depends exclusively on the
screened matrix element $g^{SCR}_{\lambda q,nmk}=1/\sqrt{2\Omega_{\lambda
q}}\bra{nk}\partial_{\lambda q}\hat{W}_{e-ion}\ket{mk-q}$ expressed in the
single-particle basis with the correct energy renormalisation. In DFPT,
$W_{e-ion}$ corresponds to the self-consistent Kohn-Sham potential. The
coupling strength which defines an overscreening error–free phonon linewidth
is

$\mathcal{G}^{\lambda q}_{mnk}=\frac{1}{2}\left\\{g^{SCR}_{\lambda
q,mnk}\left({g}^{BARE}_{\lambda q,mnk}\right)^{*}+{g}^{BARE}_{\lambda
q,mnk}\left(g^{SCR}_{\lambda q,mnk}\right)^{*}\right\\}.$ (6)

Eq.(6) is demonstrated in our companion paperPaleari and Marini (2021)222Our
companion paper Paleari and Marini (2021) contains:(i) a full derivation of
the phonon-electron self-energy without overscreening error starting from Eq.
(1), (ii) a derivation of the linewidths expressions including the requirement
of positive definitenessStrinati (1988); Stefanucci and van Leeuwen (2013)
entering Eq. (6), which was ignored in Ref. Caruso  _et al._ (2017), and
finally (iii) technical information about the code implementation. and follows
from the definition of the left and right Dyson equations for the phonon
propagator. This expression is fully real and contains a combination of both
screened and bare – with $g_{\lambda q,mnk}^{BARE}=1/\sqrt{2\Omega_{\lambda
q}}\bra{nk}\partial_{\lambda q}\hat{V}_{e-ion}\ket{mk-q}$ – matrix elements.
We have modified the Quantum Espresso codeGiannozzi  _et al._ (2017) in order
to extract the ${g}^{BARE}$ matrix elements from a DFPT calculation, and we
subsequently implemented a scheme for the calculation of Eqs. (4) and (5) in
the Yambo codeSangalli  _et al._ (2019)333The computational details needed to
reproduce our MgB2 results are the following. More information about the steps
and general structure of the calculations, as well as about the double grid
procedure, can be found in our companion paperPaleari and Marini (2021). The
input files we used for the various steps are readily available as
Supplemental Material in the interest of reproducibility. The DFT steps were
performed with PAWBlöchl (1994), scalar-relativistic pseudopotentials,
available here: https://www.quantum-espresso.org/pseudopotentials. We used a
plane-wave cutoff of $55$ Ry and a gaussian smearing of $0.025$ Ry. We
included 12 electronic bands in the calculations since they are enough to
describe transitions within one Debye energy ($\sim 100$ meV) around the Fermi
level. The electronic ground state and the derivatives of the self-consistent
potential were converged with a regular $k$-point grid equal to $9\times
9\times 9$. In order to ensure consistency between electron and phonon
calculations, we used the internally relaxed lattice parameters which are
equal to $a=3.074$ Å and $c=3.521$ Å for the in-plane and out-of-plane lattice
parameters, respectively. We then used a coarse grid of $24\times 24\times 24$
$k$-points for the electron-phonon calculations in conjunction with a fine
grid of $55296$ random $k$-points for the electronic band energies. Both these
grids were converged together and along with a broadening parameter for the
delta functions in Eqs. (4) and (5) equal to $5$ meV. Lower values generally
tended to introduce numerical noise, even for our densest coarse and fine
grids..

Application to MgB2. Magnesium diboride is a metallic layered material
composed of alternating $2$D sheets of boron and magnesium. It transitions to
phonon-mediated superconductivity at the critical temperature $T_{c}=39$
KNagamatsu  _et al._ (2001); Liu  _et al._ (2001); Bohnen  _et al._ (2001);
Goncharov  _et al._ (2001); Xi (2008). This behaviour is almost entirely due
to electron–phonon coupling relative to the boron atoms, whose electrons form
in-plane $\sigma$ and out-of-plane $\pi$ bonds. These bonds are in turn
responsible for the existence of two superconducting band gaps, with different
theoretical $T_{c}$Choi  _et al._ (2002b, a); Golubov  _et al._ (2002); Chen
_et al._ (2001); Xi (2008); Margine and Giustino (2013); Floris  _et al._
(2005). In particular, the $\sigma$ bands are considered to yield a giant
“anomalous” electron–phonon coupling due to the strong orbital overlap induced
by the in-plane optical phonon mode $E_{2g}$, as opposed to the $\pi$ bands
undergoing a weak coupling. The calculated phonon and electron dispersions are
shown in Fig. 2(a) and (b), respectively.

As they are necessary ingredients for the calculation of critical temperatures
and superconducting gapsAllen (1972); Choi  _et al._ (2002a); Campi  _et al._
(2021), the $E_{2g}$ phonon linewidths have been extensively studied, both
theoretically and experimentally, along the $\Gamma$AShukla  _et al._ (2003);
Calandra and Mauri (2005); Calandra  _et al._ (2007) and $\Gamma$MBaron  _et
al._ (2004) directions in the hexagonal Brillouin zone. In these studies, the
comparison is made between the full-width half-maximum of inelastic x-ray
scattering spectral peaks and the results of Eq. (4), showing in general
reasonable agreement, though particularly along $\Gamma$A the linewidths are
found to be larger than the theoretical results. For example, in Ref. Shukla
_et al._ (2003) a theoretical value of $20.35$ meV is found at point A, while
the experimental peak width is closer to $30$ meV444The contribution to the
theoretical value of the linewidths coming from phonon-phonon interaction
(i.e., anharmonicities) is very smallShukla  _et al._ (2003); Calandra  _et
al._ (2007) and we will not address it. Furthermore, let us mention that the
study of the linewidths of long-range phonons is instead generally conducted
with Raman spectroscopyHlinka  _et al._ (2001); Postorino  _et al._ (2001);
Bohnen  _et al._ (2001); Rafailov  _et al._ (2002); Chen  _et al._ (2001).
Results show a complex spectral structure which may be captured by considering
nonadiabatic corrections to the phonon energies and possibly a phonon self-
energy with dynamical effects beyond first orderCalandra  _et al._ (2007);
Ponosov and Streltsov (2017); Novko (2018). Therefore low-$q$ linewidths will
not be a focus of this study.

Figure 2: Electron and phonon properties of MgB2.(a) Phonon dispersion. The
phonon mode symmetries at the $\Gamma$ point are given. The modes and
$q$-regions exhibiting large coupling to the $\sigma$ and $\pi$ bands are
emphasised in blue and red, respectively, and are analysed in Fig. 3. (b) Band
structure around the Fermi level (at $0$). The $\sigma$ ($\pi$) bands are
represented in blue (red). The black horizontal dashed lines mark a region
with a width of one Debye energy ($\leavevmode\nobreak\ 100$ meV) around the
Fermi level.

Let us start the discussion from the calculation of the $E_{2g}$ mode. As can
be seen from Fig. 1, we obtain $18.2$ meV in the overscreened case, Eq. (4),
in very good agreement with the same calculation in Ref. Shukla  _et al._
(2003). The correctly screened case, Eq. (5), gives instead the value of
$28.6$ meV, showing a $57\%$ increment in the phonon linewidths. We notice
that in the case of $\sigma$ bands, around $70\%$ of the contribution comes
from “intraband” terms (i.e., elecron-hole pairs are formed within the same
$\sigma$ subband), while the remaining $30\%$ is due to “interband” terms
involving different $\sigma$ subbands. The large $57\%$ increase in the
linewidths is also the average along the full $\Gamma$A direction, as can be
seen from Fig. 3(a), while a strong increase also appears along the $\Gamma$M
direction. Along the latter, both overscreened and overscreening–free
linewidths sensibly decrease after the midpoint from $\Gamma$ to M due to a
sharp increase of the relative phonon energies. The comparison with experiment
is difficult due to the large errorbars, but overall we do obtain a better
agreement in the overscreening–free case (compare with Fig. 3 in Shukla  _et
al._ (2003) and Fig. 3 in Baron  _et al._ (2004) for the $\Gamma$A and
$\Gamma$M directions, respectively).

The $E_{2g}$ mode is not the only one undergoing large changes when
overscreening is removed. In fact, we see from Fig. 1 that also the acoustic
$A_{2u}$ mode gains a giant linewidth increase from $1$ meV (overscreened
case) to $41$ meV (overscreening–free case). By looking at Fig. 3(b) we
realise that these giant linewidths appear along the full $\Gamma$A direction,
where the acoustic $A_{2u}$ mode maintains an average linewidth of $42$ meV,
and also characterise the higher-energy, optical mode of the same symmetry.
Now, the infrared-active $A_{2u}$ modes involve out-of-plane oscillations of
the boron atoms; furthermore, we see that the linewidth is composed by purely
intraband contributions. Both these observations are consistent with the
coupling of the $A_{2u}$ modes to the $\pi$ bands. Therefore, we explain the
enormous difference between the overscreened and overscreening–free cases as
follows. The overlap of the boron $\pi$ orbitals increases because of the
$A_{2u}$ ionic oscillations; however, the effect would be weakened both by the
low charge density in the overlap volume (at equilibrium, two boron atoms on
different layers are separated by $3.5$ Å) and especially by the electronic
screening of the intervening magnesium layer. In this situation, the double
counting of the screening in the calculation of the electron–phonon matrix
elements then leads to a massive reduction of the coupling strength, which is
not found in the correct treatment.

Our results suggest that the $\pi$-band electron–phonon coupling may also be
anomalous, with consequences for the first-principles calculation of the
superconducting properties in MgB2, as well as for the role of interband
$\sigma$-$\pi$ scattering, which is considered important to stabilize the
contributions to $T_{c}$ from both superconducting gaps at $39$ KChoi  _et
al._ (2002b); Geerk  _et al._ (2005); Mazin  _et al._ (2004); Xu  _et al._
(2003); Floris  _et al._ (2005). The $\sigma$-$\pi$ carrier dynamics was
studied in two recent worksBaldini  _et al._ (2017); Novko  _et al._ (2020).

Interestingly, the acoustic $A_{2u}$ linewidths remain constant along
$\Gamma$A despite strong variations in the relative phonon energies, while the
optical $A_{2u}$ linewidths decrease from $28$ to $0.6$ meV despite the phonon
energies being roughly constant. We also point out that such a giant
overscreening effect is not limited to the coupling with the $\pi$ bands along
$\Gamma$A, but also appears – although to a lesser extent – in the linewidths
of the in-plane $E_{1u}$ modes, which couple with the $\sigma$ bands, along
$\Gamma$M (here the largest effect is on the higher-energy acoustic $E_{1u}$
mode up to $0.5\Gamma$M, where the overscreening–free linewidths rise to
$15-20$ meV).

Figure 3: Phonon linewidths of (a) the $E_{2g}$ optical modes along the
M$\Gamma$A path and (b) the $A_{2u}$ acoustic (top) and optical (bottom) modes
along $\Gamma$A. Orange: overscreened (OE) coupling. Teal: screened (no OE)
coupling.

Discussion. We have studied from first-principles the phonon linewidths of
MgB2 after removing the overscreening error in the electron-phonon coupling.
We prove that in this case, this removal entails large qualitative and
quantitative changes with respect to the state-of-the-art overscreened case.
These findings have consequences for a wide range of theoretical calculations
of phonon-related properties in crystals. We examine two such cases below.

The current theoretical approaches around phonon-mediated superconductivity
are based on Migdal-Eliahsberg theory, which is affected by the overscreening
error since the electron-phonon coupling strength enters the theory via the
Eliahsberg spectral function. We believe that strong variations in the
overscreening error–free coupling strength, like the ones presented here, will
impact the theoretical estimates of critical temperatures, quasiparticle
renormalisations and superconducting band gaps, which in some cases show a
discrepancy with the experimental results.Margine and Giustino (2013);
Giustino (2017)

The physics of out-of-equilibrium, ultrafast phenomena related to lattice
dynamics may also be very sensitive to the overscreening error. The relaxation
of the coupled electron-lattice system is governed by differential equations
in which the overscreening error is introduced at every time step, possibly
strongly affecting the resulting rates. In addition, the role of the time
evolution of the electronic screening out of equilibrium has to be
investigated carefully, along with the issues of total energy and particle
conservations, since the overscreening error–free treatment openly questions
the standard semiclassical Boltzmann approach based on Fermi’s golden
rule.Allen (1987)

Acknowledgments. A.M. and F.P. acknowledge the funding received from the
European Union projects: MaX Materials design at the eXascale
H2020-EINFRA-2015-1, Grant agreement n. 676598, and
H2020-INFRAEDI-2018-2020/H2020-INFRAEDI-2018-1, Grant agreement n. 824143;
Nanoscience Foundries and Fine Analysis - Europe H2020-INFRAIA-2014-2015,
Grant agreement n. 654360.

## References

  * W. Schäfer (2002) M. Wegener W. Schäfer,  _Semiconductor Optics and Transport Phenomena_ (Springer, 2002).
  * H. Haug (2008) A.-P. Jauho H. Haug,  _Quantum Kinetics in Transport and Optics of Semiconductors_ (Springer, 2008).
  * Giustino (2017) Feliciano Giustino, “Electron-phonon interactions from first principles,” Rev. Mod. Phys. 89, 015003 (2017).
  * van Leeuwen (2004) Robert van Leeuwen, “First-principles approach to the electron-phonon interaction,” Phys. Rev. B 69, 115110 (2004).
  * Campi  _et al._ (2021) Davide Campi, Simran Kumari, and Nicola Marzari, “Prediction of phonon-mediated superconductivity with high critical temperature in the two-dimensional topological semimetal w2n3,” Nano Letters 21, 3435–3442 (2021).
  * Choi  _et al._ (2002a) Hyoung Joon Choi, David Roundy, Hong Sun, Marvin L. Cohen, and Steven G. Louie, “The origin of the anomalous superconducting properties of MgB2,” Nature 418, 758–760 (2002a).
  * Profeta  _et al._ (2012) G. Profeta, M. Calandra, and F. Mauri, “Phonon-mediated superconductivity in graphene by lithium deposition,” Nature Physics 8 (2012), 10.1038/nphys2181.
  * Schrodi  _et al._ (2020) Fabian Schrodi, Peter M. Oppeneer, and Alex Aperis, “Full-bandwidth eliashberg theory of superconductivity beyond migdal’s approximation,” Phys. Rev. B 102, 024503 (2020).
  * Zhang  _et al._ (2021) X. Zhang, K.-H. Jin, J. Mao, M. Zhao, J. Liu, and F. Liu, “Prediction of intrinsic topological superconductivity in mn-doped gete monolayer from first-principles,” npj Computational Materials 7 (2021), 10.1038/s41524-021-00511-x.
  * Nosarzewski  _et al._ (2021) Benjamin Nosarzewski, Michael Schüler, and Thomas P. Devereaux, “Spectral properties and enhanced superconductivity in renormalized migdal-eliashberg theory,” Phys. Rev. B 103, 024520 (2021).
  * Liao  _et al._ (2015) Bolin Liao, Bo Qiu, Jiawei Zhou, Samuel Huberman, Keivan Esfarjani, and Gang Chen, “Significant reduction of lattice thermal conductivity by the electron-phonon interaction in silicon with high carrier concentrations: A first-principles study,” Phys. Rev. Lett. 114, 115901 (2015).
  * Wang  _et al._ (2016) Yan Wang, Zexi Lu, and Xiulin Ruan, “First principles calculation of lattice thermal conductivity of metals considering phonon-phonon and phonon-electron scattering,” Journal of Applied Physics 119, 225109 (2016), https://doi.org/10.1063/1.4953366 .
  * Tong  _et al._ (2019) Zhen Tong, Shouhang Li, Xiulin Ruan, and Hua Bao, “Comprehensive first-principles analysis of phonon thermal conductivity and electron-phonon coupling in different metals,” Phys. Rev. B 100, 144306 (2019).
  * Zhou  _et al._ (2020) Jiawei Zhou, Hyun D. Shin, Ke Chen, Bai Song, Ryan A. Duncan, Qian Xu, Alexei A. Maznev, Keith A. Nelson, and Gang Chen, “Direct observation of large electron–phonon interaction effect on phonon heat transport,” Nature Communications 11, 6040 (2020).
  * Gold-Parker  _et al._ (2018) Aryeh Gold-Parker, Peter M. Gehring, Jonathan M. Skelton, Ian C. Smith, Dan Parshall, Jarvist M. Frost, Hemamala I. Karunadasa, Aron Walsh, and Michael F. Toney, “Acoustic phonon lifetimes limit thermal transport in methylammonium lead iodide,” Proceedings of the National Academy of Sciences 115, 11905–11910 (2018), https://www.pnas.org/content/115/47/11905.full.pdf .
  * Pisana  _et al._ (2007) Simone Pisana, Michele Lazzeri, Cinzia Casiraghi, Kostya S. Novoselov, A. K. Geim, Andrea C. Ferrari, and Francesco Mauri, “Breakdown of the adiabatic born–oppenheimer approximation in graphene,” Nature Materials 6, 198–201 (2007).
  * Roy  _et al._ (2021) Aditya Prasad Roy, Naini Bajaj, Ranjan Mittal, Peram D. Babu, and Dipanshu Bansal, “Quasi-one-dimensional fermi surface nesting and hidden nesting enable multiple kohn anomalies in $\alpha$-uranium,” Phys. Rev. Lett. 126, 096401 (2021).
  * Nguyen  _et al._ (2020) Thanh Nguyen, Fei Han, Nina Andrejevic, Ricardo Pablo-Pedro, Anuj Apte, Yoichiro Tsurimaki, Zhiwei Ding, Kunyan Zhang, Ahmet Alatas, Ercan E. Alp, Songxue Chi, Jaime Fernandez-Baca, Masaaki Matsuda, David Alan Tennant, Yang Zhao, Zhijun Xu, Jeffrey W. Lynn, Shengxi Huang, and Mingda Li, “Topological singularity induced chiral kohn anomaly in a weyl semimetal,” Phys. Rev. Lett. 124, 236401 (2020).
  * Caruso  _et al._ (2017) Fabio Caruso, Moritz Hoesch, Philipp Achatz, Jorge Serrano, Michael Krisch, Etienne Bustarret, and Feliciano Giustino, “Nonadiabatic kohn anomaly in heavily boron-doped diamond,” Phys. Rev. Lett. 119, 017001 (2017).
  * Tong and Bernardi (2020) Xiao Tong and Marco Bernardi, “Toward precise simulations of the coupled ultrafast dynamics of electrons and atomic vibrations in materials,” (2020), arXiv:2009.07958 [cond-mat.mtrl-sci] .
  * Caruso (2021) Fabio Caruso, “Nonequilibrium lattice dynamics in monolayer mos2,” The Journal of Physical Chemistry Letters 12, 1734–1740 (2021).
  * Novko  _et al._ (2020) Dino Novko, Fabio Caruso, Claudia Draxl, and Emmanuele Cappelluti, “Ultrafast hot phonon dynamics in ${\mathrm{mgb}}_{2}$ driven by anisotropic electron-phonon coupling,” Phys. Rev. Lett. 124, 077001 (2020).
  * Reichardt and Wirtz (2020) Sven Reichardt and Ludger Wirtz, “Nonadiabatic exciton-phonon coupling in raman spectroscopy of layered materials,” Science Advances 6 (2020), 10.1126/sciadv.abb5915.
  * Shukla  _et al._ (2003) Abhay Shukla, Matteo Calandra, Matteo d’Astuto, Michele Lazzeri, Francesco Mauri, Christophe Bellin, Michael Krisch, J. Karpinski, S. M. Kazakov, J. Jun, D. Daghero, and K. Parlinski, “Phonon dispersion and lifetimes inMgB2,” Physical Review Letters 90 (2003), 10.1103/physrevlett.90.095506.
  * Ruf  _et al._ (2001) T. Ruf, J. Serrano, M. Cardona, P. Pavone, M. Pabst, M. Krisch, M. D’Astuto, T. Suski, I. Grzegory, and M. Leszczynski, “Phonon dispersion curves in wurtzite-structure gan determined by inelastic x-ray scattering,” Phys. Rev. Lett. 86, 906–909 (2001).
  * Baroni  _et al._ (1987) Stefano Baroni, Paolo Giannozzi, and Andrea Testa, “Green’s-function approach to linear response in solids,” Phys. Rev. Lett. 58, 1861–1864 (1987).
  * Baroni  _et al._ (2001) Stefano Baroni, Stefano de Gironcoli, Andrea Dal Corso, and Paolo Giannozzi, “Phonons and related crystal properties from density-functional perturbation theory,” Rev. Mod. Phys. 73, 515–562 (2001).
  * Gonze (1995a) Xavier Gonze, “Perturbation expansion of variational principles at arbitrary order,” Phys. Rev. A 52, 1086–1095 (1995a).
  * Gonze (1995b) Xavier Gonze, “Adiabatic density-functional perturbation theory,” Phys. Rev. A 52, 1096–1114 (1995b).
  * Fröhlich (1954) H. Fröhlich, “Electrons in lattice fields,” Advances in Physics 3, 325–361 (1954).
  * Keating (1968) P. N. Keating, “Dielectric screening and the phonon spectra of metallic and nonmetallic crystals,” Phys. Rev. 175, 1171–1180 (1968).
  * Calandra  _et al._ (2010) Matteo Calandra, Gianni Profeta, and Francesco Mauri, “Adiabatic and nonadiabatic phonon dispersion in a wannier function approach,” Phys. Rev. B 82, 165111 (2010).
  * Marini and Pavlyukh (2018) Andrea Marini and Yaroslav Pavlyukh, “Functional approach to the electronic and bosonic dynamics of many-body systems perturbed with an arbitrary strong electron-boson interaction,” Phys. Rev. B 98, 075105 (2018).
  * Mattuck (1992) R. D. Mattuck,  _A guide to Feynman diagrams in the many-body problem_ (McGraw-Hill, 1992).
  * Klemens and Williams (1986) P. G. Klemens and R. K. Williams, “Thermal conductivity of metals and alloys,” International Metals Reviews 31, 197–215 (1986), https://doi.org/10.1179/imtr.1986.31.1.197 .
  * Marsiglio (2020) F. Marsiglio, “Eliashberg theory: A short review,” Annals of Physics 417, 168102 (2020), eliashberg theory at 60: Strong-coupling superconductivity and beyond.
  * Migdal (1958) A B Migdal, “Interaction between electrons and lattice vibrations in a normal metal,” Soviet Physics JETP 7, 996 (1958).
  * Eliashberg (1960) G M Eliashberg, “Interactions between electrons and lattice vibrations in a superconductor,” Soviet Physics JETP 11, 696 (1960).
  * Lüders  _et al._ (2005) M. Lüders, M. A. L. Marques, N. N. Lathiotakis, A. Floris, G. Profeta, L. Fast, A. Continenza, S. Massidda, and E. K. U. Gross, “Ab initio theory of superconductivity. i. density functional formalism and approximate functionals,” Phys. Rev. B 72, 024545 (2005).
  * Marques  _et al._ (2005) M. A. L. Marques, M. Lüders, N. N. Lathiotakis, G. Profeta, A. Floris, L. Fast, A. Continenza, E. K. U. Gross, and S. Massidda, “Ab initio theory of superconductivity. ii. application to elemental metals,” Phys. Rev. B 72, 024546 (2005).
  * Allen (1987) Philip B. Allen, “Theory of thermal relaxation of electrons in metals,” Phys. Rev. Lett. 59, 1460–1463 (1987).
  * Reichardt and Wirtz (2019) Sven Reichardt and Ludger Wirtz, “Theory of resonant raman scattering: Towards a comprehensive ab initio description,” Phys. Rev. B 99, 174312 (2019).
  * Peter Y. Yu (2010) Manuel Cardona Peter Y. Yu,  _Fundamentals of Semiconductors Physics and Materials Properties_ (Springer, 2010).
  * Burkel (2002) E. Burkel, “Phonon spectroscopy by inelastic x-ray scattering,” Rep. Prog. Phys. 63 (2002), 10.1088/0034-4885/63/2/203.
  * Mahan (1990) G. D. Mahan,  _Many-Particle Physics_ (Plenum Press, 1990).
  * Note (1) In this work we use the well-known Allen approximationAllen (1972); Giustino (2017) since it has proven to work well for low-temperature MgB2Calandra and Mauri (2005). This allows us both to ease the computational burden by accelerating numerical convergence, and to properly compare with the results of Ref. Shukla  _et al._ (2003), where the Allen formula was also used.
  * Paleari and Marini (2021) Fulvio Paleari and Andrea Marini, “Electron–phonon interaction without overscreening: a strategy for first–principles modelling,” PRB joint submission XX, XXXXXX (2021).
  * Note (2) Our companion paper Paleari and Marini (2021) contains:(i) a full derivation of the phonon-electron self-energy without overscreening error starting from Eq. (1), (ii) a derivation of the linewidths expressions including the requirement of positive definitenessStrinati (1988); Stefanucci and van Leeuwen (2013) entering Eq. (6), which was ignored in Ref. Caruso  _et al._ (2017), and finally (iii) technical information about the code implementation.
  * Giannozzi  _et al._ (2017) P Giannozzi, O Andreussi, T Brumme, O Bunau, M Buongiorno Nardelli, M Calandra, R Car, C Cavazzoni, D Ceresoli, M Cococcioni, N Colonna, I Carnimeo, A Dal Corso, S de Gironcoli, P Delugas, R A DiStasio Jr, A Ferretti, A Floris, G Fratesi, G Fugallo, R Gebauer, U Gerstmann, F Giustino, T Gorni, J Jia, M Kawamura, H-Y Ko, A Kokalj, E Küçükbenli, M Lazzeri, M Marsili, N Marzari, F Mauri, N L Nguyen, H-V Nguyen, A Otero de-la Roza, L Paulatto, S Poncé, D Rocca, R Sabatini, B Santra, M Schlipf, A P Seitsonen, A Smogunov, I Timrov, T Thonhauser, P Umari, N Vast, X Wu, and S Baroni, “Advanced capabilities for materials modelling with quantum espresso,” Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Sangalli  _et al._ (2019) D Sangalli, A Ferretti, H Miranda, C Attaccalite, I Marri, E Cannuccia, P Melo, M Marsili, F Paleari, A Marrazzo, G Prandini, P Bonfà, M O Atambo, F Affinito, M Palummo, A Molina-Sánchez, C Hogan, M Grüning, D Varsano, and A Marini, “Many-body perturbation theory calculations using the yambo code,” J. Phys. Condens. Matter 31, 325902 (2019).
  * Note (3) The computational details needed to reproduce our MgB2 results are the following. More information about the steps and general structure of the calculations, as well as about the double grid procedure, can be found in our companion paperPaleari and Marini (2021). The input files we used for the various steps are readily available as Supplemental Material in the interest of reproducibility. The DFT steps were performed with PAWBlöchl (1994), scalar-relativistic pseudopotentials, available here: https://www.quantum-espresso.org/pseudopotentials. We used a plane-wave cutoff of $55$ Ry and a gaussian smearing of $0.025$ Ry. We included 12 electronic bands in the calculations since they are enough to describe transitions within one Debye energy ($\sim 100$ meV) around the Fermi level. The electronic ground state and the derivatives of the self-consistent potential were converged with a regular $k$-point grid equal to $9\times 9\times 9$. In order to ensure consistency between electron and phonon calculations, we used the internally relaxed lattice parameters which are equal to $a=3.074$ Å and $c=3.521$ Å for the in-plane and out-of-plane lattice parameters, respectively. We then used a coarse grid of $24\times 24\times 24$ $k$-points for the electron-phonon calculations in conjunction with a fine grid of $55296$ random $k$-points for the electronic band energies. Both these grids were converged together and along with a broadening parameter for the delta functions in Eqs. (4) and (5) equal to $5$ meV. Lower values generally tended to introduce numerical noise, even for our densest coarse and fine grids.
  * Nagamatsu  _et al._ (2001) Jun Nagamatsu, Norimasa Nakagawa, Takahiro Muranaka, Yuji Zenitani, and Jun Akimitsu, “Superconductivity at 39 k in magnesium diboride,” Nature 410, 63–64 (2001).
  * Liu  _et al._ (2001) Amy Y. Liu, I. I. Mazin, and Jens Kortus, “Beyond eliashberg superconductivity inMgB2: Anharmonicity, two-phonon scattering, and multiple gaps,” Physical Review Letters 87 (2001), 10.1103/physrevlett.87.087005.
  * Bohnen  _et al._ (2001) K.-P. Bohnen, R. Heid, and B. Renker, “Phonon dispersion and electron-phonon coupling inMgB2andalb2,” Physical Review Letters 86, 5771–5774 (2001).
  * Goncharov  _et al._ (2001) Alexander F. Goncharov, Viktor V. Struzhkin, Eugene Gregoryanz, Jingzhu Hu, Russell J. Hemley, Ho kwang Mao, G. Lapertot, S. L. Bud’ko, and P. C. Canfield, “Raman spectrum and lattice parameters ofMgB2as a function of pressure,” Physical Review B 64 (2001), 10.1103/physrevb.64.100509.
  * Xi (2008) X X Xi, “Two-band superconductor magnesium diboride,” Reports on Progress in Physics 71, 116501 (2008).
  * Choi  _et al._ (2002b) Hyoung Joon Choi, David Roundy, Hong Sun, Marvin L. Cohen, and Steven G. Louie, “First-principles calculation of the superconducting transition inMgB2within the anisotropic eliashberg formalism,” Physical Review B 66 (2002b), 10.1103/physrevb.66.020513.
  * Golubov  _et al._ (2002) A A Golubov, J Kortus, O V Dolgov, O Jepsen, Y Kong, O K Andersen, B J Gibson, K Ahn, and R K Kremer, “Specific heat of MgB2in a one- and a two-band model from first-principles calculations,” Journal of Physics: Condensed Matter 14, 1353–1360 (2002).
  * Chen  _et al._ (2001) X. K. Chen, M. J. Konstantinović, J. C. Irwin, D. D. Lawrie, and J. P. Franck, “Evidence for two superconducting gaps in ${\mathrm{mgb}}_{2}$,” Phys. Rev. Lett. 87, 157002 (2001).
  * Margine and Giustino (2013) E. R. Margine and F. Giustino, “Anisotropic migdal-eliashberg theory using wannier functions,” Phys. Rev. B 87, 024505 (2013).
  * Floris  _et al._ (2005) A. Floris, G. Profeta, N. N. Lathiotakis, M. Lüders, M. A. L. Marques, C. Franchini, E. K. U. Gross, A. Continenza, and S. Massidda, “Superconducting properties of ${\mathrm{m}\mathrm{g}\mathrm{b}}_{2}$ from first principles,” Phys. Rev. Lett. 94, 037004 (2005).
  * Allen (1972) Philip B. Allen, “Neutron spectroscopy of superconductors,” Phys. Rev. B 6, 2577–2579 (1972).
  * Calandra and Mauri (2005) Matteo Calandra and Francesco Mauri, “Electron-phonon coupling and phonon self-energy in ${\mathrm{mgb}}_{2}$: Interpretation of ${\mathrm{mgb}}_{2}$ raman spectra,” Phys. Rev. B 71, 064501 (2005).
  * Calandra  _et al._ (2007) Matteo Calandra, Michele Lazzeri, and Francesco Mauri, “Anharmonic and non-adiabatic effects in mgb2: Implications for the isotope effect and interpretation of raman spectra,” Physica C: Superconductivity 456, 38–44 (2007), recent Advances in MgB2 Research.
  * Baron  _et al._ (2004) A. Q. R. Baron, H. Uchiyama, Y. Tanaka, S. Tsutsui, D. Ishikawa, S. Lee, R. Heid, K.-P. Bohnen, S. Tajima, and T. Ishikawa, “Kohn anomaly in ${\mathrm{mgb}}_{2}$ by inelastic x-ray scattering,” Phys. Rev. Lett. 92, 197004 (2004).
  * Note (4) The contribution to the theoretical value of the linewidths coming from phonon-phonon interaction (i.e., anharmonicities) is very smallShukla  _et al._ (2003); Calandra  _et al._ (2007) and we will not address it. Furthermore, let us mention that the study of the linewidths of long-range phonons is instead generally conducted with Raman spectroscopyHlinka  _et al._ (2001); Postorino  _et al._ (2001); Bohnen  _et al._ (2001); Rafailov  _et al._ (2002); Chen  _et al._ (2001). Results show a complex spectral structure which may be captured by considering nonadiabatic corrections to the phonon energies and possibly a phonon self-energy with dynamical effects beyond first orderCalandra  _et al._ (2007); Ponosov and Streltsov (2017); Novko (2018). Therefore low-$q$ linewidths will not be a focus of this study.
  * Geerk  _et al._ (2005) J. Geerk, R. Schneider, G. Linker, A. G. Zaitsev, R. Heid, K.-P. Bohnen, and H. v. Löhneysen, “Observation of interband pairing interaction in a two-band superconductor: ${\mathrm{mgb}}_{2}$,” Phys. Rev. Lett. 94, 227005 (2005).
  * Mazin  _et al._ (2004) I. I. Mazin, O. K. Andersen, O. Jepsen, A. A. Golubov, O. V. Dolgov, and J. Kortus, “Comment on “first-principles calculation of the superconducting transition in ${\mathrm{mgb}}_{2}$ within the anisotropic eliashberg formalism”,” Phys. Rev. B 69, 056501 (2004).
  * Xu  _et al._ (2003) Y. Xu, M. Khafizov, L. Satrapinsky, P. Kúš, A. Plecenik, and Roman Sobolewski, “Time-resolved photoexcitation of the superconducting two-gap state in ${\mathrm{m}\mathrm{g}\mathrm{b}}_{2}$ thin films,” Phys. Rev. Lett. 91, 197004 (2003).
  * Baldini  _et al._ (2017) E. Baldini, A. Mann, L. Benfatto, E. Cappelluti, A. Acocella, V. M. Silkin, S. V. Eremeev, A. B. Kuzmenko, S. Borroni, T. Tan, X. X. Xi, F. Zerbetto, R. Merlin, and F. Carbone, “Real-time observation of phonon-mediated $\sigma\text{$-$}\pi$ interband scattering in ${\mathrm{mgb}}_{2}$,” Phys. Rev. Lett. 119, 097002 (2017).
  * Strinati (1988) Giancarlo Strinati, “Application of the green’s functions method to the study of the optical properties of semiconductors,” La Rivista del Nuovo Cimento 11, 1–86 (1988).
  * Stefanucci and van Leeuwen (2013) G. Stefanucci and R. van Leeuwen,  _Nonequilibrium Many-Body Theory of Quantum Systems_ (Cambridge University Press, 2013).
  * Blöchl (1994) P. E. Blöchl, “Projector augmented-wave method,” Phys. Rev. B 50, 17953–17979 (1994).
  * Hlinka  _et al._ (2001) J. Hlinka, I. Gregora, J. Pokorný, A. Plecenik, P. Kúš, L. Satrapinsky, and Š. Beňačka, “Phonons in MgB2 by polarized raman scattering on single crystals,” Physical Review B 64 (2001), 10.1103/physrevb.64.140503.
  * Postorino  _et al._ (2001) P. Postorino, A. Congeduti, P. Dore, A. Nucara, A. Bianconi, D. Di Castro, S. De Negri, and A. Saccone, “Effect of the al content on the optical phonon spectrum inMg1-xAlxB2,” Physical Review B 65 (2001), 10.1103/physrevb.65.020507.
  * Rafailov  _et al._ (2002) P.M Rafailov, M Dworzak, and C Thomsen, “Luminescence and raman spectroscopy on mgb2,” Solid State Communications 122, 455–458 (2002).
  * Ponosov and Streltsov (2017) Yu. S. Ponosov and S. V. Streltsov, “Raman-active ${E}_{2g}$ phonon in ${\mathrm{mgb}}_{2}$: Electron-phonon interaction and anharmonicity,” Phys. Rev. B 96, 214503 (2017).
  * Novko (2018) Dino Novko, “Nonadiabatic coupling effects in ${\mathrm{mgb}}_{2}$ reexamined,” Phys. Rev. B 98, 041112 (2018).


11institutetext: Georg-August-Universität Göttingen, Institut für Astrophysik
und Geophysik, Friedrich-Hund-Platz 1, D-37077 Göttingen, Germany email:
<EMAIL_ADDRESS>22institutetext: Nordita, KTH Royal Institute of
Technology and Stockholm University, Stockholm, Sweden

# Transition from anti-solar to solar-like differential rotation: Dependence
on Prandtl number

P. J. Käpylä 1122

###### Abstract

Context. Late-type stars such as the Sun rotate differentially due to the
interaction of turbulent convection and rotation.

Aims. The aim of the study is to investigate the effects of the effective
thermal Prandtl number, which is the ratio of kinematic viscosity to thermal
diffusivity, on the transition from anti-solar (slow equator, fast poles) to
solar-like (fast equator, slow poles) differential rotation.

Methods. Three-dimensional hydrodynamic and magnetohydrodynamic simulations in
semi-global spherical wedge geometry are used to model convection zones of
solar-like stars.

Results. The overall convective velocity amplitude increases as the Prandtl
number decreases in accordance with earlier studies. The transition from anti-
solar to solar-like differential rotation is insensitive to the Prandtl number
for Prandtl numbers below unity but for Prandtl numbers greater than unity,
solar-like differential rotation becomes significantly harder to excite.
Magnetic fields and more turbulent regimes with higher fluid and magnetic
Reynolds numbers help in achieving solar-like differential rotation in near-
transition cases where anti-solar rotation is found in more laminar
simulations. Solar-like differential rotation occurs only in cases with
radially outward angular momentum transport at the equator. The dominant
contribution to such outward transport near the equator is due to prograde
propagating thermal Rossby waves.

Conclusions. The differential rotation is sensitive to the Prandtl number only
for large Prandtl numbers in the parameter regime explored in the current
study. Magnetic fields have a greater effect on the differential rotation,
although the inferred presence of a small-scale dynamo does not lead to
drastically different results in the present study. The dominance of the
thermal Rossby waves in the simulations is puzzling given the non-detection in
the Sun. The current simulations are shown to be incompatible with the
currently prevailing mean-field theory of differential rotation.

###### Key Words.:

turbulence – convection

##  1 Introduction

The interplay of turbulent convection with the overall rotation of the Sun is
the primary cause of differential rotation observed at the solar surface and
in the interior (e.g. Rüdiger 1989; Miesch & Toomre 2009). Three-dimensional
numerical simulations solving the equations of magnetohydrodynamics (MHD)
capture the essence of this process and routinely produce solutions that are
qualitatively similar to the Sun with equatorial acceleration (e.g. Gilman
1983; Brun et al. 2004; Guerrero et al. 2013; Käpylä et al. 2014). However, it
has become increasingly clear recently that even the most sophisticated
current simulations are missing something essential. The most striking
manifestation of this is that simulations using nominal solar luminosity and
rotation rate often produce anti-solar (AS) differential rotation with
equatorial deceleration (e.g. Fan & Fang 2014; Käpylä et al. 2014; Hotta et
al. 2015), whereas solar-like (SL) differential rotation is achieved only with
significantly more rapid rotation (e.g. Viviani et al. 2018; Matilsky et al.
2020).

This is related to the convective conundrum (O’Mara et al. 2016) which is
essentially the tension between large-scale velocity amplitudes in simulations
in comparison to the Sun (e.g. Hanasoge et al. 2012, 2016; Schumacher &
Sreenivasan 2020). Until recently, the most common way to ensure SL
differential rotation in simulations with solar luminosity and rotation rate
has been to lower the convective velocities by artificially enhancing the
radiative diffusivity (e.g. Fan & Fang 2014; Käpylä et al. 2014; Hotta et al.
2016). This, however, cannot be justified based on physical grounds since
convection is thought to carry practically all of the energy flux through the
solar convection zone (CZ) with the exception of very deep layers. Another,
more plausible, effect is due to magnetic fields: is is conceivable that
sufficiently strong fields can suppress convection to a degree where the
differential rotation flips from AS to SL. Early results with relatively low
resolution simulations were mixed: Karak et al. (2015) found essentially no
dependence on magnetic field while Fan & Fang (2014) and Simitev et al. (2015)
reported more positive outcomes. Nevertheless, these simulations most probably
did not have high enough magnetic Reynolds numbers to excite a small-scale
dynamo. This was addressed by recent high-resolution simulations of Hotta &
Kusano (2021) and Hotta et al. (2022) which suggest that SL differential
rotation can indeed be achieved with the help of an efficient small-scale
dynamo.

Another important parameter is the Prandtl number, ${\rm Pr}=\nu/\chi$, where
$\nu$ is the kinematic viscosity and $\chi$ is the thermal diffusivity. A
notion that the solar convection zone is operating in a regime where the
effective Prandtl number is large, has gained popularity recently (e.g. O’Mara
et al. 2016; Bekki et al. 2017; Karak et al. 2018). While these studies
indicate that the overall velocity amplitudes are decreased in such set-ups,
the problem with the differential rotation becomes actually worse (Karak et
al. 2018). This is because it is not only the velocity amplitude that is
sensitive to ${\rm Pr}$, but also turbulent transport of angular momentum and
heat are affected (e.g. Cattaneo et al. 1991; Käpylä 2021). Furthermore,
theoretical arguments suggest that ${\rm Pr}\ll 1$ in the solar CZ (e.g.
Ossendrijver 2003; Schumacher & Sreenivasan 2020).

Prandtl numbers deviating strongly from unity are challenging numerically and
therefore most simulations are done in the ${\rm Pr}\approx 1$ regime. It is
commonly acknowledged that reaching realistic parameter regimes in terms of,
for example, Prandtl, Reynolds and Rayleigh numbers with current or
foreseeable simulations of stellar convection is infeasible (e.g. Kupka &
Muthsam 2017). The main aim of the present study is to vary the Prandtl number
within the range that is reasonably realizable with numerical simulations with
values above and below unity. The current study is also inspired by recent
results from hydrodynamic non-rotating convection in Cartesian geometry
(Käpylä 2021), where the convective energy transport and velocity statistics
were found to be sensitive to the effective Prandtl number.

##  2 The model

The simulation set-up is similar to those used in Käpylä et al. (2019) and
Käpylä et al. (2020). The simulation domain is a spherical wedge that spans
$r_{\rm in}<r<R$ in radius, where $r_{\rm in}=0.7R$ and $R$ is the radius of
the star, $\theta_{0}<\theta<\pi-\theta_{0}$ in colatitude, where
$\theta_{0}=\pi/12$, and $0<\phi<\pi/2$ in longitude. Equations of fully
compressible MHD are solved

$\displaystyle\frac{\partial{\bm{A}}}{\partial t}$ $\displaystyle=$
$\displaystyle{\bm{U}}\times{\bm{B}}-\eta\mu_{0}{\bm{J}},$ (1)
$\displaystyle\frac{{\rm D}\ln\rho}{{\rm D}t}$ $\displaystyle=$
$\displaystyle-\bm{\nabla}\bm{\cdot}{\bm{U}},$ (2) $\displaystyle\frac{{\rm
D}{\bm{U}}}{{\rm D}t}$ $\displaystyle=$
$\displaystyle{\bm{g}}-2\bm{\Omega}_{0}\times{\bm{U}}-\frac{1}{\rho}(\bm{\nabla}p-{\bm{J}}\times{\bm{B}}-\bm{\nabla}\bm{\cdot}2\nu\rho\bm{\mathsf{S}}),$
(3) $\displaystyle T\frac{{\rm D}s}{{\rm D}t}$ $\displaystyle=$
$\displaystyle\frac{1}{\rho}\left[\eta\mu_{0}{\bm{J}}^{2}-\bm{\nabla}\bm{\cdot}({\bm{F}}^{\rm
rad}+{\bm{F}}^{\rm SGS})\right]+2\nu\bm{\mathsf{S}}^{2},$ (4)

where ${\bm{A}}$ is the magnetic vector potential, ${\bm{U}}$ is the velocity,
${\bm{B}}=\bm{\nabla}\times{\bm{A}}$ is the magnetic field, $\eta$ is the
magnetic diffusivity, $\mu_{0}$ is the permeability of vacuum,
${\bm{J}}=\bm{\nabla}\times{\bm{B}}/\mu_{0}$ is the current density, ${\rm
D}/{\rm D}t=\partial/\partial t+{\bm{U}}\bm{\cdot}\bm{\nabla}$ is the
advective time derivative, $\rho$ is the density, ${\bm{g}}=-\bm{\nabla}\phi$
is the acceleration due to gravity, where $\phi=-GM/r$ is a fixed spherically
symmetric gravitational potential, with $G$ and $M$ being the universal
gravitational constant and the stellar mass, respectively.
$\bm{\Omega}_{0}=(\cos\theta,-\sin\theta,0)\Omega_{0}$ is the angular velocity
vector, where $\Omega_{0}$ is the rotation rate of the frame of reference, $p$
is the pressure, $\nu$ is the kinematic viscosity, $T$ is the temperature, and
$s$ is the specific entropy with $Ds=c_{\rm V}D\ln p-c_{\rm P}D\ln\rho$, where
$c_{\rm V}$ and $c_{\rm P}$ are the specific heat capacities in constant
volume and pressure, respectively. The gas is assumed to obey the ideal gas
law, $p=\mathcal{R}\rho T$, where $\mathcal{R}=c_{\rm P}-c_{\rm V}$ is the gas
constant. The rate of strain tensor is given by

$\displaystyle\mathsf{S}_{ij}\,=\,{\textstyle{1\over
2}}(U_{i;j}+U_{j;i})-{\textstyle{1\over
3}}\delta_{ij}\bm{\nabla}\bm{\cdot}{\bm{U}},$ (5)

where the semicolons refer to covariant derivatives (Mitra et al. 2009). The
radiative flux is given by

$\displaystyle{\bm{F}}^{\rm rad}\,=\,-\,K\bm{\nabla}T,$ (6)

where $K$ is the heat conductivity. The latter consists of two parts,
$K=K_{1}+K_{2}$, where $K_{1}=K_{1}(r)$ is a fixed function of height and
$K_{2}=K_{2}(\rho,T)$ is density- and temperature-dependent according to
Kramers opacity law (Weiss et al. 2004). The profile of $K_{1}$ is given by

$\displaystyle K_{1}=K_{\rm
top}\left[1+\tanh\left(\frac{r-R}{d_{K}}\right)\right],$ (7)

where $K_{\rm top}={2\over 3}F_{\rm bot}$, with $F_{\rm bot}=L/4\pi r_{\rm
in}^{2}$ where $L$ is the luminosity of the star, and where $d_{K}=0.015R$.
The contribution $K_{2}$ is given by

$\displaystyle
K_{2}(\rho,T)\,=\,K_{0}(\rho/\rho_{0})^{-(a+1)}(T/T_{0})^{3-b},$ (8)

where $\rho_{0}$ and $T_{0}$ are reference values of density and temperature,
and the values $a=1$ and $b=-7/2$ correspond to the Kramers opacity law. This
formulation was first used in convection simulations by Brandenburg et al.
(2000).

The subgrid scale (SGS) flux is given by

$\displaystyle{\bm{F}}^{\rm SGS}\,=\,-\,\chi_{\rm
SGS}\rho\bm{\nabla}s^{\prime},$ (9)

where $\chi_{\rm SGS}$ is the (constant) SGS diffusion coefficient for the
entropy fluctuation $s^{\prime}(r,\theta,\phi)=s-\langle
s\rangle_{\theta\phi}$, where $\langle s\rangle_{\theta\phi}$ is the
spherically symmetric part of the specific entropy. The SGS flux does not
contribute to the net radial energy transport because it is decoupled from the
mean stratification, and therefore changing $\chi_{\rm SGS}$ does not lead to
drastic changes in the boundary layer thickness near the surface.

The simulations were made using the Pencil Code111https://github.com/pencil-
code/ (Pencil Code Collaboration et al. 2021). In the present study the code
employs third-order temporal and sixth-order spatial discretisation. Advective
terms in Eqs. (1) to (4) are written as fifth-order upwinding derivatives with
a sixth-order hyperdiffusive correction where the diffusion coefficient is
flow-dependent; see Appendix B of Dobler et al. (2006).

###  2.1 System parameters and diagnostics quantities

The simulations are defined by the energy flux imposed at the bottom boundary,
$F_{\rm bot}=-(K\partial T/\partial r)|_{r=r_{\rm in}}$, the values of
$K_{0}$, $a$, $b$, $\rho_{0}$, $T_{0}$, $\Omega_{0}$, $\nu$, $\eta$,
$\chi_{\rm SGS}$, the profile of $K$, and the value of the modified Stefan-
Boltzmann constant $\sigma_{SB}$ in the upper boundary condition $\sigma_{\rm
SB}T_{\rm surf}^{4}=K\partial T/\partial r$, where $T_{\rm surf}$ is the
(unconstrained) surface temperature. The current models use a significantly
enhanced luminosity in comparison to real stars to bring the thermal and
dynamical timescales close enough to be resolved in the simulations. This
leads to correspondingly higher convective velocities and therefore the
rotation rate is increased accordingly to capture a similar rotational
influence on the flow in the simulations in comparison to real stars; see
appendix A of Käpylä et al. (2020).

The non-dimensional luminosity is given by

$\displaystyle\mathcal{L}\,=\,\frac{L_{0}}{\rho_{0}(GM)^{3/2}R^{1/2}},$ (10)

where $\rho_{0}$ is the initial density at the base of the convection zone.
The degree of luminosity enhancement is given by the ratio $L_{\rm
ratio}=\mathcal{L}/\mathcal{L}_{\odot}\approx 2.1\cdot 10^{5}$, where
$\mathcal{L}_{\odot}$ is the dimensionless solar luminosity. The initial
stratification is determined by the non-dimensional pressure scale height at
the surface

$\displaystyle\xi_{0}\,=\,\frac{\mathcal{R}T_{1}}{GM/R},$ (11)

where $T_{1}=T(R,t=0)$.

The relative strengths of viscosity, SGS diffusion, and magnetic diffusivity
are given by the SGS and magnetic Prandtl numbers

$\displaystyle{\rm Pr}_{\rm SGS}=\frac{\nu}{\chi_{\rm SGS}},\ \ {\rm
Pm}=\frac{\nu}{\eta}.$ (12)

We use ${\rm Pm}=1$ in most of the runs and vary ${\rm Pr}_{\rm SGS}$ between
0.1 and 10\. The thermal Prandtl number related to the radiative conductivity
is given by

$\displaystyle\Pr=\frac{\nu}{\chi},$ (13)

where $\chi=K/c_{\rm P}\rho$ is the radiative diffusivity, which in general
varies as a function of radius, latitude, and time. In the current simulations
$\chi_{\rm SGS}\gg\chi$ almost everywhere. The efficiency of convection is
quantified by the Rayleigh number

$\displaystyle{\rm Ra}\,=\,\frac{GM(\Delta r)^{4}}{\nu\chi_{\rm
r}R^{2}}\left(-\frac{1}{c_{\rm P}}\frac{{\rm d}s_{\rm hs}}{{\rm
d}r}\right)_{r_{\rm s}},$ (14)

where $\Delta r=0.3R$ is the depth of the layer, $s_{\rm hs}$ is the specific
entropy in a one-dimensional non-convecting hydrostatic model, evaluated near
the top of the domain at $r_{\rm s}=0.95R$, and where $\chi_{\rm s}$ is the
total thermal diffusivity $K/c_{\rm P}\rho$ from $r=r_{\rm s}$. The
hydrostatic solution is Schwarzschild-unstable only in a thin layer near the
surface (see, e.g. Barekat & Brandenburg 2014; Brandenburg 2016) which is why
the Rayleigh number is evaluated at $r_{\rm s}$. Moreover, $\chi_{\rm SGS}$
does not contribute to ${\rm Ra}$ because it only acts on deviations from the
spherically symmetric specific entropy. Additionally, a turbulent Rayleigh
number is quoted:

$\displaystyle{\rm Ra}_{\rm t}\,=\,\frac{GM(\Delta r)^{4}}{\nu\chi_{\rm
tot}R^{2}}\left(-\frac{1}{c_{\rm P}}\frac{{\rm d}\langle
s\rangle_{\theta\phi}}{{\rm d}r}\right)_{r_{\rm s}},$ (15)

where $\langle s\rangle_{\theta\phi}$ is the time- and horizontal average of
the specific entropy and $\chi_{\rm tot}=\chi_{\rm
SGS}+\langle\chi\rangle_{\theta\phi}$ is the total thermal diffusivity. ${\rm
Ra}_{\rm t}$ is always significantly smaller than ${\rm Ra}$ because
$\chi_{\rm SGS}\gg\chi$.

The magnitude of rotation is controlled by the Taylor number

$\displaystyle{\rm Ta}=\frac{4\Omega_{0}^{2}(\Delta r)^{4}}{\nu^{2}}.$ (16)

The fluid and magnetic Reynolds numbers and the Péclet number are given by

$\displaystyle{\rm Re}=\frac{u_{\rm rms}}{\nu k_{1}},\ \ {\rm Re}_{\rm
M}=\frac{u_{\rm rms}}{\eta k_{1}},\ \ {\rm Pe}=\frac{u_{\rm rms}}{\chi_{\rm
SGS}k_{1}},$ (17)

respectively, where $u_{\rm rms}=\sqrt{{\textstyle{3\over
2}}(U_{r}^{2}+U_{\theta}^{2})}$ is the time- and volume averaged rms velocity
where $U_{\phi}^{2}$ has been replaced by $(U_{r}^{2}+U_{\theta}^{2})/2$ to
avoid contributions from differential rotation. The inverse of the wavenumber
$k_{1}=2\pi/\Delta r\approx 21/R_{\odot}$ is used to characterize the radial
extent of the convection zone. Several definitions of the Coriolis number that
describes the rotational influence on the flow are discussed in Sect. 3.1.

Mean quantities are denoted by overbars are defined by the time- and azimuthal
average:

$\displaystyle\overline{f}(r,\theta)=\frac{1}{\Delta\phi\Delta
t}\int_{t_{0}}^{t_{0}+\Delta
t}\\!\\!\int_{0}^{\Delta\phi}f(r,\theta,\phi,t)d\phi dt,$ (18)

where $t_{0}$ and $\Delta t$ are the beginning and the length of the
statistically steady part of the simulation, and where $\Delta\phi=\pi/2$ is
the azimuthal extent of the simulation domain. Error estimates are obtained by
dividing the time series in three parts and computing averages over each one
of them. The largest deviation of these sub-averages from the average over the
whole time series is taken to represent the error.

Table 1: Summary of the runs.

$\begin{array}[]{p{0.08\linewidth}ccccccccccc}\hline\cr\hline\cr\vskip 3.0pt
plus 1.0pt minus 1.0pt\cr Run&{\rm Pr}_{\rm SGS}&{\rm Ta}[10^{6}]&{\rm
Co}&{\rm Co}_{\ell}&{\rm Co}_{\omega}&{\rm Co}_{\star}&{\rm Re}&{\rm Pe}&{\rm
Re}_{\rm M}&{\rm Ra}_{\rm t}[10^{6}]&{\rm DR}\\\ \hline\cr
P01-1H&0.1&2.6&1.14&1.10&0.59&0.38&36&3.6&-&0.17&\mbox{AS}\\\
P01-2H&0.1&3.1&1.25&1.15&0.63&0.42&36&3.6&-&0.19&\mbox{AS}\\\
P01-3H&0.1&3.7&1.34&1.21&0.66&0.46&37&3.7&-&0.21&\mbox{(AS)}\\\
P01-4H&0.1&4.4&1.42&1.22&0.67&0.50&37&3.7&-&0.23&\mbox{(AS)}\\\
P01-5H&0.1&5.1&1.50&1.28&0.70&0.54&38&3.8&-&0.25&\mbox{(SL)}\\\
P01-6H&0.1&5.8&1.65&1.36&0.75&0.58&37&3.7&-&0.27&\mbox{SL}\\\ \hline\cr
P01-1M&0.1&2.6&1.21&1.12&0.62&0.38&34&3.4&34&0.19&\mbox{AS}\\\
P01-1Mh&0.1&29&1.25&0.85&0.50&0.38&109&11&109&2.2&\mbox{AS}\\\
P01-2M&0.1&3.1&1.29&1.17&0.65&0.42&35&3.5&35&0.20&\mbox{AS}\\\
P01-2Mh&0.1&35&1.33&0.74&0.40&0.42&113&11&113&2.4&\mbox{(SL)}\\\
P01-3M&0.1&3.7&1.35&1.20&0.66&0.46&36&3.6&36&0.21&\mbox{AS}\\\
P01-3Mh&0.1&41&1.49&0.80&0.43&0.46&110&11&110&2.6&\mbox{SL}\\\
P01-4M&0.1&4.4&1.48&1.24&0.70&0.50&36&3.6&60&0.24&\mbox{(SL)}\\\
P01-5M&0.1&5.1&1.59&1.33&0.75&0.54&36&3.6&51&0.26&\mbox{SL}\\\
P01-6M&0.1&5.8&1.92&1.54&0.84&0.58&32&3.2&53&0.28&\mbox{SL}\\\ \hline\cr
P1-1H&1.0&2.6&1.31&1.00&0.55&0.38&31&31&-&1.4&\mbox{AS}\\\
P1-2H&1.0&3.1&1.48&1.13&0.61&0.42&30&30&-&1.5&\mbox{AS}\\\
P1-3H&1.0&3.7&1.63&1.15&0.64&0.46&30&30&-&1.7&\mbox{AS}\\\
P1-4H&1.0&4.4&1.80&1.21&0.68&0.50&30&30&-&1.8&\mbox{AS}\\\
P1-5H&1.0&5.1&1.97&1.25&0.73&0.54&29&29&-&2.0&\mbox{SL}\\\
P1-6H&1.0&5.8&2.14&1.32&0.78&0.58&29&29&-&2.1&\mbox{SL}\\\ \hline\cr
P1-1M&1.0&2.6&1.42&0.98&0.58&0.38&29&28&29&1.6&\mbox{AS}\\\
P1-2M&1.0&3.1&1.58&1.08&0.63&0.42&28&28&28&1.7&\mbox{AS}\\\
P1-2Mh&1.0&35&1.56&0.68&0.38&0.42&96&96&96&20&\mbox{(SL)}\\\
P1-3M&1.0&3.7&1.72&1.13&0.67&0.46&29&29&29&1.8&\mbox{AS}\\\
P1-3Mh&1.0&41&1.73&0.74&0.41&0.46&94&94&94&21&\mbox{SL}\\\
P1-4M&1.0&4.4&1.87&1.20&0.71&0.50&28&28&28&1.9&\mbox{SL}\\\
P1-5M&1.0&5.1&2.03&1.26&0.76&0.54&28&28&28&2.0&\mbox{SL}\\\
P1-6M&1.0&5.8&2.19&1.33&0.80&0.58&28&28&28&2.1&\mbox{SL}\\\ \hline\cr
P10-1M&10&2.6&1.75&0.99&0.58&0.38&23&233&23&13&\mbox{AS}\\\
P10-2M&10&3.1&1.98&1.06&0.64&0.42&23&227&23&14&\mbox{AS}\\\
P10-3M&10&3.7&2.16&1.16&0.70&0.46&23&227&23&15&\mbox{AS}\\\
P10-4M&10&4.4&2.31&1.26&0.74&0.50&23&229&23&15&\mbox{AS}\\\
P10-5M&10&5.1&2.53&1.31&0.80&0.54&23&226&23&16&\mbox{(AS)}\\\
P10-6M&10&5.8&2.68&1.40&0.84&0.58&23&228&23&17&\mbox{(SL)}\\\
P10-7M&10&6.6&2.80&1.46&0.88&0.62&23&234&23&17&\mbox{(AS)}\\\
P10-8M&10&7.5&3.10&1.56&0.95&0.65&22&224&22&18&\mbox{SL}\\\
P10-9M&10&8.4&3.28&1.62&1.00&0.69&22&224&22&18&\mbox{SL}\\\
\hline\cr\end{array}$

222Hydrodynamic (MHD) runs are denoted by suffix H (M). Majority of the runs
uses grid resolution $128\times 288\times 144$, the exception being runs
denoted by the suffix h where resolution $256\times 576\times 288$ was used.
The Rayleigh number in most runs is $7.0\cdot 10^{8}$, expect in the higher
resolution runs, denoted by suffix h, where it is ${\rm Ra}=2.4\cdot 10^{9}$.
${\rm Co}$, ${\rm Co}_{\ell}$, ${\rm Co}_{\omega}$, and ${\rm Co}_{\star}$
denote the definitions of the Coriolis number given by Eqs. (24), (26), (27),
and (29). ${\rm Pr}_{\rm M}=1$ in all runs except in P01-4M and P01-6M where
${\rm Pr}_{\rm M}=1.67$, and in P01-5M where ${\rm Pr}_{\rm M}=1.43$ to obtain
growing dynamos. The last column denotes the type of differential rotation
with parenthesis indicating that the result is not statistically significant.

Figure 1: Time-averaged rotation profiles from MHD runs P01-[2-5]M near the
AS-SL transition with ${\rm Pr}_{\rm SGS}=0.1$ with rotation rate increasing
from left to right. Coriolis numbers according to Eq. (24) and the mean
differential rotation at the equator $\langle\Delta\tilde{\Omega}_{\rm
eq}\rangle$ according to Eq. (23) are indicated in each panel. The white
arrows indicate the meridional flow.

###  2.2 Initial and boundary conditions

Initially the stratification is isentropic with polytropic index $n=1.5$ and
$\xi_{0}=0.02$, resulting in an initial density contrast of 30\. The value of
$K_{0}$ is chosen such that $F_{\rm rad}=F_{\rm tot}$ at the bottom of the
domain.

The radial and latitudinal boundaries are assumed impenetrable and stress-free
for the flow. On the bottom boundary, a fixed heat flux is prescribed while at
the top a black body condition is applied. On the latitudinal boundaries, the
gradients of thermodynamic quantities are set to zero; see Käpylä et al.
(2013). For the magnetic field we apply a radial field condition at the upper,
and a perfect conductor condition at the lower boundary. On the latitudinal
boundaries the field is assumed to be tangential to the boundary. These
conditions are given in terms of the magnetic vector potential by:

$\displaystyle A_{r}$ $\displaystyle=$ $\displaystyle 0,\ \ \ \frac{\partial
A_{\theta}}{\partial r}=-\frac{A_{\theta}}{r},\ \ \ \frac{\partial
A_{\phi}}{\partial r}=-\frac{A_{\phi}}{r}\ \ \ (r=R),$ (19)
$\displaystyle\frac{\partial A_{r}}{\partial r}$ $\displaystyle=$
$\displaystyle A_{\theta}=A_{\phi}=0\ \ \ (r=r_{\rm in}),$ (20) $\displaystyle
A_{r}$ $\displaystyle=$ $\displaystyle\frac{\partial
A_{\theta}}{\partial\theta}=A_{\phi}=0\ \ \
(\theta=\theta_{0},\pi-\theta_{0}).$ (21)

The azimuthal direction is periodic for all quantities. The velocity and
magnetic fields are initialized with random low-amplitude Gaussian noise
fluctuations.

##  3 Results

Three sets of simulations were done where ${\rm Pr}_{\rm SGS}=0.1$ (set P01),
$1$ (P1), and $10$ (P10), respectively. The first two sets contain
hydrodynamic and MHD runs, and a subset of the MHD runs were remeshed to
higher resolution and correspondingly higher Rayleigh, Péclet, and Reynolds
numbers; see Table 1. Only MHD variants of the P10 runs were run.

Figure 2: Same as Fig. 1 but for runs P1-[2-5]M in set P1. Figure 3: Same as
Fig. 1 but for runs P10-[6-9]M in set P10.

Figure 4: Normalized average angular velocity at the equator
$\langle\Delta\tilde{\Omega}_{\rm eq}\rangle$ for all runs as a function of
${\rm Co}$ (a), ${\rm Co}_{\ell}$ (b), ${\rm Co}_{\omega}$ (c), and ${\rm
Co}_{\star}$ (d). Circles (crosses) denote MHD (hydrodynamic) runs and the
colour of the symbols indicates the SGS Prandtl number. The sizes of the
symbols are proportional to ${\rm Re}_{\rm M}$ (${\rm Re}$) for MHD
(hydrodynamic) runs.

###  3.1 Differential rotation and meridional circulation

The main focus of the current study is to explore the effects of the Prandtl
number for the large-scale flows that develop in rotating convective systems.
The mean rotation profile is given by time- and azimuthal average:

$\displaystyle\overline{\Omega}(r,\theta)=\Omega_{0}+\frac{\overline{U}_{\phi}(r,\theta)}{r\sin\theta},$
(22)

and the meridional flow is given by $\overline{\bm{U}}_{\rm
mer}=(\overline{U}_{r},\overline{U}_{\theta},0)$. In many simulations the
latitudinal profiles of $\overline{\Omega}$ are non-monotonic such that the
rotation rate has a polar jet, a maximum at mid-latitudes or sometimes several
local minima and maxima as function of latitude. Furthermore, equatorial
asymmetries can occur, rendering the amplitude of the latitudinal shear an
unreliable diagnostic of the overall sense of differential rotation; see
representative examples in Figures 1 to 3. Therefore the classification of AS
and SL rotation profile is here based on the mean rotation profile at the
equator

$\displaystyle\langle\Delta\tilde{\Omega}_{\rm eq}\rangle=\frac{\int_{r_{\rm
in}}^{R}r^{2}[\tilde{\overline{\Omega}}(r,\theta_{\rm eq})-1]dr}{\int_{r_{\rm
in}}^{R}r^{2}dr},$ (23)

where $\theta_{\rm eq}=\pi/2$, and where the tildes refer to normalization by
the rotation rate of the frame of reference, $\Omega_{0}$. If
$\langle\Delta\tilde{\Omega}_{\rm eq}\rangle>0$
($\langle\Delta\tilde{\Omega}_{\rm eq}\rangle<0$) the run is classified as SL
(AS) rotator. This measure turns out to be a monotonic function of rotation
and it is furthermore unaffected by equatorial asymmetries or latitudinal
jets.

The current results indicate that the convective velocity increases when the
Prandtl number is decreased. This is manifested by increasing fluid Reynolds
number for decreasing SGS Prandtl number; see the eighth column of Table 1.
Naively one could then expect that achieving SL differential rotation for low
${\rm Pr}_{\rm SGS}$ would be more difficult, that is, require faster
rotation. Often the rotational influence on the flow is quantified by a simple
definition of the Coriolis number

$\displaystyle{\rm Co}\,=\,\frac{2\,\Omega_{0}}{u_{\rm rms}k_{1}},$ (24)

where the convective length scale is assumed to be unchanged by rotation.
Using this definition to characterize the results, the Coriolis number at
which the rotation profile flips from AS to SL appears to decrease
monotonically as ${\rm Pr}_{\rm SGS}$ decreases; see Fig. 4(a) where
$\langle\Delta\tilde{\Omega}_{\rm eq}\rangle$ is shown for all runs as a
function of ${\rm Co}$. That is, in the low resolution MHD runs with ${\rm
Pr}_{\rm SGS}=0.1$ ($1$) the transition occurs around ${\rm Co}\approx 1.5$
(${\rm Co}\approx 1.8$); see Figures 1 and 2, whereas for ${\rm Pr}_{\rm
SGS}=10$ the transition occurs at an even higher Coriolis number (${\rm
Co}\approx 3$); see Fig. 3.

However, the validity of this simplistic definition of the Coriolis number to
characterize the simulations can be questioned based on its very crude
estimate of the convective length scale and velocity amplitude. For example,
Gastine et al. (2014) argued that a local Rossby (inverse Coriolis) number
based on the length scale from the mean spherical harmonic degree
$\overline{\ell}_{u}$ of the $m\neq 0$ poloidal flows gives a more accurate
estimate (see also Schrinner et al. 2012). Furthermore, they showed that the
scatter near the AS-SL transition is reduced when using this definition,
essentially reducing the apparent dependence on Prandtl number significantly.
Here this was tested by computing $\overline{\ell}_{u}$ according to

$\displaystyle\overline{\ell}_{u}=\sum_{\ell}\ell\frac{\langle{\bm{U}}_{\rm
p}^{\ell}\bm{\cdot}{\bm{U}}_{\rm p}^{\ell}\rangle}{\langle{\bm{U}}_{\rm
p}\bm{\cdot}{\bm{U}}_{\rm p}\rangle},$ (25)

where ${\bm{U}}_{\rm p}$ is the non-axisymmetric poloidal flow and the
superscript $\ell$ refers to the corresponding spherical harmonic degree. Data
from a varying number of horizontal slices from near the base, at the middle,
and near the top of the CZ were analyzed for each run, and the resulting
$\overline{\ell}_{u}$ is an average over the depths and time. The number of
time slices per run varies between 7 and roughly 60. The corresponding length
scale is $d_{u}=\pi\Delta R/\overline{\ell}_{u}$, and the Coriolis number
based on this is given by

$\displaystyle{\rm Co}_{\ell}=\frac{2\Omega_{0}d_{u}}{u_{\rm rms}},$ (26)

where $u_{\rm rms}$ is defined the same way as in Eq. (17). The numbers in
Table 1, fifth column, indicate that the value of ${\rm Co}_{\ell}$ is
sensitive to the Reynolds number such that in the runs with the highest ${\rm
Re}$ the values of ${\rm Co}_{\ell}$ are roughly 30 per cent smaller than
those of the low resolution runs. Similarly, a Coriolis number based on
fluctuating vorticity can be defined as (e.g. Brun et al. 2022),

$\displaystyle{\rm Co}_{\omega}=2\Omega_{0}/\omega_{\rm rms}^{\prime},$ (27)

where ${\bm{\omega}}^{\prime}=\bm{\nabla}\times{\bm{u}}$ with
${\bm{u}}={\bm{U}}-\overline{\bm{U}}$. This quantity shows a similar
sensitivity to the fluid Reynolds number as ${\rm Co}_{\ell}$; see the sixth
column of Table 1. Both of these definitions pick up smaller length scales at
the highest Reynolds number runs resulting in lower Coriolis numbers. This is
likely an indication that the simulations are still far away from an
asymptotic regime where the results would be independent of the diffusivities.
Nevertheless, ${\rm Co}_{\ell}$ and ${\rm Co}_{\omega}$ characterize the
rotational influence on the flow more accurately than ${\rm Co}$ by being
sensitive to the actual dominant length scale. However, due to the Reynolds
number sensitivity, one should only compare results from runs with comparable
${\rm Re}$ when these definitions are used to characterize the results.

The average radial differential rotation at the equator is shown as functions
of ${\rm Co}_{\ell}$ and ${\rm Co}_{\omega}$ in Figures 4(b) and (c),
respectively. Ignoring the five runs at higher Reynolds and Péclet numbers for
the time being, these results suggest that the dependence of the differential
rotation transition as a function of the Prandtl number all but vanishes for
${\rm Pr}_{\rm SGS}<1$. However, in both cases the transition for ${\rm
Pr}_{\rm SGS}=10$ occurs still at a higher ${\rm Co}_{\ell}$ and ${\rm
Co}_{\omega}$ than for the ${\rm Pr}_{\rm SGS}=0.1$ and $1$ cases333Note that
only the two most rapidly rotating runs with ${\rm Pr}_{\rm SGS}=10$ have
statistically significant SL differential rotation; see the 12th column of
Table 1.. These results are in accordance with those reported by Karak et al.
(2018) who also found that a Prandtl number above unity promotes AS
differential rotation due to enhanced downward angular momentum transport.
Both definitions give much lower Coriolis numbers for the higher-${\rm Re}$
runs because smaller convective scales are resolved and $\overline{\ell}_{u}$
and $\omega^{\prime}$ pick these up. Furthermore, among the five higher-${\rm
Re}$ runs, the single AS model (P01-1Mh) appears to have a marginally larger
Coriolis number than the SL counterparts in both cases, although the
horizontal error estimates are large in both cases. Whether this is a real
effect or due to insufficient statistics remains open at this point. As a
side-note, the resemblance of Figures 4(b) and (c), or alternatively the
correlation between ${\rm Co}_{\omega}$ and ${\rm Co}_{\ell}$, suggests that
${\rm Co}_{\omega}$ captures the rotation dependence of the convective length
scale almost as well as ${\rm Co}_{\ell}$ without having to perform
numerically expensive spherical harmonic decomposition.

Each of the definitions of the Coriolis number discussed so far rely on
diagnostic quantities such as $u_{\rm rms}$, $\overline{\ell}_{u}$, and
$\omega^{\prime}$ that are sensitive to other details, such as the fluid
Reynolds number, of the system. Yet another alternative is to define a
Coriolis number that depends only on stellar input parameters such as the
luminosity, mass, and rotation rate of the modeled star. This can be
constructed by assuming that

$L=\rho_{\star}u_{\star}^{3}R^{2},$ (28)

where $L=4\pi r_{\rm in}^{2}F_{\rm bot}$ is the luminosity, $\rho_{\star}$ is
a reference density, and $u_{\star}$ is an estimate of an average convective
velocity. Here we assume $\rho_{\star}=\rho_{0}$444Using the average density
of the star, $\rho_{\rm av}=M/{4\over 3}\pi R^{3}$, is another option and in
that case the definition of ${\rm Co}_{\star}$ is fully determined by stellar
parameters. and construct a stellar Coriolis number

${\rm
Co}_{\star}=\frac{2\Omega_{0}R}{u_{\star}}=\frac{2\Omega_{0}R^{5/3}\rho_{0}^{1/3}}{L^{1/3}}.$
(29)

While this definition is imperfect in the sense that the actual convective
flow speed or scale do not enter, it is useful in determining whether a model
with a given rotation rate and luminosity is an AS or SL rotator. This makes
particularly sense in the homogeneous set of simulations considered here where
the stellar mass, luminosity, and radius are all fixed. The results are shown
in Fig. 4(d). It is apparent that the runs with ${\rm Pr}_{\rm SGS}=10$
require a significantly higher ${\rm Co}_{\star}$, corresponding here to
higher $\Omega_{0}$, to achieve SL differential rotation. However, the AS-SL
transition occurs at the same ${\rm Co}_{\star}$ for ${\rm Pr}_{\rm SGS}=1$
and $0.1$. Notably the MHD runs at ${\rm Pr}_{\rm SGS}=0.1$ and $1$ for ${\rm
Co}_{\star}=0.50$ are at least marginally SL in comparison to the
corresponding AS hydrodynamic cases. Similarly the higher resolution runs at
${\rm Co}_{\star}=0.46$ are SL or marginally SL for ${\rm Co}_{\star}=0.46$,
whereas the corresponding lower resolution MHD and hydrodynamic runs are all
AS. The current higher resolution runs with ${\rm Re}_{\rm M}\approx 94\ldots
113$ also have small-scale dynamos which was tested with separate runs where
the axisymmetric ($m=0$) magnetic fields were suppressed. Therefore it is
plausible that the main contribution to the earlier appearance of SL
differential rotation in these cases is due to the growing importance of the
magnetic fields. However, no corresponding higher resolution hydrodynamic runs
were made to confirm this. The current results are in accordance with the
results of Hotta & Kusano (2021) and Hotta et al. (2022) who argue in favour
of the magnetic field being the decisive factor in the transition. Finally, if
the simulations are scaled to physical units as in Käpylä et al. (2020), the
lowest value in the present study, ${\rm Co}_{\star}=0.38$, corresponds to the
case of solar rotation rate at solar luminosity.

Figure 5: Total angular momentum flux $\bm{\mathcal{T}}$ (arrows) superimposed
on the mean angular velocity profile (colour contours) from runs P1-[2-4]M and
P01-3Mh.

The conclusion of comparing the results characterized in terms of the four
variants of the Coriolis number is that SL differential rotation is
substantially more difficult to obtain for ${\rm Pr}_{\rm SGS}=10$ than for
${\rm Pr}_{\rm SGS}=1$ and $0.1$, whereas in the latter two cases there is no
clear difference. Furthermore, the use of the simplistic Coriolis number, Eq.
(24), gives misleading results and should be avoided. Contrary to the SGS
Prandtl number, the dependence on magnetic fields is clearer such that in MHD
runs SL differential rotation is easier to excite.

Figure 6: Time-averaged contributions to the radial angular momentum flux for
runs P1-2M (top row; panels a to d) and P1-4M (bottom; panels (e) to (h)). The
tildes refer to normalization by $\rho_{0}u_{\rm rms}^{2}$.

###  3.2 Angular momentum transport

When contributions from molecular viscosity can be neglected, the angular
momentum in the interior of the star is governed by the conservation equation

$\frac{\partial}{\partial
t}(\overline{\rho}{\mathcal{L}})+\bm{\nabla}\bm{\cdot}(r\sin\theta\bm{\mathcal{T}})=0,$
(30)

where ${\mathcal{L}}=r^{2}\sin^{2}\theta\overline{\Omega}$ is the specific
angular momentum and

$\displaystyle\bm{\mathcal{T}}=\overline{\rho}(\overline{{\bm{u}}u_{\phi}}+\overline{\bm{U}}\
\overline{U}_{\phi})-(\overline{{\bm{b}}b_{\phi}}+\overline{\bm{B}}\
\overline{B}_{\phi})/\mu_{0},$ (31)

is the total flux of angular momentum. The internal angular velocity profile
is thus determined by the spatial distribution of the angular momentum fluxes.
The main transporters are due to the Reynolds and Maxwell stress due to
fluctuating and mean flows and fields (cf. Rüdiger 1989; Rüdiger & Hollerbach
2004)

$\displaystyle\overline{Q}_{ij}^{\bm{u}}$ $\displaystyle=$
$\displaystyle\overline{\rho}\overline{u_{i}u_{j}},$ (32)
$\displaystyle\overline{M}_{ij}^{\bm{b}}$ $\displaystyle=$
$\displaystyle-\overline{b_{i}b_{j}}/\mu_{0},$ (33)
$\displaystyle\overline{Q}_{ij}^{\overline{\bm{U}}}$ $\displaystyle=$
$\displaystyle\overline{\rho}\overline{U}_{i}\ \overline{U}_{j},$ (34)
$\displaystyle\overline{M}_{ij}^{\overline{\bm{B}}}$ $\displaystyle=$
$\displaystyle-\overline{B}_{i}\ \overline{B}_{j}/\mu_{0},$ (35)

where ${\bm{b}}={\bm{B}}-\overline{\bm{B}}$ is the fluctuating magnetic field.

Representative results of $\bm{\mathcal{T}}$ are shown in Fig. 5 superimposed
on top of the angular velocity in Runs P1-[2-4]M and P01-3Mh. The total
angular momentum flux is nearly radially downward at high latitudes in the
bulk of the CZ. The region of downward radial flux is roughly confined inside
the tangent cylinder in all cases. Outside the tangent cylinder and near the
surface at all latitudes, $\bm{\mathcal{T}}$ is directed predominantly
equatorward irrespective whether the rotation profile is AS or SL. Outside the
tangent cylinder the angular momentum transport is increasingly axial although
the radial flux at the equator remains non-zero in all cases. The sense of the
differential rotation is related to the sign of the radial component of
$\bm{\mathcal{T}}$ at the equator: for positive (outward) flux the
differential rotation is SL (Runs P1-4M and P01-3Mh in Fig. 5) whereas it is
AS for negative (downward) flux (Run P1-2M). In the transitory case of Run
P1-3M the flux at the equator converges at the local maxima of
$\overline{\Omega}$ similarly as in the SL cases of P1-4M and P01-3Mh.
Therefore it appears sufficient to study the radial angular momentum flux to
determine the difference between AS and SL differential rotation.

Figure 6 shows all of the components of the radial angular momentum flux for
runs P1-2M (P1-4M) with AS (SL) differential rotation. In both cases the
Reynolds stress due to $m\neq 0$ contributions of the velocity is the dominant
contributors while the Reynolds stress due to meridional circulation and
Maxwell stresses are clearly subdominant. Moreover, in the AS case P1-2M the
radial transport due to the Reynolds stresses is downward almost everywhere,
while in the SL case P1-4M, both $\overline{Q}_{r\phi}^{\bm{u}}$ and
$\overline{Q}_{r\phi}^{\overline{\bm{U}}}$ are positive near the equator. The
sign of the latter near the equator is determined by the sign of
$\overline{U}_{\phi}$ because $\overline{U}_{r}>0$ in all cases; see the
meridional flow, for example, in Fig. 1.

Figure 7: Fourier-filtered and total Reynolds stress component
$\widetilde{\overline{Q}}_{r\phi}^{\bm{u}}$ from runs P1-2M (a) and P1-4M (b)
as indicated in the legend. The shaded areas indicate error estimates
according to the definition in Sect. 2.1.

To study the equatorial Reynolds stress in more detail, a number of slices of
the $m\neq 0$ velocity field ${\bm{U}}(r,\theta_{\rm eq},\phi)$ were analyzed.
Azimuthal Fourier filtering was applied to produce filtered velocity fields
${\bm{U}}_{\rm f}^{(m^{\prime}_{1},m^{\prime}_{2})}$, where azimuthal orders
ranging from $m^{\prime}_{1}$ to $m^{\prime}_{2}$ were retained555Note that
due to the $\Delta\phi=\pi/2$ azimuthal extent of the simulation domain,
$m^{\prime}$ corresponds to $4m$ in a full sphere.. These flows are used to
compute the Reynolds stress

$\overline{Q}_{r\phi}^{(m^{\prime}_{1},m^{\prime}_{2})}\\!=\\!\frac{1}{\Delta\phi\Delta
t}\int_{t_{0}}^{t_{0}+\Delta t}\\!\\!\\!\int_{0}^{\Delta\phi}\\!\\!\\!\rho
U_{r}^{(m^{\prime}_{1},m^{\prime}_{2})}U_{\phi}^{(m^{\prime}_{1},m^{\prime}_{2})}d\phi
dt.$ (36)

Here the density fluctuations are assumed to be small and no Fourier filtering
was applied to $\rho$. Representative results are shown in Fig. 7. In the AS
case P1-2M, Fig. 7(a), the total Reynolds stress is negative everywhere except
at the very base of the CZ. Contributions from the largest scale
$(m^{\prime}_{1},m^{\prime}_{2}=1,2)$ non-axisymmetric motions are
statistically almost identical with the total stress. A weak positive
contribution around the middle of the CZ is visible for
$(m^{\prime}_{1},m^{\prime}_{2}=3,5)$ but the Reynolds stress for
$(m^{\prime}_{1},m^{\prime}_{2}=1,5)$ is again very similar to the total
stress. In the SL runs the largest scales
$(m^{\prime}_{1},m^{\prime}_{2}=1,2)$ also contribute to a downward flux
whereas the main contribution to the net outward flux comes from
$(m^{\prime}_{1},m^{\prime}_{2}=3,5)$; see Fig. 7(b) for Run P1-4M. Similarly
to the AS case, the contributions from $m^{\prime}>5$ are very small. This
indicates that practically all of the Reynolds stress at the equator is due to
relatively large-scale structures which can be identified as the Busse columns
(Busse 1970a, b) which are also often referred to as banana cells. Such
features are often prominently visible in snapshots of the velocity field; see
Fig. 8 for a representative example. The Busse cells are manifestations of
non-linear prograde-propagating thermal Rossby waves. It is therefore somewhat
questionable to talk about turbulent Reynolds stress in this context since the
Busse columns are large-scale convective modes that appear essentially at a
scale corresponding to forcing of turbulence.

The mechanism by which the differential rotation is generated in the current
simulations is therefore different from that in Hotta et al. (2022) and Hotta
et al. (2022) where the small-scale Maxwell stress that is the dominant
contribution to the radial angular momentum transport. While the high-
resolution runs in the present study have small-scale dynamos and show an
increased tendency for SL differential rotation, the Reynolds stress due to
the thermal Rossby waves is still the dominant contribution to the angular
momentum flux in all of the runs considered here. In an earlier study (Käpylä
et al. 2017) the Maxwell stresses were found to be comparable to the Reynolds
stress at the highest magnetic Reynolds numbers, but in that study the modeled
stars were rotating typically three to four times faster than in the current
study. Although the Maxwell stress dominates the angular momentum transport in
the simulations of Hotta et al. (2022), large-scale Busse columns can still be
seen in the deep parts of their model; see, for example, their Fig. 6. If such
large-scale convective patters were as prominent in the Sun, they should have
been detected by helioseismology but there is no evidence currently to this
effect. Therefore it seems that although highly magnetized simulations are
more solar-like in terms of the rotation profile, the conundrum with the too
prominent large-scale structures remains.

Figure 8: Instantaneous normalized radial Reynolds stress component
$\widetilde{\overline{Q}}_{r\phi}^{\bm{u}}$ on the equatorial plane from run
P01-6M (colour contours). The $m\neq 0$ flows are indicated by the arrows, the
width of which is proportional to the local flow amplitude.

Furthermore, in hydrodynamic mean-field theories of differential rotation
(e.g. Rüdiger 1989; Kitchatinov & Rüdiger 2005; Rogachevskii & Kleeorin 2015)
the turbulence models are necessarily simplified and the large-scale
convection modes such as the Busse columns do not appear. In the most commonly
adopted approach of Kitchatinov & Rüdiger (2005), the radial angular momentum
transport is downward for slow, and vanishing for rapid rotation at the
equator. The SL differential rotation results in from a strong equatorward
transport. Numerical simulations of isothermal homogeneous anisotropic
turbulence also produce downward (slow rotation) or vanishing (rapid rotation)
radial angular momentum flux at the equator (Käpylä 2019a) in qualitative
accordance with Kitchatinov & Rüdiger (2005). Hydrodynamic mean-fields models
based on these concepts do not typically produce AS solutions unless strong
magnetic fields are present (Kitchatinov & Rüdiger 2004), although more
recently a hydrodynamic mechanism has also been discussed (Rüdiger et al.
2019). This latter process relies on poleward horizontal angular momentum flux
at slow rotation which was also found from local simulations, but which
appears to be absent in global models such as those presented here. The mean-
field theories avoid the problem of too prominent too strong thermal Rossby
waves by simply neglecting them, whereas they fail to characterize both AS and
SL cases in the current simulations. This tension is yet another facet of the
convective conundrum, the resolution of which is likely to require further
critical assessment of both theoretical and simulation approaches.

##  4 Conclusions

The transition from AS to SL differential rotation was studied as a function
of the SGS Prandtl number (${\rm Pr}_{\rm SGS}$). Four definitions of the
Coriolis number were used to quantify the exact point of transition from
simulations where the rotation of the star was varied. While this transition
occurs at a higher Coriolis number for ${\rm Pr}_{\rm SGS}=10$ than for ${\rm
Pr}_{\rm SGS}=1$ and $0.1$, no statistically relevant difference was found
between the last two cases. This suggests that the Prandtl number dependence
of the AS-SL transition is weak for ${\rm Pr}_{\rm SGS}<1$, whereas a high
Prandtl number makes it significantly more difficult to achieve SL
differential rotation.

These results are puzzling because earlier non-rotating local simulations
(Käpylä 2021) suggested that also the cases ${\rm Pr}=1$ and ${\rm Pr}=0.1$
differ significantly in many respects. A notable difference to the study of
Käpylä (2021) is that the current simulations do not include a radiative layer
below the CZ. This can explain why no subadiabatic layers develop at the base
of the CZs in the current simulations because the effects of overshooting are
absent. The latter was found to be particularly sensitive to the Prandtl
number in the local simulations Käpylä (2019b) and Käpylä (2021). The
inclusion of a radiative layer has also consequences for the dynamo solutions
(e.g. Guerrero et al. 2016; Käpylä 2022) which also couple back to
differential rotation. These aspects need to be revisited in future studies.

Many of the current simulations also included magnetic fields, albeit often in
a parameter regime where the small-scale dynamo is not excited. Thus the
influence of magnetic fields is relatively weak in most of the current runs.
Nevertheless, the magnetic fields make it easier to excite SL differential
rotation especially in the current higher resolution runs that likely also
have small-scale dynamos. However, the effects of magnetic fields are likely
to be more significant in more realistic higher-${\rm Re}_{\rm M}$ systems as
manifested by the recent results of Hotta et al. (2022). Therefore magnetism
appears to be the most promising candidate to explain the discrepancy between
solar observations and global simulations. Nevertheless, the non-detection of
thermal Rossby waves, which are still prominent in all current simulations,
from the Sun still raises questions as to the generation mechanism of solar
differential rotation. Finally, the tension between mean-field theories of
differential rotation and 3D simulation results is pointed out as another
aspect that requires further scrutiny in the future.

###### Acknowledgements.

I acknowledge the hospitality of Nordita during the program ‘The Shifting
Paradigm of Stellar Convection: From Mixing Length Concepts to Realistic
Turbulence Modelling’. The simulations were made within the Gauss Center for
Supercomputing project “Cracking the Convective Conundrum” in the Leibniz
Supercomputing Centre’s SuperMUC–NG supercomputer in Garching, Germany. This
work was supported by the Deutsche Forschungsgemeinschaft Heisenberg programme
(grant No. KA 4825/4-1).

## References

  * Barekat & Brandenburg (2014) Barekat, A. & Brandenburg, A. 2014, A&A, 571, A68
  * Bekki et al. (2017) Bekki, Y., Hotta, H., & Yokoyama, T. 2017, ApJ, 851, 74
  * Brandenburg (2016) Brandenburg, A. 2016, ApJ, 832, 6
  * Brandenburg et al. (2000) Brandenburg, A., Nordlund, A., & Stein, R. F. 2000, in Geophysical and Astrophysical Convection, Contributions from a workshop sponsored by the Geophysical Turbulence Program at the National Center for Atmospheric Research, October, 1995. Edited by Peter A. Fox and Robert M. Kerr. Published by Gordon and Breach Science Publishers, The Netherlands, 2000, p. 85-105, ed. P. A. Fox & R. M. Kerr, 85–105
  * Brun et al. (2004) Brun, A. S., Miesch, M. S., & Toomre, J. 2004, ApJ, 614, 1073
  * Brun et al. (2022) Brun, A. S., Strugarek, A., Noraz, Q., et al. 2022, ApJ, 926, 21
  * Busse (1970a) Busse, F. H. 1970a, ApJ, 159, 629
  * Busse (1970b) Busse, F. H. 1970b, Journal of Fluid Mechanics, 44, 441
  * Cattaneo et al. (1991) Cattaneo, F., Brummell, N. H., Toomre, J., Malagoli, A., & Hurlburt, N. E. 1991, ApJ, 370, 282
  * Dobler et al. (2006) Dobler, W., Stix, M., & Brandenburg, A. 2006, ApJ, 638, 336
  * Fan & Fang (2014) Fan, Y. & Fang, F. 2014, ApJ, 789, 35
  * Gastine et al. (2014) Gastine, T., Yadav, R. K., Morin, J., Reiners, A., & Wicht, J. 2014, MNRAS, 438, L76
  * Gilman (1983) Gilman, P. A. 1983, ApJS, 53, 243
  * Guerrero et al. (2016) Guerrero, G., Smolarkiewicz, P. K., de Gouveia Dal Pino, E. M., Kosovichev, A. G., & Mansour, N. N. 2016, ApJ, 819, 104
  * Guerrero et al. (2013) Guerrero, G., Smolarkiewicz, P. K., Kosovichev, A. G., & Mansour, N. N. 2013, ApJ, 779, 176
  * Hanasoge et al. (2016) Hanasoge, S., Gizon, L., & Sreenivasan, K. R. 2016, Annual Review of Fluid Mechanics, 48, 191
  * Hanasoge et al. (2012) Hanasoge, S. M., Duvall, T. L., & Sreenivasan, K. R. 2012, Proc. Natl. Acad. Sci., 109, 11928
  * Hotta & Kusano (2021) Hotta, H. & Kusano, K. 2021, Nature Astronomy, 5, 1100
  * Hotta et al. (2022) Hotta, H., Kusano, K., & Shimada, R. 2022, arXiv:2202.04183
  * Hotta et al. (2015) Hotta, H., Rempel, M., & Yokoyama, T. 2015, ApJ, 798, 51
  * Hotta et al. (2016) Hotta, H., Rempel, M., & Yokoyama, T. 2016, Science, 351, 1427
  * Käpylä (2019a) Käpylä, P. J. 2019a, A&A, 622, A195
  * Käpylä (2019b) Käpylä, P. J. 2019b, A&A, 631, A122
  * Käpylä (2021) Käpylä, P. J. 2021, A&A, 655, A78
  * Käpylä (2022) Käpylä, P. J. 2022, ApJ, 931, L17
  * Käpylä et al. (2020) Käpylä, P. J., Gent, F. A., Olspert, N., Käpylä, M. J., & Brandenburg, A. 2020, Geophysical and Astrophysical Fluid Dynamics, 114, 8
  * Käpylä et al. (2014) Käpylä, P. J., Käpylä, M. J., & Brandenburg, A. 2014, A&A, 570, A43
  * Käpylä et al. (2017) Käpylä, P. J., Käpylä, M. J., Olspert, N., Warnecke, J., & Brandenburg, A. 2017, A&A, 599, A4
  * Käpylä et al. (2013) Käpylä, P. J., Mantere, M. J., Cole, E., Warnecke, J., & Brandenburg, A. 2013, ApJ, 778, 41
  * Käpylä et al. (2019) Käpylä, P. J., Viviani, M., Käpylä, M. J., Brandenburg, A., & Spada, F. 2019, Geophysical and Astrophysical Fluid Dynamics, 113, 149
  * Karak et al. (2015) Karak, B. B., Käpylä, P. J., Käpylä, M. J., et al. 2015, A&A, 576, A26
  * Karak et al. (2018) Karak, B. B., Miesch, M., & Bekki, Y. 2018, Physics of Fluids, 30, 046602
  * Kitchatinov & Rüdiger (2004) Kitchatinov, L. L. & Rüdiger, G. 2004, Astron. Nachr., 325, 496
  * Kitchatinov & Rüdiger (2005) Kitchatinov, L. L. & Rüdiger, G. 2005, Astron. Nachr., 326, 379
  * Kupka & Muthsam (2017) Kupka, F. & Muthsam, H. J. 2017, Liv. Rev. Comp. Astrophys., 3, 1
  * Matilsky et al. (2020) Matilsky, L. I., Hindman, B. W., & Toomre, J. 2020, ApJ, 898, 111
  * Miesch & Toomre (2009) Miesch, M. S. & Toomre, J. 2009, Ann. Rev. Fluid Mech., 41, 317
  * Mitra et al. (2009) Mitra, D., Tavakol, R., Brandenburg, A., & Moss, D. 2009, ApJ, 697, 923
  * O’Mara et al. (2016) O’Mara, B., Miesch, M. S., Featherstone, N. A., & Augustson, K. C. 2016, Adv. Space Res., 58, 1475
  * Ossendrijver (2003) Ossendrijver, M. 2003, A&A Rev., 11, 287
  * Pencil Code Collaboration et al. (2021) Pencil Code Collaboration, Brandenburg, A., Johansen, A., et al. 2021, The Journal of Open Source Software, 6, 2807
  * Rogachevskii & Kleeorin (2015) Rogachevskii, I. & Kleeorin, N. 2015, Journal of Plasma Physics, 81, 395810504
  * Rüdiger (1989) Rüdiger, G. 1989, Differential Rotation and Stellar Convection. Sun and Solar-type Stars (Berlin: Akademie Verlag)
  * Rüdiger & Hollerbach (2004) Rüdiger, G. & Hollerbach, R. 2004, The Magnetic Universe: Geophysical and Astrophysical Dynamo Theory (Weinheim: Wiley-VCH)
  * Rüdiger et al. (2019) Rüdiger, G., Küker, M., Käpylä, P. J., & Strassmeier, K. G. 2019, A&A, 630, A109
  * Schrinner et al. (2012) Schrinner, M., Petitdemange, L., & Dormy, E. 2012, ApJ, 752, 121
  * Schumacher & Sreenivasan (2020) Schumacher, J. & Sreenivasan, K. R. 2020, Reviews of Modern Physics, 92, 041001
  * Simitev et al. (2015) Simitev, R. D., Kosovichev, A. G., & Busse, F. H. 2015, ApJ, 810, 80
  * Viviani et al. (2018) Viviani, M., Warnecke, J., Käpylä, M. J., et al. 2018, A&A, 616, A160
  * Weiss et al. (2004) Weiss, A., Hillebrandt, W., Thomas, H.-C., & Ritter, H. 2004, Cox and Giuli’s Principles of Stellar Structure (Cambridge, UK: Cambridge Scientific Publishers Ltd)