Alignment of a circumbinary disc around an eccentric binary with application to KH 15D

TL;DR

This study investigates the evolution and alignment of circumbinary discs around eccentric binaries using simulations and theory, with application to the KH 15D system, revealing tilt oscillations and inclination changes.

Contribution

It provides a detailed analysis of how eccentricity affects circumbinary disc alignment and applies these findings to interpret observations of KH 15D.

Findings

Eccentric binaries cause larger tilt oscillations in circumbinary discs.

Discs can increase their inclination significantly before aligning.

Simulation results agree well with linear theory.

Abstract

We analyse the evolution of a mildly inclined circumbinary disc that orbits an eccentric orbit binary by means of smoother particle hydrodynamic (SPH) simulations and linear theory. We show that the alignment process of an initially misaligned circumbinary disc around an eccentric orbit binary is significantly different than around a circular orbit binary and involves tilt oscillations. The more eccentric the binary, the larger the tilt oscillations and the longer it takes to damp these oscillations. A circumbinary disc that is only mildly inclined may increase its inclination by a factor of a few before it moves towards alignment. The results of the SPH simulations agree well with those of linear theory. We investigate the properties of the circumbinary disc/ring around KH 15D. We determine disc properties based on the observational constraints imposed by the changing binary brightness. We find that the inclination is currently at a local minimum and will increase substantially before setting to coplanarity. In addition, the nodal precession is currently near its most rapid rate. The recent observations that show a reappearance of Star B impose constraints on the thickness of the layer of obscuring material. Our results suggest that disc solids have undergone substantial inward drift and settling towards to disc midplane. For disc masses $\sim 0.001 M_\odot$, our model indicates that the level of disc turbulence is low $\alpha \ll 0.001$. Another possibility is that the disc/ring contains little gas.

Tables

Table 1. Table 1: Parameters of the initial circumbinary disc around an equal mass binary with total mass M 𝑀 M , and separation a 𝑎 a .

Binary and Disc Parameters	Symbol	Value
Mass of each binary component	$M_1/M=M_2/M$	$0.5$
Accretion radius of the masses	$r_{\mathrm{acc}}/a$	$0.25$
Initial disk mass	$M_{\mathrm{di}}/M$	$0.001$
Initial disk inner radius	$r_{\mathrm{in}}/a$	$2$
Disc viscosity parameter	$\alpha$	0.01
Disc aspect ratio	$H/r\left(r=r_{\mathrm{in}}\right)$	0.1

Table 2. Table 2: The setup of the SPH simulations which lists the eccentricity of the binary, e b subscript 𝑒 b e_{\rm b} , the initial tilt of the disc, i 0 subscript 𝑖 0 i_{0} , and the initial outer boundary of the disc, r out subscript 𝑟 out r_{\rm out} . We also list the critical inclination of a test particle derived from Equation ( 1 ). The initial tilts from each model are always below the critical to assure the disc aligns to the orbital binary plane.

Model	$e_b$	$i_0$	$i_{\mathrm{crit}}$	$r_{\mathrm{out}}/a$
Run1	$0.0$	$60°$	–	$5$
Run2	$0.3$	$50°$	$54.9°$	$5$
Run3	$0.6$	$30°$	$30.8°$	$5$
Run4	$0.8$	$15°$	$18.5°$	$5$
Run5	$0.8$	$15°$	$18.5°$	$40$
Run6	$0.3$	$10°$	$54.9°$	$5$
Run7	$0.6$	$10°$	$30.8°$	$5$
Run8	$0.8$	$10°$	$18.5°$	$5$

Table 3. Table 3: Model Parameters

Model	$r_{\mathrm{in}}/a$	$H/r\left(r_{\mathrm{in}}\right)$	$\alpha$	$p$	$s$	$e_b$
A	4	0.1	0.01	0.5	1.0	0.6
B	4	0.1	0.01	0.5	1.0	0.8
C	4	0.1	0.01	1.0	1.0	0.6
D	4	0.1	0.01	1.0	1.0	0.8
E	various	0.1	0.01	0.5	1.0	0.75
F	various	0.1	0.01	1.0	1.0	0.75

Table 4. Table 4: Period Constrained Results

Model	$r_{\mathrm{out}}/a$	$\tau$ (yr)	Max($r/{\ell }_{\mathrm{in}}\left|d\ell /dr\right|$)
A	27.6	$2.3\times {10}^5$	0.04
B	26.4	$2.8\times {10}^5$	0.06
C	37.3	$1.3\times {10}^5$	0.05
D	35.0	$1.6\times {10}^5$	0.07

Equations

$$i_{\mathrm{crit}}={\sin }^{-1}\sqrt{\frac{1-e_{\mathrm{b}}^2}{1+4e_{\mathrm{b}}^2}}$$

$$\alpha \approx \frac{{\alpha }_{\mathrm{AV}}}{10}\frac{⟨h⟩}{H},$$

$$t_{\rm c}\approx\frac{4}{(2+s)\Omega_{\rm b}h_{\rm out}}\bigg{(}\frac{r_{\rm out}}{a_{\rm b}}\bigg{)}^{3/2}$$

$$\omega_{n}(r)=k\bigg{(}\frac{a}{r}\bigg{)}^{7/2}\Omega_{\rm b}$$

$$k=-\frac{3}{4}\sqrt{1+3e_{\mathrm{b}}^2-4e_{\mathrm{b}}^4}\frac{M_1{M}_2}{M^2}.$$

$$t_{\mathrm{p},\mathrm{global}}=\frac{2(1+p)r_{\mathrm{in}}^{1+p}{r}_{\mathrm{out}}^{5/2-p}}{|k|(5-2p)a^{7/2}{\Omega }_{\mathrm{b}}},$$

$$\Sigma r^2\Omega \frac{\partial \mathbf{\ell }}{\partial t}=\frac{1}{r}\frac{\partial \mathbf{G}}{\partial r}+\mathbf{T},$$

$$\frac{\partial \mathbf{G}}{\partial t}-{\omega }_{\mathrm{a}}{\mathbf{e}}_z\times \mathbf{G}+\alpha \Omega \mathbf{G}=\frac{\mathcal{I}{r}^3{\Omega }^3}{4}\frac{\partial \mathbf{\ell }}{\partial r},$$

$${\omega }_{\mathrm{a}}(r)=\frac{3}{8}(2+3e_{\mathrm{b}}^2)\frac{M_1{M}_2}{M^2}{\left(\frac{a_{\mathrm{b}}}{r}\right)}^{7/2}{\Omega }_{\mathrm{b}}$$

$$\mathcal{I}=\int \rho z^{2}dz,$$

$$\mathbf{G}(r_{\mathrm{in}},t)=\mathbf{G}(r_{\mathrm{out}},t)=0.$$

$$\mathbf{T}=\Sigma r^2\Omega \mathbf{\tau }$$

$$\mathbf{\tau }=(a(r){\ell }_y,b(r){\ell }_x)$$

$$a(r)=(1-e_{\mathrm{b}}^2){\omega }_{\mathrm{p}}(r),$$

$$b(r)=-(1+4e_{\mathrm{b}}^2){\omega }_{\mathrm{p}}(r),$$

$${\omega }_{\mathrm{p}}(r)=\frac{3}{4}\frac{M_1{M}_2}{M^2}{\left(\frac{a_{\mathrm{b}}}{r}\right)}^{7/2}{\Omega }_{\mathrm{b}}.$$

$$i\omega \Sigma r^2\Omega \mathbf{\ell }=\frac{1}{r}\frac{d\mathbf{G}}{dr}+\Sigma r^2\Omega \mathbf{\tau }$$

$$i\omega \mathbf{G}-{\omega }_{\mathrm{a}}{\mathbf{e}}_z\times \mathbf{G}+\alpha \Omega \mathbf{G}=\frac{\mathcal{I}{r}^3{\Omega }^3}{4}\frac{d\mathbf{\ell }}{dr},$$

$$a$$

$$b$$

$$\mathbf{\ell }$$

$$\omega$$

$$\mathbf{G}$$

$$\mathbf{\tau }$$

$${\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}\Sigma r^3\Omega (i{\omega }^{(1)}{\mathbf{\ell }}^{(0)}-{\mathbf{\tau }}^{(1)})dr=0,$$

$${\mathbf{\tau }}^{(1)}(r)=(A^{(1)}(r){\ell }_y^{(0)},B^{(1)}(r){\ell }_x^{(0)}).$$

$${\omega }^{(1)}=-\frac{3}{4}\sqrt{1+3e_{\mathrm{b}}^2-4e_{\mathrm{b}}^4}\frac{M_1{M}_2}{M^2}⟨{\left(\frac{a_{\mathrm{b}}}{r}\right)}^{7/2}⟩{\Omega }_b,$$

$$⟨{\left(\frac{a_{\mathrm{b}}}{r}\right)}^{7/2}⟩=\frac{{\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}\Sigma r^3\Omega (a_{\mathrm{b}}/r{)}^{7/2}dr}{{\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}\Sigma r^3\Omega dr}.$$

$$P_{\mathrm{p}}=\frac{2\pi }{|{\omega }^{(1)}|}.$$

$${\ell }_y^{(0)}=-i\sqrt{\frac{1-e_{\mathrm{b}}^2}{1+4e_{\mathrm{b}}^4}}{\ell }_x^{(0)}.$$

$$\tan (\varphi (t))=-\frac{Re({\ell }_{\mathrm{x}}^{(0)}(t))}{Re({\ell }_{\mathrm{y}}^{(0)}(t))}.$$

$$i(t)=i_0\sqrt{\frac{2+(3-5\cos (2{\omega }^{(1)}t))e_{\mathrm{b}}^2}{2(1-e_{\mathrm{b}}^2)}}.$$

$$i_{\max }=i_0\sqrt{\frac{1+4e_{\mathrm{b}}^2}{1-e_{\mathrm{b}}^2}}.$$

Full text

††pubyear: 2019††pagerange: Alignment of a circumbinary disc around an eccentric binary with application to KH 15D–References

1 Introduction

Observations show that most stars form in relatively dense regions within stellar clusters which subsequently may be dispersed. The majority of these stars that form are members of binary star systems (Duquennoy & Mayor, 1991; Ghez et al., 1993; Duchêne & Kraus, 2013). The observed binary orbital eccentricities vary with binary orbital period (Raghavan et al., 2010; Tokovinin & Kiyaeva, 2016). For short binary orbital periods, typically less than about 10 days, the eccentricities are small, likely because the orbits are circularized by stellar tidal dissipation (Zahn, 1977). The average binary eccentricity increases as a function of binary orbital period and ranges from $0.39$ to $0.59$. In addition, there is considerable scatter in eccentricity at a given orbital period with high eccentricities $\sim 0.8$ or larger sometimes found.

Discs consisting of gas and dust likely reside within these systems at early stages. There can be multiple discs present in a binary system. A circumbinary disc orbits around the binary, while each of the binary components can be surrounded by its own disc (i.e. circumprimary and circumsecondary discs), as is found in binary GG Tau (Dutrey et al., 1994). Each of the discs may be misaligned to each other and to the binary.

Some circumbinary discs have been found to be misaligned with respect to the orbital plane of the central binary. For example, the pre-main sequence binary KH 15D has a circumbinary disc that is misaligned to the binary (Chiang & Murray-Clay, 2004; Winn et al., 2004). The circumbinary disc or ring around the binary protostar IRS 43 has a misalignment of at least $60\degree$ (Brinch et al., 2016), along with misaligned circumprimary and circumsecondary discs. The binary GG Tau A may be misaligned by $25\degree$-$30\degree$ from its circumbinary disc (Köhler, 2011; Aly et al., 2018). There is also evidence that binary 99 Herculis, with an orbital eccentricity of $0.76$, has a misaligned debris disc that is thought to be perpendicular to the orbital plane of the binary (Kennedy et al., 2012). Furthermore, there are several known circumbinary planets discovered by Kepler, two of which have a misalignment to the binary of roughly $2.5\degree$, Kepler-413b (Kostov et al., 2014) and Kepler-453b (Welsh et al., 2015). This misalignment suggests that the circumbinary disc may have been misaligned or warped during the planet formation process (Pierens & Nelson, 2018).

Misalignment between a circumbinary disc and the binary may occur through several possible mechanisms. First, turbulence in star-forming gas clouds can lead to misalignment (Offner et al., 2010; Tokuda et al., 2014; Bate, 2012). Secondly, if a young binary accretes material after its formation process, the accreted material is likely to be misaligned to the orbital binary plane (Bate et al., 2010; Bate, 2018). Finally, misalignment can occur when a binary star forms within an elongated cloud whose axes are misaligned with respect to the cloud rotation axis (e.g. Bonnell & Bastien, 1992).

The torque from binary star systems can impact the planet formation process compared to discs around single stars (Nelson, 2000; Mayer et al., 2005; Boss, 2006; Martin et al., 2014; Fu et al., 2015a, b, 2017). By understanding the structure and evolution of these discs, we can shed light on the observed characteristics of exoplanets.

Dissipation in a misaligned circumbinary disc causes tilt evolution. A disc around a circular orbit binary aligns to the orbital plane of the binary (e.g. Papaloizou & Terquem, 1995a; Lubow & Ogilvie, 2000; Nixon et al., 2011; Facchini et al., 2013; Foucart & Lai, 2014). However, for a disc around an eccentric binary, its angular momentum aligns to one of two possible orientations: alignment to the angular momentum of the binary orbit or, for sufficiently high initial inclination, alignment to the eccentricity vector of the binary (Aly et al., 2015; Martin & Lubow, 2017; Lubow & Martin, 2018; Zanazzi & Lai, 2018). The latter state is the so-called polar configuration in which the disc plane lies perpendicular to the binary orbital plane. The timescale for the polar alignment process may be shorter or longer than the lifetime of the disc depending upon the properties of the binary and the disc (Martin & Lubow, 2018).

Through SPH simulations Martin & Lubow (2017) found that an initially misaligned ($i=60\degree$) low mass circumbinary disc around an eccentric ($e_{\rm b}=0.5$) binary undergoes damped nodal oscillations and eventually evolves to a polar configuration. Martin & Lubow (2018) explored the properties of binaries and discs that lead to a final polar configuration. 1D linear models for the evolution of a low mass, nearly polar disc around an eccentric binary also show evolution to a polar configuration (Zanazzi & Lai, 2018; Lubow & Martin, 2018).

In this paper, we extend the work of Martin & Lubow (2017) and Lubow & Martin (2018) by studying the evolution of misaligned circumbinary discs around eccentric orbit binaries with lower initial inclinations that ultimately result in coplanar alignment with the binary. We apply both 3D SPH simulations and 1D linear equations for a variety of disc and binary properties.

First we examine test particle orbits around a circular and eccentric binary in Section 2. In Section 3, we use three dimensional hydrodynamical simulations of circumbinary discs to explore the evolution of aligning circumbinary discs for various values of inclination, eccentricity, and disc size. In Section 4, we apply a 1D linear model for the disc evolution. In Section 5, we apply the nearly rigid disc expansion procedure. We apply our results to the observed circumbinary disc in KH 15D in Section 6. Section 7 contains a summary.

2 Test particle orbits

In this section we consider the evolution of the orbit of an inclined test particle around a binary. For a circular orbit binary, or for a sufficiently low inclination test particle orbit around an eccentric binary, the test particle orbital angular momentum precesses about the binary angular momentum. An eccentric orbit binary generates a secular potential that is nonaxisymmetric with respect to the direction of the binary angular momentum. Consequently, the particle orbit tilt $i$ oscillates, the precession rate is nonuniform, and the precession is fully circulating. For higher inclination around an eccentric binary, the orbit precesses about the eccentricity vector of the binary and also undergoes oscillations in tilt. The particle in that case undergoes libration, rather than circulation (Verrier & Evans, 2009; Farago & Laskar, 2010; Doolin & Blundell, 2011).

We consider test particle orbits around an equal mass binary with $M_{1}=M_{2}=0.5M$, where $M$ is the mass of the binary and the semi-major axis of the binary is denoted as $a$. The particle orbits are calculated for four different binary eccentricities, $e_{\rm b}=0.0$, $0.3$, $0.6$, and $0.8$. The orbital period of the binary is given by $P_{\rm orb}=2\pi/\sqrt{G(M_{1}+M_{2})/a^{3}}$. The binary begins at periastron separation. We apply a Cartesian coordinate system $(x,y,z)$. The $x-$axis is along the binary eccentricity vector, whose direction is from the binary center of mass to the orbital pericenter. The $z-$axis is along the binary angular momentum. The test particle begins in a circular Keplerian orbit at position $(0,d,0)$ with velocity $(-\Omega_{\rm p}d\cos i_{0},0,\Omega_{\rm p}d\sin i_{0})$ where $\Omega_{\rm p}=\sqrt{G(M_{1}+M_{2})/d^{3}}$ is approximate angular frequency of a particle about the center of mass of the binary and $i_{0}$ is the initial particle orbit tilt with respect to the binary orbital plane. The longitude of the ascending node $\phi$ is measured from the $x$-axis. These initial conditions correspond to an initial longitude of the ascending node of $\phi_{0}=90\degree$.

Fig. 1 shows the test particle orbits in the $i\cos\phi$-$i\sin\phi$ phase space for binary eccentricities of $e_{\rm b}=0.0$ (upper left panel), $0.3$ (upper right panel), $0.6$ (lower left panel), and $0.8$ (lower right panel) for various initial inclinations. The test particles all begin at a separation of $d=5a$. For these test particle orbits, the separation does not affect these phase portraits, only the timescale on which the orbit precesses. Depending on the initial orbital inclination, the particle can reside on a circulating or librating orbit. The centers of the upper libration regions (for all panels except the circular case) corresponds to $i=90\degree$ and $\phi=90\degree$, while the centers for the lower librating regions correspond to $i=90\degree$ and $\phi=-90\degree$.

For higher binary eccentricity, the critical inclination angle that separates the librating solutions from circulating solutions is smaller. When the third body (in this case a test particle) is massive, the nodal libration regions shrink (see Fig. 5 in Farago & Laskar, 2010). The critical inclination for test particles that divides the librating and circulating solutions is

[TABLE]

(Farago & Laskar, 2010). For the eccentricities considered in Fig. 1 this is $i_{\rm crit}=54.9^{\circ}$ for $e_{\rm b}=0.3$, $i_{\rm crit}=30.8^{\circ}$ for $e_{\rm b}=0.6$ and $i_{\rm crit}=18.5^{\circ}$ for $e_{\rm b}=0.8$. Martin & Lubow (2018) found that the critical inclination is slightly higher for a disc than a test particle. This means that a disc is more likely to move towards coplanar alignment with the binary than a test particle. In the next section we consider the evolution of a hydrodynamic circumbinary disc and use these test particle orbits for comparison.

3 Circumbinary Disc Simulations

To model the alignment process of misaligned circumbinary discs around an eccentric binary, we use the 3D smoothed particle hydrodynamics (SPH; e.g., Price, 2012) code phantom (Lodato & Price, 2010; Price & Federrath, 2010; Price et al., 2017). phantom has been well tested and used to model misaligned accretion discs in binary systems (Nixon, 2012; Nixon et al., 2013; Martin et al., 2014; Doğan et al., 2015).

3.1 Simulation Setup

Table 1 summarises the initial conditions of the binary and disc parameters for the hydrodynamical simulations. We consider an equal mass binary with total mass $M=M_{1}+M_{2}$. The eccentric orbit of the binary lies in the $x$-$y$ plane with semi–major axis, $a$. The binary begins at time $t=0$ at apastron. The accretion radius of each binary component is $0.25a$. When a particle enters this radius, it is considered accreted and the particle’s mass and angular momentum are added to the sink particle. We consider binaries with eccentricities $e_{\rm b}=0.0$, $0.3$, $0.6$ and $0.8$. For each eccentricity, we begin with a low initial disc inclination somewhat below the critical value found from Equation (1). Table 2 summarises the setup for each simulation. For $e_{\rm b}=0.3$ we use $i=50\degree$, for $e_{\rm b}=0.6$ we use $i=30\degree$ and for $e_{\rm b}=0.8$ we use $i=15\degree$ We evolve each simulation to $3000$ binary orbits.

Each simulation has an initially low disc mass of $10^{-3}\,M$ and we ignore self–gravity. The low mass disc has a negligible dynamical affect on the orbit of the binary. Each simulation consists of $6\times 10^{5}$ equal mass gas particles that initially reside in a flat disc with an inner boundary of $2a$ and an outer boundary of $5a$. The inner boundary of the disc is chosen to be close to where the tidal torque truncates the inner edge of the disc (Artymowicz & Lubow, 1994). For misaligned discs, the tidal torque produced by the binary is much weaker allowing the disc to move closer to the binary (e.g., Lubow et al., 2015; Miranda & Lai, 2015; Nixon & Lubow, 2015; Lubow & Martin, 2018). The surface density profile is initially a power law distribution $\Sigma\propto R^{-3/2}$. We use a locally isothermal disc with sound speed $c_{\rm s}\propto R^{-3/4}$ and disc aspect ratio $H/r=0.1$ at $r=r_{\rm in}$. We take the Shakura & Sunyaev (1973) $\alpha$ to be $0.01$. From these values we derive an artificial viscosity ($\alpha_{\rm AV}$) of $0.4$ (a value of $\alpha_{\rm AV}=0.1$ represents the lower limit, below which a physical viscosity is not resolved in SPH) and set $\beta_{\rm AV}=2.0$ from the SPH description detailed in Lodato & Price (2010) which is given as

[TABLE]

where $\langle h\rangle$ is the mean smoothing length on particles in a cylindrical ring at a given radius (Price et al., 2017). With this value of $\alpha$, the disc with an initial outer radius of $5a$ is resolved with a shell-averaged smoothing length per scale height of $\langle h\rangle/H\approx 0.25$. For the simulation with a larger outer radius of $40a$, we have that $\langle h\rangle/H\approx 0.30$.

3.2 Results

In this section we describe the results of the hydrodynamical disc simulations for different values of the eccentricity of the binary orbit.

3.2.1 Circular binary with $e_{\rm b}=0.0$

The left hand panel of Fig. 2 shows the time evolution of the inclination and longitude of ascending node at a distance $3a$ (solid lines) and $5a$ (dashed lines) of a misaligned disc with an initial inclination of $60\degree$ around an circular binary (Run1 of Table 2). The inclination evolution of the disc shows that the disc is aligning to the binary orbital plane. Through viscous dissipation, the disc orbital angular momentum vector evolves towards alignment with the orbital angular momentum vector of the binary. The disc undergoes retrograde precession at a nearly constant (uniform) precession rate about the binary angular momentum vector. The disc inclination decreases monotonically. The right hand panel shows a spiral in the $i\cos\phi$-$i\sin\phi$ phase space as the disc aligns to the binary orbital plane.

3.2.2 Eccentric binary with $e_{\rm b}=0.3$

We consider a binary eccentricity of 0.3. Fig. 3 shows the time evolution of the inclination and longitude of ascending node at a distance $3a$ and $5a$ of an initially misaligned disc of $50\degree$ around the eccentric binary (Run2 of Table 2). The disc evolves towards alignment to the plane of the binary as in the circular binary case. However, during this process the disc undergoes tilt oscillations due to the eccentricity of the binary. The precession rate is nonuniform.

3.2.3 Eccentric binary with $e_{\rm b}=0.6$

The left hand panel of Fig. 4 shows the time evolution of the inclination and longitude of ascending node at a distance $3a$ and $5a$ for a misaligned disc with an initial inclination of $30\degree$ around a binary with eccentricity $e_{\rm b}=0.6$ (Run3 of Table 2). The right hand panel shows the spiral in the $i\cos\phi$-$i\sin\phi$ phase space as the disc aligns to the binary orbital plane. The precession rate is more nonuniform than in the case of $e_{\rm b}=0.3$ shown in Fig. 3 and the inclination oscillations are stronger.

3.2.4 Eccentric binary with $e_{\rm b}=0.8$

Finally, we consider a highly eccentric binary with $e_{\rm b}=0.8$. This eccentricity is at the upper end of the values for binary KH 15D determined by Johnson et al. (2004). We consider an initial misalignments of $15\degree$ (Run4 of Table 2). We show the initial orientation in the three Cartesian planes in the upper panels in Fig. 5. In the lower panels, we show the disc orientation at a time of $t=150\,P_{\rm orb}$ when the disc tilt has increased to about $50\degree$. The upper left hand panel in Fig. 6 shows the evolution of the tilt and the longitude of the ascending node. The right hand panel shows the $i\cos\phi$-$i\sin\phi$ phase space plot as the disc aligns to the binary orbital plane. As expected, as the binary eccentricity increases, the amplitude of the tilt oscillations also increases as expected from the test particle orbit case. In addition, the precession rate is more nonuniform, as seen in the lower left panel of Fig. 6.

3.2.5 Eccentric binary with a large disc

The simulations described thus far only dealt with moderately extended discs with a radial extent initially from $2a$ up to $5a$. For parameters relevant to protoplanetary discs, such discs precess in nearly solid body because the sound crossing timescale is shorter than the precession timescales. As discussed in Martin & Lubow (2018), close binaries may have a disc with a much larger radial extent relative to the binary separation. We now consider the disc evolution with a larger initial disc outer radius (Run5 of Table 2) .

We consider a disc with a radial extent initially of $40a$. Unlike previous simulations in this work, this disc has $1\times 10^{6}$ equal mass gas particles, more than in the other simulations, although the particle density and therefore spatial resolution is lower. The disc aspect ratio at the outer boundary is $0.047$. Extending the disc outer radius by a factor of $8$ increases the disc angular momentum compared to the previous simulations. We investigated whether there are significant dynamical effects that the extended disc exerts on the binary. The maximum deviation from the initial binary inclination and eccentricity is $0.0072\degree$ and $0.0104$, respectively. Thus, there are no significant dynamical effects on the binary. The initial disc setup is shown in the top panels of Fig. 7. The evolution of the tilt and longitude of ascending node are shown in Fig. 8. We show the results at three radii within the disc, $5a$, $10a$ and $25a$. For this larger disc, the sound crossing time over the radial extent of the disc is longer than the precession timescale. The inner parts of the disc begin a tilt oscillation while the outer parts of the disc remain close to their original value for longer. The lower panels of Fig. 7 show the disc at a time of $600\,P_{\rm orb}$. The outer parts of the disc have not changed much from the initial setup, while the inner parts of the disc are significantly tilted. We see evidence for disc breaking in this simulation.

To examine the behavior of the warp propagation, in Fig. 9 we show the surface density (top panel), inclination (middle panel), and longitude of the ascending node (bottom panel) as a function of radius at times $0\,\rm P_{orb}$, $10\,\rm P_{orb}$, $10^{2}\,\rm P_{orb}$, $10^{3}\,\rm P_{orb}$, and $2\times 10^{3}\,\rm P_{orb}$. The initial surface density (at $t=0$) has a profile of $\Sigma\propto r^{-3/2}$. As the disc evolves, the gas in the outer portions of the disc spreads outwards through viscosity. As time increases, the inclination of the inner portions of the disc increases due to these tilt oscillations and the wave travels outwards in time. From the $1000\,\rm P_{orb}$ curve in the middle panel, we see that the disc below a distance of about $20a$ is inclined more than the outer regions of the disc. Since the surface density at $1000\,\rm P_{orb}$ shows a dip at around $14a$, we find that the disc is broken.

Disc breaking occurs when the radial communication time-scale is larger than the is the precession time-scale, $t_{\rm c}>t_{\rm p}$. The disc is able to maintain radial communication via pressure induced bending waves that propagate at speed $c_{\rm s}/2$ for gas sound speed $c_{\rm s}$ (Papaloizou & Lin, 1995; Lubow et al., 2002). The radial communication time-scale can be approximated by

[TABLE]

(Lubow & Martin, 2018), where $h_{\rm out}$ is the disc aspect ratio at the outer edge, $s$ is related to the temperature profile of the disc ($T(r)\propto r^{-s}$), the angular frequency $\Omega_{\rm b}=\sqrt{\frac{GM}{a^{3}}}$. The nodal precession rate can be approximated by

[TABLE]

where

[TABLE]

The precession time-scale can be found by taking the inverse of the nodal precession rate. For a narrow disc we have $r_{\rm out}=r=5a$, $e_{\rm b}=0.8$, and $h_{\rm out}=0.0795$, which equates to $t_{\rm c}\approx 160\,\rm P_{orb}$ and $t_{\rm p}\approx 1317\,\rm P_{orb}$. Given that $t_{c}<t_{p}$, the narrow disc can rigidly precess. For example, we compare $t_{\rm p}$ to the numerical precession timescale $t_{\rm p,Run4}$ for simulation Run4 which is referenced in Fig. 6. We find that $t_{\rm p,Run4}\approx 1540\,\rm P_{orb}$ which is consistent with $t_{\rm p}$.

For a larger disc, $r_{\rm in}\ll r_{\rm out}$, the precession time-scale can be determined by taking the inverse of the global precession rate. The global precession rate of a disc is found by taking its angular momentum weighted average of the nodal precession rate $\omega_{n}(r)$. Therefore, the global precession time-scale is given as

[TABLE]

where $p$ is related to the initial surface density profile of the disc ($\Sigma\propto r^{-p}$), For an extended disc with $r_{\rm out}=40a$, $e_{\rm b}=0.8$ and $h_{\rm out}=0.0473$, we have $t_{\rm c}\approx 6114\,\rm P_{orb}$ and $t_{\rm p,global}\approx 2665\,\rm P_{orb}$. Since $t_{\rm c}>t_{\rm p,global}$, breaking can occur within the disc.

4 Nearly Coplanar Disc Linear Model

In this section we apply a 1D linear model to the disc evolution based on equations that assume that the level of tilt is small and that the density evolution can be ignored. The equations apply the secular torque due to an eccentric binary obtained by Farago & Laskar (2010). The advantage of using this approach is that solutions can be readily obtained over very long timescales for very large discs with far less computational effort than is required with SPH. Such an approach to modeling the circumbinary disc around KH 15D has been applied by Lodato & Facchini (2013) and Foucart & Lai (2014) for a circular orbit binary. The analysis presented in this section is similar to that of Lubow & Martin (2018) who analyzed a nearly polar disc around an eccentric orbit binary.

We consider an eccentric binary with component stars that have masses $M_{1}$ and $M_{2}$ and total mass $M=M_{1}+M_{2}$ in an orbit with semi–major axis $a_{\rm b}$ and eccentricity $e_{\rm b}$. To describe this configuration, we again apply a Cartesian coordinate system $(x,y,z)$ whose origin is at the binary center of mass and with the $z$-axis parallel to the binary angular momentum $\bm{J}_{\rm b}$ and the $x$-axis parallel to the binary eccentricity vector $\bm{e}_{\rm b}$. We consider the disc to be composed of circular rings that provide a surface density $\Sigma(r)$. The ring orientations vary with radius $r$ and time $t$ and orbit with Keplerian angular speed $\Omega(r)$. In this model, the disc surface density is taken to be fixed in time, i.e., viscous evolution of the disc density is ignored. We denote the unit vector parallel to the ring angular momentum at each radius $r$ at each time $t$ by $(\ell_{x}(r,t),\ell_{y}(r,t),\ell_{z}(r,t))$. We consider small departures of the disc from the $x-y$ plane, so that $|\ell_{x}|\ll 1,|\ell_{y}|\ll 1,$ and $\ell_{z}\approx 1$.

We apply equations (12) and (13) in Lubow & Ogilvie (2000) for the evolution of the disc 2D tilt vector $\bm{\ell}(r,t)=(\ell_{x}(r,t),\ell_{y}(r,t))$ and 2D internal torque $\bm{G}(r,t)=(G_{x},G_{y})$. The disc tilt $i$ in radians is related to the tilt vector by $i(r,t)=|\bm{\ell}(r,t)|=\ell(r,t)$. The tilt evolution equation is given by

[TABLE]

where $\bm{T}$ is the tidal torque per unit area due to the eccentric binary whose orbit lies in the $x-y$ plane. Equation (13) in Lubow & Ogilvie (2000) provides the internal torque evolution equation

[TABLE]

where $\alpha$ is the usual turbulent viscosity parameter, $\omega_{\rm a}(r)$ is the apsidal precession rate for a disc that is nearly coplanar with the binary orbital plane that is given by

[TABLE]

and

[TABLE]

for disc density $\rho(r)$. We apply boundary conditions that the internal torque vanishes at the inner and outer disc edges $r_{\rm in}$ and $r_{\rm out}$, respectively. That is,

[TABLE]

This is a natural boundary condition because the internal torque vanishes just outside the disc boundaries. Thus, any smoothly varying internal torque would need to satisfy this condition.

The torque term due to the eccentric binary follows from an application of equations (2.17) and (2.18) in Farago & Laskar (2010). The torque term is expressed as

[TABLE]

with

[TABLE]

and

[TABLE]

and

[TABLE]

where frequency $\omega_{\rm p}$ is given by

[TABLE]

We seek solutions of the form $\bm{\ell}\propto e^{i\omega t}$ and $\bm{G}\propto e^{i\omega t}$ and Equations (7) and (8) become

[TABLE]

and

[TABLE]

respectively. As usual, the physical values of $\bm{\ell}$ and $\bm{G}$ are obtained by taking their real parts.

5 Nearly Rigid Disc Expansion

5.1 Lowest order

We apply the nearly rigid tilted disc expansion procedure in Lubow & Ogilvie (2000). We expand variables in the tidal potential that is considered to be weak as follows:

[TABLE]

where $a$ and $b$ are given by Equations (14) and (15), respectively. $a$ and $b$ depend on the tidal potential and are regarded as first order quantities.

To lowest order, the disc is rigid and the tilt vector $\bm{\ell}^{(0)}=(\ell_{x}^{(0)},\ell_{y}^{(0)})$ is constant in radius. We integrate $r$ times Equation (17) over the entire disc and apply the boundary conditions given by Equation (11) to obtain

[TABLE]

where

[TABLE]

We then obtain for the disc precession rate in lowest order

[TABLE]

where the bracketed term involves the angular momentum weighted average

[TABLE]

We define the precession period as

[TABLE]

The tilt components are related by

[TABLE]

Because $|\ell^{(0)}_{x}|$ and $|\ell^{(0)}_{y}|$ differ, the disc undergoes nonuniform precession and secular tilt oscillations with tilt variations $i(t)$ with respect to the $x-y$ plane. The disc longitude of ascending node $\phi$ is related to the tilt vector by

[TABLE]

We take the initial disc longitude of ascending nodes to be $90^{\degree}$, so that 2D tilt vector $\bm{\ell}$ is initially aligned with the binary eccentricity vector. Figure 10 plots the longitude of the ascending node and the nodal precession rate as a function of time for various values of binary eccentricity. For $e_{\rm b}=0$, the precession rate is uniform and appears as the horizontal line. The precession rate becomes highly nonuniform at higher values of binary eccentricity. For $e_{\rm b}=0.8$, the precession rate varies about a factor of 10 over the precession period.

The results in the upper panel of Figure 10 for $e_{\rm b}=0.8$ are similar to those in the lower left panel of Figure 6 that are based on SPH simulations. The precession is nonuniform in both cases, with similar phase oscillations in time. One difference is that the precession period increases in time in the SPH simulations. This increase occurs because of the viscous disc density evolution that in turn changes the disc angular momentum. This effect is not included in the linear model.

The disc tilt varies in time as

[TABLE]

where $i_{0}=i(0)$ that occurs when the longitude of the ascending node is $90^{\degree}$. Figure 11 plots the tilt angle as a function of time for various values of binary eccentricity. Tilt oscillations occur because the binary potential is nonaxisymmetric around the direction of the binary angular momentum vector (the $z-$axis). For $e_{\rm b}\simeq 1$, the oscillations undergo extreme tilt variations $i(t)\propto i(0)\sqrt{(1-\cos{(2\omega^{(1)}t)})/(1-e_{\rm b})}$.

The normalised disc tilt and precession rates plotted in Figures 10 and 11 are independent of the disc properties such as its density and temperature distributions, provided that the level of disc warping is small, i.e., $\bm{\ell}(r,t)$ is nearly constant in radius.

Figure 12 plots the maximum tilt angle over time as a function of binary eccentricity implied by Equation (27) that occurs for $\cos{(2\omega^{(1)}t)})=-1$,

[TABLE]

Also plotted on the figure are the maximum inclinations for SPH simulations for models listed in Table 1 that start with $i_{0}=10^{\degree}$. As seen in the figure, the results of the SPH simulations agree well with the expected results based on linear theory. The plotted SPH results lie slightly below the expectations of linear theory, likely due to the effects of disc dissipation. Though the linear theory is valid for low inclinations, the SPH simulations that begin at higher inclinations, $i_{0}\leq 50\degree$ (Runs 2, 3, and 4) also agree quite well with the linear model.

6 Circumbinary disc of KH 15D

KH 15D is a spectroscopic binary T Tauri star in the cluster NGC 2264 and located at a distance of 760 pc (Sung et al., 1997). This system was originally thought to be a single variable star. But more in-depth observations showed this system had a stellar companion, which causes a peculiar light curve (Kearns & Herbst, 1998). The system is estimated to be an age of $3\times 10^{6}\,\rm yr$ and the total mass is roughly $1.3\,M_{\odot}$ (Hamilton et al., 2001; Johnson et al., 2004). The spectral classification of star A is K6/K7 (Hamilton et al., 2001) and star B is K1 (Capelo et al., 2012). The two stellar companions are of roughly equal mass and are on highly eccentric orbits embedded in the accretion disc which emits bipolar outflows (Hamilton et al., 2003; Deming et al., 2004; Tokunaga et al., 2004; Mundt et al., 2010). The binary has an eccentricity in the range of $e_{\rm b}=0.68$ to $0.8$ (Johnson et al., 2004), semi–major axis of $0.26\,\rm au$.

The light curve of KH 15D undergoes periodic eclipses in which the brightness drops by about $3.5$ magnitudes for a duration of roughly $24$ days with a orbital period of $48.37\,\rm days$ (Johnson et al., 2004; Winn et al., 2004; Hamilton et al., 2005). The duration of the eclipse has varied over time (see Fig.1 in Aronow et al. (2018) for the $I$-band light curve of KH 15D which shows the brightness of the system from $1951$ to $2017$). The brightness increased between 1995 and 2005 and the peak brightness decreased between 2006 and 2010 (Hamilton et al., 2001; Hamilton et al., 2005).

To understand what causes this light curve, Chiang & Murray-Clay (2004) and Winn et al. (2004) independently developed a model in which a circumbinary disc or ring that is misaligned to the orbital plane of an eccentric orbit blocks light from the binary and undergoes nodal precession. The nodal precession explains the time variations of the observed light curves. Between 1995 and 2010, the leading edge of the disc precessed across the orbit of Star A, while star B was fully occulted. During the time between 2010 and 2012, both stars A and B were only detectable through scattered light. Currently, the brightness of the system has increased as star B’s orbit has become uncovered from the the trailing edge of the precessing disc (Capelo et al., 2012; Windemuth & Herbst, 2014; Arulanantham et al., 2016; Aronow et al., 2018).

Previous hydrodynamical models for a gaseous disc in KH 15D have only modeled the binary as circular (Lodato & Facchini, 2013; Foucart & Lai, 2014). Our goal in modeling this system is to understand the properties of the disc, such as its radial extent, given the observed constraints. Based on the work by Chiang & Murray-Clay (2004) and Winn et al. (2004) we consider the disc to be observed nearly edge-on and inclined relative to the orbit of the binary. In addition, the binary eccentricity vector lies in the plane of the sky. Under these conditions, the line of ascending nodes of the disc should currently be $\phi\approx 90^{\degree}$.

We consider a model in which the disc tilt $i$ is below the critical value $i_{\rm crit}$ given by Equation (1) which implies that $20^{\degree}\la i_{\rm crit}\la 30^{\degree}$ for $0.6\la e_{\rm b}\la 0.8$ (Johnson et al., 2004). If the disc tilt is above this critical value, then the disc will evolve to a polar (perpendicular) alignment with the binary (Martin & Lubow, 2017). However, for this work, we only examine the conventional model where the disc or ring is precessing about the binary angular momentum vector.

We see from the Figures 6 and 10 that the precession rate is largest in magnitude at this phase $\phi\simeq 90^{\degree}$. For a binary eccentricity of $e_{\rm b}=0.8$, the precession rate is about $3$ times faster than the mean precession rate. The tilt at this phase is at a minimum value. At later times the retrograde precession rate $-d\phi/dt$ will be as much as an order of magnitude smaller and the tilt will be more than 3 times larger. These results are largely independent of the details of the disc/ring structure.

The observed occultation involves scattering by solid particles. Such particles would undergo differential precession of the orbits in the presence of the binary that would destroy the disc structure over time. Some mechanism is required to maintain the disc flatness. One possibility is the ring coherence is maintained by self-gravity in analogy to planetary rings (Chiang & Murray-Clay, 2004). Another possibility is that the solids are coupled to a gas disc that maintains its flatness by pressure effects (Papaloizou & Terquem, 1995b; Larwood & Papaloizou, 1997; Lubow & Ogilvie, 2000). We analyze the latter model.

To analyze the system further, we numerically solve Equations (7) and (8) subject to boundary conditions given in Equation (11) for disc modes, as is described in Lubow & Martin (2018). We analyze discs whose parameters are listed in Table 1, where $s$ and $p$ are defined by $T(r)\propto r^{-s}$ and $\Sigma(r)\propto r^{-p}$, respectively. In all cases we assume an equal mass binary $M_{1}=M_{2}$. The disc inner radii should increase somewhat with binary eccentricity, but we ignore that effect for the two values of eccentricity being considered.

6.1 Precession period constrained model

Previous disc models for this system by Lodato & Facchini (2013) and Foucart & Lai (2014) applied a constraint on the disc precession period based on the Chiang & Murray-Clay (2004) model. In that model, the precession period is approximately $3000\,\rm yr$ or about $2.09\times 10^{4}P_{\rm b}$. However, this period value is determined by considering a narrow ring and so it is not clear how well this constraint would apply to a broad disc. This model may be appropriate if the occultation is due to material in the somewhere in the middle of the radial extent of disc, rather than the outer edge. We consider an alternate model in the next subsection. We describe results for a disc period constrained model based on results from linear modes.

We adopt the disc parameters similar to those of Lodato & Facchini (2013) that are listed for Models A-D in Table 3. In addition we consider two values of binary eccentricity $e_{\rm b}=0.6$ and 0.8, while the previous models considered a circular orbit binary. Table 4 contains results for these models. The columns in the table are for the values for the disc outer radius $r_{\rm out}/a$, decay timescale of the tilt in year $\tau$, and the maximum normalized warp value across the disc Max($r/\ell_{\rm in}|\ell/dr|$). The latter is the magnitude of the logarithmic radial derivative of the tilt vector $\bm{\ell}$ divided by the magnitude of the tilt at the disc inner edge, $\ell_{\rm in}$ (see also Section 3.2 of Martin & Lubow (2018) for more details). Since this value is small, less than $H/r$, for all disc models, the disc warp is very mild and so the disc behaves quite rigidly. In addition, the linear treatment of the disc evolution is well justified for discs with small tilts.

The numerical results are similar to those in Lodato & Facchini (2013) and Foucart & Lai (2014) once slight differences in the model parameters are taken into account. For example, Table 1 in Lodato & Facchini (2013) has a value for $r_{\rm out}=26a$ for $p=0.5$, while we obtain a value of 27.6 in Model A. The small difference is likely due to binary eccentricity and the slightly different value of the binary semi-major axis adopted. In any case, as obtained previously, the disc model decays rapidly compared to the system lifetime of a few million years. The decay rate is proportional to the $\alpha$ value in the disc (for a fixed disc structure) and suggests that reductions to $\alpha\sim 10^{-3}$ are required to provide a sufficiently slow tilt decay.

The effect of binary eccentricity is to slightly decrease the required disc outer radius, as seen in comparing Models A and B and also Models C and D. In addition the decay timescale slightly increases with increasing binary eccentricity.

6.2 Velocity constrained model

There is an observational constraint on the speed of the occulting disc/ring in the plane of the sky. By comparing frames 1 and 4 in Figure 1 of Aronow et al. (2018), we estimate that the occultation occurs across distance $\simeq a(1+e_{\rm b})$ over a time $\tau_{0}$ of roughly 40 years. If we take the standard value of $a=0.26$ AU, we then have a constraint on the transverse occulting velocity $v\sim a(1+e)/\tau_{0}$, that is

[TABLE]

As discussed above, this velocity occurs for the longitude of ascending nodes that we take $\phi=90^{\degree}$. We apply this velocity constraint for various models computed from linear modes.

For a narrow ring, we determine the ring radii as a function of binary eccentricity that satisfy the velocity constraint (29) at $\phi=90^{\degree}$. The results are plotted in Figure 13. The radii agree well with the $\sim 3$ AU estimated by Chiang & Murray-Clay (2004). For larger values of binary eccentricity, the ring radius increases with eccentricity.

For a broader disc, we assume the occultation is dominated by the disc outer edge. We then apply the velocity constraint at that radius. In Figure 14 we plot the disc outer radius as a function of disc inner radius for Model E of Table 3 that has a disc with surface density parameter $p=0.5$ and assumed binary eccentricity $e_{\rm b}=0.75$. The value of $e_{\rm b}$ is close to the best fit value of 0.74 in the model of Johnson et al. (2004).

The inner radius of the circumbinary disc in KH 15D is expected to range roughly from $r=0.5{\rm AU}$ at higher viscosities $\alpha>0.01$ to $r=1{\rm AU}$ at small viscosities $\alpha<1\times 10^{-5}$ due to the balance of viscous torque with tidal torques (Artymowicz & Lubow, 1994). The disc torque increases for smaller disc inner radii and is insensitive to the disc outer radius for $r_{\rm in}\ll r_{\rm out}$. The disc angular momentum increases with the disc outer radius. For smaller disc inner radii, there is a stronger torque due to the binary that requires a larger disc outer radius to produce the same velocity at the disc outer edge. There is then an inverse relationship between the inner and outer disk radii.

In Figure 15, we plot the tilt decay timescale for Model E of Table 3 as a function of disc inner radius with parameters $s=1.0,\alpha=0.01$, and $H/r(r_{\rm in})=0.1$. In this case, the disc decay timescale is typically of order the disc lifetime of a few million years or longer. The velocity constrained model undergoes slower tilt decay than the similar models for the period constrained models of Section 6.1. In particular, no reduction of $\alpha$ below $0.01$ is required in this case to meet the requirement that the disc decay timescale exceed the disc lifetime.

We now consider the velocity constrained model with Model F in Table 3 that has the same parameters as Model E, but with $p=1$. In this case, the disc outer radius is required to be considerably larger than the $p=0.5$ case, as seen in Figure 16. We limited the plot to $r_{\rm in}\geq 1{\rm AU}$ because at smaller values of $r_{\rm in}$ the disc outer radius gets very large. The reason is that the surface density falls off faster with radius. The increased radius in the $p=1$ case is required to produce a large enough disc angular momentum that is sufficient to reduce the disc velocity at the outer edge in order the meet the velocity constraint. We find that the tilt decay timescale with $p=1$ is even longer than indicated in Figure 15. Again, no reduction in adopted $\alpha=0.01$ is required for the tilt to survive a few million years.

These models have assumed that the occultation occurs due to material at the gaseous disc outer edge. The occultation is likely due to solids (dust) that could have migrated inward somewhat from the gaseous disc outer edge. This effect would make the velocity constraint easier to satisfy. That is, the gas disc outer radius could be smaller than indicated in Figures 14 and 16 and satisfy the velocity constraint of Equation (29). The level of reduction for $r_{\rm out}$ depends on the degree to which the solids have migrated inward, as is discussed in Section 6.3.

6.3 Constraint on thickness of obscuring layer

The obscuring material likely consists of solids that form a dust embedded layer within the gaseous disc. IR observations suggest that the solids consists of 1 to 50 micron size particles (Arulanantham et al., 2016; Arulanantham et al., 2017). We define the full thickness of the obscuring layer as $2T$. The observations of KH 15D show that both stars were occulted over a time interval $\tau\sim 5$ years (see Figure  1 of Aronow et al., 2018). The disc thickness can then be expressed as

[TABLE]

where the term involving $\tau$ is due to the transverse velocity (precession) of the disc given in Equation (29) and $\tau_{0}$ is the time for the disc leading edge to precess across both stars that we estimate as $\tau_{0}\sim 40$ years, as discussed in Section 6.2. The term involving $\tau$ is then a small correction $\sim 10\%$ that we ignore. The constraint on $T$ then implies that

[TABLE]

We consider how this constraint applies to the velocity constrained model of Section 6.2. For the narrow ring case with $e_{\rm b}=0.75$ and $r_{\rm out}=3$ AU (see Figure 13), we have then $T/r_{\rm out}\sim 0.1\sin{i}$. For $e_{\rm b}=0.75$ and a circumbinary disc with $p=0.5$ that is tidally truncated by the binary at its inner radius at $r_{\rm in}\sim 1$ AU, we have from Figure 14 that $r_{\rm out}\sim 30$ AU and so $T/r_{\rm out}\sim 0.01\sin{i}$. For a circumbinary disc with the same set of parameters, but with $p=1$, we have that $T/r_{\rm out}\sim 0.001\sin{i}$. For a narrow ring, the thickness of the occulting solid layer is comparable to the thickness of the gaseous disc layer, if $\sin(i)$ is not small, which suggests that mild settling of solids has occurred. But, the broad disc $T/r_{\rm out}$ values are significantly smaller than the assumed gas disc aspect ratio $H/r\sim 0.1$, typical of protostellar disc aspect ratios. Such small $T/r_{\rm out}$ values suggest that settling of solids towards the disc midplane has occurred. Such settling suggests that the radial drift of solids might have also occurred so that the velocity constraint may be satisfied with a smaller gaseous disc outer radius, as discussed in Section 6.2.

To produce such thin layers in the broad disc cases of Section 6.2 requires that the level of disc turbulence be very low. Using equations 19 and 20 of Fromang & Nelson (2009) and setting the Schmidt number to unity, we estimate that

[TABLE]

where $t_{\rm s}$ is the stopping time for the particles given by equation 10 of Fromang & Nelson (2009). For $\sim 50$ micron particles and $H=0.1r$, we have that the upper limit to $\alpha$ is

[TABLE]

where

[TABLE]

Aronow et al. (2018) report an upper limit of the disc mass as $\simeq 1.7\times 10^{-3}M_{\odot}$ based on ALMA nondetections. For the outer parts of the velocity constrained disc in Figure 14 with $p=0.5$ and $r_{\rm in}=1$ AU and $r=r_{\rm out}=30$ AU, we obtain for a disc with $M_{\rm d}=0.001M_{\odot}$ and $e_{\rm b}=0.75$ from Equation (33) that $\alpha\sim 10^{-5}\sin^{2}{i}$. For outer parts of the velocity constrained disc in Figure 16 with $p=1.0$ and $r_{\rm in}=1$ AU and $r=r_{\rm out}=170$ AU, we obtain for a disc with $M_{\rm d}=0.001M_{\odot}$ and $e_{\rm b}=0.75$ from Equation (33) that $\alpha\sim 10^{-4}\sin^{2}{i}$ for p=1. Such levels of turbulence are extremely low. Also such thin layers suggest that the density of dust near the disc midplane is greater than the gas density. This configuration is subject to various instabilities, such as shear instability and streaming instability (Youdin & Shu, 2002; Youdin & Goodman, 2005). It is not clear that such thin layers can exist.

Less extreme values of $\alpha$ can occur if the occulting material resides at smaller radii, so that $T/H$ is larger. The smaller radii could occur due to the inward drift of solids. The velocity constraint in Equation (29) can be satisfied by the occulting solids because the precession rate is controlled by the more extended gas disc. We consider the case that $p=1$ and apply the rigid tilt approximation that assumes the disc remains flat during its evolution. In that case, the velocity constraint is satisfied provided that the outer radius of the gaseous disc satisfies

[TABLE]

where $r_{\rm s}$ is the radius of the occulting solids and $v$ is given by Equation (29). This equation holds for $r_{\rm in}=1{\rm AU}\ll r_{\rm s}<r_{\rm out}$. If the occulting occurs at $r_{\rm s}=4$ AU for the $p=1$ disc model described in the previous paragraph with $M_{\rm d}=0.001M_{\odot}$, then $T/H\sim 0.5\sin{i}$ at $r=r_{\rm s}$, then $\alpha\sim 10^{-4}\sin^{2}{i}$ by Equation (33), and $r_{\rm out}=16$ AU by Equation (35). Higher values of $\alpha\ga 0.01$ can occur for very small disc masses $M_{\rm d}\la 10^{-5}\sin^{2}{i}\,M_{\odot}$. For comparison, in the case of HL Tau, Pinte et al. (2016) found that a thin sublayer of millimeter sized grains $T/H\la 0.2$ could account for the observed properties of the system that in turn imposed an upper limit on $\alpha\sim 3\times 10^{-4}$.

The Stokes number for dust grains compares the stopping time $t_{\rm s}$ to the dynamical time. For a disc with $p=1$, its value at radius $r$ is estimated as

[TABLE]

(cf. Fromang & Nelson, 2009), where $r_{\rm g}$ is the grain size. With $r_{\rm g}=50\,\rm\micro m$, $M_{\rm d}=0.001\,\rm M_{\odot}$, $r_{\rm out}=16\,\rm AU$, and $r=r_{\rm s}=4\,\rm AU$, then ${\rm Stk}\simeq 6\times 10^{-4}$. For these parameters, the dust is well coupled to the gas. The inward radial drift velocity due to gas drag is $v_{\rm r}\sim(H/r)^{2}\,Stk\,\Omega r\sim 10^{-5}\Omega r$ with $H/r=0.1$ and $r=16\,\rm AU$ (Armitage, 2013). The drift timescale near the disc outer edge is then of order $10^{6}$ years. Its numerical value in this case is not sensitive to $p$ for $0.5\leq p\leq 1.5$ Shorter drift timescales occur for a less massive disc.

Disc warping could also influence the effective value of $T$ by making the requirements on the thickness on the solids layer even stronger (thinner layer), but we do not consider its effects here. Another possibility is that the disc does not contain significant amounts of gas with associated turbulence, but instead essentially consists of only solids. The coherence of the disc or ring against the effects of differential precession is due to the self-gravity (Chiang & Murray-Clay, 2004). For such a ring, some of the linear theory results in this paper still hold, such as those in Figures 11, 12, and 13.

7 Summary

We have analysed the behavior of a mildly tilted low mass circumbinary disc in an eccentric orbit binary star systems by means of SPH simulations and linear theory. The disc undergoes nonuniform precession and tilt oscillations due to the effects of the binary eccentricity (e.g., Figs. 6 and 10). For moderately broad discs (whose outer radii are a few times the inner radii) with typical protostellar disc parameters, the disc can precess coherently with little warping. Larger discs can undergo breaking (Fig. 7). For small initial tilts, the results of the SPH simulations agree well with linear theory (e.g., Fig.  12). The amplitude of the tilt oscillations increases with binary eccentricity. The disc tilt undergoes damped oscillation in time and ultimately approaches a coplanar alignment with the binary.

We have analyzed a model for binary KH 15D that is based on a mildly tilted precessing disc that orbits an eccentric binary. The model suggests that the disc tilt relative to the binary orbit is currently at a minimum value and that the retrogade precession rate is currently at its largest value. We considered a period constrained model for the disc, along the lines of the previous circular orbit binary studies (Lodato & Facchini, 2013; Foucart & Lai, 2014), but taking into account the binary eccentricity. We find that the large binary eccentricity changes the inferred disc outer radii by a small amount. To satisfy the disc tilt lifetime requirements, the disc $\alpha$ value must be small, less than about 0.001, as is also consistent with the earlier studies.

We then considered a model in which the outer disc edge precession velocity is constrained by the observed changes in the binary eclipse properties (e.g., Aronow et al., 2018). We determined the relation between the disc inner and outer radii subject to this constraint. We find that discs whose inner radius is tidally truncated by the binary typically have outer radii of $\sim 30-170$ AU depending on the disc density profile. The disc outer radii are reduced if there is inward radial migration of solids that are responsible for the binary occultation. Narrow disc radii are about $3$ AU, in agreement with Chiang & Murray-Clay (2004).

The recent reappearance of Star B places strong constraints on the thickness of an occulting layer of solids/dust. The most reasonable models involve a thin layer of dust that has settled towards the midplane of a low mass gaseous disc $M_{\rm d}<0.001M_{\odot}$ and has migrated considerably inward. Such thin layers suggest that the disc turbulence is very weak $\alpha\ll 0.001$. Stronger turbulence can occur for smaller mass discs. For a narrow ring, less extreme settling and levels of turbulence are required. Another possibility is that the disc/ring consists of a thin disc of solids with little gas (e.g., Chiang & Murray-Clay, 2004).

As noted in Martin & Lubow (2017), it is also possible that the disc is instead evolving to a polar (perpendicular) alignment with the binary. For this to occur, the disc tilt needs to be $i\ga 30^{\degree}$.

Acknowledgments

We much appreciate Hossam Aly for beneficial conversations and for carefully reviewing the paper. 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. SL thanks Eugene Chiang for insightful discussions. We acknowledge support from NASA through grant NNX17AB96G. Computer support was provided by UNLV’s National Supercomputing Center.