# Precession and polar alignment of accretion discs in triple (or multiple) stellar systems

Simone Ceppi, Cristiano Longarini, Giuseppe Lodato, Nicolás Cuello, Stephen Lubow

► To cite this version:

Simone Ceppi, Cristiano Longarini, Giuseppe Lodato, Nicolás Cuello, Stephen Lubow. Precession and polar alignment of accretion discs in triple (or multiple) stellar systems. Monthly Notices of the Royal Astronomical Society, 2023, 520 (4), pp.5817-5827. 10.1093/mnras/stad444. hal-04966028

HAL Id: hal-04966028

https://hal.science/hal-04966028v1

Submitted on 20 Mar 2025

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
# Precession and polar alignment of accretion discs in triple (or multiple) stellar systems

Simone Ceppi<sup>1,★</sup> Cristiano Longarini<sup>1</sup>,<sup>1</sup> Giuseppe Lodato<sup>1</sup>,<sup>1</sup> Nicolás Cuello<sup>2</sup> and Stephen H. Lubow<sup>1,3</sup>

Accepted 2023 February 7. Received 2023 February 6; in original form 2022 November 22

## ABSTRACT

We investigate the mechanism of polar alignment for accretion discs in hierarchical systems (HSs) with more than two stars. In eccentric binary systems, low-mass discs that are sufficiently tilted to the binary orbit align in a polar configuration with respect to the binary plane by aligning their angular momentum to the binary eccentricity vector. In HSs, secular evolution of the orbital parameters makes the eccentricity vector of the system precess with time. This precession undermines the stability of the polar orbit for accretion discs hosted in HSs. We analytically show that the binary criteria for polar alignment derived in the literature are necessary but not sufficient conditions for polar alignment in HSs. Then, we derive an analytical criterion for polar alignment in HSs. In general, we find that discs orbiting the innermost level of an HS can go polar. Conversely, radially extended discs orbiting the outer levels of an HS cannot polarly align and evolve as orbiting around a circular binary. We confirm our findings through detailed numerical simulations. Also, our results are compatible with the observed distribution of disc-orbit mutual inclination. Finally, we compare the observed distribution of disc inclinations in the binary and in the HS populations. Binaries host mainly coplanar discs, while HSs show a wide range of disc inclinations. We suggest that the wider range of inclinations in HSs results from the secular oscillation of their orbital parameters (such as Kozai–Lidov oscillations), rather than from a different initial condition or evolution between HSs and binaries.

**Key words:** accretion, accretion discs – hydrodynamics – planets and satellites: formation – (stars:) binaries (including multiple): close.

## 1 INTRODUCTION

Star formation takes place in clustered environments, where pre-stellar objects dynamically interact with each other, forming the seeds of multiple stellar systems (Clarke, Bonnell & Hillenbrand 2000; Reipurth et al. 2014; Offner et al. 2022). Surveys of star-forming environments (McKee & Ostriker 2007), as well as numerical simulations (Bate 2009; Bate, Lodato & Pringle 2010; Bate 2018), show that stars form from a sequence of accretion episodes, for which the angular momentum is randomly oriented. In such a scenario, at least some accretion discs are expected to be misaligned to the stellar orbital plane. Recently, observations have confirmed these theoretical predictions (Czekala et al. 2019). For example, the systems GG Tau (Köhler 2011; Aly, Lodato & Cazzoletti 2018; Keppler et al. 2020a), KH 15D (Chiang & Murray-Clay 2004; Lodato & Facchini 2013), HD142527 (Casassus et al. 2013; Price et al. 2018b), and GW Ori (Kraus et al. 2020; Bi et al. 2020; Smallwood et al. 2021) exhibit high relative misalignments.

The dynamics of misaligned circumbinary discs dynamics has been actively studied during the last few decades, both theoretically (Papaloizou & Pringle 1983; Papaloizou & Lin 1995; Ogilvie 1999;

Bate et al. 2000; Nixon et al. 2012; Nixon, King & Price 2013) and numerically (Facchini, Lodato & Price 2013; Aly & Lodato 2020). An inclined gas disc around a circular binary tends to precess around the binary angular momentum vector. If the disc sound crossing time is short enough (i.e. the inner part of the disc communicates its precession to the disc outer edge efficiently, Papaloizou & Lin 1995) the disc rigidly precesses with a frequency given by Bate et al. (2000) (see also Lodato & Facchini 2013). In addition, because of viscosity, the disc tends to align to the binary orbital plane, with an alignment time-scale that is of the order of the viscous time (Papaloizou & Terquem 1995; Larwood et al. 1996; Bate et al. 2000; Lubow & Ogilvie 2000) eventually reaching the coplanar equilibrium configuration. Also, misalignment in circumbinary discs affects dust dynamics (Longarini et al. 2021; Aly et al. 2021), creating substructures such as rings.

When the binary system is eccentric an additional stable configuration is possible. Aly et al. (2015) and Martin & Lubow (2017) showed that a circumbinary disc around an eccentric binary tends to align its angular momentum vector to the eccentricity vector of the stellar system, leading to a polar configuration. This mechanism happens if the initial inclination is above a critical value, that is a decreasing function of the eccentricity. Polar alignment is likely to occur for high initial inclination, high eccentricity (Martin & Lubow 2018), and cold and low-mass discs (Martin & Lubow 2019).

<footnote n="1"><sup>1</sup>Dipartimento di Fisica, Università degli Studi di Milano, via Celoria 16, I-20133 Milano, Italy</footnote>

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

<footnote n="3"><sup>3</sup>Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote>* E-mail: simone.ceppi@unimi.it</footnote>
Although we would expect to see several accretion discs in a polar configuration, only a few have been observed up to now. One is the disc orbiting HD98800 B (Kennedy et al. 2019). Also, 99 Her is a debris disk in a polar configuration (Kennedy et al. 2012) that likely evolved from a gaseous accretion disc (Smallwood et al. 2020). In addition, recently Kenworthy et al. (2022) showed that the light curve of an eclipse in V773 Tau multiple stellar system can be explained with the presence of a highly inclined disc. Interestingly, HD98800 B (as well as V773) is not a pure binary system (for which the polar alignment theory has been developed), but it resides inside a hierarchical quadruple stellar system. A recent study by Martin et al. (2022) finds that the high multiplicity of HD98800 does not prevent polar alignment of the circumbinary disc.

As in the case of HD98800, young multiple stellar systems often comprise more than two stars arranged in hierarchical configurations (Duchêne & Kraus 2013; Moe & Di Stefano 2017). A hierarchical configuration is made of nested binary orbits in which each semi-major axis is much larger than the semi-major axes of its subsystems. This semi-major axes hierarchy allows an analytic perturbative analysis of the orbital evolution that results in a periodic exchange of angular momentum between the orbits. As a result, hierarchical systems show periodic oscillations in the orbital shape and orientation (Naoz et al. 2013). Given this secular orbital evolution, the possibility of polar discs in systems with more than two stars needs to be further investigated.

In this work, we study whether discs in hierarchical systems can align polarly or if they fail to catch up with the secular evolution of the orbit. In Section 2, we present the analytical framework under which we derive an additional criterion for the possibility of polar alignment of circumbinary discs in hierarchical systems. In Section 3, we numerically test the analytical findings of Section 2. In Section 4, we discuss the results of sections 2 and 3. We conclude in Section 5.

# 2 POLAR CONFIGURATION STABILITY IN SYSTEMS WITH MORE THAN TWO STARS

## 2.1 Disc polar alignment in systems with more than two stars

We describe the orientation of a generic misaligned orbit around a central binary with the tilt angle, $\beta$, and the longitude of the ascending node (relative to the eccentricity vector direction), $\Omega$. The $\beta$ and $\Omega$ angles describe the direction of the specific orbital angular momentum vector ($l$) relative to the central binary angular momentum ($l_b$) and eccentricity vector ($e_b$). The frame of reference is defined by the instantaneous direction of the binary angular momentum and eccentricity vector with axes $e_b$, $e_b \times l_b$, and $l_b$. We define $\beta$ and $\Omega$ via the following relation:

$$l = (\sin \Omega \sin \beta, \cos \Omega \sin \beta, \cos \beta). \qquad (1)$$

By virtue of this definition, $\beta$ is the angle between the circumbinary orbit orbital plane and the binary orbital plane, while $\Omega$ is the angle between the direction of the binary eccentricity vector and the circumbinary orbit ascending node. Using these angles, we can describe the orientation of the orbital plane of a third body around a binary, as well as the orientation of a circumbinary accretion disc. In the latter case, $\beta$ and $\Omega$ can be functions of the disc radius $R$. In the following, we suppose the accretion disc to precess as rigid body, with a tilt and longitude of the ascending node independent of $R$.

Analytical studies of the restricted three body problem in Farago & Laskar (2010) showed that an inclined test particle orbiting an eccentric binary system undergoes two alternative kinds of precession. In the first kind, the angular momentum vector of the particle precesses

around the angular momentum vector of the binary. For a circular orbit binary, this kind of particle orbit, known as a circulating orbit, maintains a constant angle between the particle orbital plane and the binary orbital plane (i.e. a constant tilt $\beta$), while the longitude of the ascending node $\Omega$ spans $2\pi$ with a typical frequency, given for example by Lodato & Facchini (2013) (see their equation 12). For an eccentric orbit binary, the tilt angle $\beta$ and nodal precession rate vary in time by an amount that depends on the binary eccentricity (Smallwood et al. 2019). In the second kind of precession, the angular momentum of the particle precesses around the eccentricity vector of the binary. In this case, which is a librating orbit, $\beta$ oscillates around 90 degrees, while $\Omega$ oscillates around a fixed value. The orbit of the test particle belongs to the first or the second kind of orbit depending on the angles $\beta$ and $\Omega$, and on the eccentricity of the central binary $e_b$. If the tilt of the particle is larger than a critical angle $\beta_{crit}(e_b, \Omega)$, the particle will precess around the binary eccentricity vector, otherwise around the binary angular momentum vector. For an orbit co-rotating with the central system, the critical angle for polar alignment is¹ (Farago & Laskar 2010; Cuello & Giuppone 2019)

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

A gaseous disc around an eccentric binary follows a similar evolution, except that, due to viscous dissipation, the orbital plane of the disc decays towards one of the two equilibrium configurations. Orbits precessing around the binary angular momentum tend to become coplanar (aligned with the binary orbital plane), while orbits precessing around the eccentricity vector decay towards a polar configuration (perpendicular to the orbital plane). Disc polar alignment in circumbinary discs has been studied in many previous works both analytically and numerically (e.g. Aly et al. 2015; Martin & Lubow 2017; Lubow & Martin 2018; Zanazzi & Lai 2018; Cuello & Giuppone 2019).

Numerical simulations of the collapse of a molecular cloud show that the stellar and accretion disc formation process is very chaotic. In particular, we expect a randomly distributed initial misalignment between the stellar systems' orbital planes and the discs orbiting them (Bate 2018). Also, from analytical investigations, we expect no correlation between the disc and the binary initial angular momentum (Toci et al. 2018). Misaligned discs can result also from stellar flybys, which can generate an inclination up to 60° (Xiang-Gruess 2016; Cuello, Ménard & Price 2023). Thus, even if the conditions for polar alignment should be easily met, up to now we observed only one polarly aligned disc, and it orbits around the inner binary of a hierarchical quadruple system (i.e. HD 98800, Kennedy et al. 2019; Zúñiga-Fernández et al. 2021). In addition, the vast majority of observed misaligned discs (particularly highly misaligned ones, that could undergo polar alignment) are found around systems with more than two stars (Czekala et al. 2019).

In general, stable stellar systems with more than two stars are arranged in nested binary orbits in which each semi-major axis is much larger than the semi-major axes of its subsystems. These configurations are called hierarchical and are the only ones that guarantee long-lasting stability to the stellar system. For example, Fig. 1 sketches a triple system in a hierarchical configuration. Each hierarchical system level can be approximated to zeroth order by a binary system. Doing so, we could study polar alignment of

$$^1\text{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}}.$$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Hierarchical triple system sketch (not to scale). Light blue refers to the outer orbit, dark blue to the inner one. Dashed curves, crosses, and arrows are the orbit, center of mass, and eccentricity vector of the hierarchical level considered.</figcaption></figure>

discs around systems with more than two stars with the formalism developed by Lubow & Martin (2018). This crude approximation, however, completely neglects the perturbative term introduced in the gravitational potential by the presence of more than two stars. Taking into account this term means to reduce the hierarchical system to a hierarchical triple system. The exchange of angular momentum between the hierarchical levels of the system triggers further dynamical mechanisms such as Kozai–Lidov oscillations (von Zeipel 1910; Kozai 1962; Lidov 1962; Naoz 2016; Hamers 2021).

In a hierarchical triple system, we can analytically compute the evolution of the binary orbital parameters of both the inner and outer orbit. Due to their mutual torques, we expect these osculating orbital elements to evolve with time on a secular time-scale. In particular, the eccentricity vector direction, along which highly inclined discs would like to align their angular momentum, varies with time.

In coplanar hierarchical triple systems, we can analytically compute the precession frequency of the outer and inner orbit eccentricity vector. At the same time, a nearly polar aligned disc has a typical tilt oscillation frequency that is its nodal libration frequency (Lubow & Martin 2018). In the triple case, the secular evolution of the orbital parameters shifts the polar orbit position with time, possibly undermining the stability of the polar configuration. In particular, if the nodal libration frequency is not high enough, the disc will not be able to track the evolving stable polar orbit, failing to remain polarly aligned.

In a hierarchical triple system, we have two possible accretion discs orbiting a multiple stellar system: the first one orbiting the inner binary of the triple and the second one orbiting the outer orbit of the triple system. In the next sections, we discuss the time-scales involved in the two cases and the stability of the polar orbit around both the outer (Section 2.2) and the inner (Section 2.3) orbit.

## 2.2 Polar alignment in the circumtriple disc

Let us consider a coplanar hierarchical triple as in the sketch of Fig. 1. The masses of the inner binary stars are $m_a$ and $m_b$, for the primary and the secondary, respectively, while $m_c$ is the mass of the external third body. In addition, we define $q_1$, $e_1$, $a_1$, and $\Omega_1$ ($q_2$, $e_2$, $a_2$, and $\Omega_2$) the mass ratio, eccentricity, semi-major axis, and Keplerian frequency of the inner (outer) orbit. To zeroth order, the outer orbit of the triple is equivalent to the orbit of its associated binary, defined as the binary composed of $m_c$ and a single star of mass $m_a + m_b$

placed in the inner binary centre of mass. Thus, instantaneously, $m_c$ orbits the centre of mass of the inner binary along the same orbit as the one the associated binary would trace.

A nearly polar disc orbiting a pure binary undergoes nodal libration with the following theoretical frequency (Lubow & Martin 2018):

$$ \omega_{\text{pa,2}} = \frac{3\sqrt{5}}{4} e_b \sqrt{1 + 4e_b^2} \frac{M_1 M_2}{(M_1 + M_2)^2} \left\langle \left( \frac{a_b}{R} \right)^{\frac{7}{2}} \right\rangle \Omega_b, \quad (3) $$

where $e_b$ is the binary orbit eccentricity, $M_1$ and $M_2$ are the binary masses, $\Omega_b$ the binary Keplerian frequency, and angular bracket notes the average of the ratio between the binary semi-major axis $a_b$ and the disc radius $R$, weighted over the angular momentum of the disc at radius $R$.

We apply equation (3) to the outer orbit case, as instantaneously it can be approximated by its associated binary. As for the masses, $M_1 = m_c$ and $M_2 = m_a + m_b$. For the eccentricity, $e_b = e_2$ and for the orbital frequency, $\Omega_b = \Omega_2$. We also compute the average weighted over the angular momentum assuming a Keplerian velocity field and that the disc extends from $R_{\text{in}}$ to $R_{\text{out}}$ with a density profile $\Sigma(R) \propto R^{-1}$. By doing so, we obtain:

$$ \omega_{\text{pa,2}} = \frac{9\sqrt{5}}{16} e_2 \sqrt{1 + 4e_2^2} q_2 (1 - q_2) \left( \frac{a_2}{R_{\text{in}}} \right)^{\frac{7}{2}} \frac{1 - x_o^{-2}}{x_o^{3/2} - 1} \Omega_2, \quad (4) $$

where $x_o = R_{\text{out}}/R_{\text{in}}$ and $q_2 = m_c/(m_a + m_b + m_c)$ is the outer orbit mass ratio.

The librating behaviour happens due to the precession of the disc angular momentum around the eccentricity vector of the triple outer orbit. However, the triple outer orbit eccentricity vector $e_2$ precesses with time due to the perturbation of the inner binary. The precession rate for the outer orbit can be analytically computed from the Hamiltonian of the stellar system with a perturbative approach (e.g. Murray & Dermott 2000; Naoz et al. 2013):

$$ \omega_{e_2} = \frac{3}{4} a^2 q_1 (1 - q_1) \frac{1 + \frac{3}{2} e_1^2}{(1 - e_2^2)^2} \Omega_2, \quad (5) $$

where $a = a_1/a_2$ is the ratio between the inner and the outer orbit semi-major axis and $q_1 = m_b/(m_a + m_b)$ is the inner orbit mass ratio. This means that a polarly aligned disc around the outer orbit of a hierarchical triple system oscillates around a stable orbit that precesses with time with a frequency $\omega_{e_2}$. Thus, in order for the disc to stay ‘attached’ to the moving eccentricity vector, the condition for libration should be unaffected by the binary precession, meaning that $\omega_{\text{pa,2}} > \omega_{e_2}$.

To check the stability of small oscillations around the precessing polar configuration, we define the dimensionless quantity $T_2 = \omega_{e_2}/\omega_{\text{pa,2}}$: when $T_2 < 1$, polar alignment takes place; conversely, if $T_2 > 1$, it does not. The value of $T_2$ is

$$ T_2 = \frac{4}{3\sqrt{5}} a^2 \left( \frac{R_{\text{in}}}{a_2} \right)^{\frac{7}{2}} \frac{x_o^{3/2} - 1}{1 - x_o^{-2}} \frac{q_1(1 - q_1)}{q_2(1 - q_2)} \mathcal{F}_2(e_1, e_2), \quad (6) $$

where

$$ \mathcal{F}_2(e_1, e_2) = \frac{1 + \frac{3}{2} e_1^2}{(1 - e_2^2)^2 e_2 \sqrt{1 + 4e_2^2}}, \quad (7) $$

The $T_2$ factor scales quadratically with the semi-major axes ratio: thus, either widening the outer orbit or shrinking the inner one decreases $T_2$ (which favours polar alignment). The torque between the inner and the outer orbit depends on $a$. Hence, the lower the value of $a$, the lower the precession frequency (i.e. the slower the precession). In addition, $T_2$ depends on the inner and outer orbit eccentricities via the $\mathcal{F}_2$ term. This term diverges both for $e_2$ approaching unity and
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Dependency of $T_2$ (see sec. 2.2) on mass ratios, semi-major axes ratio, and eccentricities. When $T_2$ is smaller than unity (dash dotted red line), polar alignment is possible, conversely for $T_2 > 1$ it is not. On the left panel, $T_2$ is plotted as a function of the inner binary mass ratio $q_1$ (fixing $q_2$ to 0.5), while on the right one, as a function of the outer binary mass ratio $q_2$ (fixing $q_1$ to 0.5). The thicker lines refer to $e_2 = 0.5$ and $e_1 = 0$. The dotted and dashed lines show how increasing $e_1$ to 0.7 (fixing $e_2 = 0.5$) and lowering $e_2$ to 0.2 (fixing $e_1 = 0$) affect $T_2$, respectively. Different colours refer to different semi-major axes ratios.</figcaption></figure>

zero. In the former case, the precession frequency of the outer orbit diverges, making polar alignment impossible. In the latter case, the nodal libration frequency of a nearly polarly aligned disc reduces to zero. Indeed, as the outer binary becomes more circular, the process of polar alignment becomes less likely. A recent study by Lepp, Martin & Lubow (2023) examined the behaviour of circum triple test particles and found a similar criteria as given by $T_2$. Their criteria differ from $T_2$ mainly at a small $e_2$.

As for the mass ratios, polar alignment is difficult for high $q_1$ and for $q_2$ approaching zero or unity. Indeed, higher values of $q_1$ (i.e. nearly equal mass inner binaries) translate into higher precession frequencies. Meanwhile, very high and very low $q_2$ reduce the outer binary to a single star concentrating the mass in the third body or in the inner binary (respectively), thus reducing the torque on the disc to zero. Conversely, when $q_1$ tends to zero, the inner binary perturbation becomes negligible – slowing the precession. This favours polar alignment.

Fig. 2 shows how $T_2$ depends on the system orbital parameters. In this plot, we assume $R_{in} = 1.5 a_2$ (as expected in very inclined discs: Lubow, Martin & Nixon 2015; Miranda & Lai 2015) and $R_{out} = 30 a_2$, thus $x_o = 30/1.5$. In the left panel, $T_2$ is plotted as a function of the inner binary mass ratio $q_1$ (fixing $q_2$ to 0.5, i.e. equal mass outer binary), while on the right one as a function of the outer binary mass ratio $q_2$ (fixing $q_1$ to 0.5, i.e. equal mass inner binary). The thicker lines refer to a configuration in which $e_2 = 0.5$. Such outer orbit eccentricity avoids extreme values that are expected to undermine (also) the stability of the triple stellar system and it sets the easiest configuration for polar alignment to occur, being $e_2 = 0.5$ the minimum of $F_2$. Moreover, for $0.2 < e_2 < 0.8$, the impact on $F_2$ is of the order of 1. We also set $e_1 = 0$ to favour polar alignment as well. The dotted and dashed lines show how increasing $e_1$ to 0.7 (fixing $e_2$ to 0.5) and lowering $e_2$ to 0.2 (fixing $e_1$ to zero) affect $T_2$, respectively. The vast majority of the parameter space does not allow for polar alignment, as the curves are always above unity. The only exception is where the outer orbit is reduced to a two-body system,

<figure><img src="image_4.png" /><figcaption>**Figure 3.** $T_2$ factor as defined in equation (6) as a function of the ratio between the outer and the inner disc radii $x_o$. We used the same parameters of the solid line in Fig. 2 with a semi-major axis ratio of 0.1. Different colours are different inner binary mass ratios. The outer binary mass ratio is fixed to 0.5, the most favourable to polar alignment.</figcaption></figure>

that is when $q_1$ tends to zero or for extremely small semi-major axis ratios $a$.

Last, the size of the disc also impacts the possibility of polar alignment. The more extended the disc, the higher $T_2$, due to the lower nodal libration frequency of a radially extended disc. Fig. 3 shows the dependence of $T_2$ on $x_o$. The two curves refer to two different inner binary mass ratios (the outer binary mass ratio is set to 0.5, the most favourable to polar alignment). Provided $x_o$ is low enough ($\lesssim 10$), the circum triple accretion disc can become polar. Depending on the semi-major axis and mass ratios, if the triple is not
<figure><img src="image_5.png" /><figcaption>Figure 4. Ratio between the inner binary precession frequency and the disc Kozai–Lidov oscillation frequency. We use the same disc parameters as in Fig. 5 ($x_{\text{out}} = 5/1.5$ and $R_{\text{in}}/a_1 = 1.5$). The ratio is always above unity (red dash-dotted line) except for inner binary eccentricity near unity. Thus, in this configuration precession is always faster than disc Kozai–Lidov oscillation.</figcaption></figure>

isolated and the disc is truncated from the outside, $x_o$ can drop under that value, favouring the polar alignment of the disc.

## 2.3 Polar alignment in the inner circumbinary disc

The stability of the polar orbit for a disc orbiting the inner binary of a hierarchical triple system is determined by the inner binary eccentricity vector precession rate. In addition, Martin et al. (2022) showed that in this configuration also Kozai–Lidov oscillation of the disc (due to the interaction with the outer third body) could prevent circum-inner binary disc polar alignment. The disc Kozai–Lidov oscillation time-scale can be computed by the theory of rigid disks (Larwood & Papalozou 1997; Lubow & Ogilvie 2001) as done in Martin et al. (2014) (see their equation (4)). Using the same technique as in equation (3), we obtain the following frequency for the circum-inner disc Kozai–Lidov oscillation:

$$ \omega_{\text{KL}} = \frac{3}{4} \frac{q_2}{1-q_2} a^3 (1-e_2^2)^{-3/2} \left(\frac{R_{\text{in}}}{a_1}\right)^{3/2} \frac{1-x_o^{-2}}{x_o^{3/2}-1} \Omega_1. \quad (8) $$

Firstly, we compare the disc Kozai–Lidov oscillation time-scale with the other relevant ones. We compute the nodal libration frequency of a nearly polar orbit around the inner binary of a hierarchical triple system ($\omega_{\text{pa},1}$) and the precession frequency of the inner binary eccentricity vector ($\omega_{e_1}$), as we did for the outer orbit case in Section 2.2. We obtain:

$$ \omega_{\text{pa},1} = \frac{9\sqrt{5}}{16} e_1 \sqrt{1 + 4e_1^2 q_1 (1 - q_1)} \left(\frac{a_1}{R_{\text{in}}}\right)^{\frac{7}{2}} \frac{1 - x_o^{-2}}{x_o^{3/2} - 1} \Omega_1, \quad (9) $$

and

$$ \omega_{e_1} = \frac{3}{4} a^3 \frac{q_2}{1-q_2} \sqrt{\frac{(1-e_2^2)}{(1-e_2^2)^3}} \Omega_1. \quad (10) $$

Fig. 4 shows the ratio between the precession frequency of the inner binary eccentricity vector ($\omega_{e_1}$) and the disc Kozai–Lidov frequency of the circum-inner binary disc ($\omega_{\text{KL}}$). This ratio solely depends on the eccentricity of the inner orbit $e_1$ and on the disc extent,

as shown by equations (8) and (10). In Fig. 4, we use $x_{\text{out}} = 5/1.5$ and $R_{\text{in}}/a_1 = 1.5$. With such parameters, apart from $e_1 \sim 1$, precession is always faster than disc Kozai–Lidov oscillation. In general, for $\omega_{\text{KL}} > \omega_{e_1}$ we would need a combination of extremely eccentric inner binaries, small discs ($x_{\text{out}} < 5/1.5$) and large disc cavities ($R_{\text{in}}/a_1 > 2$). Thus, if polar alignment occurs because of a sufficiently slow precession of the inner binary, the disc Kozai–Lidov oscillation will not prevent polar alignment, since it is generally slower than precession. Having checked that, we go back to discussing the role of the inner binary precession. We define the corresponding $T_1 = \omega_{e_1}/\omega_{\text{pa},1}$ parameter for the inner binary as

$$ T_1 = \frac{4}{3\sqrt{5}} a^3 \left(\frac{R_{\text{in}}}{a_1}\right)^{\frac{7}{2}} \frac{x_o^{3/2} - 1}{1 - x_o^{-2}} \frac{q_2/q_1}{(1 - q_1)(1 - q_2)} \mathcal{F}_1(e_1, e_2), \quad (11) $$

with

$$ \mathcal{F}_1(e_1, e_2) = \sqrt{\frac{1 - e_1^2}{(1 - e_2^2)^3} \frac{1}{e_1 \sqrt{1 + 4e_1^2}}}, \quad (12) $$

The $T_1$ factor scales with the third power of the semi-major axes ratio. Thus, either widening the outer orbit or shrinking the inner one results in an overall decrease of $T_1$. The precession frequency lowers (i.e. the precession is slower) with lower $a$ because the torque between the inner and the outer orbit weakens. In addition, $T_1$ depends on the inner and outer orbit eccentricities via the $\mathcal{F}_1$ term. This term diverges for both $e_2$ approaching unity and $e_1$ approaching zero. In the former case, the precession frequency of the inner orbit diverges, making polar alignment impossible. In the latter case, the nodal libration frequency of the polar disc tends to zero. Again, as the inner binary becomes more circular, the process of polar alignment becomes less likely.

As for the mass ratios, polar alignment is made difficult for high values of $q_2$ and for $q_1$ approaching zero. Higher values of $q_2$ result in a higher precession frequency. Meanwhile, very low $q_1$ reduce the inner binary to a single star, reducing the torque on the disc to zero. Conversely, when $q_2$ tends to zero, the third outer body perturbation becomes negligible, which slows the precession and favours polar alignment.

Fig. 5 shows how $T_1$ depends on the system orbital parameters. In this plot, we assume $R_{\text{in}} = 1.5 a_1$ and $R_{\text{out}} = 5 a_1$, thus $x_o = 5/1.5$. In the left panel, $T_1$ is plotted as a function of the inner binary mass ratio $q_1$ (fixing $q_2$ to 0.5), while on the right one as a function of the outer binary mass ratio $q_2$ (fixing $q_1$ to 0.5). The thicker lines refer to a configuration in which $e_1 = 0.5$ and $e_2 = 0$. The dotted and dashed lines show how raising $e_2$ to 0.7 (fixing $e_1$ to 0.5) and lowering $e_1$ to 0.2 (fixing $e_2$ to zero) affect $T_1$, respectively.

A disc orbiting the innermost binary of a hierarchical system is able to go polar in the vast majority of the inner binary parameter space, contrary to an accretion disc around the outer orbit of a hierarchical system. Except for very small $q_1$ and very high $q_2$, the $T_1$ factor is always below unity no matter the eccentricities, semi-major axis ratio, or mass ratios.

In the case of HD98800, we have precise measurements of the orbital parameters of both the inner orbits and of the outer orbit of the quadruple stellar system (Zúñiga-Fernández et al. 2021). If we approximate HD98800 A to a single body, the disc around HD98800 B measures a $T_{1,\text{HD}} \approx 3 \times 10^{-4}$. This value is well below unity as expected, being the disc in a polar configuration.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Dependency of $T_1$ (see sec. 2.3) on mass ratios, semi-major axes ratio, and eccentricities. When $T_1$ is smaller than unity (dash dotted red line), polar alignment is possible, conversely for $T_1 > 1$ it is not. On the left panel, $T_1$ is plotted as a function of the inner binary mass ratio $q_1$ (fixing $q_2$ to 0.5), while on the right one as a function of the outer binary mass ratio $q_2$ (fixing $q_1$ to 0.5). The thicker lines refer to $e_2 = 0$ and $e_1 = 0.5$. The dotted and dashed lines show how lowering $e_1$ to 0.2 (fixing $e_2 = 0$) and raising $e_2$ to 0.7 (fixing $e_1 = 0.5$) affect $T_1$, respectively. Different colours refer to different semi-major axes ratios.</figcaption></figure>

## 3 NUMERICAL SIMULATIONS

We performed 3D numerical simulations using the smoothed particle hydrodynamics (SPH) code PHANTOM (Price et al. 2018a), widely used in the astrophysical community to study gas and dust in protostellar environments (e.g. Dipierro et al. 2015; Mentiplay, Price & Pinte 2019; Veronesi et al. 2020; Ragusa et al. 2020; Toci et al. 2020; Ballabio et al. 2021; Vericel et al. 2021) and for simulating the hydrodynamics of stellar systems with more than two stars embedded in accretion discs (e.g. Ragusa et al. 2017, 2021; Martin & Lubow 2017; Price et al. 2018b; Cuello et al. 2019; Poblete, Cuello & Cuadra 2019; Calcino et al. 2019; Smallwood et al. 2021; Ceppi et al. 2022). These simulations are designed to test the criteria showed in Sections 2.2 and 2.3 for polar alignment in discs orbiting a hierarchical triple outer (Section 3.1) and inner (Section 3.2) orbit.

### 3.1 Circumtriple disc simulations

Simulations S1 and S2 test the polar alignment of a disc around the outer orbit of a hierarchical triple (S1), compared to the same disc orbiting around the binary associated to the triple outer orbit (i.e. the binary with the same orbital parameters of the triple outer orbit and made of the third body and a star obtained by condensing the triple inner binary in its centre of mass, simulation S2). The triple system consists of a circular equal mass inner binary ($m_a = m_b = 0.5 M_\odot$) orbited by a coplanar third star with mass equal to the inner binary mass ($m_c = m_a + m_b = 1 M_\odot$). The semi-major axis of the inner binary is $a_1 = 1$ au. The semi-major axis and the eccentricity of the outer orbit are $a_2 = 10$ au and $e_2 = 0.5$. The outer and the inner orbital plane are coplanar. The associated binary system consists of a star of mass $m_c$ orbited by a star of mass $m_a + m_b$. The eccentricity and semi-major axis of the binary are the same as the triple outer binary orbit.

Both the binary and triple system are orbited by the same gaseous accretion disc. The disc is initially tilted to the stellar orbital plane of $\beta = 70^\circ$ and the longitude of the ascending node $\Omega$ equals $90^\circ$.

Each simulated disc is made of $3 \times 10^5$ SPH particles, resulting in a vertical resolution of $<h/H> \sim 0.24$. The total gas mass is initially the 1 per cent of the stellar mass and the discs extend from $R_{in} = 15$ au to $R_{out} = 100$ au. The gas surface density is initially distributed with a power-law profile

$$ \Sigma(R) = \Sigma_{\text{in}} \left(\frac{R}{R_{\text{in}}}\right)^{-p} \left(1 - \sqrt{\frac{R_{\text{in}}}{R}}\right), \qquad (13) $$

with $p = 1$ and $\Sigma_{\text{in}} = 22.2$ g/cm$^2$. We adopt a locally isothermal equation of state $P = c_s^2 \rho$, with

$$ c_s(R) = c_s(R_{\text{in}}) \left(\frac{R}{R_{\text{in}}}\right)^{-q}, \qquad (14) $$

with $q = 0.25$. The resulting disc aspect ratio is given by

$$ \frac{H}{R} = \frac{H_0}{R_0} \left(\frac{R}{R_{\text{in}}}\right)^{1/4}, \qquad (15) $$

with an $H_0/R_0 = 0.1$ at $R = R_{in}$. Viscosity is implemented with the artificial viscosity method (Lucy 1977; Gingold & Monaghan 1977) resulting in a Shakura & Sunyaev (1973) $\alpha$-viscosity (Lodato & Price 2010) with $\alpha_{SS} \sim 0.01$ ($\alpha_{AV} = 0.42$). In every simulation, stars are modelled as sink particles (Bate, Bonnell & Price 1995; Price et al. 2018a).

The initial tilt of the disc is well above the critical angle for polar alignment around a binary system, which is $\beta_{crit}(e = 0.5, \Omega = 90^\circ) \sim 40^\circ$ (Farago & Laskar 2010). Thus, we expect the circumbinary disc in simulation S2 to polarly align. Conversely, even if the outer orbit parameters in the hierarchical triple of simulation S1 are the same as the binary in S2, $T_2 \sim 6 > 1$. Thus, we expect the disc not to align polarly.

Fig. 6 shows radially averaged tilt ($\beta$, top panels) and longitude of the ascending node ($\Omega$, bottom panels) profiles of the disc as a function of time. The two systems behave very differently: the circumbinary disc undergoes polar alignment, since the inclination oscillates around $\beta \sim 90^\circ$ and the longitude of the ascending
<figure><img src="image_7.png" /></figure>

**Figure 6.** Radially averaged tilt (top row) and longitude of the ascending node (bottom row) of the circumbinary (left column, simulation S2) and circumtriple (right column, simulation S1) disc as a function of time. Dashed line is the longitude of the ascending node with respect to the initial semi-major axis direction, rather than the eccentricity vector. In the left column, the circumbinary disc is aligning polarly. As expected due to the mass of the disc, the polar alignment inclination is less than 90°.

In the right column, the circumtriple disc is going coplanar, even if the orbital parameters of the triple outer orbit are the same as the binary orbital parameters.

node librates, as expected for a disc going polar (Martin & Lubow 2017). On the other hand, the circumtriple disc aligns to the stellar orbital plane (i.e. the inclination decreases) and it precesses, with its longitude of the ascending node spanning $2\pi$. There are two processes that make the longitude of the ascending node to precess: (i) the eccentricity vector, to which the longitude of the ascending node refers, is precessing; (ii) the disc itself is precessing, as an inclined circumbinary disc would around a pure binary. This can be seen computing the longitude of the ascending node referring to the initial eccentricity vector position. This absolute longitude of the ascending node is plotted as a dashed curve in Fig. 6. As expected, polar alignment does not occur in this case, since the precession rate of the eccentricity vector is larger than the polar alignment one.

The disc inclination, apart from the oscillating behaviour, is exponentially decaying towards the mid-plane. We fit the exponential decay to find the time-scale of coplanar alignment, obtaining $\tau_{\text{cop}} \sim 3.8 \times 10^4$ outer orbit periods. The time-scale $\tau_{\text{cop}}$ is of the order of the disc viscous time-scale for this disc ($\sim 2 \times 10^4$ outer orbits). Bate et al. (2000) found that the alignment time-scale (given by their equation (35), see also Lubow & Ogilvie (2000)) for a tilted disc is of the order of (or slightly longer than) the disc viscous time-scale. This suggests that an accretion disc with $T_2 > 1$ orbiting the outer levels of a hierarchical system sees the central system as a circular binary. In fact, the evolution of the tilt angle of the circumtriple disc in simulation S1 is similar to a circumbinary disc orbiting a low eccentric binary system (an exponentially decaying tilt and a precessing longitude of the ascending node). This is due to the fact that the disc sees the eccentricity vector of the central stellar system averaged over the precession period.

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

**Figure 7.** Radially averaged tilt (top row) and longitude of the ascending node (bottom row) of the circum-inner binary disc (simulation S3) as a function of time. Dashed line is the longitude of the ascending node with respect to the initial semi-major axis direction, rather than the eccentricity vector. The disc is aligning to a polar configuration. The libration time-scale and the stationary inclination differ from the circumtriple disc due to the lower mass of the circum-inner binary disc.

## 3.2 Circumbinary disc simulation in a triple system

Simulation S3 studies a disc orbiting the inner binary of a hierarchical triple system. The stellar system consists of an eccentric ($e_1 = 0.5$) equal mass inner binary ($m_a = m_b = 0.5 M_\odot$) orbited by a third star with mass equal to the inner binary mass ($m_c = m_a + m_b = 1 M_\odot$) in a circular orbit. The semi-major axis of the inner and outer binary are $a_1 = 1$ au and $a_2 = 20$ au, respectively, and the inner and outer orbits are coplanar.

The inner binary is orbited by a circumbinary disc with an initial tilt of $\beta = 70^\circ$ and a longitude of the ascending node of $90^\circ$. The disc is simulated with $3 \times 10^5$ SPH particles, resulting in a vertical resolution of $\langle h/H \rangle \sim 0.23$. The total gas mass is initially the $1\%$ of the inner binary mass and the disc extends from $R_{\text{in}} = 1.5$ au to $R_{\text{out}} = 5$ au. The gas surface density profile is the same as in equation (13), with $\Sigma_{\text{in}} = 2.7 \times 10^2 \text{ g cm}^{-2}$, and the temperature and aspect ratio profiles are given by equations 14 and 15. Viscosity and stars are implemented as in previous simulations, with $\alpha_{\text{SS}} \sim 0.01$ ($\alpha_{\text{AV}} = 0.44$).

The initial disc tilt is above the critical angle for polar alignment around a pure binary system. In addition, in this configuration $T_1 \sim 0.008 < 1$, which should translate into disc polar alignment (see Section 2.3).

Fig. 7 shows radially averaged tilt ($\beta$, top panel) and longitude of the ascending node ($\Omega$, bottom panel) profiles of the disc as a function of time. The tilt plot clearly shows that the disc oscillates around $90^\circ$ with the amplitude being damped over time. The longitude of the ascending node is librating as well: thus, the disc is going polar as expected. However, the disc is librating around the eccentricity vector of the inner binary, which in turn is precessing with time. As a consequence, the disc angular momentum is precessing as well. This
can be seen computing the disc longitude of the ascending node with
respect to the initial semi-major axis direction, instead of referring
to the inner binary eccentricity vector, which is precessing. This
absolute longitude of the ascending node is plotted in the second row
of Fig. 7 with a dashed curve. The absolute longitude of the ascending
node is oscillating as the relative longitude of the ascending node,
but its mean value is decreasing with time following the precession
of the inner binary eccentricity vector.

# 4 DISCUSSION

## 4.1 Polar alignment in hierarchical stellar systems

Hierarchical systems with more than two stars are common in a young
stellar population (Duchêne & Kraus 2013; Moe & Di Stefano 2017),
and they typically host discs orbiting their hierarchical levels. A
significant fraction of misaligned stellar discs studied in the literature
are orbiting or are inside systems with more than two stars (e.g. Price
et al. 2018b; Kennedy et al. 2019; Phuong et al. 2020; Keppler et al.
2020b; Kraus et al. 2020; Bi et al. 2020). By studying the observed
misalignment distribution of discs in multiple stellar systems, we
aim at obtaining information about their physical properties, their
formation, and early stage evolution. Properties of discs (e.g. size,
viscosity) affect the evolution of the tilt angle which, in return,
provides information about them. In order to extract information
from the final observed tilt distribution, we need to understand the
way in which it evolves and, possibly, the initial conditions of its
evolution.

In particular, the fraction of discs in a polar configuration gives
us insights about the fraction of highly eccentric binary systems
with highly inclined discs. Remarkably, the time evolution towards
a polar orbit is faster than the disc lifetime. So, in principle, the final
outcome of the evolution process (i.e. the observed fraction of polar
discs) directly traces the initial conditions (i.e. the outcome of star
formation process). Indeed, the fraction of polarly aligned discs is
related to how common the critical configuration is (i.e. a binary with
eccentricity $e_b$, orbited by a disc with tilt $\beta > \beta_{crit}(e_b, \Omega)$). However,
we have to take into account that pure binary systems and systems
with more than two stars have different critical conditions for polar
alignment.

The analytical and numerical findings, reported in Sections 2 and 3, show that the secular evolution of the osculating orbital elements in hierarchical systems is a threat to disc polar alignment. In particular, a fast enough precession of the orbit eccentricity vector prevents the mechanism of polar alignment. Figs. 2 and 5 show the *T* parameter discussed in Sections 2.2 and 2.3 (respectively) as a function of the relevant parameters of a hierarchical triple (the two mass ratios *q*₁ and *q*₂, the semi-major axis ratio *a* and the inner and outer orbit eccentricities *e*₁ and *e*₂). When *T* is larger than unity, the polar orbit is unstable – no matter how high the disc–orbit misalignment is. Under this condition, the pure binary criteria for polar alignment derived in the literature cannot be applied to hierarchical systems.

Looking at the parameter space explored in Fig. 2 and at the
numerical results in section 3.1, discs orbiting the outer levels of
hierarchical systems are generally not able to go polar. Indeed, the *T*
parameter is always larger than unity, except for radially narrow discs
and where the system is reduced to an actual binary system (i.e. for
a mass ratio of the inner binary (*q*₁) approaching zero and for very
small semi-major axis ratios (*a*)). Thus, regardless of the mutual
inclination, a misaligned radially extended disc orbiting an outer
level of a hierarchical system always evolves towards a coplanar
configuration. Since the eccentricity vector is quickly precessing,

the disc evolves as orbiting a circular orbit given that the eccentricity
vector is averaged over a precession period.

In a circumbinary disc in a hierarchical system, however, polar
alignment is possible as the *T* factor is generally below unity (see
Figs. 5 and 7). As the precession of the inner binary eccentricity
vector is slow enough for the polar disc to follow it, the disc remains
locked to the eccentricity vector, precessing with it and conserving
the polar configuration. Note that the disc precession happens on
the eccentricity vector precession time-scale and not on the typical
nodal precession time-scale for a tilted circumbinary disc (Bate et al.
2000; Lodato & Facchini 2013). Superimposed to this precession,
the longitude of the ascending node is librating as well, due to the
precession of the disc angular momentum vector around the inner
binary eccentricity vector.

The evolution of the osculating elements in hierarchical systems further affects the polar alignment process in such systems. The critical angle for polar alignment significantly depends on $\Omega$, which is relative to the central orbit eccentricity vector (see equation 2). Given that this vector in hierarchical system precesses with time, $\Omega$ is constantly spanned with time. Thus, regardless of the initial $\Omega$ value, the triple configuration will precess, eventually exploring the $\Omega$ related to the lowest critical angle. Therefore, even if the initial $\Omega$ results in a configuration where $\beta < \beta_{crit}$ (so no polar alignment for a pure binary configuration), the triple configuration will polarly align as soon as the $\Omega$ precession allows it to fullfil the condition $\beta > \beta_{crit}(\Omega, e_b)$. Thus, given the eccentricity of the orbit, the critical angle for polar alignment in hierarchical systems is always the one with $\Omega = 90^\circ$ (the lower one).

Up to now, we only discussed hierarchical stellar systems with
inner and outer orbits coplanar to each other. When the inner and
the outer orbit of a hierarchical triple are misaligned, the secular
evolution of the osculating orbital parameters is even more complex.
Indeed, not only the argument of periapsis precesses with time,
but also Kozai–Lidov oscillations of the stellar orbits occur. Thus,
the eccentricity of each hierarchical level orbit and the mutual
inclination between the hierarchical orbital planes oscillate with
time. However, at least for a circular outer orbit, the Kozai–Lidov
oscillation frequency is similar to the precession frequency of the
inner binary orbit we derived in equation (10) (Antognini 2015;
Naoz 2016). Therefore, the *T* parameter is still reliable to predict
the possibility of polar alignment in the inner orbit of misaligned
hierarchical triples. However, a thorough analysis of this problem
is needed, especially regarding the possibility of polar alignment in
misaligned hierarchical triple circum-outer orbit accretion discs and
we defer this to subsequent work. An aspect of crucial importance
is that Kozai–Lidov oscillations can trigger polar alignment in
misaligned hierarchical systems. Indeed, in the oscillation, the orbit
eccentricities raise (lowering the critical angle) and the orbital
inclination oscillation makes also the tilt of the disc to change. This
can trigger polar alignment in configurations with a lower-than-unity
*T* parameter, but a low initial orbit–disc misalignment.

We thus expect polar alignment to still be possible only around the
innermost level of hierarchical systems even for misaligned multiple
systems, where it can be even promoted by Kozai–Lidov oscillations
and the $\Omega$ precession.

## 4.2 Observed tilt distribution of discs around binaries and hierarchical systems

Fig. 8, originally reported by Czekala et al. (2019), collects all
protoplanetary discs in the literature orbiting in multiple stellar
systems for which we know the inclination between the disc and
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Disc-orbit misalignment of circumbinary discs in tables 3 and 4 in Czekala et al. (2019) as a function of the central system eccentricity. Green dots are accretion discs orbiting pure binary systems, blue dots are accretion discs orbiting an outer hierarchical level of a hierarchical system, red dots are accretion discs orbiting the innermost hierarchical level of a hierarchical system (HS). The triangle represents the lower limit on the mutual inclination of R Cra. Dotted lines connect the two degenerate solutions for HD 142527, SR 24N, and GG Tau A systems.</figcaption></figure>

the system orbital plane. Each point in the plot represents a disc orbiting a pair of bound stars, either in a pure binary system or in a higher multiplicity system. Green dots refer to discs orbiting in a pure binary system, blue dots refer to discs that have more than two stars orbiting in their cavity (i.e. the outer level of hierarchical systems), and red dots refer to discs orbiting a binary with external companions outside the disc (i.e. the innermost level of hierarchical systems). The dots connected by dotted vertical lines are two possible solutions for the mutual inclination due to a 180° ambiguity in the longitude of the ascending node of the system due to lack of radial velocity information.

The three populations show different tilt distribution. Discs orbiting pure binary systems (green dots) are mostly coplanar, hence there are no polarly aligned discs. This mismatches with the theoretical expectation. If the initial tilt distribution is nearly randomly distributed, we expect a wider scatter in the distribution because the population should be slowly going coplanar on a tilt evolution time-scale comparable to the viscous time-scale (Bate et al. 2000; Lubow & Ogilvie 2000). This could be due to a lower average tilt in the initial population or to a faster than expected tilt evolution (e.g. due to an higher than expected viscosity). In addition, we lack the expected small population of polarly aligned discs (resulting from the fraction of discs initially more tilted than the critical angle for polar alignment in binaries). Hence, either the conditions for polar alignment are less populated than expected or external factors reduce the stability of the polar configuration (e.g. the interaction with the environment).

The tilt distribution of circumbinary discs in hierarchical systems is more articulated. We have a small fraction of polarly aligned discs (i.e. the disc in HD98800B and possibly in SR 24N), although there are less discs in this population than in the pure binaries one. This is in agreement with our analytical findings (polar configuration is

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Disc-orbit misalignment of circumbinary discs in tables 3 and 4 in Czekala et al. (2019) as a function of the orbit semi-major axis. Legend is the same as Fig. 8.</figcaption></figure>

stable in such systems) and suggests that the conditions for polar alignment are more likely for those systems than in pure binaries. This is also in line with what discussed in Section 4.1: phenomena that are typical of systems with more than two stars could be able to foster polar alignment. First, polar alignment for such configurations could be triggered by Kozai-Lidov oscillations of the hierarchical system. Second, by varying the longitude of the ascending node, the precession of the orbit eccentricity vector allows the disc to always reach the configuration in which the critical angle for polar alignment is minimum. We defer the study of the impact of these mechanisms on polar alignment to future works. Additionally, depending on the stellar system parameters (such as binary mass, semi-major axis, and eccentricity), general relativity could limit the possibility of polar alignment for extended discs (Lepp, Martin & Childs 2022). Thus, the presence of an additional external companion truncating the disc from the outside could facilitate polar alignment for circum-inner binary discs as well. The rest of the population is nearly coplanar. Thus, discs that do not go polar evolve as in the pure binary case.

It is worth noting that also the binary separation plays a role in the degree of alignment of circumbinary discs (see Fig. 9). Indeed, short-period binaries (no matter if isolated or with an external companion) present coplanar discs. The inclination distribution of that region of the parameter space could be affected also by the short-period binary formation mechanisms that drive initially wider binaries to shrink. Still, wider binary inclination distributions present the same trends discussed in the previous paragraph (i.e. highly misaligned discs around inner binaries and more coplanar discs around pure binaries).

Even though the statistics are low, in the population of discs orbiting more than two stars, there are no coplanar discs, neither highly misaligned (possibly polar) ones. Indeed, the highest blue point refer to the degenerate solution with a high inclination for the GG Tau A circumtriple disc. All the hydrodynamical models of GG Tau A in the literature favour the mildly inclined solution (Cazzolletti et al. 2017; Aly et al. 2018; Keppler et al. 2020a). Thus, the correct tilt for GG Tau A is around 30° the uncertainty in the astrometry and the highest blue point should be discarded. The lack of polar discs agrees
with our analytical findings. The lack of coplanar discs contrasts with our results in discs orbiting coplanar hierarchical systems. Indeed, we found that circum triple discs evolve as orbiting a circular binary (due to the fast precession of the eccentricity vector). In addition, given that they form via the same mechanisms, we expect the same initial tilt distribution for discs in pure binaries and in hierarchical systems. If the initial condition and the evolution are similar, then *why is the tilt distribution of pure binaries and systems with more than two stars so different?*

We suggest the answer lies in the orbital dynamics of systems with more than two stars. This kind of systems have access to a richer dynamical evolution compared to pure binary systems. As previously discussed their orbital parameters and orbital orientation evolve with time on a shorter time-scale compared to the disc lifetime. In particular, misaligned hierarchical systems vary their inclinations with time. Specifically, let us consider a hierarchical triple system orbited by a low-misaligned accretion disc. Looking at the pure binary population, on the long run the disc should become coplanar with the triple orbital plane. If the inner and outer orbit of the triple are misaligned, however, we expect their mutual inclination to evolve with time due to Kozai–Lidov oscillations. Thus, the inclination of the outer orbit will oscillate as well. As a consequence, the mutual inclination between the circumtriple disc and the triple outer orbit will oscillate with time, even if the disc had enough time to align to the stellar plane. Thus, the observed misalignment in the systems with more than two star population is possibly the result of these stellar orbit evolution processes, and not the outcome of accretion disc evolution. This implies that the misaligned configurations we observe in systems with more than two stars are not stable – or slowly evolving – configurations. Indeed, because we are taking a snapshot of an oscillating stellar orbital plane, we happen to be observing a tilted disc by pure chance.

## 5 CONCLUSIONS

In this work, we showed that the requirement on the parameters of a circumbinary disc for going polar are necessary but not suffi- cient when dealing with hierarchical systems. A crucial additional requirement is that the eccentricity vector precession time-scale of the system orbited by the accretion disc has to be longer than the disk libration time-scale. We derived an analytical criterion to be satisfied in order for a disc to go polar around the outer levels of a hierarchical system (equation (6)) and around the innermost hierarchical level (equation (11)).

We found that discs orbiting the outer level of a hierarchical system can hardly polarly align. Except for radially narrow discs and very small mass and semi-major axis ratios, the precession of the outer orbit eccentricity vector is always faster than the nodal libration. Conversely, discs orbiting the innermost level of a hierarchical system are able to go polar as the precession of the inner eccentricity vector is slower than polar alignment. SPH simulations confirm these results and surveys of circumbinary accretion discs are also in agreement – even though statistics are still poor.

In addition, we found that the inclination of an accretion disc orbiting a hierarchical system evolves as if it was orbiting a circular binary. The disc sees the central system eccentricity vector averaged over the precession period. This is in contrast with the fact that pure binary systems and systems with more than two stars host different disc populations. The former presents discs mostly coplanar with the stellar orbital plane with little spread in the tilt distribution. The latter consists of misaligned discs with a wider spread in tilt distribution. We suggest that the scatter in misalignment observed in system with

more than two stars is due to the secular oscillation of hierarchical system orbital parameters that continuously vary the stellar orbital plane orientation, rather than to the evolution of the accretion discs.

In conclusion, this work shows that when analysing the population of discs around binaries, it is important to separate pure binaries from binaries inside hierarchical systems. Indeed, the misalignment distributions of pure binaries and systems with more than two stars tell us different stories. The former is the result of viscous disc evolution, while the second one is related to N-body orbital parameter oscillations. This additional complexity can be a precious tool to deeply investigate the impact of multiplicity on accretion disc evolution and planet formation, and, in general, to better understand the formation of multiple stellar systems. The increasing statistics of this population will allow us also to better measure the fraction of polar discs and the distribution of mutual inclination in the two populations. This is an important test for our theory. Indeed, the final distribution have to be the result of the initial conditions we derived with models, plus the theoretical expected evolution. A mismatch between the predicted and the observed populations would highlight missing pieces in the theory of disc formation and evolution.

## ACKNOWLEDGEMENTS

The authors thank the referee for their constructive feedback and sug- gestions, which have significantly improved the original manuscript. The authors thank KAVLI Institute of Theoretical Physics in Santa Barbara for the hospitality during the completion of this work. This project and the authors have received funding from the Eu- ropean Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreements N. 823823 (DUSTBUSTERS RISE project), 896319 (SANDS), and 101042275 (project Stellar-MADE). CL acknowledges support from Fulbright Commission through VRS scholarship. This research was funded, in part, by ANR (Agence Nationale de la Recherche) of France under contract number ANR-22-ERCS-0002-01. SC and CL thank Daniel Elsender, Daniel Price, Matthew Bate, and Yan-Fei Jiang for useful discussions. We used the following Python tools and packages: NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007) and JUPYTER (Kluyver et al. 2016).

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author. The code PHANTOM used in this work is publicly available at https://github.com/danieljprice/phantom.

## REFERENCES

* Aly H., Lodato G., 2020, MNRAS, 492, 3306

* Aly H., Dehnen W., Nixon C., King A., 2015, MNRAS, 449, 65

* Aly H., Lodato G., Cazzoletti P., 2018, MNRAS, 480, 4738

* Aly H., Gonzalez J.-F., Nealon R., Longarini C., Lodato G., Price D. J., 2021, MNRAS, 508, 2743

* Antognini J. M. O., 2015, MNRAS, 452, 3610

* Ballabio G., Nealon R., Alexander R. D., Cuello N., Pinte C., Price D. J., 2021, MNRAS, 504, 888

* Bate M. R., 2009, MNRAS, 397, 232

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

* Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362

* Bate M. R., Bonnell I. A., Clarke C. J., Lubow S. H., Ogilvie G. I., Pringle J. E., Tout C. A., 2000, MNRAS, 317, 773

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

* Bi J. et al., 2020, ApJ, 895, L18
Calcino J., Price D. J., Pinte C., van der Marel N., Ragusa E., Dipierro G., Cuello N., Christiaens V., 2019, MNRAS, 490, 2579

Casassus S. et al., 2013, Nature, 493, 191

Cazzoletti P., Ricci L., Birnstiel T., Lodato G., 2017, A&A, 599, A102

Ceppi S., Cuello N., Lodato G., Clarke C., Toci C., Price D. J., 2022, MNRAS, 514, 906

Chiang E. I., Murray-Clay R. A., 2004, ApJ, 607, 913

Clarke C. J., Bonnell I. A., Hillenbrand L. A., 2000, in Mannings V., Boss A. P., Russell S. S., eds, Proc. Protostars and Planets IV, University of Arizona Press, Tucson, AZ, p. 151,

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

Cuello N. et al., 2019, MNRAS, 483, 4114

Cuello N., Ménard F., Price D. J., 2023, Eur. Phys. J. Plus, 138, 11

Czekala I., Chiang E., Andrews S. M., Jensen E. L. N., Torres G., Wilner D. J., Stassun K. G., Macintosh B., 2019, ApJ, 883, 22

Dipierro G., Price D., Laibe G., Hirsh K., Cerioli A., Lodato G., 2015, MNRAS, 453, L73

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

Facchini S., Lodato G., Price D. J., 2013, MNRAS, 433, 2142

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

Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375

Hamers A. S., 2021, MNRAS, 500, 3481

Harris C. R. et al., 2020, Nature, 585, 357

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Kennedy G. M. et al., 2012, MNRAS, 421, 2264

Kennedy G. M. et al., 2019, Nat. Astron., 3, 230

Kenworthy M. A. et al., 2022, A&A, 666, A61

Keppler M. et al., 2020a, A&A, 639, A62

Keppler M. et al., 2020b, A&A, 639, A62

Kluyver T. et al., 2016, in Loizides F., Schmidt B., eds, Proc. International Conference on Electronic Publishing, Positioning and Power in Academic Publishing: Players, Agents and Agendas. IOS Press, Amsterdam. p. 87

Köhler R., 2011, A&A, 530, A126

Kozai Y., 1962, AJ, 67, 591

Kraus S. et al., 2020, Science, 369, 1233

Larwood J. D., Papaloizou J. C. B., 1997, MNRAS, 285, 288

Larwood J. D., Nelson R. P., Papaloizou J. C. B., Terquem C., 1996, MNRAS, 282, 597

Lepp S., Martin R. G., Childs A. C., 2022, ApJ, 929, L5

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

Lidov M. L., 1962, Planet. Space Sci., 9, 719

Lodato G., Facchini S., 2013, MNRAS, 433, 2157

Lodato G., Price D. J., 2010, MNRAS, 405, 1212

Longarini C., Lodato G., Toci C., Aly H., 2021, MNRAS, 503, 4930

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

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

Lubow S. H., Ogilvie G. I., 2001, ApJ, 560, 997

Lubow S. H., Martin R. G., Nixon C., 2015, ApJ, 800, 96

Lucy L. B., 1977, AJ, 82, 1013

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

Martin R. G., Lubow S. H., 2018, MNRAS, 479, 1297

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

Martin R. G., Nixon C., Lubow S. H., Armitage P. J., Price D. J., Doğan S., King A., 2014, ApJ, 792, L33

Martin R. G., Lepp S., Lubow S. H., Kenworthy M. A., Kennedy G. M., Vallet D., 2022, ApJ, 927, L26

McKee C. F., Ostriker E. C., 2007, ARA&A, 45, 565

Mentiplay D., Price D. J., Pinte C., 2019, MNRAS, 484, L130

Miranda R., Lai D., 2015, MNRAS, 452, 2396

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

Murray C. D., Dermott S. F., 2000, Solar System Dynamics. Cambridge University Press, Cambridge

Naoz S., 2016, ARA&A, 54, 441

Naoz S., Farr W. M., Lithwick Y., Rasio F. A., Teyssandier J., 2013, MNRAS, 431, 2155

Nixon C., King A., Price D., Frank J., 2012, ApJ, 757, L24

Nixon C., King A., Price D., 2013, MNRAS, 434, 1946

Offner S. S. R., Moe M., Kratter K. M., Sadavoy S. I., Jensen E. L. N., Tobin J. J., 2022, ARA&A, preprint (arXiv:2203.10066)

Ogilvie G. I., 1999, MNRAS, 304, 557

Papaloizou J. C. B., Lin D. N. C., 1995, ApJ, 438, 841

Papaloizou J. C. B., Pringle J. E., 1983, MNRAS, 202, 1181

Papaloizou J. C. B., Terquem C., 1995, MNRAS, 274, 987

Phuong N. T. et al., 2020, aap, 635, A12

Poblete P. P., Cuello N., Cuadra J., 2019, MNRAS, 489, 2204

Price D. J. et al., 2018a, PASA, 35, e031

Price D. J. et al., 2018b, MNRAS, 477, 1270

Ragusa E., Dipierro G., Lodato G., Laibe G., Price D. J., 2017, MNRAS, 464, 1449

Ragusa E., Alexander R., Calcino J., Hirsh K., Price D. J., 2020, MNRAS, 499, 3362

Ragusa E. et al., 2021, MNRAS, 507, 1157

Reipurth B., Clarke C. J., Boss A. P., Goodwin S. P., Rodríguez L. F., Stassun K. G., Tokovinin A., Zinnecker H., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Proc. Protostars and Planets VI, University of Arizona Press, Tucson, AZ. p. 267

Shakura N. I., Sunyaev R. A., 1973, in Bradt H., Giacconi R., eds, Proc. IAU Symp. 55, X- and Gamma-Ray Astronomy. International Astronomical Union, Dordrecht, p.155

Smallwood J. L., Lubow S. H., Franchini A., Martin R. G., 2019, MNRAS, 486, 2919

Smallwood J. L., Franchini A., Chen C., Becerril E., Lubow S. H., Yang C.-C., Martin R. G., 2020, MNRAS, 494, 487

Smallwood J. L., Nealon R., Chen C., Martin R. G., Bi J., Dong R., Pinte C., 2021, MNRAS, 508, 392

Toci C., Galli D., Verdini A., Del Zanna L., Landi S., 2018, MNRAS, 474, 1288

Toci C., Lodato G., Christiaens V., Fedele D., Pinte C., Price D. J., Testi L., 2020, MNRAS, 499, 2015

Vericel A., Gonzalez J.-F., Price D. J., Laibe G., Pinte C., 2021, MNRAS, 507, 2318

Veronesi B. et al., 2020, MNRAS, 495, 1913

von Zeipel H., 1910, Astron. Nachr., 183, 345

Xiang-Gruess M., 2016, MNRAS, 455, 3086

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

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

This paper has been typeset from a *TEX/LATEX* file prepared by the author.