11institutetext: Dipartimento di Fisica, Università degli Studi di Milano, via
Celoria 16, 20133 Milano, Italy 22institutetext: Univ. Grenoble Alpes, CNRS,
IPAG, F-38000 Grenoble, France 33institutetext: School of Physics and
Astronomy, Monash University, Clayton, VIC 3800, Australia 44institutetext:
Department of Physics and Astronomy, University of Exeter, Stocker Road,
Exeter EX4 4QL, UK

# Probing initial distributions of orbital eccentricity and disc misalignment
via polar discs

Simone Ceppi<EMAIL_ADDRESS>Nicolás Cuello 22 Giuseppe Lodato 11
Cristiano Longarini 1133 Daniel J. Price 33 Daniel Elsender 44 Matthew R. Bate
44

(Accepted January 2024)

In a population of multiple protostellar systems with discs, the sub-
population of circumbinary discs whose orbital plane is highly misaligned with
respect to the binary’s orbital plane constrains the initial distribution of
orbital parameters of the whole population. We show that by measuring the
polar disc fraction and the average orbital eccentricity in the polar discs,
one can constrain the distributions of initial eccentricity and mutual
inclination in multiple stellar systems at birth.

###### Key Words.:

Accretion, accretion disks – Protoplanetary disks – Stars: protostars – Polar
alignment

##  1 Introduction

The tilt of an accretion disc orbiting an eccentric binary has two alternative
equilibrium configurations. If the mutual inclination (initial disc tilt with
respect to the stellar orbital plane) is below a critical angle for polar
alignment, (Farago & Laskar, 2010; Zanazzi & Lai, 2018; Cuello & Giuppone,
2019) the disc is expected to nodally precess, and align to the stellar
orbital plane due to viscous dissipation (Bate et al., 2000; Lubow & Ogilvie,
2000). If the initial mutual inclination is higher than the critical angle for
polar alignment, Aly et al. (2015) and Martin & Lubow (2017) showed that the
circumbinary disc angular momentum vector precesses around the binary orbit
eccentricity vector. Due to viscous dissipation, the disc evolves to a polar
configuration in which the mutual inclination is around 90∘ (depending on disc
parameters, see Martin & Lubow 2019).

Ceppi et al. (2023) showed that, in the general case of a disc in a multiple
stellar system, if the disc is orbiting more than two stars the polar
alignment mechanism is highly suppressed111Polar alignment could still occur
under the right conditions (e.g. very small semi-major axis ratios or for very
small discs). Ceppi et al. (2023) derived an analytical criterion to assess
the stability of the polar configuration in triples (see also the criterion
obtained by Lepp et al. 2023).. Conversely, if the disc is orbiting a pair of
stars with additional bodies outside the disc, polar alignment is at least as
likely as in the pure circumbinary disc case.

The properties of multiple stellar systems are already used to gain insights
into the physics driving stellar and planet formation evolution. The available
statistics on multiple stellar systems supports that the vast majority of
stars are born in a multiple stellar system and that stellar multiplicity
increases with stellar mass (Larson, 1972; Duchêne & Kraus, 2013; Offner et
al., 2022). These are crucial constraints for numerical experiments on the
collapse of molecular clouds (Bate et al., 2002; Krumholz et al., 2012; Bate,
2018, 2019; Mathew & Federrath, 2021; Mathew et al., 2023; Lebreuilly et al.,
2023). Orbital parameter distributions, such as the semi-major axis
distribution (Duquennoy & Mayor, 1991; Raghavan et al., 2010), or the binary
mass ratio distribution (Moe & Di Stefano, 2017; El-Badry et al., 2019) also
help to constrain formation mechanisms. The presence of more than two stars is
expected to change the shape of these distributions, encoding additional
information about the dependence on multiplicity (Smith et al., 1997; Ceppi et
al., 2022; Offner et al., 2022).

In this work, we show that, from the current distribution of mutual
misalignment and eccentricity — such as the one in Czekala et al. (2019) — and
with increased statistics, we are able to constrain the distributions of
eccentricity and inclination between orbital planes and discs at birth. Such
distributions, in turn, depend on which physical properties are more
significant in the formation process of these objects. In addition to
dynamical interactions between stars (Bate et al., 2002; Bate, 2018; Elsender
et al., 2023), the presence and the strength of magnetic fields (Price & Bate,
2007; Wurster et al., 2019; Zhao et al., 2020; Lebreuilly et al., 2021),
different metallicities (Elsender & Bate, 2021) and the level of turbulence in
the cloud (Bate et al., 2010; Walch et al., 2012) may play a significant role
in setting the initial properties of disc and star populations.

The paper is organised as follows: In Section 2 we present a model to compute
the fraction of polar discs and the mean eccentricity of stellar orbits
hosting polar discs. By applying this method to a population of hierarchical
systems and of pure binaries, we derive the relationship between the initial
and evolved distributions of system properties. Finally, in Section 3 we
constrain the initial conditions for both populations. In Section 4 we discuss
our results and we give our conclusions in Section 5.

##  2 Initial conditions and polar disc population

###  2.1 Polar disc fraction and mean polar systems eccentricity

As soon as the condition for polar alignment is satisfied, an accretion disc
starts oscillating around the polar configuration. The disc dissipates the
oscillation on a fraction of the viscous timescale (Lubow & Martin, 2018;
Zanazzi & Lai, 2018). Thus, it is reasonable to assume that all discs able to
go polar in the initial population will do so. Then, if we neglect the impact
of subsequent external interactions, the fraction of polarly aligned discs we
observe in a given evolved population is directly linked to the initial
conditions (eccentricity and misalignment distributions) in a forming
population.

In this section, we build a toy model to estimate the expected fraction of
polar discs and their eccentricity distribution in an evolved young stellar
population (Class II). Given the observed and theoretically predicted
preference for low mutual inclinations (e.g. Czekala et al. (2019); Elsender
et al. (2023)), we describe the initial distribution of mutual inclination
with a normalised exponential distribution:

$P_{\beta}(\beta)=\frac{1}{N_{\beta}}\exp\left(-\frac{\beta}{\sigma_{\beta}}\right),$
(1)

where $\beta$ is the mutual inclination, $\sigma_{\beta}$ is a parameter that
regulates the shape of the distribution and $N_{\beta}$ normalises the
distribution over the support considered, i.e. from $\beta=0$ to $\pi/2$.

The initial distribution of eccentricity is described by a normalised power
law distribution:

$P_{\rm e}(e)=\frac{1}{N_{\rm e}}e^{\alpha},$ (2)

where $e$ is the orbital eccentricity, $\alpha$ is a parameter regulating the
distribution shape and $N_{\rm e}$ normalises the distribution over $e=0$ to
$1$. This is in line with surveys of eccentricities (Duquennoy & Mayor, 1991;
Raghavan et al., 2010; Hwang et al., 2022) and it allows for a thermal
distribution ($P(e)\propto e$) that can be produced by repeated N-body
interactions (Jeans, 1919; Ambartsumian, 1937; Heggie, 1975).

Thus, the two-dimensional probability density function for the initial
condition is given by:

$P(e,\beta)=P_{\rm e}(e)P_{\beta}(\beta).$ (3)

For an orbit co-rotating with the central system, the critical angle for polar
alignment is222for a counterrotating orbit $\beta_{\rm crit}(e_{\rm
b},\Omega)=\pi-\arcsin{\sqrt{\frac{1-e_{\rm b}^{2}}{1-5e_{\rm
b}^{2}\cos{\Omega}^{2}+4e_{\rm b}^{2}}}}$ (Farago & Laskar, 2010; Zanazzi &
Lai, 2018; Cuello & Giuppone, 2019)

$\beta_{\rm
crit}(e,\Omega)=\arcsin{\sqrt{\frac{1-e^{2}}{1-5e^{2}\cos{\Omega}^{2}+4e^{2}}}},$
(4)

where $e$ is the orbital eccentricity and $\Omega$ is the disc longitude of
the ascending node.

Supposing that all discs above $\beta_{\rm crit}$ go polar and neglecting
subsequent external interactions, we can integrate over the initial mutual
inclination above $\beta_{\rm crit}$, to obtain the distribution of polar
discs as a function of system eccentricity $P_{\rm pol}(e)$. At this stage, we
consider $\Omega=\pi/2$ and we will fix $P_{\rm pol}(e)$ to take into account
a distribution of longitude of the ascending node in Sec. 2.2. Thus,

$P_{\rm pol}(e)=\int_{\beta_{\rm crit}(e,\frac{\pi}{2})}^{\pi/2}P(e,\beta){\rm
d}\beta.$ (5)

Integrating $P_{\rm pol}$ over the eccentricity we obtain the expected
fraction of polar discs ($F_{\rm p}$) in an evolved population, that is the
ratio between the number of polar discs and the total number of discs in the
population:

$F_{\rm p}=\int_{0}^{1}P_{\rm pol}(e){\rm d}e.$ (6)

Additionally, from $P_{\rm pol}$ we can compute the mean eccentricity of
stellar systems hosting polar discs, i.e.:

$\langle e\rangle=\int_{0}^{1}eP_{\rm pol}(e){\rm d}e.$ (7)

Integrals in Eqs. (6) and (7) are challenging to solve due to the dependencies
of $\beta_{\rm crit}$. In the next section we present two assumptions to take
into account a distribution of longitude of the ascending nodes in the case of
hierarchical systems and pure binaries.

###  2.2 Taking into account the longitude of ascending node

Figure 1: Critical angle for polar alignment as a function of the disc
longitude of the ascending node for a fixed binary eccentricity
($e=0.2,0.5,0.8$). Given a mutual inclination (e.g. $50^{\circ}$) and an
eccentricity (e.g. $e=0.5$), the dashed (dotted) line is the $\Omega$ interval
in which the disc will go perpendicular (coplanar) to the orbital plane.

Figure 1 shows how the critical angle $\beta_{\rm crit}$ in Eq. (4) depends on
the longitude of the ascending node ($\Omega$) for different orbital
eccentricities. Let us take the subpopulation of circumbinary discs with a
given mutual misalignment $\beta_{\rm sp}$ orbiting pure binaries with a given
orbital eccentricity $e_{\rm sp}$. We have a $\Omega_{\rm crit}$ so that the
critical angle for polar alignment $\beta_{\rm crit}(e_{\rm sp},\Omega_{\rm
crit})=\beta_{\rm sp}$. All discs with $\Omega<\Omega_{\rm crit}$ (dotted line
in Figure 1) will not polar align. All discs with a $\Omega>\Omega_{\rm crit}$
(dashed line) will polar align.

If we suppose a uniformly distributed $\Omega$, the fraction of discs that
will polar align is the ratio between the length of the dashed curve and the
width of the $\Omega$ interval (i.e. $\pi/2$).Thus, the fraction $f$ of discs
that will polarly align is:

$f(e_{\rm sp},\beta_{\rm sp})=1-\frac{2}{\pi}\Omega_{\rm crit}(e_{\rm
sp},\beta_{\rm sp}).$ (8)

In the case of a binary population, we have to take into account the factor
$f(e,\beta)$. To compute the distribution of polar discs with respect to
orbital eccentricity, Equation (5) becomes:

$P_{\rm pol}(e)=\int_{\beta_{\rm crit}(e)}^{\pi/2}P(e,\beta)f(e,\beta){\rm
d}\beta.$ (9)

In circumbinary discs in hierarchical systems, the precession of the
eccentricity vector drives the polar alignment process (Ceppi et al., 2023).
The orbit precesses on a shorter timescale than the timescale for coplanar
alignment. Hence, precession could lower the critical angle for polar
alignment to its minimum value because the system quickly explores the
$\Omega$ for which the critical angle is minimum. This is true only if, while
the system is exploring different longitudes of the ascending node, the tilt
of the disc does not decrease significantly. Otherwise, the polar alignment of
the disc would still be favoured compared to binary systems but the initial
longitude of the ascending node would nevertheless be relevant. If we suppose
this hypothesis to hold, for discs orbiting the inner binary of a triple what
really matters is the minimum critical angle no matter the longitude of
ascending node. Independent of $\Omega$, if the disc inclination is higher
than the minimum critical angle (the one for $\Omega=90^{\circ}$) the disc
will polar-align. Therefore, for triples $\beta_{\rm
crit}(e,\Omega)=\beta_{\rm crit}(e,\Omega=90^{\circ}).$ Thus, Eq. (5) can be
written as

$\begin{split}P_{\rm pol}(e)&=\int_{\beta_{\rm crit}(e)}^{\pi/2}P(e,\beta){\rm
d}\beta=\\\ &=-P_{\rm
e}(e)\frac{\sigma_{\beta}}{N_{\beta}}\left[\exp\left(-\frac{\pi/2}{\sigma_{\beta}}\right)-\exp\left(-\frac{\beta_{\rm
crit}(e)}{\sigma_{\beta}}\right)\right].\end{split}$ (10)

With the previous assumptions, we are left only with the eccentricity
dependence both for binaries and hierarchical systems.

###  2.3 Polar fraction and mean eccentricity for multiple systems

Figure 2: Fraction of polar discs ($F_{\rm p}$, leftmost column) and mean
orbital eccentricity of systems hosting polar discs ($\langle e\rangle$,
central column) for binary system population (top row) and triple system
population (bottom row). As shown by the rightmost column, each point of the
$\alpha-\sigma_{\beta}$ parameter space is uniquely characterized by a pair of
$F_{\rm p}$ and $\langle e\rangle$ (crossing of two contour lines, red contour
for polar fraction, black contour for mean eccentricity).

Given a pair of values $\alpha$ and $\sigma_{\beta}$ — the parameters of the
distributions of eccentricity and mutual angle, respectively — we are now able
to compute the expected polar disc fraction of an evolved population and the
mean eccentricity of systems hosting a polar disc with Equations (6) and (7),
respectively. For pure binaries and hierarchical systems, the polar disc
distributions are given by Eqs. (9) and (10), respectively.

We semi-analytically computed these integrals over the
$\alpha$–$\sigma_{\beta}$ parameter space for the initial distributions. The
parameter $\alpha$ ranges between $-1$ and $2$, while $\sigma_{\beta}$ ranges
between $0.01$ and $+\infty$. The left panels in Figure 2 show how the
fraction of expected polar discs in an evolved population $F_{\rm p}$ depends
on $\alpha$ and $\sigma_{\beta}$ for binary and hierarchical systems. In
general, the binary population is less prone to host polar discs compared to
systems with more than two stars. This is due to the different $\beta_{\rm
crit}$ we used to describe the two populations. The higher the probability of
having high eccentricity or mutual inclination, the more likely it is to find
configurations with a high mutual inclination and orbital eccentricity which
go polar more easily. Thus, $F_{\rm p}$ increases for higher $\alpha$ or
$\sigma_{\beta}$.

The central panels in Figure 2 show the mean eccentricity of orbits hosting
polar discs ($\langle e\rangle$) over the $\alpha$–$\sigma_{\beta}$ parameter
space for binaries and hierarchical systems. The binary polar population tends
to have higher $\langle e\rangle$, again due to the $f$ factor (Eq. (8)).
Indeed, the $\Omega$ interval allowing polar alignment is larger for more
eccentric systems. Thus, we have more polar discs around highly eccentric
systems. Lowering $\sigma_{\beta}$ raises $\langle e\rangle$. Indeed, higher
eccentricity values are required in regions where systems are mildly
misaligned. Conversely, lowering $\alpha$ means lowering $\langle e\rangle$.
The less likely is to have very eccentric systems in the initial population,
the less eccentric the polar population will be.

Each pair of $\alpha$–$\sigma_{\beta}$ uniquely connects to a pair of $F_{\rm
p}-\langle e\rangle$. This can be seen with contour lines in the right panel
of Figure 2. A given pair of polar fraction $F_{\rm p}$ and mean eccentricity
$\langle e\rangle$ contour lines cross only at one point of the
$\alpha$–$\sigma_{\beta}$ parameter space. Thus, we are able to numerically
invert the $\alpha$–$\sigma_{\beta}$ and $F_{\rm p}-\langle e\rangle$
coordinate systems to obtain plots in Figure 3. This plot showcases how
$\alpha$ and $\sigma_{\beta}$ depend on $F_{\rm p}$ and $\langle e\rangle$.
Doing so, once we have constrained $F_{\rm p}$ and $\langle e\rangle$ from
observations in a population we can constrain $\alpha$ and $\sigma_{\beta}$ in
the initial condition.

###  2.4 Measurement of polar fraction and mean eccentricity

Statistics on mutual inclinations between discs and stellar orbital planes are
fairly scarce at the moment. Even if it is relatively easy to measure the disc
inclination with respect to the sky plane, constraining the stellar orbital
inclination is challenging. Thus, in the most recent literature surveys
(Czekala et al., 2019; Zurlo et al., 2023) the sample size is around 15 discs.
Among these discs, only one accretion disc is confirmed in a polar
configuration: HD 98800B (Kennedy et al., 2019; Zúñiga-Fernández et al.,
2021). Additionally, two discs are likely polar333There is a degeneracy of
180∘ in the orbital longitude of the ascending node (Czekala et al., 2019).
(HD 142527, Balmer et al. (2022), and SR 24 N, Fernández-López et al. (2017)).
Another notable example is 99 Her, which is a polar debris disc that likely
evolved from a polar accretion disc (Smallwood et al., 2020). HD 98800B and SR
24N have additional companions outside the circumbinary disc. The semimajor
axis ratios for the two systems are about 0.02 and 0.01 for HD 98800 and SR
24, respectively. The presence of additional bodies could possibly affect the
process of polar alignment, for example, exciting Kozai-Lidov oscillations in
the disc. However, the outer orbits of these triples are too wide to
drastically impact the inner disc (Martin et al., 2014).

In the following, we show that the two different prescriptions for $\Omega$ in
binaries and hierarchical systems with more than two stars result in similar
values for the mean eccentricity, and in minor differences in terms of polar
fraction distributions compared to the current uncertainties we have in
surveys. Thus, here we do not distinguish between binaries and hierarchical
systems in the sample.

We estimate the mean eccentricity of systems hosting polar discs averaging
over confirmed polar systems (HD98800B), systems that are likely polar (SR24 N
and HD142527) and the polar debris disc 99 Her. We obtain an average
eccentricity of $\langle e\rangle^{\rm obs}=0.67\pm 0.11$.

To evaluate the polar disc fraction, if we take into account the single
confirmed polar system (HD98800B) we end up with a polar fraction of $0.08$.
If, however, we include also the additional two likely polar systems (SR24 N
and HD142527) the fraction raises up to $0.25$. We take the average of these
two values as a very rough estimate for the polar fraction, with their
standard deviation as an error. The result is $F_{\rm p}^{\rm obs}=0.17\pm
0.08$. This fraction should be considered as an upper limit. In a sample of
well-measured discs and orbital planes inclinations, there is a bias towards
systems in which those quantities of interest are measured properly (for
example because they are promising polar discs).

In the following, we use $F_{\rm p}^{\rm obs}$ and $\langle e\rangle^{\rm
obs}$ as references for exploring the parameter space of the initial
conditions.

##  3 Mapping observations onto the parameter space

(a)

(b)

(c)

(d)

Figure 3: $\alpha$ (left) and $\sigma_{\beta}$ (right) as a function of
$\langle e\rangle$ and $F_{\rm p}$ for binary (top) and triple (bottom)
populations. White regions are $F_{\rm p}-\langle e\rangle$ pairs that cannot
be generated by any $\alpha-\sigma_{\beta}$ pair. Purple curves are the
parameter space region in which $\alpha=2,1,0$ for solid, dotted and dashed,
respectively. Green curves are the parameter space region in which
$\sigma_{\beta}=\infty$. Black star is the solution for randomly distributed
mutual inclination and eccentricity. White box point with error bars is the
$F_{\rm p}^{\rm obs}$ and $\langle e\rangle^{\rm obs}$ observation values
derived in section 2.4. Red circle points are $F_{\rm p}-\langle e\rangle$
pairs measured in the four Bate (2019) molecular cloud collapse simulations
with different metallicities.

For a binary (triple) population, Figures 3(a) and 3(b) (3(c) and 3(d)) show
$\alpha$ and $\sigma_{\beta}$, respectively, as a function of $\langle
e\rangle$ and $F_{\rm p}$.

White regions are regions where no $\alpha$ and $\sigma_{\beta}$ pair results
in the given pair of $\langle e\rangle$ and $F_{\rm p}$. In the following, we
check which initial distribution is more likely to form the observed polar-
aligned disc population in binaries and triples.

###  3.1 Randomly distributed initial condition

The first assumption we test is flat distributions in both orbital
eccentricity and mutual inclination. This corresponds to taking the limit for
$\sigma_{\beta}$ approaching infinity and $\alpha=0$ in Eqs. (1) and (2)
respectively. The assumptions are: i) star formation at the molecular cloud
level gives no preferential orbital eccentricity; ii) there is no preferred
mutual inclination at the onset of stellar and disc formation.

Under these assumptions, we can compute analytically both the polar fraction
and the mean eccentricity of the triple polar population, solving Eqs. (6) and
(7) respectively with $P_{\rm e}=1$ and $P_{\beta}=2/\pi$. The results are
$F^{\rm tri}_{\rm p}=0.54$ and $\langle e\rangle^{\rm tri}=0.63$.

For the binary population, we numerically integrate Eq. (9) since the $f$
factor makes the integral not analytically solvable. Results for binaries are
$F^{\rm bin}_{\rm p}=0.40$ and $\langle e\rangle^{\rm bin}=0.65$. We find
these points in Fig. 3 at the crossing of the green and the dashed purple
curves, marked with a star. Indeed, green and dashed purple curves are the
limit for $\sigma_{\beta}$ approaching infinity and $\alpha=0$, respectively.
In surveys, we find a polar fraction $F_{\rm p}^{\rm obs}=0.17\pm 0.08$ and a
mean eccentricity $\langle e\rangle^{\rm obs}=0.67\pm 0.11$ (see Section 2.4).
Even if the mean eccentricity in this configuration is compatible with the
observed one, the observed polar fraction is not compatible with the expected
polar fraction for a binary/triple population.

These circumbinary discs form around young stars due to accretion of
surrounding gas and circularisation. Assuming that there is no correlation
between the forming system’s and the disc’ angular momenta, then rather than
assuming a flat distribution of mutual inclination, we expect
$P_{\beta}=\sin{\beta}$. In this case, we can take advantage of Eq. (16) in
Aly et al. (2015) which gives the fraction of configurations undergoing polar
precession for a randomly distributed $\Omega$, thus applicable to the binary
population. Integrating over the eccentricity we obtain an expected $F_{\rm
p}^{\rm bin}=1-\tanh{(2/\sqrt{5})}/\pi\approx 0.54>0.4$. This estimate is
higher than the one previously computed given that, in this case, higher
inclinations are favoured compared to coplanar configurations. Likewise, for
the triple population, we analytically solve Eq. (6) with
$P_{\beta}=\sin{\beta}$, obtaining $F_{\rm p}^{\rm tri}=(5-\sqrt{5})/4\approx
0.69>0.54$.

This result implies that non-uniform distributions — either in mutual
inclination and/or orbital eccentricity — are needed to explain the observed
polar population (or that many polar discs have been missed, which seems
unlikely).

###  3.2 Correlated orbit-disc mutual inclinations

The first hypothesis we relax is the random initial distribution of mutual
misalignment. We are still bound to move along the dashed purple curve in
Figures 3(b) and 3(d) (where $\alpha=0$). Over this restricted parameter space
region, $\langle e\rangle$ has a lower limit given by the
$\sigma_{\beta}=\infty$ case (i.e. 0.65 and 0.63 for binaries and triples,
respectively), while $F_{\rm p}$ can span values from $0$ to the
$\sigma_{\beta}=\infty$ case (i.e. 0.40 and 0.54 for binaries and triples,
respectively). The restricted parameter space region is compatible with
$\langle e\rangle^{\rm obs}$ and $F_{\rm p}^{\rm obs}$. In particular,
observations constrain $\sigma_{\beta}$ to range between $0.25$ and $0.78$ for
binaries and between $0.17$ and $0.43$ for triples which implies a narrow
$\beta$ distribution around $\beta=0$, i.e. close to aligned orbit-disc
systems should be more common at birth. Such correlation between the angular
momenta of the disc and of the binary/triple orbit is needed to describe the
observed polar population, under the assumption of randomly distributed
initial orbital eccentricities.

###  3.3 Non-flat initial conditions

We now allow $\alpha$ and $\sigma_{\beta}$ to explore the whole parameter
space to fully exploit the information contained in the measurement of
$\langle e\rangle$ and $F_{\rm p}$. The observed mean eccentricity and polar
fraction select a region of possible values of $\alpha$ and $\sigma_{\beta}$
highlighted by the error bars in Fig. 3. For binaries, we obtain $\alpha\leq
0.6$ and $0.26\leq\sigma_{\beta}\leq 1.8$. As for the triples, $\alpha\leq
0.46$ and $0.23\leq\sigma_{\beta}\leq 1.87$.

Again, the small statistics suggest the presence of a correlation between the
angular momenta of the discs and stellar systems. As for the eccentricity,
even if it is still marginally compatible with a random distribution, present
data suggest a decreasing function ($\alpha<0$). Indeed, with a flat or
slightly increasing eccentricity distribution, we would expect an higher mean
eccentricity or an higher polar fraction.

##  4 Discussion

###  4.1 Constraining initial conditions in multiple stellar systems

Our analysis suggests that, to be compatible with present data about polar
discs, initial distributions both for mutual inclinations and eccentricities
have to be non-uniform. While we expect the angular momenta of forming discs
and multiple systems to be correlated, the distribution of eccentricities we
find does not completely match with observational results from surveys of
evolved multiple stellar systems.

The observed evolved eccentricity distribution ranges from a uniform
distribution ($\alpha=0$) for orbits with semi-major axis of the order of 100
AU (Raghavan et al., 2010), up to an (increasing) thermal distribution
($\alpha=1$) for 500 AU and becomes even steeper for larger systems (Duquennoy
& Mayor, 1991; Tokovinin, 2020; Hwang et al., 2022). Conversely, we find an
$\alpha$ between $-1$ and $0.6$, pointing towards a slightly decreasing
distribution, only marginally compatible with a uniform one.

First, we have to consider that we do not observe systems hosting circum-
multiple discs with semi-major axis larger than about 100 AU (e.g. Czekala et
al., 2019) nor numerical simulations form them (e.g. Elsender et al., 2023),
thus we are interested in comparing our results with this range of semi-major
axis in which $\alpha$ is closer to our findings (systems of $\sim 200$ AU
have $\alpha\sim 0.6$). In addition, gravitational interactions tend to raise
the average eccentricity of the population. However, the thermalisation of the
distribution could be only partially explained with gravitational interaction
during cluster evolution which should act for much longer in order to push
$\alpha$ from 0 to 1 (Heggie, 1975; Weinberg et al., 1987). Thus, depending on
the thermalisation timescale, we could reconcile our findings with observed
distributions.

In addition, the orbital eccentricity can evolve through interactions between
the binary and the circumbinary disc. For example, D’Orazio & Duffell (2021)
numerically show that the hydrodynamical interaction with a coplanar disc
could lead to preferred values of binary eccentricity. Even if it is not clear
if this result is compatible with eccentricity surveys, the fact that our
model does not take into account that the eccentricity of observed systems
could be systematically different from the initial condition is also something
to investigate in the future. To the best of our knowledge, there are no
studies regarding this effect on misaligned circumbinary discs.

###  4.2 Comparison with cluster formation simulations

Nowadays, it is possible to perform detailed numerical simulations of the
collapse of molecular clouds. In such simulations, the initial conditions of
the cloud are set (e.g. amount of turbulence, strength of magnetic fields,
metallicity), and regulate the distribution of the parameters of the
population of the forming protosystems. By tuning the properties of the cloud,
one could aim at reproducing the measured distribution of inclination and
eccentricity. This would result in an indirect measurement of the molecular
cloud properties.

In this work, we applied this analysis to the molecular cloud collapse
simulations by Bate (2019). We measured $\alpha$ and $\sigma_{\beta}$ in the
newly born protostellar systems population by fitting the synthetic $P_{\rm
pol}$ (Eq. (5)) resulting from the simulation. We computed the expected
$F_{\rm p}$ and $\langle e\rangle$ with Equations (6) and (7), respectively.
Finally, we compared the values of $F_{\rm p}$ and $\langle e\rangle$ obtained
from the simulations with the ones measured in observations (see Section 2.4).

The simulation set consists of four molecular clouds collapsing with four
different metallicities (see details in Elsender & Bate 2021). We note that
the higher the metallicity, the steeper the eccentricity distribution becomes.
The steep decrease in eccentricity does not fit well with the same functional
form as the observed distributions, on which we based our parameterisation.
Thus, we are able to satisfactorily fit the three lower metallicities only.

The three red dots in Figure 3 represent the computed $F_{\rm p}$ and $\langle
e\rangle$ for three different realisations of the same cloud but with
different metallicities. The metallicity has an impact both on the
eccentricity and mutual inclination distribution, hence the scatter in Figure
3. However, it appears to have too little impact on the properties of the
population.

Regardless of the metallicity, the resulting polar disc fractions all are
generally too low compared to observations. Given that the distribution of
misalignment angles from the calculations of Bate (2019) and the observed
systems are in good agreement (Elsender et al., 2023), this mismatch must be
due to the lack of eccentric orbits in the simulations (possibly due to the
simulations being too dissipative). We note however that, as discussed in
Section 2.4, the proposed observational value has to be considered as an upper
limit. In other words, we cannot currently completely rule out such small
polar disc fractions from observations.

Finally, this analysis would benefit from numerical simulations producing
distributions of mutual inclination and eccentricity, but for different sets
of initial conditions (e.g. different amount of turbulence or magnetic field
strength) — such as Bate (2019) with metallicity. This would lead to more
robust predictions and better constraints for the polar disc population.
Provided accurate measurements of $F_{\rm p}$ and $\langle e\rangle$, the
method illustrated here constitutes a powerful way to infer the initial
conditions in molecular clouds from disc populations.

##  5 Conclusions

We showed how to measure the correlation between the inclination of accretion
discs and of forming stellar multiple systems (i.e. the distribution of mutual
inclination) and the distribution of orbital eccentricity of such stellar
systems at the onset of star and disc formation. Using our model, we were able
to compute the two fundamental parameters describing the initial distributions
of disc-orbit mutual inclination ($\sigma_{\beta}$) and orbital eccentricity
($\alpha$). The only required measurements are the fraction of polar discs in
a disc population ($F_{\rm p}$) and the mean eccentricity of systems hosting
polar discs ($\langle e\rangle$). Despite the low statistics available, we
find that:

  1. 1.

The observed disc population is not compatible with a randomly distributed
initial distribution of mutual misalignment – there must be a preference for
aligned systems;

  2. 2.

The orbital eccentricity is marginally compatible with a random distribution
as observed in field stellar systems with semi-major axis below $100$ au. The
observed increasing eccentricity for wider orbits is still compatible with
present data up to $\sim 200$ au. However, our model suggests a slight initial
preference for circular orbits. We will investigate in future works if this
discrepancy is compatible with the eccentricity evolution of young multiple
stellar systems.

The limitation of this proof-of-concept toy model lies in the simplified
distributions taken to describe the initial conditions. This simplification
allowed us to describe the distributions with only two parameters,
facilitating the computation of the relations between $\alpha-\sigma_{\beta}$
and $F_{\rm p}-\langle e\rangle$ pairs, the parameter space discussion, and
the comparison with data. This gives no degeneracy in the model and a strictly
two dimensional parameter space. A better observational constraint on the
polar disc fraction would improve the robustness of this method. Also, the
impact of the interaction with a polar disc on the eccentricity of the binary
system has to be investigated further.

In conclusion, by measuring the polar disc fraction and the distribution of
mutual inclinations, we showed that it is possible to infer the initial
eccentricity and mutual inclination distributions of binaries and triples at
birth. This will shed light on formation processes within molecular clouds
that affect the population of binary and triple stars.

## Acknowledgements

We thank the referee for the suggestions with which the manuscript
substantially improved. This project received funding from the European
Research Council (ERC) under the European Union Horizon Europe programme
(grant agreement No. 101042275, project Stellar-MADE). DJP acknowledges
Australian Research Council funding via DP220103767. SC’s visit to Monash was
funded by the EU Marie Curie RISE scheme DUSTBUSTERS (grant 823823). The ideas
discussed in this paper were partially formulated by SC, DJP, DE, and MRB
while attending the Binary22 programme at the Kavli Institute of Theoretical
Physics (KITP) at the University of California, Santa Barbara, which is
supported in part by the National Science Foundation under Grant No. NSF
PHY-1748958.

The data set consisting of the output from the calculations of Bate (2019)
that are analysed in this paper is available from the University of Exeter’s
Open Research Exeter (ORE) repository and can be accessed via the handle:
http://hdl.handle.net/10871/35993.

## References

  * Aly et al. (2015) Aly, H., Dehnen, W., Nixon, C., & King, A. 2015, MNRAS, 449, 65
  * Ambartsumian (1937) Ambartsumian, V. A. 1937, AZh, 14, 207
  * Balmer et al. (2022) Balmer, W. O., Follette, K. B., Close, L. M., et al. 2022, AJ, 164, 29
  * Bate (2018) Bate, M. R. 2018, MNRAS, 475, 5618
  * Bate (2019) Bate, M. R. 2019, MNRAS, 484, 2341
  * Bate et al. (2002) Bate, M. R., Bonnell, I. A., & Bromm, V. 2002, MNRAS, 336, 705
  * Bate et al. (2000) Bate, M. R., Bonnell, I. A., Clarke, C. J., et al. 2000, MNRAS, 317, 773
  * Bate et al. (2010) Bate, M. R., Lodato, G., & Pringle, J. E. 2010, MNRAS, 401, 1505
  * Ceppi et al. (2022) Ceppi, S., Cuello, N., Lodato, G., et al. 2022, MNRAS, 514, 906
  * Ceppi et al. (2023) Ceppi, S., Longarini, C., Lodato, G., Cuello, N., & Lubow, S. H. 2023, MNRAS[arXiv:2302.03411]
  * Cuello & Giuppone (2019) Cuello, N. & Giuppone, C. A. 2019, A&A, 628, A119
  * Czekala et al. (2019) Czekala, I., Chiang, E., Andrews, S. M., et al. 2019, ApJ, 883, 22
  * D’Orazio & Duffell (2021) D’Orazio, D. J. & Duffell, P. C. 2021, ApJ, 914, L21
  * Duchêne & Kraus (2013) Duchêne, G. & Kraus, A. 2013, ARA&A, 51, 269
  * Duquennoy & Mayor (1991) Duquennoy, A. & Mayor, M. 1991, A&A, 248, 485
  * El-Badry et al. (2019) El-Badry, K., Rix, H.-W., Tian, H., Duchêne, G., & Moe, M. 2019, MNRAS, 489, 5822
  * Elsender & Bate (2021) Elsender, D. & Bate, M. R. 2021, MNRAS, 508, 5279
  * Elsender et al. (2023) Elsender, D., Bate, M. R., Lakeland, B. S., Jensen, E. L. N., & Lubow, S. H. 2023, MNRAS, 523, 4353
  * Farago & Laskar (2010) Farago, F. & Laskar, J. 2010, MNRAS, 401, 1189
  * Fernández-López et al. (2017) Fernández-López, M., Zapata, L. A., & Gabbasov, R. 2017, ApJ, 845, 10
  * Heggie (1975) Heggie, D. C. 1975, MNRAS, 173, 729
  * Hwang et al. (2022) Hwang, H.-C., Ting, Y.-S., & Zakamska, N. L. 2022, MNRAS, 512, 3383
  * Jeans (1919) Jeans, J. H. 1919, MNRAS, 79, 408
  * Kennedy et al. (2019) Kennedy, G. M., Matrà, L., Facchini, S., et al. 2019, Nature Astronomy, 3, 230
  * Krumholz et al. (2012) Krumholz, M. R., Klein, R. I., & McKee, C. F. 2012, ApJ, 754, 71
  * Larson (1972) Larson, R. B. 1972, MNRAS, 156, 437
  * Lebreuilly et al. (2021) Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10
  * Lebreuilly et al. (2023) Lebreuilly, U., Vallucci-Goy, V., Guillet, V., Lombart, M., & Marchand, P. 2023, MNRAS, 518, 3326
  * Lepp et al. (2023) Lepp, S., Martin, R. G., & Lubow, S. H. 2023, ApJ, 943, L4
  * Lubow & Martin (2018) Lubow, S. H. & Martin, R. G. 2018, MNRAS, 473, 3733
  * Lubow & Ogilvie (2000) Lubow, S. H. & Ogilvie, G. I. 2000, ApJ, 538, 326
  * Martin & Lubow (2017) Martin, R. G. & Lubow, S. H. 2017, ApJ, 835, L28
  * Martin & Lubow (2019) Martin, R. G. & Lubow, S. H. 2019, MNRAS, 490, 1332
  * Martin et al. (2014) Martin, R. G., Nixon, C., Lubow, S. H., et al. 2014, ApJ, 792, L33
  * Mathew & Federrath (2021) Mathew, S. S. & Federrath, C. 2021, MNRAS, 507, 2448
  * Mathew et al. (2023) Mathew, S. S., Federrath, C., & Seta, A. 2023, MNRAS, 518, 5190
  * Moe & Di Stefano (2017) Moe, M. & Di Stefano, R. 2017, ApJS, 230, 15
  * Offner et al. (2022) Offner, S. S. R., Moe, M., Kratter, K. M., et al. 2022, arXiv e-prints, arXiv:2203.10066
  * Price & Bate (2007) Price, D. J. & Bate, M. R. 2007, Ap&SS, 311, 75
  * Raghavan et al. (2010) Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1
  * Smallwood et al. (2020) Smallwood, J. L., Franchini, A., Chen, C., et al. 2020, MNRAS, 494, 487
  * Smith et al. (1997) Smith, K. W., Bonnell, I. A., & Bate, M. R. 1997, MNRAS, 288, 1041
  * Tokovinin (2020) Tokovinin, A. 2020, MNRAS, 496, 987
  * Walch et al. (2012) Walch, S., Whitworth, A. P., & Girichidis, P. 2012, MNRAS, 419, 760
  * Weinberg et al. (1987) Weinberg, M. D., Shapiro, S. L., & Wasserman, I. 1987, ApJ, 312, 367
  * Wurster et al. (2019) Wurster, J., Bate, M. R., & Price, D. J. 2019, MNRAS, 489, 1719
  * Zanazzi & Lai (2018) Zanazzi, J. J. & Lai, D. 2018, MNRAS, 473, 603
  * Zhao et al. (2020) Zhao, B., Tomida, K., Hennebelle, P., et al. 2020, Space Sci. Rev., 216, 43
  * Zúñiga-Fernández et al. (2021) Zúñiga-Fernández, S., Olofsson, J., Bayo, A., et al. 2021, A&A, 655, A15
  * Zurlo et al. (2023) Zurlo, A., Gratton, R., Pérez, S., & Cieza, L. 2023, European Physical Journal Plus, 138, 411


# Short-distance HLbL contributions to the muon g-2 ∗

Johan Bijnens<EMAIL_ADDRESS>Nils Hermansson-Truedsson111Speaker,
Corresponding author<EMAIL_ADDRESS>Laetitia Laub<EMAIL_ADDRESS>Antonio Rodríguez-Sánchez<EMAIL_ADDRESS><EMAIL_ADDRESS>Department of Astronomy and Theoretical Physics,
Lund University, Sölvegatan 14A, SE 223 62 Lund, Sweden  Albert Einstein
Center for Fundamental Physics, Institute for Theoretical Physics, Universität
Bern, Sidlerstrasse 5, CH-3012 Bern, Switzerland  Present Address: Université
Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France

###### Abstract

The current $3.7\sigma$ discrepancy between the Standard Model prediction and
the experimental value of the muon anomalous magnetic moment could be a hint
for the existence of new physics. The hadronic light-by-light contribution is
one of the pieces requiring improved precision on the theory side, and an
important step is to derive short-distance constraints for this quantity
containing four electromagnetic currents. Here, we derive such short-distance
constraints for three large photon loop virtualities and the external fourth
photon in the static limit. The static photon is considered as a background
field and we construct a systematic operator product expansion in the presence
of this field. We show that the massless quark loop, i.e. the leading term, is
numerically dominant over non-perturbative contributions up to next-to-next-to
leading order, both those suppressed by quark masses and those that are not.

###### keywords:

Muon Anomalous Magnetic Moment , g-2 , Hadronic Light-by-Light , HLbL , Short-
Distance Constraints , Non-Perturbative Contributions , Operator Product
Expansion

##  1 Introduction

The experimentally measured value of the muon anomalous magnetic moment
$a_{\mu}=(g-2)_{\mu}/2$ is [1, 2]

$a^{\mathrm{exp}}_{\mu}=116\,592\,091(63)\times 10^{-11}\,.$ (1)

The Standard Model (SM) prediction on the other hand is [3]

$a_{\mu}^{\mathrm{SM}}=116\,591\,810(43)\times 10^{-11}\,,$ (2)

i.e. there is a $3.7\sigma$ discrepancy between theory and experiment. As a
consequence, the muon $g-2$ is an excellent low energy observable for the hunt
for new physics. With further improvement on precision it will be possible to
deduce the nature of this discrepancy.

The SM prediction receives contributions from several sectors, namely quantum
electrodynamics (QED), electroweak (EW) physics as well as the hadronic sector
[3]. The bulk of the value comes from QED which is known very precisely, and
the second most precise piece is the EW one. The hadronic sector is commonly
divided into two pieces, the hadronic vacuum polarisation (HVP) and the
hadronic light-by-light (HLbL). In numbers one has [3]

$\displaystyle a_{\mu}^{\mathrm{QED}}=116\,584\,718.931(104)\times
10^{-11}\,,$ (3) $\displaystyle a_{\mu}^{\mathrm{EW}}=153.6(1.0)\times
10^{-11}\,,$ (4) $\displaystyle a_{\mu}^{\mathrm{HVP}}=6845(50)\times
10^{-11}\,,$ (5) $\displaystyle a_{\mu}^{\mathrm{HLbL}}=92(18)\times
10^{-11}\,.$ (6)

The sum of these yields $a_{\mu}^{\mathrm{SM}}$. It is here clear that the
hadronic contributions dominate the uncertainty and therefore require further
consideration. Note that all higher order corrections to the HVP and HLbL here
have been included in $a_{\mu}^{\mathrm{HVP}}$ and $a_{\mu}^{\mathrm{HLbL}}$
even though we in the following will be interested in the leading order HLbL
contribution, which is usually denoted $a_{\mu}^{\mathrm{HLbL},\,\mathrm{LO}}$
and diagrammatically represented as in Fig. 1.

Figure 1:  A diagrammatic representation of the HLbL contribution to the
$g-2$. The grey blob contains hadrons and the photon momenta have been
labelled $q_{1,2,3,4}$, where the external photon is static, i.e.
$q_{4}\rightarrow 0$.

The HLbL is the time-ordered correlator of four electromagnetic (EM) currents,
and for the $g-2$ requires loop integration over three virtual photons with
momenta $q_{1}$, $q_{2}$ and $q_{3}$. The fourth, external, photon has
momentum $q_{4}$ which for the $g-2$ is in the static limit $q_{4}\rightarrow
0$. The situation is depicted in Fig. 1. In the following we will refer to the
Euclidean photon virtualities $Q_{1,2,3}^{2}=-q_{1,2,3}^{2}$. The HLbL process
can be calculated either in a data-driven approach using dispersion theory and
models, or using lattice QCD, see Ref. [3] and references therein. The data-
driven approach requires the knowledge of the HLbL tensor in different
kinematic regions, i.e. for different combinations of $Q_{i}^{2}$. The purely
short-distance (SD) region is given by $Q_{1}^{2}\sim Q_{2}^{2}\sim
Q_{3}^{2}\gg\Lambda_{\mathrm{QCD}}^{2}$. We derive SD constraints (SDCs) for
this particular kinematics by means of an operator product expansion (OPE) in
the presence of an external static field corresponding to the external photon
in Fig. 1 [4, 5]. Such SDCs provide information on omitted higher order
contributions in the dispersive approach and furthermore constrain model
calculations.

##  2 Generalities about the HLbL tensor

The HLbL tensor is a time-ordered correlation function of four EM currents.
The currents are of the form $J^{\mu}(x)=\bar{q}\,Q_{q}\gamma^{\mu}q$ where
$q=(u,d,s)$ and $Q_{q}=\textrm{diag}(e_{q})=\textrm{diag}(2/3,-1/3,-1/3)$ is
the light quark charge matrix. With this one has the tensor as

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(q_{1},q_{2},q_{3})=-i\int\frac{d^{4}q_{4}}{(2\pi)^{4}}\left(\prod_{i=1}^{4}\int
d^{4}x_{i}\,e^{-iq_{i}x_{i}}\right)$ $\displaystyle\times\langle
0|T\left(\prod_{j=1}^{4}J^{\mu_{j}}(x_{j})\right)|0\rangle\,.$ (7)

The Ward identities for $q_{1,2,3,4}$, with $q_{4}=-q_{1}-q_{2}-q_{3}$, can be
compactly written as

$q_{i,\,\mu_{i}}\,\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(q_{1},q_{2},q_{3})=0\,,$
(8)

from which it follows that we may write [6]

$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(q_{1},q_{2},q_{3})=-q_{4,\,\nu_{4}}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4,\,\mu_{4}}}(q_{1},q_{2},q_{3})\,.$ (9)

As a consequence, we may therefore calculate the HLbL by considering the
derivative in (9). Lorentz decomposing the HLbL tensor into a basis of 54
scalar functions $\Pi_{i}$ as in Refs. [7, 8] one has the HLbL contribution to
$a_{\mu}$ as

$\displaystyle
a_{\mu}^{\mathrm{HLbL}}=\frac{2\alpha^{3}}{3\pi^{2}}\int_{0}^{\infty}dQ_{1}\int_{0}^{\infty}dQ_{2}\int_{-1}^{1}d\tau\,\sqrt{1-\tau^{2}}\,$
$\displaystyle\times\,Q_{1}^{3}Q_{2}^{3}\sum_{i=1}^{12}T_{i}(Q_{1},Q_{2},\tau)\,\overline{\Pi}_{i}(Q_{1},Q_{2},\tau)\,.$
(10)

Here, $T_{i}$ are known functions and the $\overline{\Pi}_{i}$ are functions
of six linear combinations of the $\Pi_{i}$. The six linear combinations in
question are commonly denoted as $\hat{\Pi}_{1,4,7,17,39,54}$, and together
determine $a_{\mu}^{\mathrm{HLbL}}$.

As we are interested in the SD domain, we have to consider values of $Q_{i}$
greater than some $Q_{\mathrm{min}}$, which means that the integrals in (2)
have a lower cut-off. Where to choose this cut-off is a priori not clear and
is one of the sources of uncertainty in the HLbL prediction [3]. The
determination of $Q_{\mathrm{min}}$ lies beyond our scope for the moment and
we therefore let it be a variable quantity.

###  2.1 Constraints from the short-distance domain

As was explicitly shown in Ref. [4], an OPE of the HLbL tensor in (2) does not
lead to SDCs for the $g-2$ kinematics. The reason is that the systematic
expansion breaks down due to the appearance of quark propagators of the soft
momentum $q_{4}$. However, it is possible to consider the static external
field as a background in which one then can construct a systematic OPE of a
correlation function of three EM currents [4, 5]. The quantity to study then
is

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}}(q_{1},q_{2})=-\frac{1}{e}\int\frac{d^{4}q_{3}}{(2\pi)^{4}}\left(\prod_{i=1}^{3}\int
d^{4}x_{i}\,e^{-iq_{i}x_{i}}\right)$ $\displaystyle\times\langle
0|T\left(\prod_{j=1}^{3}J^{\mu_{j}}(x_{j})\right)|\gamma(q_{4})\rangle\,.$
(11)

Note that the static photon appears in the external state. The connection
between $\Pi^{\mu_{1}\mu_{2}\mu_{3}}$ and the HLbL tensor is obtained by
factoring out the external field according to [4]

$\Pi^{\mu_{1}\mu_{2}\mu_{3}}(q_{1},q_{2})=\langle
0|e_{q}F_{\nu_{4}\mu_{4}}|\gamma(q_{4})\rangle\times\frac{i}{2}\lim_{q_{4}\rightarrow
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}}{\partial
q_{4}^{\nu_{4}}}\,.$ (12)

An OPE here will therefore allow to obtain $\lim_{q_{4}\rightarrow
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}}{\partial q_{4}^{\nu_{4}}}$
which is needed to find the $\hat{\Pi}_{i}$ and from these also
$a_{\mu}^{\mathrm{HLbL}}$.

##  3 An OPE in an external EM field

We want to construct an OPE for the tensor $\Pi^{\mu_{1}\mu_{2}\mu_{3}}$, i.e.
where the external soft photon represents the background field. Such OPEs in
the presence of an external field have been considered before, e.g. for
nucleon magnetic moments in Refs. [9, 10] and also for the EW contributions to
the muon magnetic moment in Ref. [11].

An OPE is a systematic expansion of a product of operators, and when applied
to the EM currents in $\Pi^{\mu_{1}\mu_{2}\mu_{3}}$ it gives rise to a sum of
contributions depending on non-perturbative matrix elements, corresponding to
long-distance effects, multiplied by perturbative SD coefficients. In practice
this sum is obtained by expanding the time-ordered product of operators in
Dyson series and using Wick’s theorem. The non-perturbative matrix elements
come from not fully contracted terms in the Wick expansion. The resulting sum
has a systematic ordering with long- and short-distance effects separated.
Note that for each order one must analyse which operators can contribute for a
given OPE. The main difference between a vacuum OPE and ours with a background
field is that all operators which have the same quantum numbers as the
external field $F_{\mu\nu}$ can contribute.

We consider the OPE up to dimension $\alpha_{s}$ and
$(\Lambda_{\mathrm{QCD}}/Q_{\mathrm{min}})^{6}$. There are eight types of
operators which can contribute up to this order, namely

$\displaystyle S_{1,\,\mu\nu}=e\,e_{q}F_{\mu\nu}\,,\ $ (13) $\displaystyle
S_{2,\,\mu\nu}=\bar{q}\sigma_{\mu\nu}q\,,$ (14) $\displaystyle
S_{3,\,\mu\nu}=i\,\,\bar{q}G_{\mu\nu}q\,,$ (15) $\displaystyle
S_{4,\,\mu\nu}=i\,\bar{q}\bar{G}_{\mu\nu}\gamma_{5}q\,,$ (16) $\displaystyle
S_{5,\,\mu\nu}=\bar{q}q\;e\,e_{q}F_{\mu\nu}\,,$ (17) $\displaystyle
S_{6,\,\mu\nu}=\frac{\alpha_{s}}{\pi}\,G_{a}^{\alpha\beta}G^{a}_{\alpha\beta}\;e\,e_{q}F_{\mu\nu}\,,$
(18) $\displaystyle
S_{7,\,\mu\nu}=\bar{q}(G_{\mu\lambda}D_{\nu}+D_{\nu}G_{\mu\lambda})\gamma^{\lambda}q-(\mu\leftrightarrow\nu)\,,$
(19) $\displaystyle
S_{\\{8\\},\,\mu\nu}=\alpha_{s}\,(\bar{q}\,\Gamma\,q\;\bar{q}\Gamma
q)_{\mu\nu}\,.$ (20)

We have here defined $G_{\mu\nu}=ig_{S}\lambda^{a}G^{a}_{\mu\nu}$ as well as
its dual
$\bar{G}^{\mu\nu}=\frac{i}{2}\epsilon^{\mu\nu\lambda\rho}G_{\lambda\rho}$, and
$\Gamma$ is a matrix in spinor, colour and flavour space. For further
conventions see Ref. [5]. Operators $S_{\\{8\\},\,\mu\nu}$ involve flavour
mixing, as they are obtained by expanding one order in $\alpha_{s}$ and non-
contracted four-quark operators are obtained in the Wick expansion. In the
chiral limit one can obtain a basis of twelve possible four-quark operators,
as shown in Ref. [5].

The diagrams appearing in our OPE are shown in Fig. 2. This first of all shows
that the quark loop is the first term in a systematic OPE. This has for a long
time been a common assumption lacking an explicit derivation, and it has
further been assumed that it is a sufficiently good representation of the SD
HLbL behaviour [4, 5]. Note that we so far have not calculated the leading
perturbative $\alpha_{s}$ correction to the quark loop, but this is the
subject of an upcoming publication.

Figure 2:  A diagrammatic representation of the terms in the OPE up to
dimension 6 and order $\alpha_{s}$. The external field can appear either on
hard quark lines as an insertion vertex (marked by a cross) or in a non-
perturbative matrix element (in diagrams where the cross has been omitted).
The terms are ($a$) the quark loop, ($b_{1}$) the di-quark magnetic
susceptibility $\langle\bar{q}\sigma_{\mu\nu}q$, ($b_{2}$) the quark
condensate $\langle\bar{q}q\rangle$, ($c$) four-quark condensates
$\langle\bar{q}\Gamma_{1}q\,\bar{q}\Gamma_{2}q\rangle$ for colour, flavour and
Dirac matrices $\Gamma_{1,2}$, and ($d$) the gluon condensate
$\langle\alpha_{s}GG\rangle$.

###  3.1 Renormalisation

Simply following the above procedure for the OPE is not sufficient. The result
for $\Pi^{\mu_{1}\mu_{2}\mu_{3}}$ can be written

$\Pi^{\mu_{1}\mu_{2}\mu_{3}}(q_{1},q_{2})=\vec{C}^{\,T,\,\mu_{1}\mu_{2}\mu_{3}\mu_{4}\nu_{4}}(q_{1},q_{2})\,\vec{X}_{S}\,\langle
e_{q}F_{\mu_{4}\nu_{4}}\rangle\,,$ (21)

where the vector $\vec{X}_{S}$ has components $X_{S}^{i}$ defined through
$\langle S_{i,\,\mu\nu}\rangle=e\,e_{q}X_{S}^{i}\langle F_{\mu\nu}\rangle$. In
other words, the components $X_{S}^{i}$ are related to the magnetic
susceptibilities of the respective operators. The vector
$\vec{C}^{\,T,\,\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ contains the momentum
dependence. At this stage, however, long and short distances are not yet
completely separated. This can be understood since from non-contracted soft
quark lines there are contributions from the Dyson series giving rise to a
divergent series. In addition, in (21) one obtains corrections $\log
Q_{i}^{2}/m_{q}^{2}$ that lead to ill-defined series. We solve these issues by
dressing and renormalising the operators, in the $\overline{MS}$ scheme, as
explained further in Ref. [5]. After including operator mixing and defining
dressed and renormalised operators $\vec{Q}_{\overline{MS},\,\mu\nu}(\mu)$ at
scale $\mu$, the result in (21) takes the form

$\Pi^{\mu_{1}\mu_{2}\mu_{3}}(q_{1},q_{2})=\frac{1}{e}\vec{C}^{\,T,\,\mu_{1}\mu_{2}\mu_{3}\mu_{4}\nu_{4}}_{\overline{MS}}(q_{1},q_{2})\langle\vec{Q}_{\overline{MS},\mu_{4}\nu_{4}}(\mu)\rangle\,,$
(22)

where

$\langle\vec{Q}_{\overline{MS},\,\mu\nu}(\mu)\rangle=e\,\vec{X}\,\langle
e_{q}F_{\mu\nu}\rangle\,.$ (23)

The magnetic susceptibilities corresponding to the renormalised operators are
here contained in the vector. Here long and short distances have been
completely separated and no divergent mass logarithms appear. To conclude,
compared to the massless quark loop we have calculated contributions of order
$g_{s}\,\frac{\Lambda_{\mathrm{QCD}}^{4}}{Q^{4}}$, $\frac{m_{q}^{2}}{Q^{2}}$,
$g_{s}^{2}\,\frac{\Lambda_{\mathrm{QCD}}^{4}}{Q^{4}}$,
$m_{q}\frac{\Lambda_{\mathrm{QCD}}}{Q^{2}}$,
$m_{q}\frac{\Lambda_{\mathrm{QCD}}^{3}}{Q^{4}}$ and
$m_{q}^{3}\frac{\Lambda_{\mathrm{QCD}}}{Q^{4}}$.

###  3.2 Numerical estimates for the matrix elements

Having performed the renormalised OPE to the desired order, the only remaining
step to calculate $a_{\mu}^{\mathrm{HLbL}}$ is to determine the magnetic
susceptibilites in (23). First of all, $X_{5}$ and $X_{6}$ are respectively
given by the quark and gluon condensates. The quark condensate is well-studied
and its numerical value can be found in many places, e.g. Ref. [12]. The gluon
condensate was estimated in Ref. [13] to be $X_{6}\sim
0.02\,\mathrm{GeV}^{4}$. For $X_{7}$ we simply make a guess inspired from the
operator mixing, namely
$|X_{7}|\sim\frac{1}{6}\Big{\langle}\frac{\alpha_{s}}{\pi}GG\Big{\rangle}$.
For the four-quark operators only two combinations appear, which can be
estimated using large-$N_{c}$ arguments. The values coincide and are
$\overline{X}_{8,1}^{N_{c}\rightarrow\infty}=\overline{X}_{8,2}^{N_{c}\rightarrow\infty}=-2\frac{\pi\alpha_{s}}{9}X_{2}\langle\bar{q}{q}\rangle$.
The matrix element $X_{2}$ is the so-called di-quark magnetic susceptibility
of the vacuum and has been calculated on the lattice, see Ref. [14] for a
recent value. The only remaining matrix elements are $X_{3}$ and $X_{4}$ which
to our knowledge were hitherto unknown. In order to estimate them we connect
the matrix elements to vacuum QCD two-point functions and employ large-$N_{c}$
arguments. The results are
$X_{3}=-\frac{m_{0}^{2}}{6M_{\rho}^{2}}\langle\bar{q}q\rangle$ and
$X_{4}=-\frac{m_{0}^{2}}{6M_{\rho}^{2}}\langle\bar{q}q\rangle$, where
$m_{0}^{2}$ is a parameter estimated in Ref. [15]. Using the same approach for
$X_{2}$ yields $X_{2}=\frac{2}{M_{\rho}^{2}}\langle\bar{q}q\rangle$ which is
in excellent agreement with lattice determinations. For further details on the
numerical values of the matrix elements and analytic expressions from (22),
see Ref. [5].

##  4 Numerical results for $a_{\mu}^{\mathrm{HLbL}}$

Figure 3: Contributions to $a_{\mu}^{\mathrm{HLbL}}$from $X_{1,0}$,
$X_{1,m^{2}}$, $X_{2,m}$ and $X_{2,m^{3}}$.

Figure 4: Contributions to $a_{\mu}^{\mathrm{HLbL}}$from $X_{1,0}$, $X_{6}$,
$X_{7}$, $X_{8,1}$ and $X_{8,2}$.

Figure 5: Contributions to $a_{\mu}^{\mathrm{HLbL}}$from $X_{1,0}$, $X_{3}$,
$X_{4}$ and $X_{5}$.

In Figs. 4–5 we plot $a_{\mu}^{\mathrm{HLbL}}$ in (2) as a function of
$Q_{\mathrm{min}}$. We have here used $\alpha_{s}=0.33$, $m_{u}=m_{d}=5$ MeV,
$m_{s}=100$ MeV and the renormalisation scale $\mu=Q_{\mathrm{min}}$. The
massless quark loop given by $X_{1,0}$ clearly dominates with respect to the
other contributions to the OPE. This is true for both the contributions
suppressed by quark masses such as e.g. $X_{2,m}$ as well as those that are
not such as the four-quark pieces $X_{8,1}$ and $X_{8,2}$. As a consequence,
the massless quark loop seems to be a good representation of the SD HLbL
behaviour from relatively low energies. However, before such a statement can
be made, one has to study also the $\alpha_{s}$ correction to the massless
quark loop, i.e. by including two extra quark-gluon-antiquark QCD vertices in
the Dyson series. Our current preliminary evaluation shows that also this is
small compared to the massless quark loop.

##  5 Conclusions

We have developed a systematic OPE in the presence of an external EM field to
derive SDCs for the HLbL tensor. These constraints are needed to reduce the
error on the SM prediction of the muon anomalous magnetic moment. We have
considered the non-perturbative contributions up to dimension $6$ and order
$\alpha_{s}$. We have shown that the massless quark loop is the first term in
this OPE, thus putting a long-standing assumption on firm theoretical ground.
We have also shown that this leading term dominates over the non-perturbative
corrections, both those suppressed by quark masses and those that are not. The
only piece remaining to finally deduce how good a representation the massless
quark loop is of the SD HLbL tensor is the order $\alpha_{s}$ corrected
massless quark loop. Our preliminary findings show that also these corrections
are small, which will shortly be presented in an upcoming publication.

## Acknowledgements

N. H.–T. and L. L. are funded by the Albert Einstein Center for Fundamental
Physics at Universität Bern and the Swiss National Science Foundation
respectively. J. B. and A. R.–S. are supported by the Swedish Research Council
grants contract numbers 2016-05996 and 2019-03779. A. R.–S. is partially
supported by the Agence Nationale de la Recherche (ANR) under grant ANR-19-
CE31-0012 (project MORA).

## References

  * [1] G. W. Bennett, et al., Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL, Phys. Rev. D73 (2006) 072003.  arXiv:hep-ex/0602035, doi:10.1103/PhysRevD.73.072003.
  * [2] M. Tanabashi, et al., Review of Particle Physics, Phys. Rev. D98 (3) (2018) 030001\.  doi:10.1103/PhysRevD.98.030001.
  * [3] The anomalous magnetic moment of the muon in the standard model, Physics ReportsarXiv:2006.04822, doi:https://doi.org/10.1016/j.physrep.2020.07.006.
  * [4] J. Bijnens, N. Hermansson-Truedsson, A. Rodríguez-Sánchez, Short-distance constraints for the HLbL contribution to the muon anomalous magnetic moment, Phys. Lett. B 798 (2019) 134994.  arXiv:1908.03331, doi:10.1016/j.physletb.2019.134994.
  * [5] J. Bijnens, N. Hermansson-Truedsson, L. Laub, A. Rodríguez-Sánchez, Short-distance HLbL contributions to the muon anomalous magnetic moment beyond perturbation theory, JHEP 10 (2020) 203.  arXiv:2008.13487, doi:10.1007/JHEP10(2020)203.
  * [6] J. Aldins, T. Kinoshita, S. J. Brodsky, A. J. Dufner, Photon - Photon Scattering Contribution to the Sixth Order Magnetic Moments of the Muon and Electron, Phys. Rev. D1 (1970) 2378.  doi:10.1103/PhysRevD.1.2378.
  * [7] G. Colangelo, M. Hoferichter, M. Procura, P. Stoffer, Dispersion relation for hadronic light-by-light scattering: theoretical foundations, JHEP 09 (2015) 074\.  arXiv:1506.01386, doi:10.1007/JHEP09(2015)074.
  * [8] G. Colangelo, M. Hoferichter, M. Procura, P. Stoffer, Dispersion relation for hadronic light-by-light scattering: two-pion contributions, JHEP 04 (2017) 161\.  arXiv:1702.07347, doi:10.1007/JHEP04(2017)161.
  * [9] I. I. Balitsky, A. V. Yung, Proton and Neutron Magnetic Moments from QCD Sum Rules, Phys. Lett. 129B (1983) 328–334.  doi:10.1016/0370-2693(83)90676-7.
  * [10] B. L. Ioffe, A. V. Smilga, Nucleon Magnetic Moments and Magnetic Properties of Vacuum in QCD, Nucl. Phys. B232 (1984) 109–142.  doi:10.1016/0550-3213(84)90364-X.
  * [11] A. Czarnecki, W. J. Marciano, A. Vainshtein, Refinements in electroweak contributions to the muon anomalous magnetic moment, Phys. Rev. D67 (2003) 073006, [Erratum: Phys. Rev.D73,119901(2006)].  arXiv:hep-ph/0212229, doi:10.1103/PhysRevD.67.073006,10.1103/PhysRevD.73.119901.
  * [12] S. Aoki, et al., FLAG Review 2019, Eur. Phys. J. C80 (2) (2020) 113.  arXiv:1902.08191, doi:10.1140/epjc/s10052-019-7354-7.
  * [13] M. A. Shifman, A. I. Vainshtein, V. I. Zakharov, QCD and Resonance Physics. Theoretical Foundations, Nucl. Phys. B147 (1979) 385–447.  doi:10.1016/0550-3213(79)90022-1.
  * [14] G. S. Bali, G. Endrődi, S. Piemonte, Magnetic susceptibility of QCD matter and its decomposition from the lattice, JHEP 07 (2020) 183.  arXiv:2004.08778, doi:10.1007/JHEP07(2020)183.
  * [15] V. Belyaev, B. Ioffe, Determination of Baryon and Baryonic Resonance Masses from QCD Sum Rules. 1. Nonstrange Baryons, Sov. Phys. JETP 56 (1982) 493–501.