Warps and breaks in circumbinary discs

ABSTRACT

Disc warping, and possibly disc breaking, has been observed in protoplanetary discs around both single and multiple stars. Large warps can break the disc, producing multiple observational signatures. In this work, we use comparisons of disc time-scales to derive updated formulae for disc breaking, with better predictions as to when and where a disc is expected to break and how many breaks could occur. Disc breaking is more likely for discs with small inner cavities, cooler temperatures, and steeper power-law profiles, such that thin, polar-aligning discs are more likely to break. We test our analytical formulae using three-dimensional grid-based simulations of protoplanetary discs warped by the gravitational torque of an inner binary. We reproduce the expected warp behaviours in different viscosity regimes and observe disc breaking at locations in agreement with our derived equations. As our simulations only show disc breaking when disc viscosity is low, we also consider a viscous criterion for disc breaking, where rapid alignment to the precession vector can prevent a break by reducing the maximum misalignment between neighbouring rings. We apply these results to the GW Orionis circumtriple disc and find that the precession induced from the central stars can break the disc if it is relatively thin. We expect repeated or multiple disc breaking to occur for discs with sufficiently steep power-law profiles. We simulate a polar-aligning disc around an eccentric binary with steep power-law profiles and observe two separate breaking events at locations in rough agreement with our analytical predictions.

1 INTRODUCTION

Many astrophysical discs are observed to not lie entirely within a single plane. These discs are generally referred to as ‘warped’ discs, changing their orientation in three-dimensional space with radius. Warped discs have been found at all different astrophysical scales through direct observation, from planetary rings to the discs of galaxies (Burke 1957; Kerr, Hindman & Carpenter 1957; Shu, Cuzzi & Lissauer 1983; Burrows et al. 1995). Many other warps have been inferred from their effects on their nearby surroundings, including maser emission in the vicinity of supermassive black holes (Miyoshi et al. 1995; Martin 2008), to explain long-term periodicity in close X-ray binaries (Katz 1973; Petterson, Rothschild & Gruber 1991; Scott, Leahy & Wilson 2000), and shadows cast in protoplanetary discs (Marino, Perez & Casassus 2015; Debes et al. 2017). The ubiquity of such a structure warrants careful study of its behaviour and evolution.

The warping of discs is a well-studied phenomenon, and a large amount of work has gone into deriving the analytical behaviour of the evolution and testing this behaviour via numerical methods (see Papaloizou & Pringle 1983; Pringle 1992; Papaloizou & Lin 1995; Ogilvie 1999; Lodato & Pringle 2007; Lodato & Price 2010, as well as Nixon & King 2016 for an overall review). In a warped disc, the angular momentum vector of the disc (⁠$\mathit{L}$⁠) changes with radius (r). The overall evolution of a warped disc can be effectively described with two different behaviours, depending on the relative sizes of the disc aspect ratio of scale height to radius, $h/r$⁠, to its $\alpha$-viscosity coefficient (Shakura & Sunyaev 1973). Thin, high-viscosity discs with $\alpha >h/r$ are said to be in the diffusive or viscous regime, where the warp flattens diffusively on time-scales inversely proportional to a second viscosity coefficient ${\nu }_2$⁠. Thicker, low-viscosity discs with $\alpha <h/r$ are instead considered to be in the bending-wave or wave-like regime, in which the warp launches a wave that propagates through the disc.

In the diffusive regime, the evolution of a warped disc may be described as a diffusion equation in $\mathit{L}$ with terms proportional to the standard viscosity ${\nu }_1$ and a warp viscosity ${\nu }_2$⁠. Terms including ${\nu }_1$ describe the standard evolution of a flat disc (i.e. Lynden-Bell & Pringle 1974), while terms including ${\nu }_2$ describe the evolution of the warp in the disc (Pringle 1992). There are also additional terms present, including a third viscosity term ${\nu }_3$⁠, which describe an induced precession from the warp as neighbouring annuli induce a torque perpendicular to the direction of the warp (Ogilvie 1999). Subsequent works have also derived evolution equations in the bending-wave regime in which the effects of pressure dominate over viscosity (Papaloizou & Lin 1995; Demianski & Ivanov 1997; Lubow & Ogilvie 2000), as well as more general sets of equations that encompass both regimes (Martin et al. 2019; Dullemond, Kimmig & Zanazzi 2022).

The relative sizes of the viscosity coefficients ${\nu }_1$⁠, ${\nu }_2$⁠, and ${\nu }_3$ also depend on the warp amplitude, characterized by the dimensionless parameter $\psi$⁠, which is defined as

where $\mathit{l}$ is the unit angular momentum vector at r. The dependence of the coefficients on $\psi$ for large, non-linear warps and the disc viscosity $\alpha$ is characterized in Ogilvie (1999).

When warps become very large, the disc may be unable to communicate across the warp effectively. The warp may become unstable, and disc breaking can occur, splitting the disc into separate rings. This phenomenon has been studied analytically (Doǧan et al. 2018; Raj, Nixon & Doğan 2021), and previously observed in simulations of discs in binary systems and around single and binary black holes (Larwood et al. 1996; Nixon et al. 2012; Nixon, King & Price 2013; Nealon, Price & Nixon 2015). Disc breaking is also thought to have occurred in some discs as a possible explanation for the observed geometry or accretion kinematics (Casassus et al. 2015; Facchini, Juhász & Lodato 2018; Zhu 2019; Kraus et al. 2020; Nealon et al. 2022). However, the exact conditions that lead a disc to break are still unclear, as well as where the location of the break (known as the breaking radius) will occur.

In this paper, we attempt to improve the currently known methods of determining a disc breaking event by studying the evolution of a disc under the influence of an external precession torque. In our context, we consider a misaligned disc around a binary star system, analogous to protoplanetary discs forming in young star clusters. We will see that under the influence of the binary, the disc can develop warps and breaks depending on the internal conditions of the disc.

The previous studies of disc breaking in binaries have been largely based on the use of smoothed particle hydrodynamics (SPH) codes.With this Lagrangian code, the resolution is determined by the particle density and is therefore a concern near disc breaks. The process of the initial breaking event, as well as the subsequent evolution of the break and interaction between the misaligned inner and outer discs, involves sparse regions that may be poorly characterized in an SPH simulation unless large numbers of particles are used. In this paper, we apply an Eulerian code whose resolution is instead limited by the grid size. In addition, we take advantage of the fact that disc breaking occurs more easily for circumbinary discs that are evolving to a polar orientation, compared to evolving to a coplanar orientation, around eccentric orbit binaries.

Our paper is organized as follows. In Section 2, we discuss the time-scales relevant to disc breaking and describe our new method of calculating the location of a break. Section 3 outlines the numerical methods we use to test our new method for calculating the breaking radius. We present our results in Section 4. We discuss these results in Section 5 and conclude in Section 6.

2 TIME-SCALES OF DISC EVOLUTION

In this section, we examine the different evolutionary time-scales that affect a warped disc and consider how these time-scales determine the breaking radius of the disc. We consider the situation of a circumbinary disc surrounding a central binary. The binary consists of two stars with masses $M_1$ and $M_2$⁠, with total mass labelled $M_{\mathrm{t}\mathrm{o}\mathrm{t}}=M_1+M_2$⁠, orbital semimajor axis $a_{\mathrm{b}}$⁠, and orbital eccentricity $e_{\mathrm{b}}$⁠. We choose radial power-law profiles for mid-plane density, temperature, and surface density, using the indices d, s, and p given by¹

and

For a steady-state Keplerian disc, these power-law indices are related through the equation

since $\rho =\mathrm{\Sigma }/H$ (e.g. Pringle 1981) and the sound speed is given by $c_{\mathrm{s}}=h\mathrm{\Omega }\propto T^{1/2}$⁠.

2.1 Disc breaking in the wave-like regime

The torque from the binary potential causes a ring of material to nodally precess at a frequency (e.g. Farago & Laskar 2010)

where ${\mathrm{\Omega }}_{\mathrm{b}}$ is the orbital frequency of the binary, r is the radius of the ring, and k is a constant that depends on the binary masses, eccentricity, and whether the disc is close to a coplanar or polar configuration (Aly et al. 2015; Martin & Lubow 2017, 2018; Lubow & Martin 2018; Zanazzi & Lai 2018). For a circumbinary disc, the inner regions feel a stronger torque from the binary than the outer regions. In the absence of communication between the rings of the disc, the inner rings precess faster than the outer rings. Therefore, depending upon the communication time-scale between the rings of the disc, the disc may become warped.

If the radial communication time-scale is short compared with the precession time-scale, then the disc can precess like a solid body. For a disc with inner radius $r_{\mathrm{i}\mathrm{n}}$ and outer radius $r_{\mathrm{o}\mathrm{u}\mathrm{t}}$⁠, the global precession time-scale $t_{\mathrm{p}}$ of the disc as a solid body is given by

The radial communication time-scale in the wave-like regime $t_{\mathrm{c}}$ is given by (equations 33 and 31 of Lubow & Martin 2018)

where $h_{\mathrm{o}\mathrm{u}\mathrm{t}}$ is the scale height at the outer edge of the disc.

We expect disc breaking to occur when the global precession time-scale is shorter than the communication time-scale, $t_{\mathrm{p}}<t_{\mathrm{c}}$⁠. Analytical estimates of the disc breaking radius have been given in many previous works (e.g. Nixon et al. 2013; Lubow & Martin 2018). Combining equations (7) and (8) and solving for $r_{\mathrm{o}\mathrm{u}\mathrm{t}}$ as the breaking radius gives

The time-scales given in equations (7) and (8) assume that $r_{\mathrm{i}\mathrm{n}}\ll r_{\mathrm{o}\mathrm{u}\mathrm{t}}$⁠; i.e. the discs have a large radial extent. During breaking events, newly created discs may break off as radially thin rings, so a more complete description of these equations is required.

Here, we attempt to derive a more refined estimate of $t_{\mathrm{p}}$ and $t_{\mathrm{c}}$ by revisiting these equations. The disc nodal precession frequency is given in equations (16) and (17) of Lubow & Martin (2018) as

where $\mathrm{\Omega }$ is the orbital frequency and k is a general constant depending on the parameters of the binary. Equation (16) of Lubow & Martin (2018) and equation (5) of Smallwood et al. (2019) give the value of k as

The above expressions for k are of simple analytical forms for nearly coplanar and polar orientations. The more general case is discussed in Section 5.1.

The fraction containing the binary masses, $M_1{M}_2/M_{\mathrm{t}\mathrm{o}\mathrm{t}}^2$⁠, may also be written as $q/(1+q{)}^2$⁠, where $q=M_2/M_1$ is the binary mass ratio. Equation (11) implies that the value of k is largest (and thus the precession rate quickest) for equal-mass binaries. Low-eccentricity binaries produce larger k values for coplanar discs, while higher eccentricity binaries give larger values for polar discs. When considering disc breaking, larger values of k reduce the precession time-scale, and thus are more favourable for producing a break. The sign of k also determines the direction of precession; a positive value of k signifies counterclockwise precession (in the opposite direction of disc rotation), while a negative k denotes clockwise precession.

We define the disc precession time to be $t_{\mathrm{p}}=1/{\omega }_{\mathrm{p}}$⁠. Evaluating the integrals in equation (10) for a disc with the surface density profile in equation (4) and Keplerian rotation profile $\mathrm{\Omega }={\mathrm{\Omega }}_{\mathrm{K}}\propto r^{-3/2}$⁠, we find for the disc precession time

In the limit that $r_{\mathrm{o}\mathrm{u}\mathrm{t}}\gg r_{\mathrm{i}\mathrm{n}}$⁠, this tends to equation (7).

To calculate the communication time of a warp across the disc, we consider the travel speed of bending waves starting from the inner edge. These waves travel at half the local sound speed (Papaloizou & Lin 1995). Assuming a radial power-law profile for the disc temperature as in equation (3), we calculate the radial communication time as

An example of these time-scales is shown as a function of disc radius in Fig. 1. For geometrically thin rings $(r_{\mathrm{i}\mathrm{n}}\approx r_{\mathrm{o}\mathrm{u}\mathrm{t}})$⁠, equations (12) and (13) reproduce the correct limiting behaviour $(t_{\mathrm{p}}\approx 1/{\omega }_{\mathrm{p}}$ and $t_{\mathrm{c}}\approx 0)$⁠, which equations (7) and (8) do not.

Using equations (12) and (13), we can derive an updated location for the breaking radius. In regions where $t_{\mathrm{p}}>t_{\mathrm{c}}$⁠, a precession-induced warp can be communicated via bending waves and allow the disc to precess rigidly. Taking the breaking radius $r_{\mathrm{b}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{k}}$ to be the value of $r_{\mathrm{o}\mathrm{u}\mathrm{t}}$ at which $t_{\mathrm{p}}=t_{\mathrm{c}}$⁠, for $r_{\mathrm{o}\mathrm{u}\mathrm{t}}>r_{\mathrm{b}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{k}}$ we have $t_{\mathrm{p}}<t_{\mathrm{c}}$ and the disc will precess faster than it is able to communicate the precession to its outer edge, causing the disc to break. This criterion approximately says that the breaking radius occurs at the distance that a bending wave can travel over a disc precession time-scale, starting from the disc inner edge.

Whether the curves of $t_{\mathrm{p}}$ and $t_{\mathrm{c}}$ intersect is determined by the disc geometry$(r_{\mathrm{i}\mathrm{n}},r_{\mathrm{o}\mathrm{u}\mathrm{t}},h/r)$⁠, disc structure (power-law slopes d, s, and p), and inner binary arrangement (q, $a_{\mathrm{b}}$⁠, and $e_{\mathrm{b}}$⁠, which determine the constant k). The slope of equation (12) is determined by the surface density power-law slope p and the ratio $r_{\mathrm{o}\mathrm{u}\mathrm{t}}/r_{\mathrm{i}\mathrm{n}}$⁠, whereas the slope of equation (13) is determined by the temperature power-law slope s. The location of the inner radius, $r_{\mathrm{i}\mathrm{n}}$⁠, sets the fastest precession rate of the entire disc and the initial height of the $t_{\mathrm{p}}$ curve, and is a strong factor in determining whether or not a disc will break. Because of this sensitive dependence on $r_{\mathrm{i}\mathrm{n}}$⁠, nearly polar-aligned discs are more likely to have breaks (Sections 3.1 and 3.2), since the reduced strength of the binary torque allows them to maintain smaller inner cavities (Miranda & Lai 2015; Franchini, Lubow & Martin 2019) and therefore faster precession rates. In addition, for fixed binary and disc parameters, the precession rate normalized by binary frequency of a nearly polar disc is faster than for a nearly coplanar disc for binary eccentricity $e_{\mathrm{b}}>0.41$⁠, as follows from equation (11).

Once the disc breaks into separate rings, the outer disc will have its time-scales reset as its inner radius has now moved to the location of the break. These time-scales can grow again starting from the inner edge of the outer disc and, depending on the conditions of the outer disc, they may intersect again and cause the outer disc to break into yet another set of rings. Thus, being applied repeatedly, our breaking condition could predict the breaking radii of multiple breaks. In Section 2.4, we explore the possibility of multiple disc breaking in more detail, examining which disc parameters are required to cause multiple disc breaking.

2.2 A criterion for viscous disc breaking

Discs align towards their precession vector (⁠${\mathit{L}}_{\mathrm{b}}$ for coplanar orientations, or ${\mathit{e}}_{\mathrm{b}}$ for polar orientations) in both the wave-like and viscous regimes on a time-scale that is roughly inversely proportional to $\alpha$ (King et al. 2013; Lubow & Martin 2018, equation 29):

If the alignment time-scale is significantly faster than the precession rate, the disc will align to the precession vector before it can precess a significant amount. This can suppress large warps generated by azimuthal ‘twisting’ motions (Raj et al. 2021), which may prevent disc breaking even if the other disc parameters would normally allow a break to occur. We can express this relation by finding the ratio of $t_{\mathrm{p}}$ to $t_{\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{n}}$

When this ratio is less than 1, $t_{\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{n}}<t_{\mathrm{p}}$ and alignment occurs before disc precession, preventing disc breaking. This criterion is primarily applicable to discs in the diffusive regime, although the factor of $(h/r{)}^2/\alpha$ suggests that it may apply to some discs in the wave-like regime as well. Equation (15) is a simple estimate that only considers the disc scale height at one particular radius; other approaches may take into account how the scale height changes radially across the disc. Various disc breaking criteria have been suggested in the past. Equation (13) of King et al. (2013) considers $t_{{\nu }_2}{\omega }_{\mathrm{p}}<1$ as the criterion to prevent large disc warping in the viscous regime, where $t_{{\nu }_2}$ is the viscous time-scale associated with the ${\nu }_2$ viscosity. We note that in the linear regime of small warps, this criterion and the equation (15) give opposite stability results because $({\omega }_{\mathrm{p}}{t}_{\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{n}}{)}^{-1}\sim {\omega }_{\mathrm{p}}{t}_{\nu 2}$⁠. A short (long) alignment time-scale occurs with a long (short) warp viscous time-scale. However, the ${\nu }_2$ value for small warps may not apply for cases of disc breaking.

2.3 An application to GW Orionis

GW Orionis (GW Ori) is a hierarchical triple star system surrounded by a large protoplanetary disc (Mathieu, Adams & Latham 1991; Berger et al. 2011). Resolved images of the system (Bi et al. 2020; Kraus et al. 2020) show a series of dust rings in differing orientations around the central star system as well as distorted velocity maps in ${}^{12}\mathrm{C}\mathrm{O}$⁠, suggesting that the disc is warped. Observations indicate that the outer disc is inclined with respect to the outer binary orbital plane at an angle of roughly $i={38}^{\circ }$⁠.

The GW Ori system has been simulated in previous works, but simulations have shown differing results as to whether the disc will break solely under the influence of the inner star system. In their observational paper, Kraus et al. (2020) use SPH simulations to explain the warped geometry of the system. They simulate the GW Ori disc with an initial radial range of 40–$200$ au, using power-law slopes of $p=0.2$ and $s=1.0$ and a very thin disc with $h/r=0.02$⁠. The disc expands inwards until the inner radius reaches about 30 au, where a thin ring breaks off from the inner edge of the disc and precesses independently, eventually aligning with the binary plane. Kraus et al. (2020) interpret this inner ring as a potential origin of the R3 dust ring seen in observations.

Bi et al. (2020) and Smallwood et al. (2021) ran a similar suite of SPH simulations for the GW Ori system, using a disc with the same initial inner edge but a slightly thicker aspect ratio of $h/r=0.05$ and power-law slopes of $p=1.5$ and $s=1.0$⁠. Their simulations find a disc that warps strongly but does not break through disc precession alone. They suggest that a giant planet embedded in the disc is required to carve a gap first in order for an inner ring to break off. They also find that a disc with a smaller inner radius, $r_{\mathrm{i}\mathrm{n}}=20\mathrm{a}\mathrm{u}$⁠, can spontaneously generate a break, with a breaking radius of roughly 75 au. However, note that this is a result of the initial surface density profile, which is a power law in radius with a sharp cut-off at the initial disc inner edge. The disc is not in equilibrium and over time material is cleared from this region as a result of the binary torque. The set-up is therefore somewhat artificial. Despite these problems, they determine that the geometry of the disc is an important factor in determining whether a disc is able to break.

By applying equations (12) and (13) to the GW Ori system, it is possible to explain the appearances of both SPH simulations. To model this triple star system with our analytical equations, we assume that the triple star orbit can be approximated by the outer binary orbit of the triple. This assumption is justified for a low binary eccentricity (e.g. Smallwood et al. 2021; Lepp, Martin & Lubow 2023). The components of the outer binary are composed of two masses with $M_1=0.742M_{\mathrm{t}\mathrm{o}\mathrm{t}}$ and $M_2=0.258M_{\mathrm{t}\mathrm{o}\mathrm{t}}$ and an eccentricity of $e_{\mathrm{b}}=0.379$⁠. Including the effect of the inner binary will cause apsidal precession in the orbit of the outer binary (Lepp et al. 2023), but we do not consider this effect when considering the precession of the disc. For the circumtriple disc, we consider the two disc models with the parameters listed above. We use the parameters from two disc models listed above, with the size of $r_{\mathrm{i}\mathrm{n}}$ varying from 40 to 20 au to simulate the inward drift of the disc.

For the Kraus et al. (2020) simulation, the initial conditions have $t_{\mathrm{p}}>t_{\mathrm{c}}$ everywhere, indicating a disc with no breaks. However, the two curves approach as the inner edge of the disc drifts inwards. Once the inner edge of the disc reaches about 21 au, the two curves intersect and a break spontaneously appears at $r_{\mathrm{b}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{k}}=49\mathrm{a}\mathrm{u}$⁠, producing a thin secondary ring along the inner edge. By the time the inner cavity has shrunk to $20\mathrm{a}\mathrm{u}$⁠, the breaking radius has also moved inwards to a distance of $r_{\mathrm{b}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{k}}=35\mathrm{a}\mathrm{u}$⁠, a distance that is consistent with the observed location of the ring seen in the SPH simulations. Fig. 2 shows the evolution of $t_{\mathrm{p}}$ and $t_{\mathrm{c}}$ in the Kraus disc model. If spontaneous disc breaking occurs in this manner, the material in the inner ring can be shuttled inwards into a tight ring similar to what is seen in the GW Ori observations.

The discs used by Bi et al. (2020) and Smallwood et al. (2021) are thicker and have a higher disc temperature, creating smaller values for $t_{\mathrm{c}}$⁠. Although the two curves begin to intersect at nearly the same value of $r_{\mathrm{i}\mathrm{n}}$⁠, the change in the shape of the curves causes the value of $r_{\mathrm{b}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{k}}$ to start from the outer edge of the disc and sweep inwards, instead of starting in the middle like the Kraus et al. (2020) disc. The inner cavity must shrink even more, down to $r_{\mathrm{i}\mathrm{n}}$  = 15 au, in order for the break to produce a similarly confined inner ring. This model is unlikely to break the disc using the binary torque alone, making a planet-driven explanation more likely for these disc parameters.

By comparing these simulations, it is clear that the details of the simulation parameters are important in determining the location of a disc break. Observations of the GW Ori system are currently unable to put strong constraints on the disc parameters, so the GW Ori disc may be able to break under the torque of the binary if it is closer to the Kraus et al. (2020) simulation, i.e. cooler and with lower values of $h/r$⁠. In Section 3.2, we describe our simulations for the GW Ori system.

2.4 Multiple disc breaking

After a disc breaks into two, the inner and outer discs will begin to precess independently. The inner disc is guaranteed to precess rigidly, by the conditions set by the break, but the outer disc may once again precess differentially with its own $t_{\mathrm{c}}$ and $t_{\mathrm{p}}$ starting from the breaking radius, which is its new inner edge. If the breaking radius of this new outer disc is within the radial extent of the outer disc, then the outer disc may break again. This process of multiple disc breaking can repeat and break the disc into multiple rings for as long as there is disc material, or until the estimates for $t_{\mathrm{c}}$ and $t_{\mathrm{p}}$ above are no longer valid. Some SPH simulations have observed this phenomenon previously, in the context of black hole binaries and the Bardeen–Petterson effect (Nixon et al. 2013; Nealon et al. 2015).

To get a better understanding of how the various disc parameters affect the number of disc breaks, in Fig. 3 we map out the number of disc breaks for different initial values of $r_{\mathrm{i}\mathrm{n}}$⁠, $r_{\mathrm{o}\mathrm{u}\mathrm{t}}$⁠, and $(h/r{)}_0$ (the disc geometry) and as a function of the power-law slopes p and s (the disc structure). In this figure, the disc is placed around an equal-mass $(M_1=M_2=0.5),e_{\mathrm{b}}=0.5$ central binary. Each panel shows the number of expected disc breaking events for the given initial disc geometry as a function of p and s. In each panel, the horizontal line at $s=1$ represents discs with a constant $h/r$ throughout; points below this line represent flaring discs, with $h/r$ increasing with radius. In general, discs that are thin have small inner cavities, and steep profiles are more likely to have multiple breaks.

The effect of changing the inner radius can be seen by comparing the first and second panels. Doubling the inner radius reduces the number of breaks in each region by roughly 1 each. This behaviour is consistent with the disc breaking seen in the simulations of Smallwood et al. (2021), who also noticed that disc breaking was dependent on the disc inner radius. Fig. 3 also shows that increasing the disc scale height increases the width of each ‘stripe’, requiring a steeper temperature gradient to achieve multiple breaks. For typical protoplanetary disc parameters, disc breaking will usually be limited to a single break unless the disc is extremely thin. However, multiple breaking may still occur naturally when considering other sources of precession, such as around black holes. Tidal interactions from the binary component can truncate the disc edges, providing limits on the values of $r_{\mathrm{i}\mathrm{n}}$ and $r_{\mathrm{o}\mathrm{u}\mathrm{t}}$ (Larwood et al. 1996; Lubow, Martin & Nixon 2015; Miranda & Lai 2015). Miranda & Lai (2015) study tidal truncation of circumbinary discs due to Lindblad torques and find that for highly inclined discs, the inner radius is truncated close to the 1:3 and 1:4 commensurabilities, approximately from 2.0a_b to 2.5$a_{\mathrm{b}}$⁠.

3 METHODS

We use the grid-based code athena++ (Stone et al. 2020) to simulate an inclined disc around a central binary. To study the various effects of disc warping and breaking, we run three different sets of simulations, each with slightly different disc and binary set-ups. Each set is described in the following sections. The disc and binary parameters used for each simulation set are listed in Table 1. Below, we describe features common to all of the simulations performed in this paper.

We choose spherical polar coordinates $(r,\theta ,\varphi )$ for our simulations. Our simulation domain extends from $5$ to ${175}^{\circ }$ in $\theta$⁠, and 0 to $2\pi$ in $\varphi$⁠, and is divided into a grid of $176\times 384$ cells in $(\theta ,\varphi )$⁠. We use logarithmically spaced cells in the radial domain. The extent and number of cells in the radial domain r vary depending on the simulation set. At the simulation boundaries, we set a one-way outflow boundary condition in the radial direction, and reflecting boundary conditions in the polar (⁠$\theta$⁠) direction. The choice of boundary conditions in the radial direction allows bending waves in the disc to dissipate without reflection at the disc edges. Though the boundary conditions in the polar direction are reflecting, the large extent of the polar domain and settling of the disc to the simulation mid-plane keeps the disc far away from these boundaries.

We initialize the density profile of the disc by numerically integrating the density at each grid cell to establish hydrostatic equilibrium according to the power-law profile in equation (2), and set the vertical profile by numerically integrating the density at each grid cell to establish hydrostatic equilibrium according to the disc scale height $h=c_{\mathrm{s}}/{\mathrm{\Omega }}_{\mathrm{K}}$⁠. Here, $r_0=a_{\mathrm{b}}$⁠, ${\rho }_0=1$ at $r_0$⁠, $c_{\mathrm{s}}$ is the local isothermal sound speed $\sqrt{P/\rho }$ at r, and ${\mathrm{\Omega }}_{\mathrm{K}}$ is the Keplerian orbital frequency at r. The disc temperature is initialized using the power-law profile in equation (3), where $T_0=1$ at $r_0$⁠. We truncate the edges of the disc using an exponential cut-off of the form $\exp [(r-{\mu }_r)/{\sigma }_r]$⁠, where ${\mu }_r$ is the location of the disc edge and ${\sigma }_r$ is the relative scale length of the cut-off. We also use a spherically symmetric density floor with a value of ${\rho }_{\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{o}\mathrm{r}}={10}^{-4}{\rho }_0$ at $r=1$ and a power-law slope of d, identical to that of the density profile. We initialize the disc with an initial scale height of $(h/r{)}_0$ at $r=r_0$⁠. We use the orbital cooling scheme outlined in equation (5) of Zhu et al. (2015), using a dimensionless cooling time of $t_{\mathrm{c}\mathrm{o}\mathrm{o}\mathrm{l}}=0.01{\mathrm{\Omega }}_{\mathrm{K}}^{-1}$⁠.

The binary components are simulated as gravitational bodies with masses $M_1$ and $M_2$⁠, which orbit with semimajor axis $a_{\mathrm{b}}$ and eccentricity $e_{\mathrm{b}}$⁠. We integrate the motion of the binary components using a second-order leapfrog integrator, but we do not model the change in the orbits due to interaction with the disc material.

3.1 Polar-aligning warped discs

Our first set of simulations studies the effects of disc warping in a polar-aligning disc around an eccentric binary. As described in Section 2, polar-aligning discs benefit from faster precession times due to smaller central cavities and the response from the central binary, making them ideal for studying disc breaking. Our set-up is similar to that used in Rabago et al. (2023), which we review in brief here. The binary components are initialized as equal-mass particles $(M_1=M_2=0.5)$ with $a_{\mathrm{b}}=0.28$ and $e_{\mathrm{b}}=0.5$⁠, placing the binary as close to the inner radial domain as possible. We simulate the region from $0.3=1.07a_{\mathrm{b}}$ to $10.0=35.7a_{\mathrm{b}}$ using 216 cells in the radial domain. The initial scale height is set at $(h/r{)}_0=0.103$⁠, giving the disc a scale height of $h/r=0.075$ at $r=3.5a_{\mathrm{b}}$⁠. We use power-law slopes of $d=2.25$ and $s=1.5$$(p=1.5)$ and truncate the inner edge of the disc using ${\mu }_r=2a_{\mathrm{b}}$ and ${\sigma }_r=0.35a_{\mathrm{b}}$⁠. The outer edge of the disc extends to the outer radial domain. We vary the $\alpha$-viscosity parameter from the values $\alpha ={10}^{-1},{10}^{-2},{10}^{-3}$⁠, and ${10}^{-5}$⁠. These values correspond to discs with $h/r<\alpha$ (diffusive regime), $h/r\sim \alpha$ (intermediate case), $h/r>\alpha$ (wave-like regime), and a case for the inviscid limit. We place the binary along the simulation $xz$-plane and initialize the disc with a ${60}^{\circ }$ inclination to the binary orbital plane. We run these simulations for 1000 binary orbits.

3.2 GW Orionis

We create a suite of simulations that replicates the GW Ori system. We use the known binary parameters from Kraus et al. (2020), using the AB-C binary separation of $a_{\mathrm{b}}=8.89$ au as the system scale. Previous simulations of the GW Ori system from Smallwood et al. (2021) found that the motion of the inner AB binary provides smaller effects to the disc compared to the larger motion of the AB-C binary, so we model the GW Ori triple system as a binary to simplify calculations as in Section 2.3. Our binary model consists of the AB binary as one mass $M_{\mathrm{AB}}=0.742M_{\mathrm{t}\mathrm{o}\mathrm{t}}$ and the outer C component as the second mass $M_{\mathrm{C}}=0.258M_{\mathrm{t}\mathrm{o}\mathrm{t}}$⁠, placed in an orbit with semimajor axis $a_{\mathrm{b}}=8.89$ au and eccentricity $e_{\mathrm{b}}=0.379$⁠.

We simulate the region from $1.5a_{\mathrm{b}}$ to $50a_{\mathrm{b}}$⁠, covering a range of roughly 13–450 au. We use 192 cells across the radial domain. The disc is inclined at an angle of ${38}^{\circ }$ relative to the binary. We use power-law exponents of $d=2.5$ and $s=1.0$⁠. This gives the disc a surface density profile of $p=1.5$ and a constant $h/r$ throughout, which we choose to be $(h/r{)}_0=0.075$⁠. We truncate the inner edge of the disc at ${\mu }_r=24$ au = $2.7a_{\mathrm{b}}$ and use a disc viscosity of $\alpha ={10}^{-5}$⁠. We expect this simulation to be similar to the Smallwood et al. (2021) simulation with a small inner radius, generating a break at around $10a_{\mathrm{b}}$⁠. Although the disc viscosity is lower than their simulations, we expect the inner edge of the disc to be truncated at a similar radius due to the torque of the binary and, with $\alpha <h/r$⁠, the disc remains in the wave-like regime and should break in a similar fashion. We run these simulations for 5000 binary orbits, enough time for the inner cavity to settle and the warp to propagate through the disc. We do not attempt to reproduce the disc from Kraus et al. (2020), due to the thin disc scale height of $h/r=0.02$ used in their work. Resolving this disc at the required inclination would require a large increase in the overall grid resolution, which we consider prohibitively expensive.

3.3 Multiple breaks

As described in Section 2.4, the equations for $t_{\mathrm{p}}$ and $t_{\mathrm{c}}$ suggest that a disc can undergo multiple breaks if the power-law slope of $t_{\mathrm{p}}$ is less than that of $t_{\mathrm{c}}$⁠. This can occur given the right combination of disc geometry and density profiles. To examine this possibility, we run an additional simulation with a disc that is predicted to undergo multiple breaks by equations (7) and (8). To maximize the chance of multiple breaks, we choose a polar-aligning disc with steep density and temperature profiles.

We modify our set-up for the polar disc simulations, changing the radial domain so that the inner domain remains at $r_{\mathrm{i}\mathrm{n}}=1.07a_{\mathrm{b}}$ but extending the outer radial domain to $r_{\mathrm{o}\mathrm{u}\mathrm{t}}=100a_{\mathrm{b}}$ using 272 radial cells across the radial domain. We use power-law profiles of $d=4.0$ and $s=2.5$$(p=3.75)$⁠. We choose the disc scale height to be $(h/r{)}_{r=5a}=0.075$⁠. Combined with the temperature profile, the disc scale height ranges from $(h/r{)}_{\mathrm{i}\mathrm{n}}=0.251$ at the disc inner edge to $(h/r{)}_{\mathrm{o}\mathrm{u}\mathrm{t}}=0.013$ at $r=50a_{\mathrm{b}}$⁠. To ensure that the disc is resolved as the scale height decreases, we use a single level of mesh refinement to refine the regions between $r=[8.89a_{\mathrm{b}},52.1a_{\mathrm{b}}]$ and $\theta =[{50.8}^{\circ },{129.1}^{\circ }]$⁠. The disc viscosity is set to $\alpha ={10}^{-5}$⁠, as in the nearly inviscid cases.

A disc with these parameters is predicted to undergo either two or three breaking events, depending on the exact location of the inner radius. The location of this set-up is marked on each panel of Fig. 3 with a ‘$\times$’. For the left panel, with $r_{\mathrm{i}\mathrm{n}}=2.5a_{\mathrm{b}}$⁠, the breaks are predicted to occur at distances of 4.0a_b, 8.6a_b, and 29.1$a_{\mathrm{b}}$⁠. We evolve this disc for 1500 binary orbits, enough time for the outer parts of the disc to evolve for a few orbital time-scales.

4 RESULTS

4.1 Disc warping and breaking

Our polar-aligning disc simulations exhibit both diffusive and wave-like behaviours, depending on the choice of $\alpha$⁠. Fig. 4 shows the various discs at $t=1000$ binary orbits, displaying varied behaviours in their evolution. For the $\alpha ={10}^{-1}$ simulation, the disc develops a strong warp that evolves diffusively, while for the $\alpha ={10}^{-2},{10}^{-3},$ and ${10}^{-5}$ simulations, the disc breaks and exhibits wave-like dispersion. Fig. 5 shows the time evolution for both regimes, following the evolution of the z-component of the unit angular momentum vector. In all simulations, an initial warp is excited in the inner disc by the binary precession in the first $\sim 200$ binary orbits. For the diffusive case, the inner regions of the disc quickly align towards polar due to the short alignment time-scale, and the initial warp spreads outwards in the disc diffusively. For the wave-like cases, the initial warp develops into a sharp discontinuity, characteristic of a disc breaking event. After the disc breaks, a smaller disturbance propagates through the outer disc in a wave-like manner.

To determine the location for the warp and follow its movement within the disc, we track the position of maximum $\psi$ as the disc evolves. In the case of a disc that breaks, two large peaks in $\psi$ are present, one at the breaking radius and another at the bending wave in the outer disc. We identify the radial location of these features by searching for the largest two maxima in $\psi$ to determine the locations of the disc breaking radius and the warp in the outer disc. The warp location and propagation speed for our warp simulations are shown in Fig. 6.

For simulations in the diffusive regime, the inner disc generates a large warp as the disc angular momentum aligns to the precession axis. The warp creates a large peak in $\psi$ that travels outwards at roughly 0.2 times the local sound speed. The left panels of Figs 5 and 6 show that the warp starts as a very steep and local feature, but broadens over time, encompassing a range of approximately $10a_{\mathrm{b}}$ by the end of the simulation. A small wave is launched ahead of the main warp (visible in Fig. 6 as a series of small secondary peaks between 10a_b and 25$a_{\mathrm{b}}$⁠, beneath the larger curves drawn at later times). This wave does not appear with a similar feature in Fig. 5. The highly diffusive nature of this disc quickly removes this feature, causing it to disappear by $t=500T_{\mathrm{b}}$⁠.

Simulations in the wave-like regime show the initial disc breaking at $r_{\mathrm{b}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{k}}\sim 4a_{\mathrm{b}}$⁠, independent of the small disc viscosity. After the initial breaking event, the breaking radius moves outwards, faster for the intermediate case $\alpha =0.01$⁠. The bottom panel of Fig. 6 shows that the outward warp decreases in velocity as it moves outwards in the disc, the propagation speed following close to one-half of the local sound speed in the disc.

We examine the drift of the breaking radius in Fig. 7. Each panel shows the location of the breaking radius for a simulation with a particular $\alpha$-viscosity, with curves plotting the normalized surface density and vertical dashed lines marking the location of the breaking radius at each time. Outward drift of the break is seen for the $\alpha ={10}^{-2}$ and ${10}^{-5}$ simulations, but not for the $\alpha ={10}^{-3}$ simulation. For the $\alpha ={10}^{-5}$ simulation, the break forms at roughly $r_{\mathrm{b}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{k}}=4.25a_{\mathrm{b}}$ and slowly drifts outwards until $t=800T_{\mathrm{b}}$⁠, where the break appears to stabilize at a final location of $5.1a_{\mathrm{b}}$⁠. The $\alpha ={10}^{-2}$ simulation breaks the disc slightly farther out, at $r_{\mathrm{b}\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{k}}\sim 4.65a_{\mathrm{b}}$⁠, and drifts outwards at a rate nearly twice as fast as the $\alpha ={10}^{-5}$ simulation. By the end of the simulation, the break has nearly reached $7a_{\mathrm{b}}$⁠, and does not show any signs of slowing down. In both simulations, the vertical lines representing the location of breaking radius are unevenly spaced, suggesting that the breaking radius may not be growing constantly. Large growth events may be related to the times when the inner and outer discs are momentarily aligned, allowing for increased accretion between the two discs.

Our analytical equations for disc breaking show close agreement with the location of the break observed in the numerical simulations. In Fig. 8, we compare the location of the break in our $\alpha ={10}^{-5}$ simulation to the location predicted by equations (12) and (13). The breaking radius is visible as a sharp drop in the disc surface density. We show the disc surface density at $t=100$ binary orbits, just after the initial break has formed. Comparing the location of the break at early times is important for accurately determining the location of the breaking radius, since the break slowly drifts outwards in all of our simulations.

To determine the value of $r_{\mathrm{i}\mathrm{n}}$ for our analytical equations, we compare the density-weighted angular momentum of the disc $\int \mathrm{\Sigma }{r}^3\mathrm{\Omega }\mathrm{d}r$ to the analytical value expected for a disc with the same surface density power-law profile. For this calculation, we choose a snapshot shortly after the initialization, when the disc has had time to equalize but before a break has had time to open up. We consider $r_{\mathrm{i}\mathrm{n}}$ to be the location at which these two values diverge. For the polar-aligning discs, this gives a value of $r_{\mathrm{i}\mathrm{n}}=2.2a_{\mathrm{b}}$⁠.

We also observe substructures in the individual discs created after the breaking event. Fig. 9 shows the $\alpha ={10}^{-5}$ simulation as it would appear ‘flattened’ on to a single plane, as well as the disc eccentricity and vorticity. Simulations with $\alpha ={10}^{-3}$ and $\alpha ={10}^{-2}$ show similar details in the gap regions, but no substructures in the inner disc. We calculate the disc eccentricity and vorticity in a similar manner to that of Rabago et al. (2023). With low disc viscosities, the inner disc develops a large Rossby wave instability (RWI) vortex (Lovelace et al. 1999; Li et al. 2000), as well as a localized overdensity and a single-armed spiral. This feature is identical to the vortices seen in the polar discs of Rabago et al. (2023). We observe less features in the outer disc; some faint two-armed spirals are visible at certain times, and no vortices are present. The outer disc has undergone less dynamical times, and so may not have had enough time to develop a large vortex at its inner edge. The gap edge between the inner and outer discs may not be as steep as the edge of the inner disc truncated by the central binary, which may also inhibit the outer disc from forming a strong vortensity minimum and the growth of RWI vortices (Bae, Hartmann & Zhu 2015).

4.2 GW Orionis

Fig. 10 shows the disc surface density and warp strength of the simulation representing the GW Ori system (Section 3.2). The disc develops a warp within the first 1000 binary orbits. However, no disc breaking event occurs, even after the next several thousand orbits of the simulation. A light warp persists in the disc throughout the simulation time, but is mostly confined to the inner regions of the disc, and never grows large enough to become unstable and evolve into a full break. We run tests with different $\alpha$ viscosities, simulating out to at least 1000 binary orbits, but none of the simulations are able to create a break in the disc.

The arrangement of the GW Ori binary changes the strength of the induced precession. The terms in equation (11) suggest that the precession induced by the binary is strongest for polar discs around high-eccentricity, equal-mass binaries. The GW Ori system has a moderate binary eccentricity of $e_{\mathrm{b}}=0.379$ and a mass ratio, $q=M_{\mathrm{C}}/M_{\mathrm{AB}}=0.348$⁠, making the precession rate normalized by ${\mathrm{\Omega }}_{\mathrm{b}}$ of this system about one-fifth as strong as that induced by the polar discs in Section 4.1.

In our simulations, the binary torque truncates the inner edge of the disc at roughly $r_{\mathrm{i}\mathrm{n}}=3.5a_{\mathrm{b}}$⁠. This is a larger value of $r_{\mathrm{i}\mathrm{n}}$ than observed in the Smallwood et al. (2021) simulations, and larger than the value we initially used for calculating the breaking radius above. The size of the inner cavity remains larger than the Smallwood et al. (2021) simulations, even for runs with the same $\alpha$-viscosity. Recalculating the location of the breaking radius with the value of $r_{\mathrm{i}\mathrm{n}}$ seen in the simulations changes the disc time-scales so that no break occurs in the disc.

Although we try to reproduce the simulations of Smallwood et al. (2021) with a small inner cavity as best as possible, using nearly identical parameters for both the disc and binary, we are unable to reproduce the observed disc breaking using our grid-based simulations. Our use of a slightly larger scale height due to computational limitations decreases the value of $t_{\mathrm{c}}$⁠, but the analytical equations we present in Section 2 predict that we should still observe a break as long as the inner edge of the disc is truncated in the same location as the Smallwood et al. (2021) simulation. As described in Section 2.3, the initial surface density profile of the Smallwood et al. (2021) simulation is artificial; too much material is initialized close to the binary with a sharp cut-off, allowing the inner disc to precess and break before the inner cavity is fully cleared. We were unable to replicate their initial density profile exactly in our grid-based simulations, as the inner cavity of our disc is a gradual truncation starting at $r_{\mathrm{i}\mathrm{n}}=24$ au in order to ensure simulation stability. This slightly increases the value of $t_{\mathrm{p}}$ and may be enough to prevent breaking. It is also possible that differences in the code formulation may play a role in how the initial stages of the disc evolve, when the cavity depletes and the disc is most likely to break.

4.3 Multiple breaks

A visualization of our multiple break simulation (Section 3.3) is shown in Fig. 11. Two disc breaking events are visible, separating the disc into three distinct rings. At the end of the simulation, the disc breaks are located at radial distances of 7a_b and $28a_{\mathrm{b}}$⁠, with the outer edge of the disc reaching equilibrium at roughly $r_{\mathrm{o}\mathrm{u}\mathrm{t}}=60a_{\mathrm{b}}$⁠. The third ring shows a strong warp at roughly $35a_{\mathrm{b}}$⁠. This warp is strong enough for the disc to appear broken due to how thin the disc is at this radius, but subsequent analysis will show that the disc is still communicating across this region. Similar to the polar discs in Section 4.1, we see that RWI vortices form in the innermost disc, but not in the second or third discs.

Fig. 12 compares the observed locations of each break to the locations predicted by the time-scale equations. We find that it is difficult to determine a suitable inner radius for this disc. Choosing an inner boundary of $r_{\mathrm{i}\mathrm{n}}=2.5a_{\mathrm{b}}$ our analytical equations predict three breaks, with the second and third breaks in rough agreement with the locations of the breaks observed in the simulation. The first break, expected at $r=4.25a_{\mathrm{b}}$⁠, is not observed in the innermost disc. Instead, the observed breaks may be roughly fitted by choosing $r_{\mathrm{i}\mathrm{n}}=4.0a_{\mathrm{b}}$⁠, but this places the ‘inner edge’ of the disc much farther out, in a region where the local surface density is higher than the initial condition.

To better identify the regions in the disc that are communicating, we calculate the azimuthal angle of the shell-averaged angular momentum vector, defined as

In our multiple break simulation, the z direction is along the binary eccentricity vector and the binary orbit is in the $x-z$ plane, so ${\varphi }_L$ can be viewed as the longitude of ascending node of the disc in the $x-y$ plane as it circulates around the binary eccentricity vector. Regions of the disc in radial communication will precess together, and thus have roughly the same value of ${\varphi }_L$⁠. In Fig. 13, we plot the value of ${\varphi }_L$ across the disc as a function of simulation time. As the disc evolves and breaks off rings, distinct horizontal bands appear, with each band corresponding to a separate disc in the final simulation. Disc breaking events are visible as sharp boundaries between two horizontal bands, as the value of $\mathit{L}$ changes sharply from one disc to another.

In Fig. 13, it can be seen that the innermost disc breaks off almost immediately, in about 100 binary orbits. The fast precession of the inner disc causes ${\varphi }_L$ to cycle several times during the rest of the simulation, and creates the first band at the bottom of the figure. The outer disc then develops a warp that gradually extends outwards to about $25a_{\mathrm{b}}$⁠, before breaking off a second disc at roughly $1000T_{\mathrm{b}}$⁠. After the disc breaks at $25a_{\mathrm{b}}$⁠, the warp between the first break and the second break gradually vanishes and this disc region becomes coplanar. This is another indication that the disc loses the communication at the breaking radius and the outer disc cannot torque the inner disc any more. A faint band can be seen extending from the second disc to about $35a_{\mathrm{b}}$ in the outermost disc, signifying a strong warp. However, since a sharp discontinuity has not yet formed, this band is not yet considered a third breaking event.

The multiple break simulation showcases the difficulty in using equations (12) and (13) to predict the exact location of the breaking radii when multiple breaks are packed together in a small radial range. Breaking events require a few local dynamical time-scales to fully develop, so the innermost break will have time to drift outwards while subsequent breaks form and develop. The bottom panel of Fig. 12 shows that the surface density, normalized to an initial power-law distribution, is quite disturbed even in regions far away from the breaks, suggesting that the disc no longer follows a power-law distribution even though the analytical equations assume one when calculating the location of a break. We discuss this behaviour further in Section 5.1.

5 DISCUSSION

For the polar discs used in Section 3.1, we show the effects of disc viscosity by plotting $\tau ={\mathit{t}}_{\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{n}}/{\mathit{t}}_{\mathrm{p}}$ as a function of radius in Fig. 14, with the assumption that while calculating $t_{\mathrm{p}}$ all material interior to r is able to precess coherently. We see that for our selected disc parameters, the discs with $\alpha \lesssim 0.01$ maintain $\tau >1$ for the entirety of the disc and thus are able to break under the induced precession of the binary.

We can estimate ${\nu }_2$ and ${\nu }_3$ for our $\alpha =0.1$ simulation using figs 4 and 5 of Ogilvie (1999). From our Fig. 6, the maximum warp amplitude is roughly ${\psi }_{\mathrm{m}\mathrm{a}\mathrm{x}}=0.5-0.75$⁠, at a distance of $r=10a_{\mathrm{b}}$⁠. This corresponds to warp viscosities of ${\nu }_2=0.019a_{{}_{\mathrm{b}}}^2{\mathrm{\Omega }}_{\mathrm{b}}$⁠, ${\nu }_3=0.0013a_{{}_{\mathrm{b}}}^2{\mathrm{\Omega }}_{\mathrm{b}}$⁠, and $t_{{\nu }_2}=5000{\mathrm{\Omega }}_{\mathrm{b}}^{-1}$⁠.

We note that, for the simulations in Section 4.1, there is a change in the overall disc orientation of about 5–10° from its initial orientation. This change occurs faster than the local precession time-scale induced by the binary, and instead occurs as the warp (diffusive or wave-like) propagates through the disc. This effect is partially visible in Fig. 4, since the disc, initialized in alignment with the simulation x-axis, is now misaligned at the end of the simulation. We consider this to be an effect of the ${\nu }_3$ warp viscosity, which has caused the entire disc to precess.

5.1 Efficacy of the disc breaking equations

The location and number of breaking radii determined by equations (12) and (13) are a sensitive function of the disc parameters, particularly the inner disc radius $r_{\mathrm{i}\mathrm{n}}$⁠. Because of this sensitivity, small adjustments in the value of $r_{\mathrm{i}\mathrm{n}}$ result in large changes in the breaking radius. Conversely, small adjustments in the location of the breaking radii only require minor changes in the required value of $r_{\mathrm{i}\mathrm{n}}$⁠. An important caveat of this behaviour is that it easily allows our description of disc breaking to remain consistent with observed discs, but it is harder for the same equations to act as predictive tools. For example, equations (12) and (13) can be made consistent with simulations or observations by finely adjusting the disc parameters to match the observed location(s) of the breaking radius, but it is more difficult to predict the location of emergent breaks starting from the initial disc conditions.

When examining our polar disc simulations in Section 4.1, we initially find that the expected breaking radii are different from our initial guess as the inner radius shifts and reaches a quasi-steady state. Once we adjust the value of $r_{\mathrm{i}\mathrm{n}}$ to account for this change, the predicted breaks are in much closer agreement. For the multiple break simulation (Section 4.3), we have a harder time trying to match both the observed breaking locations and the location of the inner radius simultaneously, especially when multiple breaks are packed in a small radial range.

The nature of this sensitivity leads to some difficulty in using equations (12) and (13) to predict the location of disc breaks. There is no clear location at which the gas at the inner radius of the disc acts as an ‘edge’ and stops precessing as a solid body with the rest of the disc. Thus, different methods of calculating the inner radius such as using surface density cut-offs, precession rates, angular momentum deviations, or other methods will result in slightly differing values of $r_{\mathrm{i}\mathrm{n}}$⁠, leading to greatly differing estimates of the breaking radii in the disc. Differences in $r_{\mathrm{i}\mathrm{n}}$ caused by code formulation, such as those present in the GW Ori simulations in Section 4.2, can further complicate this issue. This behaviour was also noted in the simulations by Young et al. (2023), who comment that the wide dependence on disc and binary parameters makes the creation of a disc breaking equation a challenging matter.

Our analytical equations are also restricted by their lack of full coverage of the possible parameters. Previous SPH simulations have observed a dependence of disc breaking with inclination to the precession vector (Facchini, Lodato & Price 2013; Nealon et al. 2015), suggesting that disc breaks are more likely for highly inclined or retrograde discs, and may not occur if the mutual inclination between the disc and precession vectors (prograde or retrograde) is small. Analyses by Doǧan et al. (2018) and Raj et al. (2021) suggest that a warp may continue to grow unstably if the warp amplitude $\psi$ is beyond a critical limit ${\psi }_{\mathrm{c}}$⁠, dependent on the disc viscosity. The current analysis is relatively simple and is strictly analytical, using the initial power-law profiles set at the beginning of the simulation to calculate the breaking radius. As the binary–disc system evolves, the gas may settle towards different power-law profiles than those set by the initial conditions, or towards a radial distribution not approximated by a power-law profile, as seen in Fig. 12. An example of this is shown in Fig. 15, where we compare the analytical, power-law disc time-scales in our multiple break simulation to the time-scales calculated by directly integrating equation (10) (since our simulations are nearly isothermal, direct integration for $t_{\mathrm{c}}$ is identical to the analytical result of equation 13). There is a slight outward shift in the location of the breaking radius in this case, but the new location is similar to the location predicted by the analytical equations. However, this adjustment may lead to cascading outward shifts for the subsequent breaks farther out in the disc.

We can estimate the effects of the disc inclination on our results by examining the analytical equations of Farago & Laskar (2010) describing the orbits of test particles around eccentric binaries. We determine the correction factor for the particle precession period due to inclination away from nearly polar or coplanar orbits. This same correction factor also applies to a flat disc, since the factor is independent of radius. For a circumbinary test particle, the precession period of the particle is given in equation (2.32) of Farago & Laskar (2010) as

where $n_1$ is the mean motion of the binary (⁠$n_1=1$ in our set-up), and the orbital elements $a_{\mathrm{p}}$ and $e_{\mathrm{p}}$ refer to the semimajor axis and eccentricity of the orbiting particle, respectively. The constant $h=-4e_{\mathrm{b}}^2+{\sin }^{2}i(1+4e_{\mathrm{b}}^2)$ describes the orientation of the test particle’s orbit, and $K(k^2)$ is an elliptic integral of the first kind integrated over the azimuthal angle $\varphi$⁠, defined as

where $k^2=\frac{5e_{\mathrm{b}}^2}{1-e_{\mathrm{b}}^2}{\cot }^{2}i$ and ${\varphi }_0=\arcsin \left(1/k\right)$⁠. The dependence of K and h on the inclination angle causes the value of T to vary with changing disc inclination.

The effects of the change in precession time are shown as a function of disc inclination in Fig. 16 for the polar-aligning set-up in Section 4.1 and the GW Ori set-up in Section 4.2. The precession time lengthens considerably only in the vicinity of the critical inclination, and so we consider the effect on our simulations to be relatively small. The precession period diverges logarithmically with the difference in inclination from the critical inclination. Importantly, Fig. 16 shows that for discs within a few degrees of the critical inclination the disc will precess very slowly and is therefore unlikely to break. Simultaneously, inclinations close to 0 or ${90}^{\circ }$ are also unlikely to break, since they are aligned close to their precession vectors and will not produce a strong warp amplitude. Thus, certain inclinations between $0^{\circ }$ and $i_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}$ or $i_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}$ and ${90}^{\circ }$ are the most likely orientations for disc breaking.

The dependence on inclination derived here is different from that expected from Nixon et al. (2013), who predict that breaking is most efficient for inclination angles ${45}^{\circ }<i<{135}^{\circ }$ and is strongest at $45$ and ${135}^{\circ }$⁠. This criterion is made by considering the cancellation of angular momentum from counter-aligned rings produced as the inner disc precesses. Although this range of inclination angles is correct when considering circular binaries and Lense–Thirring precession induced from black holes, the addition of librating orbits for eccentric central binaries produces new regions at potentially moderate inclinations where disc breaking is unlikely.

5.2 Observational signatures

Our simulations of GW Ori in Section 4.2 were conducted using a binary–disc inclination of ${38}^{\circ }$⁠, as listed in previous works (Bi et al. 2020; Kraus et al. 2020; Smallwood et al. 2021). After the simulations were conducted, we reanalysed the observational data collected in Kraus et al. (2020) and found incorrectly listed values for the mutual inclination. We recalculate these values using the data listed in their tables S5 and S6 and find the mutual inclinations of the outer AB-C binary to the R1, R2, and R3 dust rings (largest to smallest) to be ${\theta }_{\mathrm{R}1}={28.8}^{\circ }$⁠, ${\theta }_{\mathrm{R}2}={28.2}^{\circ }$⁠, and ${\theta }_{\mathrm{R}3}={24.9}^{\circ }$⁠. These values suggest a slightly different structure of the GW Ori disc than has previously been reported in Bi et al. (2020). In this model, the inclination of the disc does not change much with respect to the AB-C binary orbital plane, but the longitude of the ascending node $\mathrm{\Omega }$ changes significantly (measured in the plane of the AB-C binary) between the R2 and R3 rings. The reduction of the assumed misalignment angle keeps the GW Ori disc in the regime for nodal circulation and coplanar alignment, so we expect the previous and current simulations of the GW Ori disc to remain broadly consistent. This disc model is consistent with a disc in a mildly non-coplanar orientation that has been warped by the precession around the outer binary, with possible breaking near the location of the R3 dust ring. Our new disc model better supports the possibility of disc breaking from azimuthal twisting motions rather than a changing binary–disc inclination as one moves inwards from R1 to R3. The same misalignment at the inner R3 ring may imply that the alignment time-scale is long relative to the disc precession time-scale; from equation (15), this implies a value of $\tau <1$ for GW Ori, allowing the disc to break without impedance from viscous forces.

The changing orientation of warped discs is known to produce specific observable signatures. Warping changes the projection of the disc on to the sky plane and produces twisted or S-shaped contours when viewed in molecular lines. This feature has been observed in the GW Ori circumtriple system (Bi et al. 2020; Kraus et al. 2020) as well as within the cavities of some transition discs such as HD 142527 (Casassus et al. 2015), and can be reproduced using simulations of warped discs (Juhász & Facchini 2017; Smallwood et al. 2021). Channel maps can also explore details about the larger surrounding structures around a warped disc, and reveal whether an observed warp is caused by a central source such as a binary or whether the warp is related to infalling matter from the surrounding molecular cloud. In the case of GW Ori, the identification of at least one streamer connected to the disc indicates that both effects may be at play (Czekala et al. 2017; Fang et al. 2017). High-resolution kinematics of warped discs may be used to examine radial flows between broken discs and determine the strength of a disc warp, including whether an observed warp is also a break.

Broken discs have similar observational signatures. If the inner and outer discs are in different orientations, each disc will have differently oriented and possibly separated butterfly patterns in CO channel maps (Zhu 2019). The inner disc will cast shadows as a pair of dark lanes on to the outer disc. The orientation of these lanes may not be symmetric around the disc, and their pattern speed around the disc may not follow the precession rate of the inner disc, depending on the relative orientations of the inner disc, outer disc, and binary (Facchini et al. 2018; Zhu 2019). Discs with multiple breaks may have a complex series of shadowing features, as each disc shadows all of the discs external to it. However, Facchini et al. (2018) find that the flaring geometry of the disc is an important feature for producing asymmetric illumination of the outer disc when the inner and outer discs have a low relative inclination. Since multiple disc breaking is less favoured for flaring discs, the shadowing patterns from these discs may appear more symmetric. The local temperature changes caused by these shadows can generate spiral arms in the outer disc (Montesinos et al. 2016), so the thermal effects of these lanes may be important when attempting to reproduce observations.

The long-term movement of these shadows can give insights into the orientation of the central disc relative to the binary. From equation (11), discs in coplanar and polar orientations have opposite signs of k, and will precess in opposite directions compared to their precession axis. Long-term monitoring of discs, combined with rotational information from kinematic observations, can be used to determine the direction of the shadow’s movement, and thus the precession direction of the inner disc. This can be used to determine whether the inner disc is in a coplanar or polar configuration, and can place constraints on the arrangement of the central binary.

Once a break is formed in the disc, it may slowly drift outwards over time. From Fig. 7, the average drift rate is roughly ${10}^{-3}{a}_{\mathrm{b}}/T_{\mathrm{b}}$ for the inviscid simulations, and roughly twice as fast for the intermediate $\alpha ={10}^{-2}$ simulation. It is unclear whether this behaviour is part of the disc settling towards a new equilibrium state, or whether it is part of the long-term evolution of the disc after it has broken. In our bending-wave simulations (Section 4.1), we observe a range of behaviours, from breaks that drift outwards for thousands of binary orbits to breaks that do not drift at all, with a relationship that is not monotonically dependent on the disc viscosity.

For the break observed in the GW Ori system, the drift rates listed above translate to an outward movement of around $0.7-1.4\times {10}^{-3}\mathrm{a}\mathrm{u}\mathrm{y}{\mathrm{r}}^{-1}$⁠. If the break continues outwards at a constant rate, it will reach the edge of the observed disc (500 au; Kraus et al. 2020) in roughly ${10}^5$ yr. This may be relatively short compared to the disc lifetime of 1–10 Myr (Andrews et al. 2009), and may imply that broken discs are transient features. On the other hand, the precession of the binary is constantly warping the disc, and should constantly generate breaks in the disc as long as it is still misaligned, similar to how other disc substructures such as gaps and spirals may be continuously generated by companions or orbiting planets. In this way, it may be possible for a disc to generate multiple breaks even though the disc conditions are only enough to satisfy a single break, as breaks move outwards in the disc and new ones are generated at the breaking radius.

6 CONCLUSIONS

We have used grid-based hydrodynamic simulations to study the behaviour of warped discs in the context of circumbinary discs created during early star formation. We find that the previously known behaviours for warped discs are reproduced well with the grid-based code, and that disc breaking is readily achieved for discs satisfying $\alpha \lesssim h/r$⁠. From this, we propose a viscous criterion for disc breaking, where the disc must undergo significant precession before aligning perpendicular to the precession axis in order to break.

We also derive new formulations of the disc time-scales in order to predict the location of a disc break. Our formulation for disc breaking suggests that breaking events are more likely when the disc is thinner, the inner cavity is smaller, and the disc power-law profiles are steeper. These criteria are well supported by our simulation results, and the predicted location for the breaking radius is in agreement with the breaks observed in our simulations. We also show that repeated disc breaking, previously seen in simulations of active galactic nucleus discs, is predicted by our analytical formulae, though it is unlikely to occur in a protoplanetary disc. When compared against our simulations, the analytical equations produce breaking radii that are consistent with the observed locations. However, the sensitivity of these equations makes it difficult to use them in a predictive manner and precisely determine the location of the breaking radius. We have also applied our disc breaking formulation against the GW Ori system to explain the discrepancies between previous simulations and their observed breaks, but better observational constraints are required to accurately assess whether the disc will break on its own.

We have explored only a small area of the total parameter space for disc breaking in this paper. Our simulations primarily focus on discs around moderately eccentric $(e\sim 0.5)$ binaries, with equal or moderate mass ratios. Many characteristics of the overall parameter space, such as the regions where a break forms in the middle of the disc (Section 2.3), have yet to be constrained. In the future, well-constrained disc parameters may be combined with our formulation to determine the breaking radii of warped discs.

ACKNOWLEDGEMENTS

We thank our reviewer Alison Young for pointing out the different values for the misalignment angle in the GW Ori system, and for the many suggestions that helped improve the text of this manuscript. The simulations in this work were conducted using the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center. This material is based upon work supported in part by the National Aeronautics and Space Administration under grant no. 80NSSC20M0043 issued through the Nevada NASA Space Grant Consortium. ZZ acknowledges support from NASA award 80NSSC22K1413. Figures in this paper were made with the help of matplotlib (Hunter 2007), numpy (Harris et al. 2020), and visit (Childs et al. 2012), as well as the NAS Visualization Team. RGM and SL acknowledge support from NASA through grants 80NSSC21K0395 and 80NSSC19K0443.

DATA AVAILABILITY

The data used in this paper are available upon request to the corresponding author.

Footnotes

REFERENCES