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

Simone Ceppi, Nicolás Cuello, Giuseppe Lodato, Cristiano Longarini, Daniel

J. Price, Daniel Elsender, Matthew R. Bate

► To cite this version:

Simone Ceppi, Nicolás Cuello, Giuseppe Lodato, Cristiano Longarini, Daniel J. Price, et al. Prob-
ing initial distributions of orbital eccentricity and disc misalignment via polar discs. Astronomy &
Astrophysics - A&A, 2024, 682, 10.1051/0004-6361/202348375. insu-04836904

HAL Id: insu-04836904

https://insu.hal.science/insu-04836904v1

Submitted on 13 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Probing initial distributions of orbital eccentricity and disc misalignment via polar discs

Simone Ceppi<fnref n="1" />,<fnref n="3" />✉, Nicolás Cuello<fnref n="2" />✉, Giuseppe Lodato<fnref n="1" />, Cristiano Longarini<fnref n="1" />,<fnref n="3" />, Daniel J. Price<fnref n="3" />✉,
Daniel Elsender<fnref n="4" />, and Matthew R. Bate<fnref n="4" />✉

Received 24 October 2023 / Accepted 15 December 2023

## ABSTRACT

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 – hydrodynamics – methods: analytical – methods: statistical – protoplanetary disks – binaries: close

## 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 suppressed<fnref n="1" />. 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. 2023). 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. 2023).

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 Sect. 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 Sect. 3 we constrain the initial conditions for both populations. In Sect. 4, we discuss our results and we give our conclusions in Sect. 5.

<footnote n="1">¹ Dipartimento di Fisica, Università degli Studi di Milano, via Celoria 16, 20133 Milano, Italy
e-mail: simone.ceppi@unimi.it</footnote>

<footnote n="2">² Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="3">³ School of Physics and Astronomy, Monash University, Clayton, VIC 3800, Australia</footnote>

<footnote n="4">⁴ Department of Physics and Astronomy, University of Exeter, Stocker Road, Exeter EX4 4QL, UK</footnote>

<footnote n="1">¹ Polar 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).</footnote>
## 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), \qquad (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_e(e) = \frac{1}{N_e} e^{\alpha}, \qquad (2)$$

where $e$ is the orbital eccentricity, $\alpha$ is a parameter regulating the distribution shape and $N_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_e(e)P_{\beta}(\beta). \qquad (3)$$

For an orbit co-rotating with the central system, the critical angle for polar alignment is<fnref n="2" /> (Farago & Laskar 2010; Zanazzi & Lai 2018; Cuello & Giuppone 2019)

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

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

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

<figure><img src="image_2.png" /><figcaption>Fig. 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°) 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.</figcaption></figure>

a distribution of longitude of the ascending node in Sect. 2.2. Thus,

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

Integrating $P_{\text{pol}}$ over the eccentricity we obtain the expected fraction of polar discs ($F_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_p = \int_0^1 P_{\text{pol}}(e)de. \qquad (6)$$

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

$$\langle e \rangle = \int_0^1 e P_{\text{pol}}(e)de. \qquad (7)$$

Integrals in Eqs. (6) and (7) are challenging to solve due to the dependencies of $\beta_{\text{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 shows how the critical angle $\beta_{\text{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_{\text{sp}}$ orbiting pure binaries with a given orbital eccentricity $e_{\text{sp}}$. We have a $\Omega_{\text{crit}}$ so that the critical angle for polar alignment $\beta_{\text{crit}}(e_{\text{sp}}, \Omega_{\text{crit}}) = \beta_{\text{sp}}$. All discs with $\Omega < \Omega_{\text{crit}}$ (dotted line in Fig. 1) will not polar align. All discs with a $\Omega > \Omega_{\text{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_{\text{sp}}, \beta_{\text{sp}}) = 1 - \frac{2}{\pi} \Omega_{\text{crit}}(e_{\text{sp}}, \beta_{\text{sp}}). \qquad (8)$$

<footnote n="2">² for a counterrotating orbit $\beta_{\text{crit}}(e_b, \Omega) = \pi - \arcsin \sqrt{\frac{1-e_b^2}{1-5e_b^2 \cos \Omega^2 + 4e_b^2}}$.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Fraction of polar discs ($F_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_p$ and $\langle e \rangle$ (crossing of two contour lines, red contour for polar fraction, black contour for mean eccentricity).</figcaption></figure>

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, Eq. (5) becomes:

$$P_{\text{pol}}(e) = \int_{\beta_{\text{crit}}(e)}^{\pi/2} P(e, \beta) f(e, \beta) d\beta. \quad (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_{\text{crit}}(e, \Omega) = \beta_{\text{crit}}(e, \Omega = 90^\circ)$. Thus, Eq. (5) can be written as

$$ \begin{aligned} P_{\text{pol}}(e) &= \int_{\beta_{\text{crit}}(e)}^{\pi/2} P(e, \beta) d\beta \\ &= -P_e(e) \frac{\sigma_\beta}{N_\beta} \left[ \exp\left(-\frac{\pi/2}{\sigma_\beta}\right) - \exp\left(-\frac{\beta_{\text{crit}}(e)}{\sigma_\beta}\right) \right]. \end{aligned} \quad (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

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 Eqs. (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 Fig. 2 show how the fraction of expected polar discs in an evolved population $F_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_{\text{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_p$ increases for higher $\alpha$ or $\sigma_\beta$.

The central panels in Fig. 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_p-\langle e \rangle$. This can be seen with contour lines in the right panel of Fig. 2. A given pair of polar fraction $F_p$ and mean eccentricity $\langle e \rangle$ contour lines cross only at one point of the $\alpha-\sigma_\beta$ parameter space. Thus,
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** $\alpha$ (left) and $\sigma_{\beta}$ (right) as a function of $\langle e \rangle$ and $F_{p}$ for binary (top) and triple (bottom) populations. White regions are $F_{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_{p}^{\text{obs}}$ and $\langle e \rangle^{\text{obs}}$ observation values derived in Sect. 2.4. Red circle points are $F_{p}-\langle e \rangle$ pairs measured in the four Bate (2019) molecular cloud collapse simulations with different metallicities.</figcaption></figure>

we are able to numerically invert the $\alpha-\sigma_{\beta}$ and $F_{p}-\langle e \rangle$ coordinate systems to obtain plots in Fig. 3. This plot showcases how $\alpha$ and $\sigma_{\beta}$ depend on $F_{p}$ and $\langle e \rangle$. Doing so, once we have constrained $F_{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 polar<fnref n="3" /> (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 (HD 98800B), systems that are likely polar (SR24 N and HD 142527) and the

<footnote n="3">³ There is a degeneracy of 180° in the orbital longitude of the ascending node (Czekala et al. 2019).</footnote>
polar debris disc 99 Her. We obtain an average eccentricity of $\langle e \rangle^{\text{obs}} = 0.67 \pm 0.11$.

To evaluate the polar disc fraction, if we take into account the single confirmed polar system (HD 98800B) we end up with a polar fraction of 0.08. If, however, we include also the additional two likely polar systems (SR24 N and HD 142527) 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_p^{\text{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_p^{\text{obs}}$ and $\langle e \rangle^{\text{obs}}$ as references for exploring the parameter space of the initial conditions.

# 3. Mapping observations onto the parameter space

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

White regions are regions where no $\alpha$ and $\sigma_{\beta}$ pair results in the given pair of $\langle e \rangle$ and $F_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_c = 1$ and $P_{\beta} = 2/\pi$. The results are $F_p^{\text{tri}} = 0.54$ and $\langle e \rangle^{\text{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_p^{\text{bin}} = 0.40$ and $\langle e \rangle^{\text{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_p^{\text{obs}} = 0.17 \pm 0.08$ and a mean eccentricity $\langle e \rangle^{\text{obs}} = 0.67 \pm 0.11$ (see Sect. 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_p^{\text{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_p^{\text{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 Figs. 3b and 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_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^{\text{obs}}$ and $F_p^{\text{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_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 \le 0.6$ and $0.26 \le \sigma_{\beta} \le 1.8$. As for the triples, $\alpha \le 0.46$ and $0.23 \le \sigma_{\beta} \le 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 simula-
tions 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 α is closer to our findings (systems of ~200 AU have
α ~ 0.6). In addition, gravitational interactions tend to raise the
average eccentricity of the population. However, the thermali-
sation of the distribution could be only partially explained with
gravitational interaction during cluster evolution which should
act for much longer in order to push α from 0 to 1 (Heggie
1975; Weinberg et al. 1987). Thus, depending on the thermali-
sation timescale, we could reconcile our findings with observed
distributions.

In addition, the orbital eccentricity can evolve through inter-
actions 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 ini-
tial 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 simula-
tions of the collapse of molecular clouds. In such simulations,
the initial conditions of the cloud are set (e.g. amount of turbu-
lence, 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 α and σβ
in the newly born protostellar systems population by fitting the
synthetic Ppol (Eq. (5)) resulting from the simulation. We com-
puted the expected Fp and <e> with Eqs. (6) and (7), respectively.
Finally, we compared the values of Fp and <e> obtained from
the simulations with the ones measured in observations (see
Sect. 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 Fig. 3 represent the computed F<sub>p</sub> and ⟨e⟩ 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 Fig. 3. However, it appears to have too little impact on the properties of the population.

Regardless of the metallicity, the resulting polar disc frac-
tions all are generally too low compared to observations. Given
that the distribution of misalignment angles from the calcu-
lations 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 Sect. 2.4, the proposed observational value has to
be considered as an upper limit. In other words, we cannot cur-
rently completely rule out such small polar disc fractions from
observations.

Finally, this analysis would benefit from numerical simula-
tions producing distributions of mutual inclination and eccen-
tricity, 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 pre-
dictions and better constraints for the polar disc population.
Provided accurate measurements of F<sub>p</sub> and ⟨e⟩, 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 incli-
nation 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 (σβ) and orbital eccentricity (α).
The only required measurements are the fraction of polar discs
in a disc population (Fp) and the mean eccentricity of systems
hosting polar discs (⟨e⟩). Despite the low statistics available, we
find that:

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

2. The orbital eccentricity is marginally compatible with a ran-
dom 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 ~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 α-σβ and Fp-⟨e⟩ 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 distri-
butions 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). D.J.P. acknowledges Australian Research Council funding via DP220103767. S.C.'s visit to Monash was funded by the EU Marie Curie RISE scheme DUST-BUSTERS (grant 823823). The ideas discussed in this paper were partially formulated by S.C., D.J.P., D.E., and M.R.B. 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, H., Dehnen, W., Nixon, C., & King, A. 2015, MNRAS, 449, 65

Ambartsumian, V. A. 1937, *AZh*, **14**, 207

Balmer, W. O., Follette, K. B., Close, L. M., et al. 2022, AJ, 164, 29

Bate, M. R. 2018, MNRAS, 475, 5618

Bate, M. R. 2019, MNRAS, 484, 2341

Bate, M. R., Bonnell, I. A., Clarke, C. J., et al. 2000, MNRAS, 317, 773

Bate, M. R., Bonnell, I. A., & Bromm, V. 2002, MNRAS, 336, 705

Bate, M. R., Lodato, G., & Pringle, J. E. 2010, MNRAS, 401, 1505

Ceppi, S., Cuello, N., Lodato, G., et al. 2022, MNRAS, 514, 906

Ceppi, S., Longarini, C., Lodato, G., Cuello, N., & Lubow, S. H. 2023, MNRAS, 520, 5817

Cuello, N., & Giuppone, C. A. 2019, A&A, 628, A119

Czekala, I., Chiang, E., Andrews, S. M., et al. 2019, ApJ, 883, 22

D'Orazio, D. J., & Duffell, P. C. 2021, ApJ, 914, L21

Duchêne, G., & Kraus, A. 2013, ARA&A, 51, 269

Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485

El-Badry, K., Rix, H.-W., Tian, H., Duchêne, G., & Moe, M. 2019, MNRAS, 489, 5822

Elsender, D., & Bate, M. R. 2021, MNRAS, 508, 5279

Elsender, D., Bate, M. R., Lakeland, B. S., Jensen, E. L. N., & Lubow, S. H. 2023, MNRAS, 523, 4353

Farago, F., & Laskar, J. 2010, MNRAS, 401, 1189

Fernández-López, M., Zapata, L. A., & Gabbasov, R. 2017, ApJ, 845, 10

Heggie, D. C. 1975, MNRAS, 173, 729

Hwang, H.-C., Ting, Y.-S., & Zakamska, N. L. 2022, MNRAS, 512, 3383

Jeans, J. H. 1919, MNRAS, 79, 408

Kennedy, G. M., Matrà, L., Facchini, S., et al. 2019, Nat. Astron., 3, 230

Krumholz, M. R., Klein, R. I., & McKee, C. F. 2012, ApJ, 754, 71

Larson, R. B. 1972, MNRAS, 156, 437

Lebreuilly, U., Hennebelle, P., Colman, T., et al. 2021, ApJ, 917, L10

Lebreuilly, U., Vallucci-Goy, V., Guillet, V., Lombart, M., & Marchand, P. 2023, MNRAS, 518, 3326

Lepp, S., Martin, R. G., & Lubow, S. H. 2023, ApJ, 943, L4

Lubow, S. H., & Ogilvie, G. I. 2000, ApJ, 538, 326

Lubow, S. H., & Martin, R. G. 2018, MNRAS, 473, 3733

Martin, R. G., & Lubow, S. H. 2017, ApJ, 835, L28

Martin, R. G., & Lubow, S. H. 2019, MNRAS, 490, 1332

Martin, R. G., Nixon, C., Lubow, S. H., et al. 2014, ApJ, 792, L33

Mathew, S. S., & Federrath, C. 2021, MNRAS, 507, 2448

Mathew, S. S., Federrath, C., & Seta, A. 2023, MNRAS, 518, 5190

Moe, M., & Di Stefano, R. 2017, ApJS, 230, 15

Offner, S. S. R., Moe, M., Kratter, K. M., et al. 2023, ASP Conf. Ser., 534, 275

Price, D. J., & Bate, M. R. 2007, Ap&SS, 311, 75

Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1

Smallwood, J. L., Franchini, A., Chen, C., et al. 2020, MNRAS, 494, 487

Smith, K. W., Bonnell, I. A., & Bate, M. R. 1997, MNRAS, 288, 1041

Tokovinin, A. 2020, MNRAS, 496, 987

Walch, S., Whitworth, A. P., & Girichidis, P. 2012, MNRAS, 419, 760

Weinberg, M. D., Shapiro, S. L., & Wasserman, I. 1987, ApJ, 312, 367

Wurster, J., Bate, M. R., & Price, D. J. 2019, MNRAS, 489, 1719

Zanazzi, J. J., & Lai, D. 2018, MNRAS, 473, 603

Zhao, B., Tomida, K., Hennebelle, P., et al. 2020, Space Sci. Rev., 216, 43

Zúñiga-Fernández, S., Olofsson, J., Bayo, A., et al. 2021, A&A, 655, A15

Zurlo, A., Gratton, R., Pérez, S., & Cieza, L. 2023, Eur. Phys. J. Plus, 138, 411