Formation of the polar debris disc around 99 Herculis

ABSTRACT

We investigate the formation mechanism for the observed nearly polar aligned (perpendicular to the binary orbital plane) debris ring around the eccentric orbit binary 99 Herculis. An initially inclined non-polar debris ring or disc will not remain flat and will not evolve to a polar configuration, due to the effects of differential nodal precession that alter its flat structure. However, a gas disc with embedded well coupled solids around the eccentric binary may evolve to a polar configuration as a result of pressure forces that maintain the disc flatness and as a result of viscous dissipation that allows the disc to increase its tilt. Once the gas disc disperses, the debris disc is in a polar aligned state in which there is little precession. We use three-dimensional hydrodynamical simulations, linear theory, and particle dynamics to study the evolution of a misaligned circumbinary gas disc and explore the effects of the initial disc tilt, mass, and size. We find that for a wide range of parameter space, the polar alignment time-scale is shorter than the lifetime of the gas disc. Using the observed level of alignment of 3° from polar, we place an upper limit on the mass of the gas disc of about |$0.014 \, \mathrm{M}_\odot$| at the time of gas dispersal. We conclude that the polar debris disc around 99 Her can be explained as the result of an initially moderately inclined gas disc with embedded solids. Such a disc may provide an environment for the formation of polar planets.

accretion, accretion discs, hydrodynamics, planets and satellites: formation, binaries: general

1 INTRODUCTION

The majority of stars that form within dense regions of stellar clusters are formed as binary systems (Duquennoy & Mayor 1991; Ghez, Neugebauer & Matthews 1993; Duchêne & Kraus 2013), which are most likely accompanied by circumstellar and circumbinary discs (e.g. Dutrey, Guilloteau & Simon 1994; Beust & Dutrey 2005). The evolution of circumbinary disc structure and orientation has been studied extensively. An initially slightly misaligned circumstellar or circumbinary disc involving a circular orbit binary precesses about the binary angular momentum vector and evolves towards alignment with it due to viscous dissipation in the disc. As a result, the disc becomes coplanar with the binary (Papaloizou & Terquem 1995; Lubow & Ogilvie 2000; Nixon, King & Pringle 2011; Facchini, Lodato & Price 2013; Foucart & Lai 2014). If the binary orbit is eccentric, a low-mass circumbinary disc with a large enough initial inclination can precess around the eccentricity vector (semimajor axis) of the binary. The disc’s angular momentum vector eventually aligns with the eccentricity vector. This means that the disc angular momentum is aligned polar (perpendicular) with respect to the binary angular momentum (Aly et al. 2015; Martin & Lubow 2017, 2018; Lubow & Martin 2018; Zanazzi & Lai 2018). The disc then lies perpendicular to the orbital plane of the binary. A massive disc aligns to a generalized polar state at lower misalignment to the binary orbital plane (Zanazzi & Lai 2018; Martin & Lubow 2019).

Disc misalignment can occur in various phases of stellar evolution. Misaligned discs around binaries can arise from chaotic accretion of turbulent molecular clouds in star-forming regions (Offner et al. 2010; Bate 2012; Tokuda et al. 2014) or whenever a young binary system accretes material post-formation (Bate, Lodato & Pringle 2010; Bate 2018). Furthermore, misalignments can be present if the binary forms within an elongated cloud, where the binary axis is misaligned with respect to the cloud rotation axis (e.g. Bonnell & Bastien 1992).

There are currently a number of observed systems with misaligned circumbinary discs. KH 15D is an eccentric spectroscopic binary T Tauri star with a misaligned circumbinary disc (Chiang & Murray-Clay 2004; Winn et al. 2004; Lodato & Facchini 2013; Smallwood et al. 2019). High misalignment has been observed in the binary protostar IRS 43, where the tilt of the disc is at least 60° with respect to the orbital plane of the binary (Brinch et al. 2016). The binary GG Tau consists of T Tauri stars with a circumbinary disc misaligned by 25°–30° along with misaligned discs around each of the binary components (Dutrey et al. 1994; Köhler 2011; Cazzoletti et al. 2017; Aly, Lodato & Cazzoletti 2018). The binary system HD 98800 BaBb shows evidence of having a nearly polar circumbinary gas disc (Kennedy et al. 2019). Lastly, misalignment can be observed also after the gas disc has been dispersed. The binary 99 Herculis (99 Her) has a misaligned circumbinary debris disc that is almost perpendicular to the binary orbital plane (Kennedy et al. 2012).

The lifetimes of discs around single stars are observed to be around 1–|$10\, \rm Myr$| (Haisch, Lada & Lada 2001; Hernández et al. 2007, 2008; Mamajek 2009; Ribas, Bouy & Merín 2015). Mass accretion rates through circumbinary discs may be inhibited due to the tidal torques exerted by the binary, resulting in extended disc lifetimes (e.g. Alexander 2012). There is observational evidence for extended disc lifetimes for circumbinary discs. For example, the circumbinary gas discs HD 98800 B, V4046 Sgr, and AK Sco have disc ages of |$10\pm 3$|⁠, |$23\pm 3$|⁠, and |$18\pm 1\, \rm Myr$|⁠, respectively (Soderblom et al. 1998; Mamajek & Bell 2014; Czekala et al. 2015).

The lifetime of protoplanetary discs is fundamentally linked to their dispersal mechanisms. The main processes that remove mass and/or angular momentum from the disc include viscous evolution of the disc (e.g. Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974; Hartmann et al. 1998), photoevaporation by stellar radiation (e.g. Shu, Johnstone & Hollenbach 1993; Hollenbach et al. 1994), magnetically launched jets and winds (e.g. Königl & Salmeron 2011), and the interaction with newborn planets (e.g. Kley & Nelson 2012). However, the |$\sim \!\rm Myr$| time-scale for planet formation (Pollack et al. 1996) suggests that this is not a major factor for rapid disc dispersal required by observations (Alexander et al. 2014). It has also been shown that planets likely account for |$\lesssim 1{{\ \rm per\ cent}}$| of the initial disc mass budget (Mayor et al. 2011; Wright et al. 2011)

After the gaseous protoplanetary disc is dispersed, the remnant planetesimals produce a second generation of dust through collisions which leads to the formation of a gas-poor, less massive disc called a debris disc. These debris discs are much cooler in temperature and are analogous to the Solar system Kuiper belt (e.g. Wyatt 2008; Hughes, Duchêne & Matthews 2018).

The larger solids in a debris disc can be modelled as set of particles on nearly circular ballistic circumbinary orbits. The particles are not interacting except during close encounters or collisions. If such a disc is initially inclined with respect to the binary by some arbitrary tilt angle, the orbits will undergo differential nodal precession. As a consequence, the disc will not maintain its flat form and initially nearby orbits may undergo violent collisions (e.g. Nesvold et al. 2016). A low mass polar (or coplanar) debris disc is an exception to this rule because it does not undergo nodal precession.

Kennedy et al. (2012) suggested that the polar aligned debris disc in 99 Her could be the result of the capture of material or a stellar exchange which leaves the circumbinary debris in a polar orbit. Although such a scenario is possible, it involves fine-tuning of the conditions that result in a polar state, as they point out. Martin & Lubow (2017) instead suggested that the 99 Her debris disc began as a somewhat inclined gas disc with embedded solids that later evolved to a polar configuration. Since gas discs can evolve to a polar configuration, small embedded solids in the disc should follow the gas to a polar configuration. During this tilt evolution, dust (and perhaps somewhat larger solid bodies) can be well coupled dynamically to the gas that evolves as a nearly flat disc due to gas pressure communication. The disc tilt evolves as a consequence of viscous dissipation in the gas. Once in the polar configuration, the disc evolved to become a debris disc in the usual manner. This process might then operate over a wide range of initial conditions without fine-tuning. The purpose of this paper is to explore the viability of this scenario.

Circumbinary dust particles experience various degrees of coupling to the gas depending on their Stokes number (e.g. Birnstiel, Dullemond & Brauer 2010). Over time, small, initially well coupled, dust grains grow to higher Stokes number and gradually decouple from the gas disc. If significant decoupling occurs in a misaligned circumbinary disc, the dust particle orbits may evolve independently of the gas disc (e.g. Aly & Lodato 2020). Thus the gas disc around 99 Her must evolve to a polar configuration before the dust decouples because a polar debris ring is observed in 99 Her. According to the dust coagulation model, dust grains remain coupled to the gas with |$\rm St \lt 1$| (Birnstiel, Klahr & Ercolano 2012). Large decimetre-sized dust grains have large relative velocities which induce fragmentation rather than sticking during collisions. This fragmentation barrier limits dust growth which in-turn prevents the dust from decoupling from the gas (Blum & Wurm 2008; Brauer, Dullemond & Henning 2008; Birnstiel et al. 2010). Nevertheless, planet formation may still proceed when the dust particles are concentrated to high density and gravitationally collapse into planetesimals, e.g. by the streaming instability (e.g. Johansen, Youdin & Lithwick 2012; Yang, Johansen & Carrera 2017). Therefore, if the fragmentation barrier dominates in the disc of 99 Her, a significant fraction of the solid materials may still remain coupled to the gas while the gas disc evolves to a polar state.

The structure of this paper is as follows. In Section 2, we discuss the orbital parameters of 99 Her and the observational evidence for a polar aligned debris disc. In Section 3, we show the orbital evolution of a circumbinary particle in the potential of the binary system 99 Her. In Section 4, we present the set-up for the SPH simulations and discuss the results. In Section 5, we apply linear theory for the disc evolution. Finally, we discuss the possible implications of this work in Section 6 and we draw our conclusions in Section 7.

2 OBSERVED PROPERTIES OF 99 HERCULIS

99 Her is a particularly interesting system since it hosts a resolved polar debris disc and the binary orbit is well characterized. The binary has been observed since the latter half of the 1800’s (e.g. Burnham 1878; Flammarion 1879; Gore 1890). It consists of an F7V primary and a K4V secondary. The primary star has an estimated age of 6–10 Gyr, consistent with being on the main sequence (Nordström et al. 2004; Takeda 2007). There has been an abundance of observations of this system in recent years. Kennedy et al. (2012) performed a more precise derivation of the system orbital parameters. The orbit was found by fitting position angles (PAs) and separations from the Washington Double Star (WDS) Catalog (Mason et al. 2018). The semimajor axis, eccentricity, inclination, and orbital period are |$a=16.5\, \rm au$|⁠, e_b = 0.766, i = 39°, and |$P_{\rm orb}=56.3\, \rm yr$|⁠, respectively. The longitude of the ascending node and longitude of pericentre are Ω = 41° and ω = 116°, respectively. The longitude of pericentre is measured anticlockwise from the ascending node, and projected on to the sky plane which has a PA of 163°. The reason that the sum of Ω and ω do not equal the sky plane PA is due to the fact that the binary orbit is inclined. The total mass of the binary inferred from the data is |$M=1.4 \, \rm M_{\odot }$| using a distance of |$15.64\, \rm pc$| (van Leeuwen 2008). Using the spectroscopic mass function, a mass ratio of 0.49 was inferred. Therefore, the primary mass is |$M_1=0.94\, \rm M_{\odot }$| while the mass of the secondary is |$M_2=0.46 \, \rm M_{\odot }$|⁠.

Debris discs around binary systems are as common as around single star systems (Trilling et al. 2007). The debris around 99 Her was first detected using 100 and |$160\, \rm \mu m$| data from Herschel Photodetector and Array Camera and Spectrometer (PACS; Griffin et al. 2010; Poglitsch et al. 2010). The debris was first detected but not resolved by SPIRE at 250 and |$350\, \rm \mu m$|⁠. Kennedy et al. (2012) provided resolved PACS images of the disc and were able to estimate the debris structure, inclination, and PA using two-dimensional Gaussian models. To best fit the observations, they invoked a debris ring model located in a small range of radii near |$120 \, \rm au$|⁠. The grain properties and size distributions cannot accurately be constrained due to the low number of measured disc emission wavelengths.

The projection of the binary pericentre direction on the plane of the sky has a PA of 163° ± 2°, and a line perpendicular to this has a PA of 73° ± 2°. Since the observed debris disc has a PA of 72°, Kennedy et al. (2012) concluded that this is consistent with the disc being at 87° with respect to the binary pericentre direction, or 3° away from polar alignment. However, they found that the ring could be mirrored in the sky plane causing the observed misalignment to be 30°.

Kennedy et al. (2012) used circumbinary test particles around the eccentric binary 99 Her to model each of these possible ring misalignments. The test particles undergo secular perturbations which lead to nodal precession, as well as inclination oscillations. A debris ring with a tilt of 30° will eventually spread out into a broader, non-ring like structure due to the effects of differential nodal precession. By assuming the largest fragments are at least |$1\, \rm mm$|⁠, the secular precession period is estimated to be |$0.5\, \rm Myr$| (Kennedy et al. 2012). By their calculations, within 10 precession periods (⁠|$5\, \rm Myr$|⁠), the ring structure will be erased. Thus, the model that best fits the PACS images with considerable agreement is a ring in a nearly polar configuration.

The polar particle orbits in their model are stable over the stellar lifetime, which means that the observed dust could be the steady-state collision products of the polar planetesimal ring. Though a thin ring of debris best fit their data, the initial extent of a gas disc must have been at least as large as this. Radial drift can reduce the outer radius of the debris disc due to the gas drag force (Adachi, Hayashi & Nakazawa 1976; Weidenschilling 1977; Takeuchi & Lin 2002). This process can lead to the outer radius of the dust being less than the outer radius of the gas.

3 INCLINED CIRCUMBINARY PARTICLE ORBITS AROUND 99 HER

We first investigate the evolution of inclined circumbinary particle orbits in the potential of the binary system 99 Her, following Chen et al. (2019). For a circular orbit binary, the orbital angular momentum of a circumbinary test (massless) particle always precesses about the binary angular momentum vector. We adopt a Cartesian coordinate system (x, y, z), where the x-axis is along the initial binary eccentricity vector and the z-axis is along the initial binary angular momentum vector. The origin of the coordinate system is taken to be the initial centre of mass of the entire system. An eccentric orbit binary generates a secular potential that is non-axisymmetric about the z-axis. For test particle simulations, the binary orbit is fixed. However, for a massive particle the binary eccentricity vector precesses. In addition, the particle angular momentum vector oscillates.

We consider nearly circular circumbinary particle orbits with an initial semimajor axis of |$d = 82.5\, \rm au$| (five times the binary semimajor axis) with inclination i₀ with respect to the binary orbital plane. We consider test particles with various inclinations that have stable orbits. We define the precession angle, ϕ_pb, of the particle angular momentum vector relative to the binary eccentricity vector as

$$\begin{eqnarray*} \phi _{\rm pb} = \phi - \phi _{\rm b}. \end{eqnarray*}$$

(1)

The longitude of the ascending node of the particle relative to the x-axis (the initial binary eccentricity vector) is given by

$$\begin{eqnarray*} \phi = \tan ^{-1} \bigg (\frac{l_{\rm py}}{l_{\rm px}} \bigg) + \pi /2, \end{eqnarray*}$$

(2)

where l_px, l_py are the x and y components of the particle angular momentum vector l_p. The azimuthal angle of the eccentricity vector is defined as

$$\begin{eqnarray*} \phi _{\rm b} = \tan ^{-1} \bigg (\frac{e_{\rm by}}{e_{\rm bx}} \bigg) + \pi /2, \end{eqnarray*}$$

(3)

where e_bx, e_by are the x and y components of the eccentricity vector e_b. The initial precession angle of the particle is ϕ_pb = ϕ₀ = 90° in all of the three body simulations we consider. We define the angular momentum ratio of the particle to the binary as

$$\begin{eqnarray*} j=\frac{J_{\rm p}}{J_{\rm b}}. \end{eqnarray*}$$

(4)

We consider three different values of j = 0, 0.1, and 0.51 by varying the mass of the particle to 0, 0.006 45, and 0.034 09 M_⊙, while keeping the particle separation fixed.

The upper panel of Fig. 1 shows the circumbinary massless test particle orbits in the icos ϕ_pb–isin ϕ_pb phase space, where i is the inclination between the angular momentum of the particle and the binary. For a massless (test) particle, we have j = 0. For the eccentric binary orbit, the orbit of a particle with sufficiently low-inclination precesses about the binary angular momentum (this is the circulating phase). However, the evolution is different to the circular orbit binary case because the particle orbit displays tilt oscillations during this process (Smallwood et al. 2019). If the test particle starts with a higher initial inclination, its orbit instead precesses about the eccentricity vector of the binary. In this configuration, the longitude of the ascending node of the particle orbit undergoes oscillations about ϕ_pb = 90° (or ϕ_pb = 270°) while the tilt undergoes oscillations (this is the librating regime) (Verrier & Evans 2009; Farago & Laskar 2010; Doolin & Blundell 2011). The librating orbits are identified by the red curves while the green curves represent prograde circulating orbits and the blue curves denote retrograde circulating orbits. For a test particle (j = 0), the critical inclination angle between the circulating and librating orbits is given by

$$\begin{eqnarray*} i_{\rm crit} = \sin ^{-1}\sqrt{\frac{1-e_{\rm b}^2}{1+4e_{\rm b}^2}} \end{eqnarray*}$$

(5)

(e.g. Farago & Laskar 2010). For the eccentricity of 99 Her (e_b = 0.766), the critical angle is i_crit = 20.6°, in agreement with the upper panel of Fig. 1.

Figure 1.

The icos ϕ_pb – isin ϕ_pb phase portrait for an initially circular particle orbits around the 99 Her binary system for different initial inclination, i, and initial precession angle, ϕ_pb = ϕ₀ = 90°. The initial separation of the test particle is |$d = 82.5\, \rm au$| in all panels. The green lines show prograde circulating solutions. The red lines show librating solutions that have initial inclination i < i_s, the stationary inclination. The cyan lines show librating solutions that have initial inclination i > i_s and the blue lines show the retrograde circulating solutions. The black arrows represent the initial position of the particle in the phase portrait and its orbital direction. Upper panel: test particle orbits with j = 0.0. Middle panel: particle orbits with j = 0.1. Lower panel: particle orbits with j = 0.51.

We define the ‘stationary inclination’, i_s, as the inclination corresponding to the centre of the librating region. We show librating orbits that begin with i < i_s in red, and those that begin with i > i_s in cyan. For the test particle orbits, the binary orbit is fixed and so the orbit is the same no matter where it begins within the librating region. Thus, in this case the red and cyan lines are exactly the same.

The middle panel of Fig. 1 shows particle orbits for j = 0.1 and the bottom panel shows particle orbits for j = 0.51. As the angular momentum ratio of the particle to the binary increases, the critical inclination between prograde circulating and librating orbits is higher. For j > 0, we apply equations (31) and (37) of Martin & Lubow (2019) to obtain for j = 0.1 that i_crit = 24.6° and for j = 0.51 that i_crit = 40.7°. These values agree with the three-body simulations shown in the figure.

For the higher mass particle, the binary orbit eccentricity oscillates and precesses. Consider a particle in a librating orbit that begins at ϕ₀ = 90° with i < i_s. When it reaches ϕ = 90° again, but with i > i_s, the eccentricity of the binary has decreased. This means that the red and the cyan lines are no longer the same for a massive third body. We do not plot any circulating retrograde orbits for the high-mass case because they are unstable (Chen et al. 2019).

While these three-body simulations are on stable repeating orbits, the presence of viscosity in a circumbinary accretion disc will act to damp its tilt oscillations. Thus, the final inclination of the disc is either aligned (or counter aligned) to the binary orbit or polar aligned. Thus, the inclination at the centre of the librating region is an important parameter for the disc simulations. For the test particle case, j = 0, the centre of the librating region corresponds to i_s = 90° and ϕ_s = 90°. For the particle cases with j = 0.1 and j = 0.51, the centre of the librating regions corresponds to i_s = 87.9° and i_s = 80.4° (calculated from equation 15 in Martin & Lubow 2019) and ϕ_s = 90°, respectively, in agreement with Fig. 1.

4 CIRCUMBINARY DISC SIMULATIONS

The time-scale for polar alignment may be shorter or longer than the lifetime of the gas disc, depending on the binary and disc parameters (Martin & Lubow 2018). Therefore, we explore the parameter space for the polar alignment of a circumbinary gas disc in 99 Her by varying the mass, inclination, and outer radius of the disc. The gas disc needs to evolve to a polar configuration within its lifetime for a polar debris disc to remain after the gas disc is dispersed.

4.1 Simulation set-up

We use the three-dimensional smoothed particle hydrodynamics (SPH; e.g. Price 2012) code phantom (Lodato & Price 2010; Price & Federrath 2010; Price et al. 2018). phantom has been well tested and used to model misaligned accretion discs in binary systems (e.g. Nixon, King & Price 2013; Martin et al. 2014; Franchini, Martin & Lubow 2019a). The disc simulations are in the so-called bending waves regime where the disc aspect ratio H/R is larger than the viscosity coefficient α (Shakura & Sunyaev 1973). In this regime, the warp induced in the disc by the binary torque propagates as a pressure wave with speed c_s/2 (Papaloizou & Pringle 1983; Papaloizou & Lin 1995).

Each simulation consists of N equal mass particles initially distributed from the inner disc radius, R_in, to the outer disc radius, R_out. The inner disc radius is chosen to be |$R_{\rm in}=2a=33\, \rm au$|⁠, which is close to the radius where tidal torque truncation is important (Artymowicz & Lubow 1994). However, for a misaligned disc the tidal torque produced by the binary is much weaker, allowing the disc to survive closer to the binary (e.g. Lubow, Martin & Nixon 2015; Miranda & Lai 2015; Nixon & Lubow 2015; Lubow & Martin 2018). For simulations with outer disc radius |$R_{\rm out} = 120\, \rm au$|⁠, we take N = 300 000 particles, and with |$R_{\rm out} = 200\, \rm au$| we take N = 500 000. The binary begins at apastron with an eccentricity of 0.766. The accretion radius of each binary component is |$R_{\rm acc}=4\, \rm au$|⁠. Particles within this radius are accreted and their mass and angular momentum are added to the star. We ignore the effect of self-gravity since it has no effect on the nodal precession rate of flat circumbinary discs. For the narrow discs, the simulation lifetime is |$1000\, \rm P_{orb}$| or |$56\,000\, \rm yr$|⁠, and for the extended discs it is |$1500\, \rm P_{orb}$| or |$90\,000\, \rm yr$|⁠.

We chose to model the physical disc viscosity by using the artificial viscosity α^AV, implemented in phantom (Lodato & Price 2010). The surface density profile of the disc is initially set as a power-law distribution Σ ∝ R^−3/2. The disc is locally isothermal with sound speed c_s ∝ R^−3/4 and H/R = 0.1 at R = R_in. With this prescription, the disc is uniformly resolved meaning that 〈h〉/H and therefore the disc viscosity parameter α is constant over the radial extent of the disc (Lodato & Pringle 2007). We take the Shakura & Sunyaev (1973) α_SS to be 0.01. There exists a lower limit for the artificial viscosity in this type of simulation below which a physical viscosity is not resolved: α^AV = 0.1.

In our simulations, the disc is resolved with shell-averaged smoothing length per scale height 〈h〉/H ≈ 0.29, independent of the disc size since we increase the number of particles for the larger disc radius. We show the initial parameters of the circumbinary disc for all simulations in Table 1. We analyse the data from the SPH simulations by dividing the disc into 300 bins in spherical radius, R, which range from the radius of the innermost particle to |$170\, \rm au$| for narrow discs and |$250\, \rm au$| for extended discs. Within each bin we calculate the mean properties of the particles such as the surface density, inclination, longitude of ascending node, and eccentricity. The longitude of the ascending node of the disc is calculated by using the prescription from equation (2). In the next Section, we describe how varying the initial tilt, the initial disc mass, and the initial disc size affect the evolution.

Table 1.

Parameters of the initial circumbinary disc for all simulations. The first column is the model name. The second column is the figure number in which the model appears. The third column is the initial number of particles. The fourth column is the initial disc tilt. The fifth column is the initial outer radius of the disc. The sixth column is the initial disc mass. The seventh column is the initial ratio of the angular momentum of the disc to the angular momentum of the binary. The eighth column describes whether the disc is undergoing circulation (C) or libration (L). The ninth column in the tilt decay time-scale calculated from equation (6). The tenth column is the approximate polar alignment time-scale. The eleventh column describes whether the disc breaks. Finally, the twelfth column denotes the final disc inclination.

Model	Fig.	N	Tilt	Rout	Mdisc	j0	C/L	τ	tpolar	Disc breaking	ifinal
			(°)	|$(\rm au)$|	(M⊙)			(yr)	(yr)		(°)
A	–	300 000	20	120	0.001	0.01	C	–	–	No	–
B	2 (upper)	300 000	30	120	0.001	0.01	L	11 500	49 500	Yes	89
C	2 (middle)	300 000	40	120	0.001	0.01	L	11 700	30 000	Yes	89
D	2 (lower), 3	300 000	60	120	0.001	0.01	L	17 000	26 000	No	89
E	4	300 000	40	120	0.01	0.10	L	15 000	36 000	Yes	87
F	–	300 000	40	120	0.05	0.51	C	–	–	No	–
G	5 (upper)	300 000	60	120	0.01	0.10	L	18 000	31 000	No	87
H	5 (lower)	300 000	60	120	0.05	0.51	L	67 000	–	No	74.5
I	7 (upper), 8	500 000	60	200	0.001	0.012	L	–	67 000	No	89
J	7 (lower)	500 000	60	200	0.01	0.12	L	–	45 000	No	85

4.2 Alignment time-scale calculations

In this section, we describe the process for calculating the tilt decay time-scale and for estimating the polar alignment time-scale from simulations. As a disc evolves to a polar configuration, it undergoes librations of the longitude of the ascending node and tilt. Dissipation causes tilt oscillations that damp towards a polar configuration. The tilt decay time-scale, τ, to the polar state (that may correspond to an increase in tilt relative to the binary) is the time of exponential decay of the damped oscillations. The polar alignment time-scale, t_polar, is the time at which the disc is nearly fully orientated in a polar fashion. Several tilt decay time-scales can be used as an estimate for the polar alignment time-scale (Martin & Lubow 2018).

We determine the times of the first two local maxima of the inclination of the disc relative to the binary orbital plane to be t₁ and t₂. The tilt decay time-scale towards the polar orientation is given by

$$\begin{eqnarray*} \tau = \frac{t_2 - t_1}{\log [(i_1-i_{\rm final})/(i_2-i_{\rm final})]} , \end{eqnarray*}$$

(6)

(Martin & Lubow 2018), where i₁ and i₂ are the local maxima of inclination at t₁ and t₂, respectively, and i_final is the final inclination of the disc. As can be seen in the following sections, some of the simulations exhibit strongly warped discs. Consequently, we a use the density-weighted average of the disc tilt. In all simulations, the resolution decreases with time. Thus, the tilt decay time-scale in this manner needs to be calculated from early on in the simulations.

Where possible, we define the polar alignment time-scale, t_polar, of the disc to be the time after which the magnitude of the difference in the density-weighted average of the longitude of the ascending node of the disc and the azimuthal angle of the binary eccentricity vector does not vary by more than 1°. Estimating the polar alignment time-scale from the inclination of the disc is difficult because the angular momentum of the disc evolves in time. The stationary inclination, i_s, increases as the disc loses mass (Martin & Lubow 2019). The tilt decay time-scale and the polar alignment time-scale from all simulation models are listed in the ninth and tenth columns of Table 1, respectively.

4.3 Effect of initial disc tilt for a low-mass disc

We first vary the initial tilt of the disc while keeping the initial disc mass and initial disc size constant (models A, B, C, and D in Table 1). The initial outer radius of the disc is |$R_{\rm out}=120\, \rm au$| and the initial disc mass is |$0.001\, \rm \mathrm{M}_\odot$|⁠. Note that this initial value of R_out is the same as the ring radius obtained by Kennedy et al. (2012). The disc outer radius increases as the disc evolves. Consequently, the disc could produce debris at |$120\, \rm au$|⁠. We consider the effects of larger initial disc outer radii in Section 4.5. The critical angle between librating and circulating solutions for a circumbinary test particle around 99 Her found in Section 3 is i = 20.57°. In these simulations, the disc mass is very low (⁠|$M_{\rm d}=0.001\, \mathrm{ M}_{\odot }$|⁠) so the critical inclination angle that separates librating and circulating solutions is consistent with the test particle calculations. The disc with initial misalignment of 20° (model A) aligns to the binary orbital plane rather than evolving towards a polar configuration. We do not show this result in the paper since we are interested in investigating discs that evolve to polar alignment.

Fig. 2 shows the time evolution of the disc tilt, i, and the longitude of the ascending node, ϕ, at two different radii, |$50$| (dashed) and |$120\, \rm au$| (dotted) inside the disc for the three simulations that go to polar alignment. We also show the density weighted average over the disc for both i and ϕ, given by the sold lines. The initial misalignment is i₀ = 30° (model B), 40° (model C), and 60° (model D) for the top, middle, and bottom panel, respectively. The blue line represents the azimuthal angle of the eccentricity vector given by equation (3).

Figure 2.

Evolution of the inclination, i, and longitude of the ascending node, ϕ, both as a function of time with varying initial disc tilt, i₀, evaluated at two radii within the disc, 50 (dashed) and |$120\, \rm au$| (dotted). The solid lines represent the density weighted averages over the entire disc. The azimuthal angle of the eccentricity vector of the binary is shown by the blue line. The vertical grey line marks the polar alignment time-scale. Top panel: i₀ = 30°, (model B from Table 1). Middle panel: i₀ = 40° (model C). Bottom panel: i₀ = 60° (model D).

The tilt decay time-scales, calculated from equation (6), are |$11\, 500$|⁠, |$11\, 700$|⁠, and |$17\, 000\, \rm yr$| for discs with initial tilts of i₀ = 30°, 40°, and 60°, respectively. For all three of these models, the final inclination, i_final, is 89°. The tilt decay time-scale decreases when the initial disc tilt becomes farther from the stationary inclination because the disc becomes more warped, leading to stronger dissipation and a faster decay time-scale. The polar alignment time-scale for the three models are |$\sim \!49\, 500$|⁠, |$\sim \!30\, 000$|⁠, and |$\sim \!26\, 000\, \rm yr$|⁠, respectively. The estimated polar alignment time-scales are a few tilt decay time-scales. The polar alignment time-scale is shown by the grey vertical lines in Fig. 2. For the low-mass disc, if its initial tilt is closer to the stationary inclination angle, then its evolution to polar alignment occurs on a shorter time-scale. However, in all cases, the disc aligns to polar on a time-scale much shorter than the expected disc lifetime.

The evolution of a circumbinary disc whose initial tilt is close to the critical angle (i = 30° and i = 40°) is somewhat different to the case with i₀ = 60°. A protoplanetary disc precesses nearly as a solid body if the radial communication time-scale is shorter than the precession time-scale (Papaloizou & Terquem 1995; Larwood & Papaloizou 1997). For a disc with a lower initial inclination, the librating region is larger which causes the disc to break into disjointed rings (e.g. Nixon et al. 2013). A disc with an initial tilt much greater than the critical inclination and closer to polar does not undergo this breaking and smoothly transitions to polar.

Fig. 3 shows the initial (upper panel) and final (lower panel) disc–binary system for the disc that is initially tilted by 30° (model B). The disc centre initially lies at the centre of mass of the binary and the binary angular momentum vector initially lies along the z-axis and the binary eccentricity vector initially along the x-axis. At a time of |$56\, 000\, \rm yr$|⁠, the disc is in the polar configuration.

Figure 3.

Low-mass circumbinary disc with an initial tilt of 30° (model B). The binary components are shown by the red circles, with the primary positioned to the left and the secondary to the right. The size of the circle is the accretion radius of the sink. Upper panel: initial disc setup for the SPH simulation in the x–z plane. Lower panel: the disc in a polar configuration at a time of |$t = 56\,000\, \rm yr$|⁠. The colour denotes the gas density with yellow regions being about two orders of magnitude larger than the blue.

While we have not run any simulations with inclination that is closer to retrograde than prograde, the alignment time-scale of such a disc would be shorter than the prograde simulations shown here. A disc that is closer to retrograde feels a weaker binary torque and thus is able to extend closer to the binary (Miranda & Lai 2015; Nixon & Lubow 2015). The precession time-scale is shorter closer to the binary and therefore the alignment is faster (Martin & Lubow 2018; Cuello & Giuppone 2019).

4.4 Effect of the disc mass

We now additionally consider two higher values of initial disc mass, |$M_{\rm d}=0.01,\, 0.05\, \rm M_{\odot }$| for initial tilts of |$i_0=40^{\circ },\, 60 ^{\circ }$|⁠. The ratio of the angular momentum of the disc to the angular momentum of the binary is initially j = 0.01, 0.10, and 0.51 for the initial disc masses M_d = 0.001, 0.01, and |$0.05\, \mathrm{ M}_{ \odot}$|⁠, respectively.

Fig. 4 shows the evolution of the inclination and the longitude of the ascending node for an initial disc tilt of i₀ = 40° with an initial disc mass |$M_{\rm d}=0.01\, \rm M_{\odot }$| (model E), which corresponds to j = 0.1 initially. The critical inclination between librating and circulating solutions found from the three body problem in Section 3 is i_crit = 24.6°. Thus, this moderate mass disc evolves towards a polar configuration. Comparing to the low-mass disc simulation with the same initial tilt (middle panel of Fig. 2), the mass of the disc has not qualitatively changed the behaviour. The precession rate of the binary eccentricity vector is much higher for the higher mass disc as expected. However, the alignment time-scale does not change significantly compared to the low-mass disc.

Figure 4.

Same as Fig. 2 but for a circumbinary disc with an initial tilt of i₀ = 40° and initial disc mass of |$M_{\rm d}=0.01\, \rm M_{\odot }$| (model E).

We also considered a higher disc mass of |$M_{\rm d}=0.05\, \rm M_{\odot }$| for a disc inclined by 40° (model F). This disc has angular momentum ratio of the disc to the binary of j = 0.51 initially. In this case, the disc does not evolve to polar alignment but instead aligns to the binary orbital plane. For a particle with angular momentum ratio 0.51, we can estimate the critical inclination between librating and oscillating solutions to be i_crit = 43.5° (with equation 37 in Martin & Lubow 2019). Thus, the simulations agree with the three body problem solutions in Fig. 1. The higher the disc mass, the more likely it is that the disc will evolve to alignment with the binary orbital plane rather than polar alignment.

Fig. 5 shows the evolution of the inclination and longitude of the ascending node for different initial disc masses for a disc with initial inclination of i₀ = 60°. Even for high-disc mass, providing that the initial inclination is sufficiently high, the disc still evolves to polar alignment. We estimate the stationary inclination angle for each simulation by taking the final inclination. The more massive the disc, the lower the stationary tilt inclination. This is consistent with our massive particle calculations in Section 3 where we show that the centre of the librating region, i.e. the stationary inclination angle, decreases as the mass of the particle increases. The stationary inclination angles inferred from our suite of SPH simulations are |$i_{\rm s}\approx 89^{\circ },\, 87^{\circ },\, 74.5^{\circ }$| for discs with initial mass |$M_{\rm d} = 0.001,\, 0.01,\, 0.05\, \rm M_{\odot }$|⁠, respectively. Since these disc masses correspond to j = 0.01, 0.10, and 0.51 initially, these correspond to the same angular momentum ratios as the particle orbits in Fig. 1. In comparison, the particle model described in Section 2 predicts i_s ≈ 90.0°, 87.9°, 80.4°. Consequently, the departures from polar alignment increase with mass in both the disc simulations and particle model cases. They are in rough quantitative agreement.

Figure 5.

Same as Fig. 2 but with a circumbinary disc with i₀ = 60°. Top panel: |$M_{\rm d}=0.01\, \rm M_{\odot }$| (model G). Bottom panel: |$M_{\rm d}=0.05\, \rm M_{\odot }$| (model H).

The discrepancy between the particle and disc results for the largest j value is likely due in part to the effects of accretion torques on the binary that are not accounted for in the particle model. The orbital evolution of the binary is dependent on both the accretion of angular momentum from the disc and the gravitational torques from the disc. Fig. 6 shows the time evolution of the binary eccentricity for the three circumbinary disc masses. The most massive disc may eventually cause the binary to circularize. There is however, limited resolution in the inner gap which leads to uncertainties in the accretional torque (Martin & Lubow 2019). A lower binary eccentricity leads to a lower stationary inclination, as seen in the simulations. For a constant eccentricity, the stationary angle is shown as a function of the angular momentum ratio in the lower panel of fig. 10 in Martin & Lubow (2019) for a very similar eccentricity to that of 99 Her. Thus, in the following section, where we discuss the effects of disc size, we omit the high-disc mass case due to the uncertain effects of accretional torques on the binary eccentricity evolution.

Figure 6.

The time evolution of the binary eccentricity, e_b, for three circumbinary disc masses: |$M_{\rm d} = 0.001$| (solid), |$M_{\rm d} = 0.01$| (dashed), and |$M_{\rm d} = 0.05\, \rm M_{\odot }$| (dotted).

The models with an initial disc mass of |$M_{\rm d}=0.01\, \rm M_{\odot }$| and initial tilts of |$i_0 = 40^{\circ }, 60^{\circ }$|⁠, have tilt decay time-scales of |$15\, 000$| and |$18\, 000\, \rm yr$|⁠, respectively. The polar alignment time-scale for the two models are estimated to be |$36\, 000$| and |$31\, 000\, \rm yr$|⁠, respectively. As expected, a massive disc that has an initial tilt closer to the stationary inclination angle aligns to that state on a shorter time-scale.

The models with a disc mass of |$M_{\rm d}=0.01\, \rm M_{\odot }$| and initial tilts of |$i_0 = 40^{\circ }, 60^{\circ }$| undergo polar alignment within |$\sim \!23\, 000$| and |$\sim \!12\, 000\, \rm yr$|⁠, respectively. A disc that has an initial tilt closer to the stationary inclination angle polar aligns on a faster time-scale. The disc with a mass of |$0.05\, \rm M_{\odot }$| and an initial tilt of 60° has not fully aligned polar within the simulation time. Therefore, we can only estimate the tilt decay time-scale, which is |$67\, 000\, \rm yr$|⁠. The large tilt decay time-scale implies that the polar alignment time-scale for this massive disc would be of the order of |$10^5\, \rm yr$|⁠. The comparison between discs with different masses (⁠|$M_{\rm d}=0.01, 0.05\, \rm M_{\odot }$|⁠) with the same initial inclination (60°) shows that a more massive disc aligns more slowly than a less massive disc. However, the alignment time-scale is still shorter than the expected lifetime of a disc.

4.5 Effect of the disc size

Finally, we investigate how the disc size affects the polar alignment process in the context of the binary 99 Her. The models that have been described so far have investigated the behaviour of a disc that extends from |$R_{\rm in}=33$| to |$R_{\rm out}= 120\, \rm au$|⁠. We now increase the disc outer radius to |$R_{\rm out} = 200\, \rm au$| and increase the simulation time to |$1500\, \rm P_{orb}$| or |$90\, 000\, \rm yr$|⁠. The initial disc tilt is set at i₀ = 60°. Because of uncertainties in the high-disc mass model described in the previous subsection, we consider the two different lower initial disc masses, |$M_{\rm d}=0.001$| and |$M_{\rm d}=0.01\, \rm M_{\odot }$| (see models I and J from Table 1). The initial ratio of the angular momentum of the disc to the initial angular momentum of the binary is 0.0122 and 0.122, respectively. We also increase the total initial number of equal mass particles to |$N=500\, 000$| in order to uniformly resolve the disc as in the previous sections with a shell-averaged smoothing length per scale height of 〈h〉/H ≈ 0.29.

Fig. 7 shows the evolution of the tilt and the longitude of ascending node. The top panel represents the extended disc model with a mass of |$M_{\rm d}=0.001\, \rm M_{\odot }$| and the bottom panel represents the model with a mass of |$M_{\rm d}=0.01\, \rm M_{\odot }$|⁠. The inner regions of the disc undergo tilt oscillations on a shorter time-scale compared to the outer regions, which slowly increase their tilt as the warp wave propagates outwards. A wider disc is less likely to precess as a rigid body because the radial communication time-scale may be longer than the precession time-scale, (see section 4.5 in Lubow & Martin 2018, for details). Therefore the disc can become warped or it can even break (e.g. Nixon & King 2012). Fig. 7 shows a warp in the inclination and a twist in the precession angle that moves outwards in time for both disc masses.

Figure 7.

Evolution of the inclination, i, and longitude of the ascending node, ϕ, both as a function of time for a circumbinary disc with an initial tilt of i₀ = 60° and an outer radius of |$r_{\rm out} = 200\, \rm au$|⁠. Measurements are evaluated at two radii within the disc, |$50$| (dashed) and |$200\, \rm au$| (dotted). The solid lines represent the density weighted averages over the entire disc. The azimuthal angle of the eccentricity vector of the binary is shown by the blue line. The vertical grey line marks the polar alignment time-scale. Top panel: disc mass of |$M_{\rm d}=0.001\, \rm M_{\odot }$| (model I). Bottom panel: disc mass of |$M_{\rm d}=0.01\, \rm M_{\odot }$| (model J).

Due to the intense warping in these larger discs, the tilt decay time-scale cannot be calculated because it is no longer close to the linear regime. Though, the polar alignment time-scale can still be estimated. The less massive disc (⁠|$M_{\rm d}=0.001\, \rm M_{\odot }$|⁠) polar aligns after |$t_{\rm polar} \sim 67\, 000\, \rm yr$|⁠, while the more massive disc (⁠|$M_{\rm d}=0.01\, \rm M_{\odot }$|⁠) has a polar alignment time-scale of |$t_{\rm polar} \sim 45\, 000\, \rm yr$|⁠. Therefore increasing the disc radial extent still results in a significantly shorter polar alignment time-scale compared to the expected lifetime of a gas disc.

To investigate the local behaviour of the wider accretion disc around 99 Her, in Fig. 8 we show the surface density (top panel), inclination (middle panel), and longitude of the ascending node (bottom panel) as a function of radius at |$t=0,\, 670, \, 6700, 67\, 000\, \rm yr$| for a disc with initial mass |$M_{\rm d}=0.001\, \rm M_{\odot }$|⁠. The initial surface density (at |$t = 0\, \rm yr$|⁠, identified by the solid line) is a power law Σ ∝ r^−3/2. As the disc evolves, the initial outer radius of the disc spreads outwards because of the presence of the disc viscosity. The inclination of the inner portions of the disc increase towards polar alignment as the wave travels outwards in time. Looking at the inclination of the disc as a function of radius at |$6700\, \rm yr$| (dotted line) in the middle panel, we see that the disc inner regions inside about |$100\, \rm au$| have larger misalignment compared to the outer regions. At the same time, the surface density profile shows a dip at around |$100\, \rm au$|⁠. The disc is strongly warped but a higher resolution is required in order to properly resolve if the disc actually breaks (Nixon et al. 2013). Eventually, after roughly |$t = 67\, 000\, \rm yr$|⁠, the whole disc aligns polar as shown by the dot–dashed line in the middle panel of Fig. 7.

Figure 8.

Surface density (top panel), tilt (middle panel), and longitude of the ascending node (bottom panel) as a function of radius at different times. The solid, dashed, dotted, and dash–dotted curves correspond to |$t = 0,\, 6.7 \times 10^2, \, 6.7 \times 10^3, \, 6.7 \times 10^4 \, \rm yr$|⁠, respectively. The initial conditions of the circumbinary disc are for model I, which is an extended disc with a disc mass |$M_{\rm d}=0.001\, \rm M_{\odot }$|⁠.

5 ANALYTIC ESTIMATES

In this section, we first compare the alignment time-scale calculated from linear theory to the results of the hydrodynamical simulations. We then use analytic results from the three body problem to constrain the mass of the gas disc at the time of disc dispersal based on the observed debris disc inclination.

5.1 Alignment time-scale

The main purpose of this work is to understand whether a primordial circumbinary accretion disc in 99 Her can undergo polar alignment within its lifetime, under a wide range of initial conditions, in order to explain the observations showing a polar debris disc around the binary star (Kennedy et al. 2012). We can see from the results of the simulations presented in the previous sections that the time-scale on which the disc evolves to polar alignment is of the order of a few tens of thousands of years. The mechanism that leads to alignment between the accretion disc and the binary angular momenta is the natural presence of the disc viscosity that acts to damp the inclination oscillations on a time-scale ∝ α⁻¹ dissipating the warp (Papaloizou & Terquem 1995; Lubow & Ogilvie 2000; Lubow & Martin 2018). The linear theory of warped discs describes the evolution of a disc that remains nearly flat and it has recently been applied to discs that are close to a polar aligned state (Lubow & Martin 2018; Zanazzi & Lai 2018).

The warp is dissipated and the disc aligns on a time-scale

$$\begin{eqnarray*} t_{\rm align} = \frac{1}{\alpha } \left(\frac{H}{R}\right)^2 \frac{\Omega _{\rm b}}{\Omega _{\rm p}^2} \end{eqnarray*}$$

(7)

(e.g. Lubow & Martin 2018), where |$\Omega _{\rm b}=\sqrt{G(M_1+M_2)/a^3}$|⁠. For a disc that is aligning to polar, the global disc precession frequency is given by

$$\begin{eqnarray*} \Omega _{\rm p}=\frac{3\sqrt{5}}{4}e_{\rm b}\sqrt{1+4e_{\rm b}^2}\frac{M_1M_2}{M^2}\left\langle \left(\frac{a}{R}\right)^{7/2}\right\rangle \end{eqnarray*}$$

(8)

(see equation 16 in Lubow & Martin 2018), where

$$\begin{eqnarray*} \left\langle \left(\frac{a}{R}\right)^{7/2}\right\rangle =\frac{\int _{R_{\rm in}}^{R_{\rm out}}\Sigma R^3 \Omega (a/R)^{7/2}\, \mathrm{ d}R }{\int _{R_{\rm in}}^{R_{\rm out}} \Sigma R^3 \Omega \, \mathrm{ d}r}, \end{eqnarray*}$$

(9)

where |$\Omega =\sqrt{G(M_1+M_2)/R^3}$| is the Keplerian angular frequency. Equation (7) is a rough estimate because it neglects the variation of quantities such as α and disc aspect ratio H/R with radius R. Precise predictions of linear theory are made by computing the complex eigenfrequency for the linear tilt evolution equations (e.g. Lubow & Martin 2018). But we do not carry out such calculations here.

The linear model we apply does not account for the evolution of the density distribution, as discussed in Martin & Lubow (2019). It takes as input, the disc surface density profile that is approximated as being fixed in time. Since the density evolves over time, the profiles used in the calculation are taken to be representative at some intermediate time. We first take the surface density profile from the SPH simulations and calculate the alignment time-scale. Fig. 9 shows the surface density profile taken from model D at a time of |$t=10\, 000\, \rm yr$| (black line) and from model I at a time of |$t=30\, 000\, \rm yr$|⁠. Fig. 10 shows the alignment time-scale as a function of disc aspect ratio for the two surface density profiles with α = 0.01 and varying disc aspect ratio. The alignment time-scale for the disc initially truncated at |$120\, \rm au$| (model D with H/R ≈ 0.07 at r_out), |$t_{\rm polar} = 26\, 000\, \rm yr$|⁠, is in rough agreement with the alignment time-scale predicted by the black line in Fig. 10, which is |$\sim \!50\, 000\, \rm yr$|⁠. For the larger initial truncation radius of |$200\, \rm au$| (model I with H/R ≈ 0.06 at r_out), the alignment time-scale predicted here, |$\sim \!140\, 000\, \rm yr$|⁠, is reasonable compared with |$t_{\rm polar} = 67\, 000\, \rm yr$| from the simulation. For both initial outer disc radii, the linear theory is off by a factor of two compared with simulations.

Figure 9.

The surface density profile taken from model D at a time of |$t=10\, 000\, \rm yr$| (black line) and from model I at a time of |$t=30\, 000\, \rm yr$|⁠.

Figure 10.

Alignment time-scale calculated with equation (7) as a function of H/R for the surface density profiles taken from Fig. 9. The lower black solid line corresponds to the surface density profile taken from model D with initial truncation radius |$R_{\rm out}=120\, \rm au$|⁠, while the blue dashed line corresponds to the surface density profile taken from model I with initial truncation radius |$R_{\rm out}=200\, \rm au$|⁠.

Next, we consider a power-law surface density with Σ ∝ R^−3/2 extending from |$R_{\rm in}=1.6\, a$| (see Franchini, Lubow & Martin 2019b, for a discussion of the inner disc radius of a polar aligned disc) up to R_out. We take α = 0.01 and consider different values for H/R. The results are shown in Fig. 11 as a function of the disc outer radius. A radially wider accretion disc is expected to align on a longer time-scale compared to a more narrow disc with the same mass and initial inclination. This is consistent with the comparison between the results for model D (shown in the bottom panel of Fig. 2) and model I, shown in the upper panel of Fig. 7. Even if the gas disc is extended significantly father out than our simulations, the alignment time-scale is still shorter than the expected disc lifetime. However, this holds for α = 0.01 used in the simulations. For much smaller |$\alpha \lesssim10^{-4}$|⁠, equation (7) predicts that the polar alignment time-scale could become longer than the disc lifetime.

Figure 11.

Alignment time-scale calculated with equation (7) as a function of the disc outer radius for a power-law density profile Σ ∝ R^−3/2, |$R_{\rm in}=1.6\, a$|⁠, and α = 0.01. The solid black line has H/R = 0.05, the dashed blue line has H/R = 0.075 and the dotted red line has H/R = 0.1.

Dust grains are not included in our simulations but the time at which they decouple from the disc provides a more stringent constraint on the required polar alignment time-scale than the lifetime of the disc. The dust particles must be coupled to the gas disc until it reaches polar alignment. Otherwise, a misaligned disc of debris will undergo differentiated precession and lead to spherical distribution rather than a disc (e.g. Nesvold et al. 2016). Therefore, the time it takes the dust particles to grow and decouple from the gas must be longer than the polar alignment time-scale of the gas disc.

The strong disc warping can be seen in Fig. 12, where the disc tilt is given as a function of radius with <h > /H ≤ 0.5. the black line shows the narrow disc (model D) at a time of |$t = 5000\, \rm yr$| and the blue shows the wider disc (model I) for a time of |$t = 10\, 000\, \rm yr$|⁠. Strong non-linear disc warping is beyond the regime of applicability for linear theory. Strong warping may lead to additional dissipation and reduce the alignment time-scale to below the values we estimate.

Figure 12.

The inclination profile taken from model D at a time of |$t=5000\, \rm yr$| (black line) and from model I at a time of |$t=10\, 000\, \rm yr$| (blue line).

We see from Figs 2 and 5 that less massive discs tend to align polar on a slightly shorter time-scale, regardless of the initial misalignment while more massive discs tend to align on a longer time-scale. The linear theory we apply assumes that the disc mass/angular momentum is very small compared to the binary mass/angular momentum. However, in general we find the theoretical estimates to be consistent with the time-scale for the disc to evolve to polar alignment inferred from our SPH simulations.

5.2 Stationary inclination

Since the time-scale for polar alignment is shorter than the disc lifetime for a broad range of disc parameters, the debris disc should be at an inclination given by the generalized polar state. In the limit of a zero angular momentum disc, the generalized polar state is a disc aligned to the binary eccentricity vector and inclined to the binary angular momentum vector by 90°. A more massive disc has a generalized polar alignment with a lower level of misalignment.

We follow Martin & Lubow (2019) and calculate the maximum possible disc gas mass at the time of decoupling of gas and solids given that the observed debris disc is within 3° of being polar aligned. The decoupling may occur when the gas disc disperses. We assume that the disc has aligned to a stationary inclination that depends upon binary eccentricity and disc angular momentum given by

$$\begin{eqnarray*} \cos i_{\rm s} = \frac{ -\left(1+ 4e_{\rm b}^2\right)+ \sqrt{ \left(1+ 4e_{\rm b}^2 \right)^2+60 \left(1-e_{\rm b}^2\right) j^2} }{10 j} \end{eqnarray*}$$

(10)

(see equation 17 in Martin & Lubow 2019). This was derived using the secular theory of a circumbinary particle with mass (Farago & Laskar 2010). Since we know the binary eccentricity and we set i_s = 87°, we can find an upper limit on the angular momentum ratio of

$$\begin{eqnarray*} j=\frac{J_{\rm d}}{J_{\rm b}}=0.14, \end{eqnarray*}$$

(11)

where J_d is the angular momentum of the disc. We assume that the surface density of the disc is distributed as a power Σ ∝ R^−q between |$R_{\rm in}=1.6\, a=26.4\, \rm au$| and R_out. We parametrize the disc angular momentum with

$$\begin{eqnarray*} J_{\rm d} = k M_{\rm d} \, a^2 \Omega _{\rm b}, \end{eqnarray*}$$

(12)

where k depends upon the density profile. For |$R_{\rm out}=120\, \rm au$|⁠, k varies from 2.2 to 2.0 as q varies from 0 to 1.5. For |$R_{\rm out}=200\, \rm au$|⁠, k varies from 2.8 to 2.4 as q varies from 0 to 1.5. The angular momentum of the binary is

$$\begin{eqnarray*} J_{\rm b}= 0.14 M a_{\rm b}^2 \Omega _{\rm b}. \end{eqnarray*}$$

(13)

For |$R_{\rm out}=120\, \rm au$| and q = 1.5, with equation (11) we find the disc mass to be |$M_{\rm d} \approx 0.01 M=0.014\, \rm \mathrm{M}_\odot$|⁠. But since the 99 Her debris ring has a radius of about |$120\, \rm au$|⁠, we must consider larger discs. For |$R_{\rm out}=200\, \rm au$| and q = 1.5, we find the disc mass to be |$M_{\rm d} \approx 0.0085 M=0.012\, \rm \mathrm{M}_\odot$|⁠. If the debris disc has an inclination of less than 3° to polar, the disc mass may be lower than these estimates. These masses are in agreement with SPH simulations. Models E and G that begin with a disc mass of |$0.01\, M$| both have a final inclination of 87°. This disc mass is also quite reasonable for a protostellar disc. This mass represents the gas disc mass at the time that the gas and the dust-producing solids within the disc decouple. In principle, this imposes further constraints that we do not examine here.

6 DISCUSSION

Hundreds of debris discs around single stars have been detected to date in far-IR surveys. Since the dust has a lifetime shorter than the age of the host star these observations imply that the dust contained in the debris disc is likely not primordial, at least for main-sequence stars with age |$\gt 10\, {\rm Myr}$| (Wyatt 2018). The dust content is therefore thought to be replenished through planetesimal collisions. Therefore, in principle, debris discs provide information on the planet formation process showing the location and properties of the planetesimals and also their collision velocities.

The detected extrasolar debris discs have shown a variety of morphologies (Booth et al. 2013), including evidence of multiple components similar to the outer Solar system configuration (Chen et al. 2014; Kennedy & Wyatt 2014). This provides only circumstantial support to the idea of the presence of planets within debris disc systems since they can remove dust from certain regions of the disc.

All currently known imaged planets have been observed in systems with debris discs (Ricci et al. 2015; Bowler 2016; Marino et al. 2016). In principle, it is then possible to study imaged planets and discs in the same system while this is difficult in protoplanetary discs due to their high optical depths (Hughes et al. 2018). Several efforts are currently being made to find a possible correlation between debris discs properties and imaged planets. Since the planets and the debris disc typical locations are separated by tens of au, it is not clear how the planets influence the disc detectability.

There are currently 21 confirmed circumbinary planets, of which 10 discovered through transit within 9 binary systems (Doyle et al. 2011; Welsh et al. 2012; Orosz et al. 2012a, b; Kostov et al. 2013, 2014; Schwamb et al. 2013; Welsh et al. 2015; Kostov et al. 2016; Orosz et al. 2019). For these systems, the planet orbital inclination with respect to the binary orbital plane is very low (Kostov et al. 2014; Welsh et al. 2015). This is likely the result of the small orbital period of all the observed binaries leading to tidal circularization of the binary and therefore alignment of the disc (Czekala et al. 2019). We showed that a sufficiently misaligned and not too massive accretion disc around an eccentric binary can undergo polar alignment within its lifetime and therefore form a polar debris disc. This implies that if planets formed in this disc they would be on polar orbits around the binary. However, there are no confirmed circumpolar planets detected so far. Planets with polar orbits would be harder to detect than the nearly coplanar planets found by recurrent transits with Kepler. Polar planets may be detectable as non-recurring transits of the binary or eclipse timing variations of the binary (Zhang & Fabrycky 2019). The best-fitting model to reproduce the 99 Her observations is a polar debris ring (Kennedy et al. 2012). A polar ring could be produced by the presence of polar planets. It is known that shepherding planets can cause a debris disc to form ring-like structures (Rodigas, Malhotra & Hinz 2014). In addition, the debris disc may contain or evolve to contain polar planets.

7 CONCLUSIONS

We have explored parameters of a circumbinary gas disc around the eccentric binary system 99 Her that result in polar alignment of the disc. We investigated the effects of the initial disc inclination, mass, and size on the time-scale over which the disc aligns polar to the binary orbital plane. Since the eccentricity of 99 Her is very high (e_b = 0.766), the initial disc misalignment does not have to be very large for the disc to evolve to polar alignment, as suggested by Martin & Lubow (2017). The critical angle depends on the mass of the disc and is larger for a more massive disc (see Fig. 1).

Since the observed inclination of the debris disc relative to the binary is only a few degrees away from 90°, the initial mass of the gas disc can be constrained to be M_di ≲ 0.01. For low-mass discs, an initially mild inclination (i₀ > 20°) results in polar alignment, a more massive disc is able to align polar only if the initial misalignment is moderate (e.g. i₀ > 40°). Gas discs that are initially inclined by more than 60° undergo polar alignment regardless of the initial disc mass. Wider gas discs are still able to align polar even though they might warp and break if the sound crossing time-scale is longer than the precession time-scale induced by the binary torque.

All simulations that evolved to a polar state, did so on time-scales of order of tens of thousands of years. The polar alignment time-scale for these simulations is therefore much shorter than the lifetime of the gas disc, which has been estimated for circumstellar gas discs to be a few Myr (Haisch et al. 2001) or longer for circumbinary gas discs, as discussed in the Section 1. We note that the time at which the dust grains decouple from the gas disc may provide a more stringent constraint than the lifetime of the gas disc. Moreover, our simulations have a turbulent viscosity parameter α = 0.01, while the linear theory suggests that for significantly smaller disc viscosity |$\alpha \lesssim10^{-4}$| the alignment time-scale becomes longer than the disc lifetime. Mindful of these caveats, we conclude that the presence of a remnant polar debris disc around 99 Her can be explained by an initially inclined gas disc evolving into a polar configuration before dispersing the gas under a wide range of initial conditions and disc parameters. Such a disc may provide an environment for the formation of polar planets.

ACKNOWLEDGEMENTS

We thank the anonymous referee for helpful suggestions that positively impacted the work. We thank Daniel Price for providing the phantom code for SPH simulations and acknowledge the use of splash (Price 2007) for the rendering of the figures. Computer support was provided by UNLV’s National Supercomputing Center. We acknowledge support from NASA through grants NNX17AB96G and 80NSSC19K0443.