# Warps and Breaks in Circumbinary Discs

Ian Rabago1,2 , Zhaohuan Zhu1,2 , Stephen Lubow3 , and Rebecca G. Martin1,2
1Department of Physics and Astronomy, University of Nevada, Las Vegas 4505 S.
Maryland Parkway Las Vegas, NV 89154, USA
2Nevada Center for Astrophysics, University of Nevada, Las Vegas 4505 S.
Maryland Parkway Las Vegas, NV 89154, USA
3Space Telescope Science Institute, Baltimore, MD 21218
0000-0001-5008-2794 0000-0003-3616-6822 0000-0002-4636-7348
0000-0003-2401-7168

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Disc warping, and possibly disc breaking, has been observed in protoplanetary
discs around both single and multiple stars. Large disc warps can break the
disc, producing observational signatures in kinematics and shadowing along the
outer disc. In this work, we use comparisons of disc timescales 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. Thus, thin, polar-aligned discs
are more likely to break. We test our analytic formulae using 3D grid-based
simulations of protoplanetary discs warped by the gravitational torque of an
inner binary. We reproduce the expected warp behaviors in different viscosity
regimes and observe disc breaking at locations in agreement with our derived
equations. As disc breaking is observed only for discs with low viscosity, 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 neighboring 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 the disc is relatively thin. We expect repeated or
multiple disc breaking to occur for discs with sufficiently steep power law
profiles. We simulate a disc with steep power-law profiles, and observe two
separate breaking events at locations in rough agreement with our analytical
predictions.

###### keywords:

accretion, accretion discs – protoplanetary discs – hydrodynamics – methods:
numerical – binaries:general – stars:individual:GW Orionis

††pubyear: 20XX††pagerange: Warps and Breaks in Circumbinary Discs–Warps and
Breaks in Circumbinary Discs

##  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 et al., 1957; Shu
et al., 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
et al., 1991; Scott et al., 2000), and shadows cast in protoplanetary discs
(Marino et al., 2015; Debes et al., 2017). The ubiquity of such a structure
warrants careful study of its behavior and evolution.

The warping of discs is a well-studied phenomenon, and a large amount of work
has gone into deriving the analytic behavior of the evolution and testing this
behavior 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 $\bm{L}$ changes with radius, $r$. The
overall evolution of a warped disc can be effectively described with two
different behaviors, 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 timescales 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 which propagates through the disc.

In the diffusive regime, the evolution of a warped disc may be described as a
diffusion equation in $\bm{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 neighboring 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 et al., 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

$\psi=r\left|\frac{\partial\bm{l}(r)}{\partial r}\right|.$ (1)

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 et al., 2021), and previously observed
in simulations of discs around single and binary black holes (Nixon et al.,
2012, 2013; Nealon et al., 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 et al., 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 current 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 the SPH code. With this Lagrangian code the resolution is
determined by the particle density and is therefore a concern near disc
breaks. 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 timescales
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 Timescales of Disc Evolution

In this section, we examine the different evolutionary timescales that affect
a warped disc, and consider how these timescales 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 labeled $M_{\rm tot}=M_{1}+M_{2}$, orbital semi-major
axis $a_{\rm b}$, and orbital eccentricity $e_{\rm b}$. We choose radial
power-law profiles for midplane density, temperature, and surface density,
using the indices $d$, $s$, and $p$ given by 111Previous works commonly use
$p$ for the midplane density index and $q$ for either the temperature or sound
speed index. There is also some inconsistency when considering the sign of the
index, i.e. if the negative sign of the exponent is included in the variable
(e.g. Takeuchi & Lin (2002)). We choose positive indices and our variable
names to remain consistent with the variables in Lubow & Martin (2018), and to
avoid confusion with the variable used for the binary mass ratio.

$\rho(r)=\rho_{0}\left(\frac{r}{r_{0}}\right)^{-d},$ (2)
$T(r)=T_{0}\left(\frac{r}{r_{0}}\right)^{-s}$ (3)

and

$\Sigma(r)=\Sigma_{0}\left(\frac{r}{r_{0}}\right)^{-p}.$ (4)

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

$2d+s-2p=3$ (5)

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

###  2.1 Disc breaking in the wave-like regime

Figure 1: Disc timescales $t_{\rm p}$ and $t_{\rm c}$ as a function of disc
radius on a logarithmic scale. The red vertical line indicates the expected
location of a disc breaking event based on Equations (12) and (13). For this
figure, we consider a nearly coplanar disc around an equal mass binary with
$e_{\rm b}=0.5$, extending from $r_{\rm{in}}=2a_{\rm b}$ to
$r_{\rm{out}}=50a_{\rm b}$. The disc has a constant scale height $h/r=0.1$,
with power law slopes of $p=1.5$ and $s=1.0$.

The torque from the binary potential causes a ring of material to nodally
precess at a frequency

$\omega_{\rm p}(r)=k\left(\frac{a_{\rm b}}{r}\right)^{7/2}\Omega_{\rm b}$ (6)

(e.g. Farago & Laskar, 2010), where $\Omega_{\rm b}$ is the orbital frequency
of the binary, $r$ is the radius of the ring, and $k$ is a constant which
depends on the binary masses, eccentricity, and whether the disc is close to a
coplanar or polar configuration (Martin & Lubow, 2017; Zanazzi & Lai, 2018;
Lubow & Martin, 2018; Martin & Lubow, 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 timescale between the rings of the disc, the disc may become
warped.

If the radial communication timescale is short compared with the precession
timescale then the disc can precess like a solid body. For a disc with inner
radius $r_{\rm{in}}$ and outer radius $r_{\rm{out}}$, the global precession
timescale $t_{\rm p}$ of the disc as a solid body is given by

$t_{\rm p}=\frac{2(1+p)r_{\rm{in}}^{1+p}r_{\rm{out}}^{5/2-p}}{|k|(5-2p)a_{\rm
b}^{7/2}\Omega_{\rm b}}.$ (7)

The radial communication timescale in the wave-like regime $t_{\rm c}$ is
given by

$t_{\rm c}=\frac{4}{(2+s)\Omega_{\rm
b}h_{\rm{out}}}\left(\frac{r_{\rm{out}}}{a_{\rm b}}\right)^{3/2}$ (8)

where $h_{\rm out}$ is the scale height at the outer edge of the disc.

We expect disc breaking to occur when the global precession timescale is
shorter than the communication timescale, $t_{\rm p}<t_{\rm c}$. Analytic
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_{\rm out}$ as the breaking radius gives

$r_{\rm
break}=\left[\frac{(p+1)(s+2)h_{\rm{out}}}{2(5-2p)|k|}\left(\frac{r_{\rm{in}}}{a_{\rm
b}}\right)^{p+1}\right]^{\frac{1}{p-1}}a_{\rm b}.$ (9)

The timescales given in Equations (7) and (8) assume that $r_{\rm{in}}\ll
r_{\rm{out}}$, i.e. that 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_{\rm p}$ and $t_{\rm
c}$ by revisiting these equations. The disc nodal precession frequency is
given in Equations (16)-(17) of Lubow & Martin (2018) as

$\omega_{\rm p}=k\frac{\int_{r_{\rm{in}}}^{r_{\rm{out}}}\Sigma
r^{3}\Omega\left(a_{\rm b}/r\right)^{7/2}\
dr}{\int_{r_{\rm{in}}}^{r_{\rm{out}}}\Sigma r^{3}\Omega\ dr}\Omega_{\rm b}$
(10)

where $\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

$k=\begin{cases}-\frac{3}{4}\sqrt{1+3e_{\rm b}^{2}-4e_{\rm
b}^{4}}\frac{M_{1}M_{2}}{M_{\rm tot}^{2}}&\rm{for\ a\ nearly\ coplanar\
disc,}\\\ \frac{3}{4}\sqrt{5}e_{\rm b}\sqrt{1+4e_{\rm
b}^{2}}\frac{M_{1}M_{2}}{M_{\rm tot}^{2}}&\rm{for\ a\ nearly\ polar\ disc.}\\\
\end{cases}$ (11)

The fraction containing the binary masses $M_{1}M_{2}/M_{\rm tot}^{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 are better for polar discs.

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

$t_{\rm p}=\frac{-2(1+p)}{k(5-2p)\Omega_{\rm
b}}\left(\frac{r_{\rm{in}}}{a_{\rm
b}}\right)^{7/2}\left[\frac{\left(r_{\rm{out}}/r_{\rm{in}}\right)^{5/2-p}-1}{\left(r_{\rm{out}}/r_{\rm{in}}\right)^{-(1+p)}-1}\right].$
(12)

In the limit that $r_{\rm out}\gg r_{\rm in}$, 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

$t_{\rm c}=\int_{r_{\rm{in}}}^{r_{\rm{out}}}\frac{2}{c_{\rm
s}(r)}dr=\frac{4}{s+2}\left[r_{\rm{out}}^{s/2+1}-r_{\rm{in}}^{s/2+1}\right].$
(13)

An example of these timescales are shown as a function of disc radius in
Figure 1. For geometrically thin rings $(r_{\rm{in}}\approx r_{\rm{out}})$,
Equations (12) and (13) reproduce the correct limiting behavior $(t_{\rm
p}\approx 1/\omega_{\rm p}$ and $t_{\rm 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_{\rm p}>t_{\rm c}$, a precession induced
warp can be communicated via bending waves and allow the disc to precess
rigidly. Taking the breaking radius $r_{\rm break}$ to be the value of $r_{\rm
out}$ for which $t_{\rm p}=t_{\rm c}$; for $r_{\rm out}>r_{\rm break}$ we have
that $t_{\rm p}<t_{\rm 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 timescale, starting from
the disc inner edge.

Whether the curves of $t_{\rm p}$ and $t_{\rm c}$ intersect is determined by
the disc _geometry_ $(r_{\rm{in}},r_{\rm{out}},h/r)$, disc _structure_ (power-
law slopes $d$, $s$, and $p$), and inner binary arrangement ($q$, $a_{\rm b}$
, and $e_{\rm 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_{\rm{out}}/r_{\rm{in}}$, whereas the slope of Equation (13) is determined
by the temperature power-law slope $s$. The location of the inner radius,
$r_{\rm{in}}$, sets the fastest precession rate of the entire disc and the
initial height of the $t_{p}$ curve, and is a strong factor in determining
whether or not a disc will break. Because of this sensitive dependence on
$r_{\rm{in}}$, nearly polar aligned discs are more likely to have breaks
(Section 3.1 and 3.2), since they are able to maintain smaller inner cavities
(Franchini et al., 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_{\rm b}>0.41$, as follows from Equation (11).

Once the disc breaks into separate rings, the outer disc will have its
timescales reset as its inner radius has now moved to the location of the
break. These timescales 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.

Figure 2: Predicted disc evolution of the Kraus et al. (2020) parameters, as
inferred from our timescale equations. In each panel the lines are the same as
Fig. 1. The shaded region indicates the inner cavity of the disc. _Top:_ Early
on, the inner edge of the disc is far enough out that a break does not occur.
_Middle:_ As the disc evolves, the inner edge drifts inwards until
$r_{\rm{in}}=21.25\,$au, where a break spontaneously occurs in the disc at
$49\,$au. _Bottom:_ As the inner edge continues to drift inwards, so does the
expected breaking radius. Once the inner disc reaches an inner radius of 20
au, the expected breaking radius has moved to a radius of roughly 35 au.

###  2.2 A Criterion for Viscous Disc Breaking

Figure 3: Predicted number of breaks for a disc with a given initial geometry,
shown as a function of the surface density and temperature power-law exponents
$p$ and $s$. The radial extent of the disc is given as $r_{\rm
disc}=\left[r_{\rm{in}},r_{\rm{out}}\right]$ and the initial scale height is
given at a distance of $r=5a_{\rm b}$. The central binary is an equal mass
binary with eccentricity $e_{\rm b}=0.5$. _Left:_ Initial inner radius and
scale height. _Middle:_ Effect of a larger inner disc radius. _Right:_ Effect
of a larger disc aspect ratio. The ‘$\times$’ in each panel marks the disc
parameters used in our multiple break simulation (Section 3.3.)

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

$t_{\rm align}=\frac{(h/r)^{2}\Omega_{\rm b}}{\alpha(\omega_{\rm p})^{2}}.$
(14)

If the disc aligns significantly faster than the precession rate, the disc
will become coplanar 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 for a break to occur. We can
express this relation by finding the ratio of $t_{\rm p}$ to $t_{\rm align}$

$\tau(r)=\frac{t_{\rm align}}{t_{\rm p}}=\frac{(h/r)^{2}}{\alpha}t_{\rm
p}\Omega_{\rm b}.$ (15)

When this ratio is less than 1, $t_{\rm align}<t_{\rm 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 it may apply to some discs in the wave-like regime
as well. Equation (15) is a simple estimate which 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_{\rm p}<1$ as the criterion to prevent large disc
warping in the viscous regime, where $t_{\nu_{2}}$ is the viscous timecale
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_{\rm p}\,t_{\rm align})^{-1}\sim\omega_{\rm
p}\,t_{\nu 2}$. A short (long) alignment timescale occurs with a long (short)
warp viscous timescale. But the $\nu_{2}$ value for small warps may not apply
for cases of disc breaking.

###  2.3 An Application to GW Orionis

GW Orionis is a hierarchical triple star system surrounded by a large
protoplanetary disc (Mathieu et al., 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}\rm{CO}$, 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 Orionis 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 to $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
30au, 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. However, they find that a disc with a smaller inner radius,
$r_{\rm{in}}=20\rm{au}$, can spontaneously generate a break, with a breaking
radius of roughly 75 au. From this, they determine that the geometry of the
disc is an important factor in determining if a disc is able to break.

Using Equations (12) and (13), it is possible to explain the appearances of
both SPH simulations. To model this triple star system with our analytic
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. Lepp et al., 2023). The binary is composed of
$M_{1}=0.742M_{\rm tot}$ and $M_{2}=0.258M_{\rm tot}$ with an eccentricity of
$e_{\rm b}=0.379$. The effect of the inner binary is that it causes apsidal
precession of the outer binary (Lepp et al., 2023), but we do not consider
this effect when considering the precession of the disc.

For the Kraus et al. (2020) simulation, the initial conditions have $t_{\rm
p}>t_{\rm c}$ everywhere, indicating a disc with no breaks. However, the two
curves approach as the inner edge of the disc drifts inward. Once the inner
edge of the disc reaches about 21au, the two curves intersect and a break
spontaneously appears at $r_{\rm break}=49\rm{au}$, producing a thin secondary
ring along the inner edge. By the time the inner radius has decreased to
$20\rm{au}$, the breaking radius has also moved inwards to a distance of
$r_{\rm break}=35\rm{au}$, a distance that is consistent with the observed
location of the ring seen in the SPH simulations. Figure 2 shows the evolution
of $t_{\rm p}$ and $t_{\rm c}$ in the Kraus disc model. If spontaneous disc
breaking occurs in this manner, the material in the inner ring can be shuttled
inward 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_{\rm c}$. This
prevents the timescale curves from intersecting unless the inner radius is
driven to very small values. Additionally, the change in the shape of the
curves causes the value of $r_{\rm break}$ 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. This disc model is unlikely to break using the binary torque
alone, and thus a planet-driven explanation is 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 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 disc will begin to precess
independently. The inner disc is guaranteed to precesss rigidly, by the
conditions set by the break, but the outer disc may once again precess
differentially with its own $t_{\rm c}$ and $t_{\rm 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_{\rm c}$ and $t_{\rm 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 Figure 3 we map out the number of disc breaks for
different initial values of $r_{\rm{in}}$, $r_{\rm{out}}$, 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_{\rm 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 behavior 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. Figure 3 also shows that
increasing the disc scale height increases the width of each "stripe",
requiring a steeper temperature gradient to achieve multiple breaks.

##  3 Methods

Figure 4: Density contours of the polar-aligning discs at $t=1000$ binary
orbits. Horizontal slice shows disc density, along the plane of the binary.
Discs in the $\alpha=0.01,0.001$, and $10^{-5}$ cases are broken; the inner
disc in the $\alpha=0.01$ simulation has momentarily reconnected with the
outer disc as it precesses about the binary eccentricity vector.

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 setups. Each set is described in the
following subsections. Below we describe features common to all of the
simulations performed in this paper.

We choose spherical-polar coordinates $(r,\theta,\phi)$ for our simulations.
Our simulation domain extends from $5\degree$ to $175\degree$ in $\theta$, and
$0$ to $2\pi$ in $\phi$, and is divided into a grid of $176\times 384$ cells
in $(\theta,\phi)$. We use logarithmically spaced cells in the radial domain.
The extent and number of cells in the radial domain $r$ varies depending on
the simulation set.

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_{\rm s}/\Omega_{\rm K}$.
Here $r_{0}=a_{\rm b}$, $\rho_{0}=1$ at $r_{0}$, $c_{\rm s}$ is the local
isothermal sound speed $\sqrt{P/\rho}$ at $r$, and $\Omega_{\rm 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 cutoff of the form
$\exp{\left[(r-\mu_{r})/\sigma_{r}\right]}$, where $\mu_{r}$ is the location
of the disc edge and $\sigma_{r}$ the relative scale length of the cutoff. We
also use a spherically symmetric density floor with a value of $\rho_{\rm
floor}=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_{\rm
cool}=0.01\Omega_{\rm K}^{-1}$. At the simulation boundaries, we set a one-way
outflow boundary condition in the radial direction, and reflecting boundary
conditions in the polar direction.

The binary components are simulated as gravitational bodies with masses
$M_{1}$ and $M_{2}$, which orbit with semi-major axis $a_{\rm b}$ and
eccentricity $e_{\rm b}$. We integrate the motion of the binary components
using a 2nd-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-
aligned 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 setup 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_{\rm b}=0.28$ and $e_{\rm
b}=0.5$, placing the binary as close to the inner radial domain as possible.
We simulate the region from $0.3=1.07a_{\rm b}$ to $10.0=35.7a_{\rm 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_{\rm 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_{\rm b}$ and
$\sigma_{r}=0.35a_{\rm b}$. The outer edge of the disc extends to the outer
radial domain. We vary the $\alpha$-viscosity parameter between
$\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 Orionis system. We use
the known binary parameters from Kraus et al. (2020), using the AB-C binary
separation of $a_{\rm 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_{AB}=0.742M_{\rm tot}$ and the outer C component as the
second mass $M_{C}=0.258M_{\rm tot}$, placed in an orbit with semi-major axis
$a_{\rm b}=8.89$ au and eccentricity $e_{\rm b}=0.379$.

We simulate the region from $1.5a_{\rm b}$ to $50a_{\rm b}$, covering a range
of roughly 13 au to 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 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_{\rm 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_{\rm p}$ and $t_{\rm c}$
suggest a disc can undergo multiple breaks if the power-law slope of $t_{\rm
p}$ is less than that of $t_{\rm 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 setup for the polar disc simulations, changing the radial domain
so the inner domain remains at $r_{\rm{in}}=1.07a_{\rm b}$ but extending the
outer radial domain to $r_{\rm{out}}=100a_{\rm 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 a scale height ranges from
$(h/r)_{\rm in}=0.251$ at the disc inner edge to $(h/r)_{\rm out}=0.013$ at
$r=50a_{\rm b}$. To ensure 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_{\rm b},52.1a_{\rm 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 setup is marked on each panel of Figure 3 with an ‘$\times$’.
For the left panel, with $r_{\rm{in}}=2.5a_{\rm b}$, the breaks are predicted
to occur at distances of 4.0, 8.6, and 29.1 $a_{\rm b}$. We evolve this disc
for 1500 binary orbits, enough time for the outer parts of the disc to evolve
for a few orbital timescales.

##  4 Results

###  4.1 Disc Warping and Breaking

Figure 5: Evolution of the warp profile for the viscous ($\alpha=0.1$, left)
and inviscid cases ($\alpha=10^{-5}$, right). Each curve represents the
$z$-component of the unit angular momentum vector $\mathbf{l_{z}}$ as a
function of radius. Curves are plotted every 100 binary orbits. Figure 6:
Motion of the warp for the viscous ($\alpha=0.1$, left) and inviscid cases
($\alpha=10^{-5}$, right). _Top row:_ The $\psi$ parameter plotted against the
disc radius. Curves are plotted every 100 binary orbits as in Figure 5.
_Middle row:_ The radial location of the warp, determined as the local
location of maximum $\psi$. For wave-like propagation, both the initial break
and secondary wave propagation are plotted. _Bottom row:_ The outward
propagation velocity of the warp. The dotted line represents the local sound
speed velocity of $c_{\rm s}$ at the warp radius.  Figure 7: Inner regions of
the $\alpha=10^{-5}$ polar-aligning simulation, projected back onto a single
plane. _Left:_ Midplane density. The inner disc features an RWI vortex.
_Center:_ Disc eccentricity. The inner disc is somewhat eccentric. _Right:_
Disc vorticity. Note the strong vorticity minimum at the location of the
vortex.  Figure 8: Comparison of the simulated breaking radius to analytical
predictions. _Top:_ Analytical timescales as predicted from Equations (12) and
(13). Here, $r_{\rm{in}}$ is calculated to be $2.2a_{\rm b}$ (see text for
details). _Bottom:_ Surface density of the $\alpha=10^{-5}$ simulation at
$t=100$ binary orbits, in units of the initial surface density profile
$\Sigma_{0}$. The breaking radius is denoted by the sharp dip in
$\Sigma/\Sigma_{0}$. In both plots, vertical dashed lines indicate the
predicted locations of the breaking radius.

Our polar-aligning simulations exhibit both diffusive and wave-like behavior,
depending on our choice of $\alpha$. Figure 4 shows the various discs at
$t=1000$ binary orbits, displaying varied behaviors in their evolution. For
the $\alpha=10^{-1}$ simulation, the disc develops a strong warp which evolves
diffusively, while for the $\alpha=10^{-2},10^{-3},$ and $10^{-5}$ simulations
the disc breaks and exhibits wave-like dispersion. Figure 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 timescale, 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 peaks in $\psi$ are present, one at the breaking
radius and one in the warped outer disc. We find the radial location of both
maxima 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 is shown in Figure 6.

For simulations in the diffusive regime, the inner disc generates a large warp
as the disc 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.
Figures 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_{\rm b}$ by
the end of the simulation. A small wave is launched ahead of the main warp
(visible in Figure 6 as a series of secondary purple peaks between 10 to 25
$a_{\rm b}$, beneath the yellow curves drawn at later times). This wave does
not appear with a similar feature in Figure 5. The highly diffusive nature of
this disc quickly removes this feature, causing it to disappear by
$t=500T_{\rm b}$.

Simulations in the wave-like regime show the initial disc breaking at $r_{\rm
break}\sim 4a_{\rm 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 Figure 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.

Our analytical equations for disc breaking show close agreement with the
location of the break observed in the numerical simulations. In Figure 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_{\rm{in}}$ for our analytic equations, we compare
the density-weighted angular momentum of the disc $\int\Sigma r^{3}\Omega dr$
to the analytic value expected for a disc with the same surface density power
law profile. For this calculation, we choose a snapshot shortly after the
initalization, when the disc has had time to equalize but before a break has
had time to open up. We consider $r_{\rm{in}}$ to be the location at which
these two values diverge. For the polar-aligning discs, this gives a value of
$r_{\rm{in}}=2.2a_{\rm b}$.

We also observe substructures in the individual discs created after the
breaking event. Figure 7 shows the $\alpha=10^{-5}$ simulation as it would
appear “flattened” onto 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 disk. 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 disc 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 et al., 2015).

Figure 9: Surface density and warp profile of the GW Ori simulation at
$t=0,2500$, and $5000T_{\rm b}$.

###  4.2 GW Orionis

Figure 9 shows the disc surface density and warp strength of the simulation
representing the GW Orionis system (Sec.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 Orionis 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 Orionis system has a moderate binary eccentricity of $e_{\rm
b}=0.379$ and a mass ratio, $q=M_{C}/M_{AB}=0.348$, making the precession rate
normalized by $\Omega_{\rm b}$ of this system about one-fifth as strong as
that induced by the polar disks in Section 4.1.

The binary torque truncates the inner edge of the disc at roughly
$r_{\rm{in}}=3.5a_{\rm b}$. This is a larger value of $r_{\rm{in}}$ than
observed in the Smallwood et al. (2021) simulations, and larger than the value
we initially used for calculating the breaking radius above. Recalculating the
location of the breaking radius with the empirical value of $r_{\rm{in}}$
changes the disc timescales so that no break occurs in the disc.

###  4.3 Multiple Breaks

Figure 10: Rendering of the multiple break simulation at $t=1200$ binary
orbits. The directions of $\bm{L}_{\rm b}$ and ${\bm{e}}_{\rm b}$ are
indicated in the top right, with $\bm{L}_{\rm b}$ pointing into the page. The
first two breaks at $7a_{\rm b}$ and $28a_{\rm b}$ are visible. There is a
large warp at $35a_{\rm b}$ that appears similar to a break, but the disc is
not broken at this distance. See text for details. An animation of the disc
evolution is also available as an online video.  Figure 11: Similar to Figure
8, but for the multiple break simulation. In this figure, each panel shows
different choices of the inner radius, $r_{\rm{in}}=2.5a_{\rm b}$ (top) and
$r_{\rm{in}}=4.0a_{\rm b}$ (middle). Both choices produce outer breaking
locations that are roughly consistent with the observed breaks in the
simulation. The bottom panel is from a simulation at a time of 1500 $T_{\rm
b}$.  Figure 12: Azimuthal rotation angle plotted as a function of radial
location and time. Areas in the same connected disc precess together, and
create a band of the same color. By the end of the simulation, three distinct
bands are visible, corresponding to the three discs in the simulation. The
surface density profile from Figure 11 is also shown at right to compare the
locations of the observed gaps. A detailed time evolution of the disc surface
density and warp profile is also available as an online video.

A visualization of our multiple break simulation (Sec. 3.3) is shown in Figure
10. 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
roughly at radial distances of $7$ and $28a_{\rm b}$, with the outer edge of
the disc reaching equilibrium at about $r_{\rm{out}}=60a_{\rm b}$. The third
ring shows a strong warp at roughly $35a_{\rm 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 RWI vortices form in
the innermost disc, but not in the second or third discs.

Figure 11 compares the observed locations of each break to the locations
predicted by the timescale equations. We find it is difficult to determine a
suitable inner radius for this disc. Choosing an inner boundary of
$r_{\rm{in}}=2.5a_{\rm b}$ our analytical equations predict 3 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_{\rm b}$, is
not observed in the innermost disc. Instead the observed breaks may be roughly
fit by choosing $r_{\rm{in}}=4.0a_{\rm 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

$\phi_{L}=\arctan{\left(\frac{L_{y}}{L_{x}}\right)},$ (16)

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 $\phi_{L}$
can be viewed as the longitude of ascending node of the disk 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 $\phi_{L}$. In Figure 12, we plot the value of $\phi_{L}$ across
the disc as a function of simulation time. As the disc evolves and breaks off
rings, distinct horizontal bands of color appear, 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 $\bm{L}$
changes sharply across this region.

In this figure, it can be seen that the innermost disc breaks off almost
immediately, in about 100 binary orbits. The outer disc then develops a warp
that gradually extends outwards to about $25a_{\rm b}$, before breaking off a
second disc at roughly $1000T_{\rm b}$. After the disc breaks at $25a_{\rm
b}$, the warp between the first break and the second break gradually vanishes
and this disc region becomes co-planar. This is another indication that the
disc loses the communication at the breaking radius and the outer disc cannot
torque the inner disc anymore. A faint band can be seen extending from the
second disc to about $35a_{\rm 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 timescales to fully develop, so the innermost break will
have time to drift outwards while subsequent breaks form and develop. The
bottom panel of Figure 11 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 disk no longer follows a power-law
distribution even though the analytic equations assume one when calculating
the location of a break. We discuss this behavior further in Section 5.1.

##  5 Discussion

Figure 13: $\tau$ as calculated by Equation (15) for different values of
$\alpha$. The vertical line indicates the rough location of the break as seen
in our low viscosity simulations.

For the polar discs used in Section 3.1, we show the effects of disc viscosity
on $\tau$ by plotting $\tau$ as a function of radius in Figure 13, with the
assumption that while calculating $t_{\rm 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
Figures 4 and 5 of Ogilvie (1999). From our Figure 6, the maximum warp
amplitude is roughly $\psi_{\rm max}=0.5-0.75$, at a distance of $r=10a_{\rm
b}$. This corresponds to warp viscosities of
$\nu_{2}=0.019a_{\rm{}_{b}}^{2}\Omega_{\rm b}$,
$\nu_{3}=0.0013a_{\rm{}_{b}}^{2}\Omega_{\rm b}$, and
$t_{\nu_{2}}=5000\Omega_{\rm b}^{-1}$.

We note that, for the simulations in Section 4.1, there is a change in the
overall disc orientation of about 5 to 10 degrees from its initial
orientation. This change occurs faster than the local precession timescale
induced by the binary, and instead occurs as the warp (diffusive or wave-like)
propagates through the disc. This effect is partially visible in Figure 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) is a sensitive function of the disc parameters, particularly the inner
disc radius $r_{\rm{in}}$. Because of this sensitivity, small adjustments in
the value of $r_{\rm{in}}$ results in large changes in the breaking radius.
Conversely, small adjustments in the location of the breaking radii produce
relatively minor changes in the required value of $r_{\rm{in}}$. Thus,
Equations (12) and (13) can easily 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
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_{\rm{in}}$ to account for this change, the predicted breaks are in much
closer agreement. For the multiple break simulation (Section 12), 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 definitions of $r_{\rm{in}}$ using surface density cutoffs,
precession rates, or other methods will calculate slightly different values,
which will lead to greatly differing estimates of the breaking radii. Our
simulations of GW Ori in Section 4.2 show that the differences in code
formulation (SPH vs. grid-based) can also change the location of
$r_{\rm{in}}$, further complicating this issue. This behavior was also noted
in the simulations by Young et al. (2023), who comment that the wide
dependence on disc and binary parameters make the creation of a disk breaking
equation a challenging matter.

Our analytic equations are also restricted by their lack of full coverage of
the possible parameters. Previous SPH simulations have observed a dependence
of disc breaking and inclination to the precession vector (Facchini et al.,
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. Analysis 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_{c}$, dependent on the disc viscosity. Our current
equations do not include any factors related to these parameters, and so do
not capture any related effects when determining a breaking radius. The
current analysis is relatively simple and is strictly analytic, 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 Figure 11. Future analysis may improve this method by directly
integrating the disc radial profiles at later times to provide better
calculations of the disc timescales.

###  5.2 Observational Signatures

The changing orientation of warped discs is known to produce specific
observable signatures. Warping changes the projection of the disc onto the sky
plane, and produces twisted or S-shaped contours when viewed in molecular
lines. This feature has been observed in the GW Orionis 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). High resolution kinematics of warped discs may be
able to determine the strength of a disc warp, and may aid in determining if
an observed warp is also a break, as well as the flows in the region between
the individual discs.

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 onto 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 of
the disc is an important geometrical feature for producing asymmetric
shadowing when the inner and outer disc are at low relative angles. Since
multiple disc breaking is less favored for flaring discs, the shadowing
pattern may appear different. 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.

Once a break is formed, it slowly drifts outwards over time. From the polar-
aligning discs in Section 4.1, the drift rate is roughly $10^{-3}a_{\rm
b}/T_{\rm b}$ for the inviscid simulations, and roughly twice as fast for the
intermediate $\alpha=10^{-2}$ simulation. For the break observed in the GW
Orionis system, this translates to a drift rate of at least $7.7\times
10^{-4}\rm{au/yr}$. Assuming a constant outward drift rate, the break will
reach the edge of the observed disc (500 au, Kraus et al. 2020) in roughly
$6\times 10^{5}$ years. 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 outward in the disc and new ones are generated at
the breaking radius.

##  6 Conclusions

We have used grid-based hydrodynamic simulations to study the behavior of
warped discs in the context of circumbinary discs created during early star
formation. We find that the previously known behaviors 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 to the precession axis in order to break.

We also derive new formulations of the disc timescales 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, predicted by our formulations, is
possible.

When compared against our simulations, the analytic 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 breaking formulation against the GW Orionis 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. In the future, well-constrained disc parameters
may be combined with our formulation to determine the breaking radii of warped
discs.

## Acknowledgements

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. Z. Z. 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).
RGM and SHL acknowledge support from NASA through grants 80NSSC21K0395 and
80NSSC19K0443.

## Data Availability

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

## References

  * Andrews et al. (2009) Andrews S. M., Wilner D. J., Hughes A. M., Qi C., Dullemond C. P., 2009, ApJ, 700, 1502
  * Bae et al. (2015) Bae J., Hartmann L., Zhu Z., 2015, ApJ, 805, 15
  * Berger et al. (2011) Berger J. P., et al., 2011, A&A, 529, L1
  * Bi et al. (2020) Bi J., et al., 2020, ApJ, 895, L18
  * Burke (1957) Burke B. F., 1957, AJ, 62, 90
  * Burrows et al. (1995) Burrows C. J., Krist J. E., Stapelfeldt K. R., WFPC2 Investigation Definition Team 1995, in American Astronomical Society Meeting Abstracts. p. 32.05
  * Casassus et al. (2015) Casassus S., et al., 2015, ApJ, 811, 92
  * Childs et al. (2012) Childs H., et al., 2012, in , High Performance Visualization–Enabling Extreme-Scale Scientific Insight. pp 357–372
  * Debes et al. (2017) Debes J. H., et al., 2017, ApJ, 835, 205
  * Demianski & Ivanov (1997) Demianski M., Ivanov P. B., 1997, A&A, 324, 829
  * Doǧan et al. (2018) Doǧan S., Nixon C. J., King A. R., Pringle J. E., 2018, MNRAS, 476, 1519
  * Dullemond et al. (2022) Dullemond C. P., Kimmig C. N., Zanazzi J. J., 2022, MNRAS, 511, 2925
  * Facchini et al. (2013) Facchini S., Lodato G., Price D. J., 2013, MNRAS, 433, 2142
  * Facchini et al. (2018) Facchini S., Juhász A., Lodato G., 2018, MNRAS, 473, 4459
  * Farago & Laskar (2010) Farago F., Laskar J., 2010, MNRAS, 401, 1189
  * Franchini et al. (2019) Franchini A., Lubow S. H., Martin R. G., 2019, ApJ, 880, L18
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Juhász & Facchini (2017) Juhász A., Facchini S., 2017, MNRAS, 466, 4053
  * Katz (1973) Katz J. I., 1973, Nature Physical Science, 246, 87
  * Kerr et al. (1957) Kerr F. J., Hindman J. V., Carpenter M. S., 1957, Nature, 180, 677
  * King et al. (2013) King A. R., Livio M., Lubow S. H., Pringle J. E., 2013, MNRAS, 431, 2655
  * Kraus et al. (2020) Kraus S., et al., 2020, Science, 369, 1233
  * Lepp et al. (2023) Lepp S., Martin R. G., Lubow S. H., 2023, ApJ, 943, L4
  * Li et al. (2000) Li H., Finn J. M., Lovelace R. V. E., Colgate S. A., 2000, ApJ, 533, 1023
  * Lodato & Price (2010) Lodato G., Price D. J., 2010, MNRAS, 405, 1212
  * Lodato & Pringle (2007) Lodato G., Pringle J. E., 2007, MNRAS, 381, 1287
  * Lovelace et al. (1999) Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805
  * Lubow & Martin (2018) Lubow S. H., Martin R. G., 2018, MNRAS, 473, 3733
  * Lubow & Ogilvie (2000) Lubow S. H., Ogilvie G. I., 2000, ApJ, 538, 326
  * Lynden-Bell & Pringle (1974) Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603
  * Marino et al. (2015) Marino S., Perez S., Casassus S., 2015, ApJ, 798, L44
  * Martin (2008) Martin R. G., 2008, MNRAS, 387, 830
  * Martin & Lubow (2017) Martin R. G., Lubow S. H., 2017, ApJ, 835, L28
  * Martin & Lubow (2018) Martin R. G., Lubow S. H., 2018, MNRAS, 479, 1297
  * Martin et al. (2019) Martin R. G., et al., 2019, ApJ, 875, 5
  * Mathieu et al. (1991) Mathieu R. D., Adams F. C., Latham D. W., 1991, AJ, 101, 2184
  * Miyoshi et al. (1995) Miyoshi M., Moran J., Herrnstein J., Greenhill L., Nakai N., Diamond P., Inoue M., 1995, Nature, 373, 127
  * Montesinos et al. (2016) Montesinos M., Perez S., Casassus S., Marino S., Cuadra J., Christiaens V., 2016, ApJ, 823, L8
  * Nealon et al. (2015) Nealon R., Price D. J., Nixon C. J., 2015, MNRAS, 448, 1526
  * Nealon et al. (2022) Nealon R., Ragusa E., Gerosa D., Rosotti G., Barbieri R., 2022, MNRAS, 509, 5608
  * Nixon & King (2016) Nixon C., King A., 2016, in Haardt F., Gorini V., Moschella U., Treves A., Colpi M., eds, , Vol. 905, Lecture Notes in Physics, Berlin Springer Verlag. p. 45, doi:10.1007/978-3-319-19416-5_2
  * Nixon et al. (2012) Nixon C., King A., Price D., Frank J., 2012, ApJ, 757, L24
  * Nixon et al. (2013) Nixon C., King A., Price D., 2013, MNRAS, 434, 1946
  * Ogilvie (1999) Ogilvie G. I., 1999, MNRAS, 304, 557
  * Papaloizou & Lin (1995) Papaloizou J. C. B., Lin D. N. C., 1995, ApJ, 438, 841
  * Papaloizou & Pringle (1983) Papaloizou J. C. B., Pringle J. E., 1983, MNRAS, 202, 1181
  * Petterson et al. (1991) Petterson J. A., Rothschild R. E., Gruber D. E., 1991, ApJ, 378, 696
  * Pringle (1981) Pringle J. E., 1981, ARA&A, 19, 137
  * Pringle (1992) Pringle J. E., 1992, MNRAS, 258, 811
  * Rabago et al. (2023) Rabago I., Zhu Z., Martin R. G., Lubow S. H., 2023, MNRAS,
  * Raj et al. (2021) Raj A., Nixon C. J., Doğan S., 2021, ApJ, 909, 81
  * Scott et al. (2000) Scott D. M., Leahy D. A., Wilson R. B., 2000, ApJ, 539, 392
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33
  * Shu et al. (1983) Shu F. H., Cuzzi J. N., Lissauer J. J., 1983, Icarus, 53, 185
  * Smallwood et al. (2019) Smallwood J. L., Lubow S. H., Franchini A., Martin R. G., 2019, MNRAS, 486, 2919
  * Smallwood et al. (2021) Smallwood J. L., Nealon R., Chen C., Martin R. G., Bi J., Dong R., Pinte C., 2021, MNRAS, 508, 392
  * Stone et al. (2020) Stone J. M., Tomida K., White C. J., Felker K. G., 2020, arXiv e-prints, p. arXiv:2005.06651
  * Takeuchi & Lin (2002) Takeuchi T., Lin D. N. C., 2002, ApJ, 581, 1344
  * Young et al. (2023) Young A. K., Stevenson S., Nixon C. J., Rice K., 2023, arXiv e-prints, p. arXiv:2306.11809
  * Zanazzi & Lai (2018) Zanazzi J. J., Lai D., 2018, MNRAS, 473, 603
  * Zhu (2019) Zhu Z., 2019, MNRAS, 483, 4221
  * Zhu et al. (2015) Zhu Z., Dong R., Stone J. M., Rafikov R. R., 2015, ApJ, 813, 88