© The Authors 2024

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¹. 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.

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 β is the mutual inclination, σ_β is a parameter that regulates the shape of the distribution and Ν_β normalises the distribution over the support considered, i.e. from β = 0 to π/2.

The initial distribution of eccentricity is described by a normalised power law distribution: $$P_{\text{e}}(e)=\frac{1}{N_{\text{e}}}{e}^{\alpha },$$(2)

where e is the orbital eccentricity, α 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) ∝ e) that can be produced by repeated Ν -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_{\text{e}}(e)P_{\beta }(\beta ).$$(3)

For an orbit co-rotating with the central system, the critical angle for polar alignment is² (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}},$$(4)

where e is the orbital eccentricity and Ω is the disc longitude of the ascending node.

Supposing that all discs above β_crit go polar and neglecting subsequent external interactions, we can integrate over the initial mutual inclination above β_crit, to obtain the distribution of polar discs as a function of system eccentricity P_pol(e). At this stage, we consider Ω = π/2 and we will fix P_pol(e) to take into account a distribution of longitude of the ascending node in Sect. 2.2. Thus, $$P_{\text{pol}}(e)={\displaystyle {\int }_{{\beta }_{\text{crit}}\left(e,\frac{\pi }{2}\right)}^{\pi /2}P}(e,\beta )\text{d}\beta .$$(5)

Integrating P_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_{\text{p}}={\displaystyle {\int }_0^1{P}_{\text{pol}}}(e)\text{d}e.$$(6)

Additionally, from P_pol we can compute the mean eccentricity of stellar systems hosting polar discs, i.e.: $$\langle e\rangle ={\displaystyle {\int }_0^{1}e}{P}_{\text{pol}}(e)\text{d}e.$$(7)

Integrals in Eqs. (6) and (7) are challenging to solve due to the dependencies of β_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.

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 Ω interval in which the disc will go perpendicular (coplanar) to the orbital plane.

2.2 Taking into account the longitude of ascending node

Figure 1 shows how the critical angle β_crit in Eq. (4) depends on the longitude of the ascending node (Ω) for different orbital eccentricities. Let us take the subpopulation of circumbinary discs with a given mutual misalignment β_sp orbiting pure binaries with a given orbital eccentricity e_sp. We have a so that the critical angle for polar alignment β_crit(e_sp, Ω_crit) = β_sp. All discs with Ω < Ω_crit (dotted line in Fig. 1) will not polar align. All discs with a Ω > Ω_crit (dashed line) will polar align.

If we suppose a uniformly distributed Ω, the fraction of discs that will polar align is the ratio between the length of the dashed curve and the width of the Ω interval (i.e. π/2).Thus, the fraction f of discs that will polarly align is: $$f\left(e_{\text{sp}},{\beta }_{\text{sp}}\right)=1-\frac{2}{\pi }{\Omega }_{\text{crit}}\left(e_{\text{sp}},{\beta }_{\text{sp}}\right).$$(8)

In the case of a binary population, we have to take into account the factor f(e, β). To compute the distribution of polar discs with respect to orbital eccentricity, Eq. (5) becomes: $$P_{\text{pol}}(e)={\displaystyle {\int }_{{\beta }_{\text{crit }}(e)}^{\pi /2}P}(e,\beta )f(e,\beta )\text{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 Ω 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 Ω, if the disc inclination is higher than the minimum critical angle (the one for Ω = 90°) the disc will polar-align. Therefore, for triples β_crit(e, Ω) = β_crit(e, Ω = 90°). Thus, Eq. (5) can be written as $$\begin{array}{l}{P}_{\text{pol}}(e)={\displaystyle {\int }_{{\beta }_{\text{crit}}(e)}^{\pi /2}P}(e,\beta )\text{d}\beta \\ =-P_{\text{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{array}$$(10)

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

Fig. 2 Fraction of polar discs (Fp, leftmost column) and mean orbital eccentricity of systems hosting polar discs (〈e〉, central column) for binary system population (top row) and triple system population (bottom row). As shown by the rightmost column, each point of the α−σβ parameter space is uniquely characterized by a pair of Fp and 〈e〉 (crossing of two contour lines, red contour for polar fraction, black contour for mean eccentricity).

2.3 Polar fraction and mean eccentricity for multiple systems

Given a pair of values α and σ_β – 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 α−σ_β parameter space for the initial distributions. The parameter α ranges between −1 and 2, while σ_β ranges between 0.01 and +∞. The left panels in Fig. 2 show how the fraction of expected polar discs in an evolved population F_p depends on α and σ_β 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 β_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 α or σ_β.

The central panels in Fig. 2 show the mean eccentricity of orbits hosting polar discs (〈e〉) over the α−σ_β parameter space for binaries and hierarchical systems. The binary polar population tends to have higher 〈e〉, again due to the f factor (Eq. (8)). Indeed, the Ω interval allowing polar alignment is larger for more eccentric systems. Thus, we have more polar discs around highly eccentric systems. Lowering σ_β raises 〈e〉. Indeed, higher eccentricity values are required in regions where systems are mildly misaligned. Conversely, lowering α means lowering 〈e〉. The less likely is to have very eccentric systems in the initial population, the less eccentric the polar population will be.

Each pair of α−σ_β uniquely connects to a pair of F_p−〈e〉. 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 〈e〉 contour lines cross only at one point of the α−σ_β parameter space. Thus, we are able to numerically invert the α−σ_β and F_p−〈e〉 coordinate systems to obtain plots in Fig. 3. This plot showcases how α and depend on F_p and 〈e〉. Doing so, once we have constrained F_p and 〈e〉 from observations in a population we can constrain α and σ_β in the initial condition.

Fig. 3 α (left) and σβ (right) as a function of 〈e〉 and Fp for binary (top) and triple (bottom) populations. White regions are Fp−〈e〉 pairs that cannot be generated by any α−σβ pair. Purple curves are the parameter space region in which α = 2, 1, 0 for solid, dotted and dashed, respectively. Green curves are the parameter space region in which σβ = ∞. Black star is the solution for randomly distributed mutual inclination and eccentricity. White box point with error bars is the $F_{\text{p}}^{\text{obs}}$ and 〈e〉obs observation values derived in Sect. 2.4. Red circle points are Fp−〈e〉 pairs measured in the four Bate (2019) molecular cloud collapse simulations with different metallicities.

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³ (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 Ω 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 polar debris disc 99 Her. We obtain an average eccentricity of 〈e〉^obs = 0.67 ± 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_{\text{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_{\text{p}}^{\text{obs}}$ and 〈e〉^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 α and σ_β, respectively, as a function of 〈e〉 and F_p.

White regions are regions where no α and pair results in the given pair of 〈e〉 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 σ_β approaching infinity and α = 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_e = 1 and P_β = 2/π. The results are $F_{\text{p}}^{\mathrm{tri}}=0.54$ and 〈e〉^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_{\text{p}}^{\mathrm{tri}}=0.40$ and 〈e〉^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 σ_β approaching infinity and α = 0, respectively. In surveys, we find a polar fraction $F_{\text{p}}^{\text{obs}}=0.17\pm 0.08$ and a mean eccentricity 〈e〉^obs = 0.67 ± 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_β = sin β 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 Ω, thus applicable to the binary population. Integrating over the eccentricity we obtain an expected $F_{\text{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_β = sin β obtaining $F_{\text{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 α = 0). Over this restricted parameter space region, 〈e〉 has a lower limit given by the σ_β = ∞ case (i.e. 0.65 and 0.63 for binaries and triples, respectively), while F_p can span values from 0 to the σ_β = ∞ case (i.e. 0.40 and 0.54 for binaries and triples, respectively). The restricted parameter space region is compatible with 〈e〉^obs and $F_{\text{p}}^{\text{obs}}$. In particular, observations constrain σ_β to range between 0.25 and 0.78 for binaries and between 0.17 and 0.43 for triples which implies a narrow β distribution around β = 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 α and σ_β to explore the whole parameter space to fully exploit the information contained in the measurement of 〈e〉 and F_p. The observed mean eccentricity and polar fraction select a region of possible values of α and σ_β highlighted by the error bars in Fig. 3. For binaries, we obtain α ≤ 0.6 and 0.26 ≤ σ_β ≤ 1.8. As for the triples, α ≤ 0.46 and 0.23 ≤ σ_β ≤ 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 (α < 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 (α = 0) for orbits with semi-major axis of the order of 100 AU (Raghavan et al. 2010), up to an (increasing) thermal distribution (α = 1) for 500 AU and becomes even steeper for larger systems (Duquennoy & Mayor 1991; Tokovinin 2020; Hwang et al. 2022). Conversely, we find an α 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 α 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 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 α 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 α and in the newly born protostellar systems population by fitting the synthetic P_pol (Eq. (5)) resulting from the simulation. We computed the expected F_p and 〈e〉 with Eqs. (6) and (7), respectively. Finally, we compared the values of F_p 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_p 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 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 Sect. 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_p 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 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 (σ_β) and orbital eccentricity (α). The only required measurements are the fraction of polar discs in a disc population (F_p) 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 randomly distributed initial distribution of mutual misalignment – there must be a preference for aligned systems;

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 ~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 F_p−〈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 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). 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/18871/35993.

References

All Figures

	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 Ω interval in which the disc will go perpendicular (coplanar) to the orbital plane.
In the text

Fig. 2 Fraction of polar discs (Fp, leftmost column) and mean orbital eccentricity of systems hosting polar discs (〈e〉, central column) for binary system population (top row) and triple system population (bottom row). As shown by the rightmost column, each point of the α−σβ parameter space is uniquely characterized by a pair of Fp and 〈e〉 (crossing of two contour lines, red contour for polar fraction, black contour for mean eccentricity).

In the text

Fig. 3 α (left) and σβ (right) as a function of 〈e〉 and Fp for binary (top) and triple (bottom) populations. White regions are Fp−〈e〉 pairs that cannot be generated by any α−σβ pair. Purple curves are the parameter space region in which α = 2, 1, 0 for solid, dotted and dashed, respectively. Green curves are the parameter space region in which σβ = ∞. Black star is the solution for randomly distributed mutual inclination and eccentricity. White box point with error bars is the $F_{\text{p}}^{\text{obs}}$ and 〈e〉obs observation values derived in Sect. 2.4. Red circle points are Fp−〈e〉 pairs measured in the four Bate (2019) molecular cloud collapse simulations with different metallicities.

In the text