# The Bardeen-Petterson effect in accreting supermassive black-hole binaries:
disc breaking and critical obliquity

Rebecca Nealon1,2, Enrico Ragusa3,4, Davide Gerosa5,6,7, Giovanni Rosotti3,8,
Riccardo Barbieri9
1Centre for Exoplanets and Habitability, University of Warwick, Coventry CV4
7AL, UK
2Department of Physics, University of Warwick, Coventry CV4 7AL, UK
3School of Physics and Astronomy, University of Leicester, University Road,
Leicester LE1 7RH, UK
4Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique
de Lyon UMR5574,
F-69230, Saint-Genis-Laval, France
5Dipartimento di Fisica “G. Occhialini”, Universitá degli Studi di Milano-
Bicocca, Piazza della Scienza 3, 20126 Milano, Italy
6INFN, Sezione di Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy
7School of Physics and Astronomy & Institute for Gravitational Wave Astronomy,
University of Birmingham,
Birmingham, B15 2TT, UK
8Leiden Observatory, Leiden University, P.O. Box 9513, NL-2300 RA Leiden, The
Netherlands
9Max Planck Institute for Gravitational Physics (Albert Einstein Institute),
Am Mühlenberg 1, Potsdam 14476, Germany

###### Abstract

The inspiral of supermassive black-hole binaries in gas-rich environment is
driven by the presence of an accretion disc and viscous interactions tend to
align the spin of the black holes with the orbital angular momentum of the
disc. Recent work introduced a new iterative approach to describe the
alignment process and the resulting non-linear evolution of the surrounding
warped accretion disc. Their model predicted that black-hole spins reach
either full alignment or a ‘critical obliquity’ where solutions to the warp
equations cease to exist. In this paper, we show that this critical region
corresponds to the disc breaking phenomenon, where the disc is disrupted into
two or more discrete sections. We use 3D hydrodynamical simulations to (i)
recover the predictions of the semi-analytic model and (ii) unveil a richer
phenomenology where the disc exhibits either unsuccessful, single and multiple
breaks. We additionally identify hydrodynamic effects such as spiral arms that
are able to stabilise the disc against breaking beyond criticality. Our
results show that when disc breaking occurs, the ability of black holes and
disc to align is compromised and in some cases even prevented as the binary
inspirals.

###### keywords:

accretion, accretion discs — black-hole mergers — gravitational waves —
hydrodynamics

††pubyear: 2021††pagerange: The Bardeen-Petterson effect in accreting
supermassive black-hole binaries: disc breaking and critical
obliquity–LABEL:lastpage

##  1 Introduction

Accretion discs play a pivotal role in a variety of astrophysical processes,
ranging from planet formation to interacting binary stars, and active-galactic
nuclei (AGN) (Pringle, 1981; Frank et al., 2002). While accretion onto a
single Newtonian object results into a planar disc configuration, the presence
of external torques might induce a distorted, or ‘warped’, structure. Known
processes that can excite disc warps include higher-order harmonics of the
central gravitational potential (Tremaine et al., 2009), the presence of a
binary companions (Lubow & Ogilvie, 2000), embedding in a stellar clusters
(Bregman & Alexander, 2012), radiation pressure from the central object
(Pringle, 1996), magnetic fields (Lai, 1999), as well as relativistic effects
(Bardeen & Petterson, 1975).

General-relativistic frame dragging is at the heart of the so-called ‘Bardeen-
Petterson effect’. For a disc surrounding a spinning black hole (BH), Lense-
Thirring precession preferentially dissipates angular momentum in a direction
perpendicular to the BH spin, thus acting towards aligning the disc with the
equatorial plane of the BH (Lense & Thirring, 1918). Crucially, the relevant
precession frequency decreases rather steeply with the distance from the BH
($\Omega_{\rm LT}\propto 1/R^{3}$, where $R$ a radial coordinate; e.g. Kumar &
Pringle 1985), implying that momentum can be efficiently re-distributed only
for gas rings that are sufficiently close to the BH. While the inner disc —up
to the so-called ‘warp radius’— aligns with the BH spin, the outer disc
maintains its generically misaligned orientation. However, most of the angular
momentum resides in this outer disc that thus reacts by pulling the BH spin
toward a fully planar configuration (Rees, 1976). For typical AGN-disc
parameters, this process takes place on a timescale of $1-10$ Myr (Natarajan &
Pringle, 1998; Lodato & Gerosa, 2013).

The Bardeen-Petterson effect has been invoked to explain misaligned jets in
both AGNs (Caproni et al., 2006, 2007; Falceta-Gonçalves et al., 2010) and
microquasars (Maccarone, 2002; Martin et al., 2008), as well as quasi-periodic
oscillations in X-ray binaries (Fragile et al., 2001) and the light curves of
some tidal disruption events (Lei et al., 2013). For discs surrounding
supermassive BH binaries, gas-driven spin alignment is thought to be a key
process to prevent the ejections of BH-merger remnants from their host
galaxies following relativistic recoils (Bogdanović et al., 2007; Miller &
Krolik, 2013; Dotti et al., 2010; Gerosa et al., 2015). The upcoming
gravitational-wave mission LISA (Amaro-Seoane et al., 2017) has the potential
of directly measuring the spin directions of several of these systems, thus
providing a complementary probe to further test the occurrence of the Bardeen-
Petterson effect (Berti & Volonteri, 2008; Sesana et al., 2014; Sayeb et al.,
2021).

Building on earlier explorations by Scheuer & Feiler (1996), Martin et al.
(2007, 2009), and Tremaine & Davis (2014), some of the authors recently
presented a systematic investigation of the Bardeen-Petterson effect in
accreting supermassive BH binaries (Gerosa et al., 2020). They put forward a
one-dimensional (1D) numerical scheme that takes into account, in a consistent
fashion, both (i) the non-linear character of the fluid viscosities in warped
configurations (Ogilvie, 1999; Ogilvie & Latter, 2013) and (ii) the combined
effect of the Lense-Thirring and companion torques. Their study highlighted
the occurrence of a ‘critical obliquity’ —a specific region in the parameter
space where solutions to the underlying 1D boundary-value problem cease to
exist. Hints of this behaviour were previously reported by Tremaine & Davis
(2014) with a different numerical scheme. Gerosa et al. (2020) conjectured
that their numerical divergences correspond to a physical scenario where the
disk breaks into disconnected regions, hindering the subsequent spin-alignment
process.

In the context of BH accretion, disc breaking and tearing has almost
exclusively been explored using numerical simulations. Nelson & Papaloizou
(2000) conducted the first three-dimensional (3D) simulations of a misaligned
accretion disc around a BH and found that the inner disc aligns with the BH
spin and the outer disc maintains its original misalignment, in broad
agreement with the theoretical expectations of the Bardeen-Petterson effect.
Additionally, for large initial misalignments, they found that the transition
between the inner and outer disc plane was no longer continuous so that the
disc was ‘close to breaking into two discrete pieces’ (Nelson & Papaloizou,
2000). The concept of disc breaking was later expanded upon by Nixon et al.
(2012), who demonstrated that discs could tear into more than just two pieces,
and under certain circumstances could break into precessing rings of gas.
Breaking was further confirmed in the wave-like regime by Nealon et al.
(2015), suggesting that disc breaking is an inevitable consequence of moderate
to strongly misaligned flows accreting onto rotating BHs. Disc breaking has
also been confirmed using both a grid based magneto-hydrodynamic treatment
(Liska et al., 2021) and analytic arguments for discs subjected to a non-
Keplerian potentials (Doǧan & Nixon, 2020).

In this paper, we investigate the interplay between disc criticality and disc
breaking. Guided by the 1D predictions of Gerosa et al. (2020), we present a
large suite of 3D smoothed particle hydrodynamics (SPH) simulations of
misaligned accretion discs surrounding spinning BHs in binary systems. We
confirm that (i) the occurrence of a critical obliquity corresponds to disc
breaking and (ii) its importance increases as the influence of the BH
companion increases. Our simulations further allow us to unveil a richer
phenomenology which includes (iii) disc breaking into both single and multiple
precessing rings as well as (iv) the stabilizing effect of spiral arms in the
disc. A future publication will make use of these prescriptions to investigate
the spin directions of large populations of supermassive BHs and their
relevance to the LISA mission.

This paper is organised as follows. In Sec. 2 we briefly summarise the physics
of warped accretion discs and the breaking conditions we employ. In Sec. 3, we
present our numerical implementation. In Sec. 4, we illustrate our main
results in terms of both disc morphology and BH spin alignment. In Sec. 5 and
6 we present our conclusions and highlight prospects for future work in this
area.

##  2 Disc breaking

In this section we summarise the relevant analytic framework of warped discs
around BHs that we will make use of in this work. Consistent with most of the
previous literature (but see Raj & Nixon 2021), we use the terms ‘breaking’
and ‘tearing’ interchangeably to both describe discs that separate into a
discontinuous structure.

###  2.1 Defining a break

To define where and when the disc tears in our simulations we consider both
the mass surface density and warp profiles. First, we require that the surface
density profile $\Sigma(R)$ present a sustained local minimum $\Sigma_{\rm
min}$. This straightforward definition was also used by both Nixon et al.
(2012) and Nealon et al. (2015). Second, we consider the gradient of the
angular momentum profile (also refereed to as ‘warp profile’, e.g. Lodato &
Price 2010)

$\displaystyle\psi(R)=R\absolutevalue{\frac{\partial\hat{\textbf{{L}}}(R)}{\partial
R}},$ (1)

where $\hat{\textbf{{L}}}$ is the unit vector pointing in the direction of the
disc angular momentum. We require an increasing local maximum $\psi_{\rm max}$
to determine that disc breaking is occurring. This is similar to previous work
by both Nelson & Papaloizou (2000), who considered a steep increase in the
radial warp profile to identify a potential break, and Raj & Nixon (2021).

Figure 1: The different disc structure in a simulations with warping,
unsuccessful breaking, successful breaking/tearing with one ring and
successful breaking/tearing with multiple rings (left to right). The
renderings in the top row show the $x$-$y$ plane, with the BH oriented such
that $\bm{J}=(\sin\theta,0,\cos\theta)$. The second and third rows show the
surface density $\Sigma$ and the warp profile $\psi$. The simulations used
here, from left to right, have parameters $\kappa=3.9\times 10^{-5}$ and
$\theta=60\degr$, $\kappa=1.2$ and $\theta=40\degr$, $\kappa=1.4\times
10^{-4}$ and $\theta=80\degr$, $\kappa=2.6\times 10^{-1}$ and
$\theta=60\degr$. All simulations are shown at $t=4.3$ companion orbits.

With these two metrics in mind, we find that the outcome of our simulations
can be classified into the following outcomes. These four scenarios are
showcased in Fig. 1 and a description of how we discriminate these can be
found in Sec. 3.4.

  1. 1.

Warping: A disc that is not flat but does not show signs of breaking. Here the
angular momentum profile varies as a function of radius, i.e. $\psi(R)\neq 0$,
but there is no indication of minima in the surface density profile.

  2. 2.

Unsuccessful breaking: A disc that shows the characteristics of breaking but
does not actually separate into distinct smaller discs. The warp profile shows
a local maximum with a corresponding local minimum in the surface density
profile. However, as the disc evolves, $\psi_{\rm max}$ may increase but
$\Sigma_{\rm min}$ does not continue to decrease.

  3. 3.

Successful breaking (single): A disc that shows the characteristics of
breaking with an increasing $\psi_{\rm max}$ and decreasing $\Sigma_{\rm
min}$, successfully separating into two smaller discs. In these simulations,
$\Sigma_{\rm min}$ continues to decrease until $\Sigma_{\rm min}\approx 0$.

  4. 4.

Successful breaking (multiple): A disc that shows the characteristics of
breaking at many radii simultaneously. The warp profile has local maxima and
the surface density approaches zero at multiple locations, leading to several
rings tearing off the disc.

As shown in Fig. 1, the condition $\psi_{\rm max}$ and $\Sigma_{\rm min}$ are
strongly correlated. In the case that the disc is stable, there is no or a
weak $\psi_{\rm max}$ and no $\Sigma_{\rm min}$. As this disc evolves
$\psi_{\rm max}$ decreases, making the disc more stable against breaking (see
Sec. 5.4 for further discussion). In the case that the disc is unstable, the
disk starts to break when both the surface density decreases and the warp
profile increases at the same radial location. The disc then actually
separates if this process continues until $\Sigma_{\rm min}\approx 0$.
Previous 1D analyses (Doǧan et al., 2015, 2018; Gerosa et al., 2020) could
only perform a coarser distinction between connected and disconnected discs.
On the other hand, the 3D hydrodynamical simulations presented in this paper
allow us to delineate between unsuccessful, single, and multiple breaking.

###  2.2 Influence of the binary companion

Let us consider a BH of mass $M$ and spin $J=GM^{2}\chi/c$ (where
$\chi\in[0,1]$ is the Kerr parameter) embedded in a disc with Shakura &
Sunyaev (1973) viscosity $\alpha$. In the absence of a binary companion, the
disc angular-momentum profile with radius is self-similar: BHs with different
masses and spins will be surrounded by re-scaled versions of the same gaseous
structure (Scheuer & Feiler, 1996; Martin et al., 2007).

If instead the BH is part of a binary system, the disc is subjected to both
the Lense-Thirring torque at small radii as well as the tidal torque from the
companion at large radii. In this case the dynamics depends on the interaction
between the two external torques in the disc. Each torque has an associated
radius where it creates the largest warp; at $R_{\rm LT}$ the Lense-Thirring
torque most strongly affects the warp profile while at $R_{\rm tid}$ the tidal
torque does the same. Following Gerosa et al. (2020), we parameterise the
interaction of the binary by introducing the non-dimensional ‘companion
parameter’ $\kappa$ with

$\displaystyle\kappa=\left(\frac{R_{\rm tid}}{R_{\rm LT}}\right)^{-7/2}\,.$
(2)

We set $R_{\rm LT}$ to be the reference radius such that $R_{\rm LT}=R_{0}$.
The advantage of this parameterisation is that it encodes information on the
companion’s mass $M_{\star}$ and orbital separation $R_{\star}$. In the
notation of Gerosa et al. (2020), this ‘companion parameter’ reads

$\displaystyle\kappa$ $\displaystyle\simeq
0.66\left(\frac{M}{10^{7}M_{\odot}}\right)^{2}\left(\frac{\chi}{0.5}\right)^{2}\left(\frac{M_{\star}}{10^{7}M_{\odot}}\right)\left(\frac{R_{\star}}{0.1{\rm
pc}}\right)^{-3}$
$\displaystyle\;\;{\times\left(\frac{H_{0}/R_{0}}{0.002}\right)^{-6}\left(\frac{\alpha}{0.2}\right)^{-3}\left[\frac{\zeta}{1/(2\\!\times\\!0.2^{2})}\right]^{-3}\,,}$
(3)

where

$\displaystyle\zeta=\frac{2(1+7\alpha^{2})}{\alpha^{2}(4+\alpha^{2})}\,.$ (4)

where in the limit of $\alpha\rightarrow 0$, $\zeta\rightarrow 1/(2\alpha)$
and $\kappa\propto\alpha^{3}$. Here $H_{0}/R_{0}$ is the aspect ratio set at
radius $R_{0}$ where the Lense-Thirring torque most strongly warps the disc.
From Eqs. (16) and (20) in Gerosa et al. (2020), this radius is equivalent to

$\displaystyle
R_{0}=\frac{GM}{c^{2}}\left(\frac{H_{0}}{R_{0}}\right)^{-4/3}\left(\frac{4\chi}{\alpha\zeta}\right)^{2/3}\,.$
(5)

The case $\kappa=0$ corresponds to the self-similar case of a single BH, while
larger values of $\kappa$ correspond to configurations where the companion BH
strongly perturbs the disc evolution. This could be due to the binary being
massive (large $M,M_{\star}$), the BH being more rapidly rotating (large
$\chi$), the orbit being tight (small $R_{\star}$), the disc being thin (small
$H_{0}/R_{0}$), a low viscosity (small $\alpha$), or some combination of
these. For binaries inspiralling under gas-assisted migration, the parameter
$\kappa$ increases with time and may eventually lead the system to
criticality.

###  2.3 Break radius

If the disc breaks, this is likely to happen at a radius $R_{\rm break}$ that
maximizes the warp profile $\psi$. Martin et al. (2009) estimated this
location by matching the external torques due to Lense-Thirring precession and
the companion to find

$\displaystyle R_{\rm
break}=\left(\frac{8G^{1/2}M^{1/2}JR_{\star}^{3}}{3c^{2}M_{\star}}\right)^{2/9}\,.$
(6)

In particular, Eq. (6) is independent of the disc viscosity $\alpha$.
Therefore, we do not expect significant differences in the location of the
breaking radius if the disc is in the diffusive or wave-like regime. However,
we note that this prediction relies solely on the balance between Lense-
Thirring and companion tidal torques and neglects, by definition, the response
of the disc to the warp propagation. Because of this, Eq. (6) does not
(somewhat unphysically) depend on the relative inclination between the BH and
the companion’s orbit. As described below, the relative inclination does play
a role in the breaking dynamics (see also Fig. 1).

###  2.4 Disc backreaction

As the BH warps the disc, the disc reacts by aligning the BH with its own
angular momentum (Rees, 1976). The evolution of the angle $\theta$ between the
BH spin J and the binary’s orbital angular momentum $\textbf{{L}}_{*}$ is
given by

$\displaystyle\frac{\rm d\cos\theta}{\rm dt}=\frac{\rm
d\hat{\textbf{{J}}}}{\rm dt}\cdot\hat{\textbf{{L}}}_{*}\,.$ (7)

The change of angular momentum of the BH $\rm d\textbf{{J}}/dt$ can be
calculated from the integral of the torque exerted by the disc onto the BH
(Gerosa et al., 2020)

$\displaystyle\frac{\rm d\textbf{{J}}}{dt}=-\int_{R_{\rm min}}^{R_{\rm
max}}\frac{2G}{c^{2}}\frac{\textbf{{J}}\times\textbf{{L}}}{R^{3}}2\pi R\rm
dR\,,$ (8)

where L is the angular momentum of the disc,
$\absolutevalue{\textbf{{L}}}=\Sigma\sqrt{GMR}$ and $R_{\rm min}$ and $R_{\rm
max}$ are the inner and outer extent of the disc.

Ideally, one would like to fully take into account the back-reaction of the
disc onto the BH and evolve the system self-consistently. This was possible
with the 1D scheme of Gerosa et al. (2020) but they had to rely on a quasi-
adiabatic treatment and could only follow the evolution up to the critical
obliquity. Although tracking the binary inspiral is prohibitive for our 3D
simulations, one can still use the above description to investigate how the
spin aligns on short timescales while the disc reaches its steady state. This
is measured from our simulations by assuming azimuthal asymmetry to evaluate
Eq. (8) and (7) for each simulation snapshot (see Appendix A for a detailed
explanation).

##  3 Numerical simulations

Our simulations are performed with the 3D SPH code Phantom (Price et al.,
2018). This code has been used extensively to model inclined discs around BHs
(Nixon et al., 2013; Nealon et al., 2015), to examine disc breaking (Nixon et
al., 2012; Doǧan et al., 2015), and for comparison with analytical predictions
of warped discs (Lodato & Price, 2010).

We conduct a total of 143 simulations considering multiple viscosities, binary
separations, disc aspect ratios and inclinations. Here we detail the relevant
aspects of the code for our application and the initial conditions used in all
of our runs. Additionally, to the best of our knowledge a portion of our
parameter suite constitutes the thinnest discs around BHs simulated to date
with SPH (albeit marginally). Simulations are run until they show successful
breaking or 15 binary orbits. Although this is a short time compared to the
viscous time of the disc, the entirety of 143 simulations required $\sim 2.5$
million CPU hours (we discuss this limitation in Section 5.4).

The vast majority of the discs in our parameter suite are comfortably in the
‘diffusive’ regime, where $\alpha\gtrsim H/R$ and the warp propagates
diffusively (Papaloizou & Pringle, 1983). For the discs that have $\alpha\sim
H/R$, we may expect tilt oscillations to occur as seen by Ivanov & Illarionov
(1997); Lubow et al. (2002) and Nealon et al. (2015). However, our discs are
often not simulated until they reach a steady state and, as we will discuss
next, we set a slightly larger outer boundary than was considered in these
previous works. We thus do not expect (nor do we recover) evidence of tilt
oscillations here.

###  3.1 Black-hole modeling

In Phantom the rotating BH is modelled using a fixed potential. To do this we
use a post-Newtonian approximation to model the potential of the rotating BH,
achieved with a first order in $v/c$ correction in the momentum equation
(Nelson & Papaloizou, 2000)

$\displaystyle\frac{\text{d}\bm{v}}{\text{d}t}=-\frac{1}{\rho}\nabla
P+\bm{v}\times\bm{h}-\nabla\Phi+S_{\rm visc}\,,$ (9)

where $\bm{v}$, $\rho$ and $P$ are the gas velocity, density and pressure
respectively and $S_{\rm visc}$ is the viscous force per unit mass. The term
$\bm{v}\times\bm{h}$ represents the gravitomagnetic force per unit mass, where

$\displaystyle\bm{h}=\frac{2\chi
G^{2}M^{2}}{R^{3}c^{3}}\left(\bm{\hat{J}}-3\frac{(\bm{\hat{J}}\cdot\bm{r})\bm{r}}{R^{2}}\right)\,,$
(10)

r is the spherical coordinate vector, and $R$ is the distance to the primary
BH. We use the modified potential of Nelson & Papaloizou (2000)

$\displaystyle\Phi(R)=-\frac{GM}{R}\left(1+\frac{3R_{\rm g}}{R}\right)\,,$
(11)

where $R_{\rm g}=GM/c^{2}$ is the gravitational radius of the BH. While
preventing the gravitational force from tending to infinity close to the BH,
this expression also accurately recovers both apsidal precession frequency at
large radii and the sign of the nodal precession frequency. We note that the
modified potential looses accuracy for $R\lesssim 10R_{\rm g}$ but this is
within our numerical accretion radius.

We note that using a fixed potential limits our simulations, in that
physically it is equivalent to assuming that the primary BH is located at the
centre of mass of the binary system. To respect this assumption, we only
consider systems with a mass ratio $M_{\star}/M\approx 0.01$ and a low disc
mass such that the centre of mass of the system is within the accretion radius
of the BH (see Sec. 3.2). We vary the other parameters entering Eq. (3) to
span a wide range of $\kappa$ values, above and below criticality.

###  3.2 Initial conditions

Each of our simulations is initialised with a disc and binary companion in
orbit around the primary BH. Here we detail the properties of the initial
conditions, noting that simulations are non-dimensionalised by using $\kappa$
from Eq. (3). Our parameters and the corresponding values of $\kappa$ are
summarised in Table 1.

$\kappa$ | $R_{\star}/R_{\rm g}$ | $R_{\rm out}/R_{\rm g}$ | $\alpha$ | $H_{\rm ref}/R_{\rm ref}$ | $R_{0}/R_{\rm g}$ | $H_{0}/R_{0}$ | $\theta$ | $N$
---|---|---|---|---|---|---|---|---
$3.9\\!\times\\!10^{-5}$ | 398 | 250 | 0.10 | 0.08 | 19.2 | 0.089 | 20$\degr$–160$\degr$ | 14
$1.4\\!\times\\!10^{-4}$ | 398 | 250 | 0.05 | 0.05 | 25.9 | 0.052 | 20$\degr$–160$\degr$ | 12
$1.7\\!\times\\!10^{-4}$ | 398 | 250 | 0.15 | 0.08 | 26.8 | 0.082 | 20$\degr$–160$\degr$ | 15
$5.9\\!\times\\!10^{-4}$ | 250 | 150 | 0.05 | 0.05 | 25.9 | 0.052 | 10$\degr$–160$\degr$ | 10
$2.7\\!\times\\!10^{-3}$ | 398 | 150 | 0.10 | 0.05 | 49.3 | 0.044 | 20$\degr$–160$\degr$ | 13
$1.2\\!\times\\!10^{-2}$ | 398 | 250 | 0.15 | 0.05 | 68.5 | 0.041 | 20$\degr$–160$\degr$ | 11
$2.8\\!\times\\!10^{-2}$ | 398 | 250 | 0.20 | 0.05 | 83.0 | 0.039 | 20$\degr$–160$\degr$ | 10
$2.6\\!\times\\!10^{-1}$ | 398 | 250 | 0.10 | 0.03 | 136.8 | 0.021 | 20$\degr$–160$\degr$ | 14
$1.2$ | 398 | 250 | 0.15 | 0.03 | 190.3 | 0.019 | 20$\degr$–160$\degr$ | 14
$2.1\\!\times\\!10^{4}$ | 250 | 150 | 0.10 | 0.01 | $1.2\times 10^{3}$ | 0.0040 | 20$\degr$–160$\degr$ | 11
$9.2\\!\times\\!10^{4}$ | 250 | 150 | 0.15 | 0.01 | $1.7\times 10^{3}$ | 0.0036 | 20$\degr$–160$\degr$ | 10
$2.2\\!\times\\!10^{5}$ | 250 | 150 | 0.20 | 0.01 | $2.1\times 10^{3}$ | 0.0034 | 20$\degr$–160$\degr$ | 9
Table 1: Summary of the initial conditions of our simulations, ordered by
increasing $\kappa$. Here $\kappa$ is a measure of the relative influence of
the secondary [cf. Eq. (3)], $R_{\star}$ is the binary separation, $R_{\rm
out}$ is the outer radius of the disc, $\alpha$ is the disc viscosity, $H_{\rm
ref}/R_{\rm ref}$ the disc aspect ratio at $R_{\rm ref}$, $R_{0}$ is the
Lense-Thirring radius [cf. Eq. (5)], $H_{0}/R_{0}$ is the disc aspect ratio at
$R_{0}$, $\theta$ is the initial inclination between the disc and BH, and $N$
is the number of simulations we have completed in the listed $\theta$ range.
The simulations in the upper (lower) part of the table are presented in Fig. 3
(Fig. 4).

The disc is initialised as a flat disc in the $x$-$y$ plane (i.e., in the
plane of the binary orbit) with a mass of $10^{-6}M$. This disc mass is
deliberately low to respect our assumption of using a fixed potential and to
avoid any back-reaction effect on the properties of the binary companion. The
surface density profile is given by

$\displaystyle\Sigma(R)=\Sigma_{\rm ref}\left(\frac{R}{R_{\rm
ref}}\right)^{-1}\left(1-\sqrt{\frac{R_{\rm in}}{R}}\right)\,,$ (12)

where the normalisation $\Sigma_{\rm ref}$ is determined from the disc mass,
$R_{\rm ref}=30R_{\rm g}$ is the reference radius, and $R_{\rm in}=15R_{\rm
g}$ is the inner edge of the disc. The outer radius of the disc $R_{\rm out}$
is set to either $150R_{\rm g}$ or $250R_{\rm g}$ depending on the orbit of
the companion (cf. Table 1). Our results are reported in orbits of the binary
companion, with all simulations running for a minimum of 3 orbits of the
binary or $\sim$150 orbits at $R_{\rm ref}$.

We assume that the disc is vertically isothermal, such that the sound speed in
the disc can be described by $c_{\rm s}(R)=c_{\rm s,ref}(R/R_{\rm ref})^{-q}$
with $q=3/4$. Here $c_{\rm s,ref}$ is determined by the disc thickness (aspect
ratio), with $(H_{\rm ref}/R_{\rm ref})=0.01,0.03,0.05$ and $0.08$ set at the
reference radius $R_{\rm ref}$ (cf. Table 1). While these values are
relatively large for AGN discs (Haiman et al., 2009), we are limited by our
numerical resolution to $H/R\gtrsim 0.01$. We stress that the dynamics only
depend on the companion parameter $\kappa$: results obtained with larger
aspect ratio will still be robust for discs with a lower $H/R$ but the same
value of $\kappa$.

We model the viscosity in the disc using the Shakura & Sunyaev (1973)
prescription, with $\alpha=0.05,0.10,0.15$ and $0.20$. This is implemented in
Phantom using the shock viscosity term described by Price et al. (2018) (but
see also Artymowicz & Lubow, 1994; Murray, 1996; Lodato & Price, 2010). The
artificial viscosity coefficient $\alpha_{\rm AV}$ is related to the physical
viscosity by

$\displaystyle\alpha\approx\frac{1}{10}\alpha_{\rm AV}\frac{\langle
h\rangle}{H}\,,$ (13)

where $\langle h\rangle$ is the shell-averaged smoothing length. For a given
resolution determined by $\langle h\rangle/H$, we set $\alpha_{\rm AV}$ to
give the targeted $\alpha$. Price et al. (2018) suggests that $\alpha_{\rm
AV}\gtrsim 0.1$ is necessary to resolve the physical viscosity and for all of
our simulations we have $\alpha_{\rm AV}\gtrsim 0.85$. Our choice of surface
density and sound speed power-law profiles implies that $\alpha$ varies with
radius throughout the disc following a power-law with index $-1/6$ (Lodato &
Price, 2010). Once simulations begin evolving, we find that $\langle
h\rangle/H$ becomes roughly constant across the vast majority of the disc
(except where breaking occurs), essentially removing any radial variation in
$\alpha$.

In the simulations the BH has a spin of $\chi=0.9$ and the mass is set to
$M=1$ (alongside $G=c=1$), noting that for our figures we rescale it to
$M=10^{7}M_{\odot}$. To accommodate a relative inclination between the disc
and BH, the spin angular momentum vector of the BH is set to
$\bm{\hat{J}}=(\sin\theta,0,\cos\theta)$ where $\theta$ is the relative
inclination. Thus a BH with $\theta=0$ would have spin along the $z$ axis
relative to a disc that is initialised in the $x$-$y$ plane. Material that
falls inside $R_{\rm in}$ is accreted without further checks.

We restrict the mass of the binary companion to be relatively low to
accommodate our use of a fixed BH potential (see Sec. 3.1). For all of our
simulations, the mass of the binary companion is set to $M_{\star}=10^{-2}M$
and the semi-major axis $R_{\star}=250,398R_{\rm g}$ as in Table 1. The binary
companion is initially set on a circular orbit neglecting the disc mass in the
$x$-$y$ plane. The accretion radius of the companion is set to $0.25$ of the
Hill radius, equivalent to $14.6R_{\rm g}$ for $R_{\star}=398R_{\rm g}$ and
$9.2R_{\rm g}$ for $R_{\star}=250R_{\rm g}$, in line with previous guidance by
Nealon et al. (2018). Although the orbit of the companion is free to evolve
and feels the back-reaction from the disc, this effect is negligible due to
the low mass ratio ($M_{\rm disc}/M_{\star}=10^{-4}$).

Accurately resolving the disc is crucial to recovering disc breaking (Nealon
et al., 2015). We use $N=5\times 10^{6}$ particles for each simulation, which
ensures that the discs are well resolved. To check this we measure the average
smoothing length to disc scale height ratio $\langle h\rangle/H$ as in Lodato
& Price (2010). The disc is discretised into 300 radial annuli and the
particle properties in each annuli averaged to produce radial profiles. To
take into account the warping of the disc, we additionally measure the disc
scale height from the instantaneous warped mid-plane. We meet the resolution
criteria $\langle h\rangle/H<0.5$ in all but the innermost region of the disc
across all of our simulations.

###  3.3 Evaluating the companion parameter $\kappa$

In order to scale our simulations consistently with the description of Gerosa
et al. (2020), we need to evaluate $\kappa$ which in turn depends on both
$R_{0}$ and $H_{0}/R_{0}$. To connect these parameters, we use the analytic
expression of the aspect ratio

$\displaystyle\frac{H}{R}=\frac{H_{\rm ref}}{R_{\rm ref}}\left(\frac{R}{R_{\rm
ref}}\right)^{1/2-q}\,,$ (14)

and Eq. (5), solving for $R_{0}$ in terms of $R_{\rm g}$. This yields

$\displaystyle R_{0}=\left(\frac{4\chi}{\alpha\zeta}\right)\left(\frac{H_{\rm
ref}}{R_{\rm ref}}\right)^{-2}\left(\frac{R_{\rm
ref}}{R_{g}}\right)^{-1/2}R_{\rm g}\,,$ (15)

where $(H_{\rm ref}/R_{\rm ref})$ is the aspect ratio at the reference radius
in our simulation initial conditions. From $R_{0}$, one can then evaluate
$H_{0}/R_{0}$ from Eq. (14) and thus $\kappa$ from Eq. (3). The resulting
values are reported in Table 1.

###  3.4 Quantifying a breaking disc

Given the size of our parameter suite, we desire an automated process to
identify if and where any given simulation exhibits tearing. Based on the
criteria highlighted in Sec. 2.1, we consider either unsuccessful or
successful tearing to occur when we find a radially correlated local minimum
in $\Sigma$ and local maximum in $\psi$ —that is, where a given $\Sigma_{\rm
min}$ and $\psi_{\rm max}$ occur at the same $R$. In particular, we consider
these to be at the same radii when they are within three radial bins of each
other, corresponding to 1% of the full radial domain. Successful tearing
occurs when additionally $\Sigma_{\rm min}\approx 0$, which we accept when
$\Sigma$ drops below 10% of the maximum of the initial surface density
profile.

Our discs feature spiral arms at large radii due to the tidal interaction with
the binary. Crucially, these can be picked up as false positives by the above
criteria. We thus also require that $\psi_{\rm max}>1.0$ to confirm tearing,
as most the spirals tend to be associated with values of $\psi_{\rm max}$ that
are considerably lower. All borderline cases were also visually inspected. We
use this procedure across our simulations for all time-steps to identify and
locate the breaking radius.

Figure 2: Density renderings for a representative set of our suite of 143
simulations showing examples of warping, unsuccessful breaking and successful
breaking. Common to all simulations are the spiral arms at the outer disc edge
driven by the tidal interaction with the BH companion (shown in green). Top to
bottom the view alternates between the $x$-$y$ and the $x$-$z$ plane. The
colour scale shows column density in $\rm kg/\rm m^{2}$ and all panels are
shown on the same spatial scale. Simulations are shown at their end time
(different for each). A summary of the full suite is shown in Fig. 3.

##  4 Results

###  4.1 Qualitative behaviour

Figure 2 shows a demonstrative selection of the discs from our simulated
suite. Broadly speaking, all of our simulations evolve similarly, with a warp
developing in the inner regions while the outer disc is shaped by the tidal
interaction with the binary companion.

Within the first couple of orbits, the binary companion induces two spiral
arms which tends to be more pronounced in the simulations with the smaller
aspect ratio. As expected, the structure of the spiral arms is not affected by
the relative inclination of the primary BH and are sustained throughout the
duration of the simulations (irrespective of whether the disc tears or not).
Additionally, the orbit of the binary does become inclined slightly, but this
is a small effect ($\lesssim 1\degr$).

The inner disc evolution depends on whether the disc warps, breaks into two
sections, or breaks into multiple rings (cf. Fig. 1). In all cases, the inner
region ($R\lesssim 50R_{\rm g}$) shows evidence of a warp within the first
couple of orbits of the binary companion. As the disc continues to evolve,
numerous discs in our suite show visual evidence of break. For those that show
multiple rings, as in previous work we find that the rings tear off
successively from the inner region outwards (e.g. Nixon et al., 2013). The
thickness of the sections that tear off correspond to the disc thickness, with
thinner rings forming for those discs with the smaller aspect ratio. The
broken components precess differentially resulting in a range of misalignments
relative to the outer disc. In a handful of simulations, the broken ring or
disc also exhibits some local asymmetric perturbations but these are short
lived.

Figure 3: Comparison of our 3D hydrodynamical simulations to the 1D semi-
analytic prediction by Gerosa et al. (2020) as a function of $\kappa$ (Eq. 3)
and the initial BH misalignment $\theta$. Here the behaviour of the disc in
our simulations is coded with warped discs as squares, unsuccessful breaking
as pentagons, successful breaking into two discs with diamonds, and tearing
into many rings as stars. The shaded regions indicate the critical region
identified by Gerosa et al. (2020) where 1D disc solutions cannot be found.
Colours represent the viscosity of the discs.

By the end of our simulations the entire disc has tilted away from its initial
orientation in the $x$-$y$ plane, but the final inclination depends on whether
the disc has broken or not. In the cases where the disc has not broken, the
outer edge of the disc is almost co-planar with the orbit of the binary and
the inner region is warped. Alternatively, when the disc tears, the outer disc
has an inclination that lies between that of the binary companion and the
equatorial plane of the BH, with the inner broken regions more strongly
misaligned (for example, see the lower row of Fig. 2).

###  4.2 Testing the semi-analytic predictions

Figure 3 summarises our results for the different $\kappa$ values and relative
inclinations $\theta$. Overplotted on Fig. 3 are the semi-analytic predictions
by Gerosa et al. (2020) for where disc breaking should occur for each of the
$\alpha$ values considered. Here we have used their $\beta=3/2$ solutions (in
their notation, $\beta$ is the power-law index of the viscosity profile) but
note that this choice does not alter our conclusions. The solid lines indicate
the locations of the critical obliquity: for each value of $\alpha$, 1D
solutions cannot be found in the central shaded region of the plot bounded by
the two solid lines.

Overall, we find very good agreement between the 1D analytic prediction from
Gerosa et al. (2020) and the results of our 3D hydrodynamic simulations. In
the region where Gerosa et al. (2020) predicts breaking, we additionally
distinguish cases of unsuccessful tearing, successful breaking and breaking
with single vs. multiple rings. Far from criticality (i.e. above and below the
solid lines), simulated discs tend to warp without breaking (squares in Fig.
3). The semi-analytic approach correctly describes the transition to a
different regime characterised by either unsuccessful (pentagons) or breaking
(diamonds or stars). Once the inclination is greater than criticality, we find
that whether discs break into two or many discs depends on their disc
properties and not their relative inclination.

Overall, we thus confirm that the critical obliquity does indeed correspond to
disc breaking.

####  4.2.1 Large aspect ratios, $H/R$

We do identify some cases where the 3D and 1D results differ substantially.
Across our full parameter suite we do not observe breaking in any of the discs
with $H/R=0.08$. This suggests disagreement with the semi-analytic model, in
particular for the discs with $\kappa=3.9\times 10^{-5}$ and $1.7\times
10^{-4}$ and lower inclinations. This discrepancy is unlikely to be due to the
$\kappa$ and $\alpha$ values used, as the $\kappa=1.4\times 10^{-4}$ series
has a similar $\kappa$ and a lower $\alpha$ but good agreement. We also
dismiss any resolution effects because these discs are thicker than the rest
of our suite and thus present slightly better numerical resolution than the
thinner discs. For these disc parameters, the simulations are also comfortably
in the diffusive regime so we can safely disregard any potential issues due to
wave-like behaviour. We thus conclude that the semi-analytic calculation,
which is inherently designed to a model a disc-like structure, has an
additional limitation in assuming that the disc is sufficiently thin. For the
other parameters that we have held constant (i.e. $R_{*}$), this limitation
corresponds to roughly $H/R\sim 0.08$. In support of this we find good
agreement with the semi-analytic model for our $H/R=0.05$ discs and excellent
agreement when $H/R=0.03$.

####  4.2.2 Large companion parameters, $\kappa$

We additionally run three sets of simulations at very large values of
$\kappa>10^{4}$, summarised in Fig. 4. Here the influence of the binary
companion is much larger and our corresponding $R_{0}$ values are well outside
the outer radius of the disc (c.f. Table 1). Such large $R_{0}$ values are
problematic when comparing to the semi-analytic model, as the current
implementation of Gerosa et al. (2020) assumes that $R_{0}$ is located within
the disc. For the runs presented earlier ($\kappa\lesssim 1.2$), the
combination of the dimensionless scaling and radial range resulted in an inner
boundary that was consistent with the inner boundary in our 3D simulations. If
$\kappa$ is orders of magnitudes larger, however, their inner boundary falls
well outside the inner boundary we have adopted. For example, with
$\kappa=2.1\times 10^{4}$ the inner boundary assumed by Gerosa et al. (2020)
sits at $123R_{\rm g}$ which is almost at our simulated outer disc edge of
$150R_{\rm g}$. A direct comparison between the two approaches is thus not
possible for these large-$\kappa$, large-$R_{0}$ sets of simulations.

Despite this, at high $\kappa$ we find that Fig. 4 still shows that for
moderate to large inclinations we can expect breaking and for small
inclinations, the disc to remain stable. We thus expect our predictions in
relation to the BH disc alignment to be indicative even at higher $\kappa$
(see Section 4.5).

####  4.2.3 High viscosity, $\alpha$

We additionally have two sets of simulations that have $\alpha=0.2$, which at
low $\kappa$ do not necessarily show agreement with the semi-analytic model.
For the case with $\kappa=2.8\times 10^{-2}$ and $\alpha=0.2$ this is
particularly surprising because the semi-analytic model should work quite well
at low $\kappa$ and high $\alpha$. As before, for these discs we can rule out
any differences due to our numerical implementation as the cause of this
discrepancy. However, consideration of the surface density profiles suggests
that the inner regions — where we expect breaking to occur — is rapidly
accreted. Figure 5 compares two $\Sigma$ profiles for simulations with
$\alpha=0.05$ and $\alpha=0.2$. In the low $\alpha$ case the disc shows clear
breaking at $\approx 50R_{\rm g}$ and there is still a significant amount of
material inside this radius making up the inner ring. The high $\alpha$ case
has accreted much more material from the inner edge, potentially prohibiting
the ring from breaking off.

Figure 4: Same as Fig. 3 but for additional simulation sets with higher
$\kappa$ which violates the assumption of the semi-analytic 1D model. Even at
these large $\kappa$ values We still find evidence of disc breaking in our 3D
simulations. Figure 5: Surface density profiles showing that at high $\alpha$
the material where we might expect to see disc breaking is rapidly accreted.
Here the simulations both have $\theta=80\degr$ with $\alpha=0.05$,
$\kappa=1.4\times 10^{-4}$ and $\alpha=0.2$, $\kappa=2.8\times 10^{-2}$. Both
simulations are shown at 5.2 orbits.

###  4.3 Location of the break radius

Figure 6: Location of the breaking radius measured in the simulations scaled
by its predicted location, Eq. (6). We show our simulations with unsuccessful
and successfully broken discs with $\kappa=2.6\times 10^{-1}$. In the early
stages of these simulations there appears to be a dependence on the
inclination. Retrograde discs (stars) consistently break at larger radii than
their prograde counterparts (circles).

Figure 6 shows a representative sample of the location of the breaking radius
for our simulations. We find that, when breaking starts, the radius at which
it occurs depends on the inclination but the relationship is unclear. The
measured break radius is also systematically smaller than the prediction of
Eq. (6). In most cases, the breaking radius decreases with time, although the
opposite happens for a few runs (e.g. $\theta=40\degr$, $140\degr$ in Fig. 6).
Both of these observations are likely a consequence of the short duration of
our simulations, as we only aim to demonstrate stability against breaking.
Overall, we find that the disc breaks at a location that is up to $\sim 5$
times smaller than the prediction by Martin et al. (2009).

Across our parameter suite we consistently find that the break radius for
retrograde discs is larger than for their prograde counterparts. Figure 7
shows an example prograde/retrograde pairing (with $\kappa=1.45\times 10^{-2}$
and $\theta=90\degr\pm 40\degr$). The disc structure only differs in the inner
regions with the retrograde case breaking at a slightly larger radius. Gerosa
et al. (2020) finds that these cases are perfectly symmetric, with an
identical dynamics. However, the direction of the spirals relative to the BH
spin does depend on whether the BH is prograde or retrograde and thus breaks
this assumption of symmetry. This an exclusively 3D effect that cannot be
captured with simpler 1D models. We confirm this behaviour for all of our
prograde/retrograde pairings and note that, even though the details of the
breaking are slightly different, the prediction of whether the disc will break
or not is robust (Fig. 3).

Figure 7: Comparing the density evolution of a prograde/retrograde pair of
simulations after three companion orbits. While the 1D treatment implies
perfect symmetry, in 3D simulations this is broken by the spirals induced by
the binary companion. The prediction of whether the disc will break or not
remains robust, see Fig. 3, although these discs break at slightly different
locations. The simulations shown here have $\kappa=1.45\times 10^{-2}$ and
$\theta=90\degr\pm 40\degr$.

The predicted value of $R_{\rm break}$ in Eq. (6) is independent of viscosity.
In Fig. 8, we compare two of our simulations with $\alpha=0.05$, $H/R=0.05$
and $\kappa=1.4\times 10^{-4}$ vs $\alpha=0.2$, $H/R=0.01$ and
$\kappa=2.2\times 10^{5}$. Here we hold $\theta=40\degr$ and $R_{*}$ (and thus
the external torques) as constant. The lower $\alpha$ case demonstrates
tearing while the higher $\alpha$ case shows evidence of trying to tear but
ultimately is unsuccessful, cf. Fig. 3. This demonstrates that our simulations
are consistent with predictions from both Martin et al. (2009) and Gerosa et
al. (2020): while successful disc tearing depends strongly on the disc
viscosity, the radius where it tears does not.

Figure 8: Location of signatures of tearing for simulations with the same BH
and binary torque but different viscosities. Here the simulation with
$\alpha=0.05$ has $\kappa=1.4\times 10^{-4}$ and $\alpha=0.2$ has
$\kappa=2.2\times 10^{5}$. Although only the lower $\alpha$ disc successfully
tears, both show indications of tearing at the same radius, simultaneously
confirming the predictions by Martin et al. (2009) and Gerosa et al. (2020).

###  4.4 Spiral arms can prevent disc tearing

Here we further investigate our simulation with $\kappa=1.4\times 10^{-4}$ and
$\theta=40\degr$, which we identify as unsuccessfully breaking. Figure 9 shows
the evolution of the surface density and warp profiles throughout this
simulation (purple curves). As the warp propagates from the inner region, it
grows in amplitude with a peak at $R\sim 130R_{\rm g}$, until it reaches the
region where the spiral arms are located. Upon meeting the spiral arms just
before $R\sim 200R_{\rm g}$, the warp sharply decreases in amplitude and
breaking is halted. This simulation thus suggests that the spiral arms can
stabilise the disc against breaking when the disc is set to break near or
beyond them. Such an effect is likely due to the local variations in the disc
viscosity: spiral arms increase the local viscosity and a larger viscosity
makes it harder to break the disc (e.g. Fig. 3). We note that this effect can
be captured only using hydrodynamical simulations and is thus absent in the
predictions of Gerosa et al. (2020).

To further investigate this behaviour, Fig. 9 also shows a second simulation
(green curves) where the companion is placed at a larger radius
$R_{\star}=597R_{\rm g}$ resulting in $\kappa=4.3\times 10^{-5}$. With the
increased separation, the companion has a weaker effect on the disc, the
spiral arms are weaker, and one would anticipate that it is harder to break
the disc. However, Fig. 9 shows that the disc breaking is actually more
successful ($\psi\sim 1.5$) in the case with the larger separation and weaker
spiral arms. Figure 10 displays the density rendering of these comparison
simulations at the same time as the final time-step shown on Fig. 9. In the
case with the strong spiral arms ($R_{\star}=398R_{\rm g}$) we find only
unsuccessful breaking, whilst in the case with weak spiral arms
($R_{\star}=597R_{\rm g}$) there is a break forming in the disc — indeed
roughly at the location of the spiral arms.

This point illustrates that, while the parameter $\kappa$ captures the
qualitative occurrence of the breaking, a more complex dynamics is present and
can only be captured with detailed simulations.

Figure 9: Surface density ($\Sigma$, top panel) and warp profiles ($\psi$,
bottom panel) for discs with $\theta=40\degr$, $\kappa=1.4\times 10^{-4}$
(‘Strong spirals’, purple) and $\theta=40\degr$, $4.3\times 10^{-5}$ (‘Weak
spirals’, green). Solid (dashed) lines are computed at $t=4$ orbits ($t=10$
orbits) of the closer companion. The spirals are located at $R\sim 180R_{\rm
g}$. Disc breaking is more successful in the case where the spirals are weaker
as shown by the greater warp amplitude. Figure 10: Breaking can be prevented
by local hydrodynamical structures, like the spirals driven by the interaction
with the binary companion. Upper and lower panel show a mass rendering of the
simulations with $\theta=40\degr$ and $\kappa=1.4\times 10^{-4},4.3\times
10^{-5}$, respectively, at the same timestep of the dotted lines in Fig. 9.
Increasing the binary separation reduces the strength of the spirals, making
it easier to break the disc even though $\kappa$ decreases.

###  4.5 Disc-black hole alignment

We now calculate $d\cos\theta/dt$ (cf. Sec. 2.4 and Appendix A) to evaluate
how disc breaking affects alignment between the disc and BH. In Fig. 11 we
show two representative sub-sets of our simulations with all discs prograde
and $\kappa$ equal to either $1.4\times 10^{-4}$ or $2.6\times 10^{-1}$. While
both sets contain simulations that have both warping and successful breaking,
the discs in the low $\kappa$ set only include discs with a single break while
the high $\kappa$ only presents multiple breaks.

When the disc remains warped with no signs of breaking (dotted lines, Fig.
11), $d\cos\theta/dt$ is an increasing function that appears to asymptote to a
value that is in fair agreement from the prediction from Gerosa et al. (2020).
As soon as a disc shows signs of successful breaking (solid lines), our
measure for $d\cos\theta/dt$ develops strong oscillations on a timescale set
by the precession of the inner broken disc, with $d\cos\theta/dt$ at its
lowest when the inner disc(s) most strongly oppose the outer disc.

The amplitude of the oscillations depends on whether the disc has broken into
two smaller discs (upper panel, low $\kappa$) or multiple rings (lower panel,
high $\kappa$). We can understand this by considering Eq. (8), where we
integrate from the inner edge to the outer edge of the disc. When the inner
disc breaks and precesses, it can develop opposing angular momentum to the
outer disc. Thus when we evaluate Eq. (8), it is possible that the
contribution of the inner disc cancels out part of the contribution of the
outer disc in the integral — particularly when we take into account the
$R^{-3}$ dependence of the Lense-Thirring torque.

In the case that there is a single break, the inner disc has a relatively
large radial extent and thus holds an appreciable fraction of the total disc
angular momentum. Our estimation of $d\cos\theta/dt$ for these discs is mostly
governed by the orientation of the inner disc and has a large amplitude. If
instead there are multiple rings, each of these is radially narrower and so
hold a smaller fraction of the total disc angular momentum compared to the
single-break case. Additionally, they precess with a rate determined by their
radius and so have a range of orientations. Configurations where they oppose
each other are also possible. This leads to an oscillatory modulation on the
$d\cos\theta/dt$ profile, with the oscillations reflecting the opposing
angular momenta from the sum of the rings.

This behaviour is replicated across our parameter suite. For discs with one
break we find that $d\cos\theta/dt$ is slowed and alignment between the disc
and BH may be prevented. For discs with multiple rings, alignment between the
outer disc and BH is hindered but not necessarily prevented.

Figure 11: Predicted BH-disc alignment for simulations with one break (upper
panel, $\kappa=1.4\times 10^{-4}$) and multiple breaks (lower panel,
$\kappa=2.6\times 10^{-1}$). In each panel discs that have not broken are
shown with a dotted line. The oscillatory profiles reflect the precession of
the broken disc segments. Here we only show the prograde cases for clarity but
note that the retrograde cases are qualitatively similar.

##  5 Discussion

###  5.1 Limitations of the 1D model

Our comparison in Fig. 3 highlights the limitations of the semi-analytic model
of Gerosa et al. (2020). As expected, at large values of $\kappa$ the boundary
assumptions of the 1D model mean we cannot compare directly. Additionally, we
may expect non-axisymmetric effects to play a major role when the companion
has such a strong influence. Here our simulations still recover combinations
of stable and broken discs, with the broken discs at higher inclinations. Our
results thus suggest that even when the companion has a strong influence, disc
breaking can inhibit the alignment between the disc and the BH. At this
extreme we do note that our assumption of a fixed potential may become
problematic (but see Sec. 5.3).

We do not recover the predicted breaking for very thick ($H/R=0.08$) or large
viscosity discs ($\alpha=0.2$). In the former case, this is likely to be
because the disc is thick enough that it has violated the assumption
inherently made in the semi-analytic model that the disc is thin, even though
the calculation of $\kappa$ depends on the aspect ratio. In the latter case
this is likely due to the inner part of the disc accreting rapidly, preventing
a ring from successfully breaking off.

###  5.2 Numerical viscosity at low $\alpha$

A potential cause for any discrepancy between the 1D semi-analytic model and
the results of our simulations could be our modelled value of $\alpha$. While
Lodato & Price (2010) demonstrates that the $\alpha$ prescription in our
simulations is appropriate (down to smaller $\alpha$ than we have used here),
this has not yet been confirmed in presence of a break. Prior to breaking, the
disc is continuous and resolved and so the viscosity treatment is robust.
After the break has occurred, at the location of the break we have locally
poor resolution which can naturally lead to higher local viscosity with our
chosen viscosity implementation. In practice, this may mean that the
low-$\alpha$ cases might appear to be more viscous than we expect, increasing
the angle of criticality and moving points to larger $\kappa$ in Fig. 3. We
speculate this is the cause of the slight discrepancy observed for
$\alpha=0.05$ and $\kappa=1.4\times 10^{-4}$, but note that this does not
appear to be an issue with any of our other sets.

###  5.3 Assumption of a fixed potential

As introduced in Sec. 3.1, we use a fixed potential to model the primary
rotating BH. Using a fixed potential in this manner is common; for example in
other SPH codes (Nelson & Papaloizou, 2000), in grid codes (Dyda & Reynolds,
2020) and applications other than accretion discs (Bonnerot et al., 2016). A
consequence of this approach is that it is equivalent to assuming the primary
BH is the centre of mass of the binary system. For the parameters we have
chosen (see Table 1), the centre of mass of our simulation is within
$3.8R_{\rm g}$ of the primary BH. This is well within our numerical accretion
radius of $15R_{\rm g}$ and so will not significantly impact the simulated
disc evolution.

###  5.4 Long term behaviour

With (unfortunately) finite computational time at our disposal, we have chosen
to run more simulations with different parameters for a shorter time rather
than fewer simulations for longer times. This means that our simulations run
for a timescale that is much shorter compared to both the viscous time of the
disc or the binary inspiral time. However, for the vast majority of our
simulations we do not require long time-scales to confirm if the disc is
stable to breaking or not.

In the instance that the disc is unstable (i.e. in the region of the parameter
space well beyond criticality), breaking occurs within a few orbits and we do
not have to simulate further. Confirmation that breaking occurs on a few
dynamical time-scales rather than the viscous time-scale can be inferred from
Fig. 8, where we see signatures of breaking for discs with high and low
viscosities at the same number of companion orbits. If instead the disc is
stable, a low amplitude $\psi(R)$ profile (e.g. Fig. 1) develops without local
maxima. As the disc evolves, the $\psi(R)$ profile decreases in amplitude,
moving the disc further towards stability (e.g. Doǧan et al., 2018). With no
changing external torques to influence the disc, this profile continues to
decrease and it becomes increasingly stable against breaking. The difference
in the evolution of the $\psi(R)$ profile for the stable/unstable cases is the
primary feature we test for when we summarise our results in Fig. 3.

The only discs in our suite that would benefit from longer simulation times
are those that show unsuccessful breaking. We identify numerous cases of
unsuccessful breaking where the disc is moving towards breaking, stalls and
then subsequently stabilises. While we have made every effort to rule out that
these discs are simply in the early stages of a successful break, the
intriguing subtleties of this behaviour deserve further attention.

###  5.5 Isothermal equation of state

The semi-analytic model of Gerosa et al. (2020) assumed a locally isothermal
equation of state - that is, the temperature varies as a function of radius
but does not allow for any disc heating. In order to compare with our 3D
numerical simulations, we then assume the same equation of state which results
in the disc thickness being kept constant at a given radius. While this is
valid up until the point in our simulation where the disc breaks, it is not
clear that this holds once the disc has broken. The broken disc geometry
naturally leads to large relative velocities, and when these are of order of
the Keplerian velocity these could easily lead to strong shocks. Such shocks
may heat up the disc, altering the disc geometry and properties significantly.

Although assuming a particular equation of state is unlikely to be correct
after the disc has broken, to date it is the most widely adopted approach
(e.g. Nelson & Papaloizou, 2000; Ogilvie & Latter, 2013; Tremaine & Davis,
2014; Nealon et al., 2015; Liska et al., 2021; Doǧan & Nixon, 2020, amongst
others). Simulations that include radiative transfer will be able to account
for heating due to shocks as well as any viscous heating or cooling that may
occur in the broken disc. While this is not necessary for our comparison with
the semi-analytic model and does not affect our assessment of whether the disc
breaks or not, it should be taken into account when considering subsequent
detailed evolution of broken discs.

##  6 Conclusions

In this work we have considered the structure of a misaligned accretion disc
surrounding BHs in binary systems. With our suite of 143 SPH simulations we
have shown that:

  1. 1.

The ‘critical obliquity’ first identified by Tremaine & Davis (2014) and
explored at length by Gerosa et al. (2020) does indeed correspond to disc
breaking/tearing, where the disc separates into distinct segments.

  2. 2.

Our numerical simulations recover the qualitative predictions of the 1D semi-
analytic model, with the mutual inclination that causes disc breaking
decreasing with increasing $\kappa$ and decreasing $\alpha$.

  3. 3.

At the same time, 3D hydrodynamics allows us to unveil a richer phenomenology.
Disc breaking hinders (and in some cases can prevents) alignment between the
disc and the BH. The difference in this behaviour depends on the whether the
disc breaks into two discs or multiple precessing rings.

  4. 4.

Hydrodynamic effects not taken into account by the semi-analytic model such as
spiral arms are able to stabilise the disc against breaking.

Our results have strong implications for inspiralling binary BHs. Most current
models are based on the idea that the lighter binary member accretes more than
its heavier companion (a.k.a. ‘differential accretion’, e.g. Gerosa et al.
2015; Siwek et al. 2020; Sayeb et al. 2021), predicting that BH binaries
should reach their merger phase with the primary’s (secondary’s) spins aligned
(misaligned) with the orbital angular momentum of the binary. Once the
critical obliquity is included in this picture, disc breaking implies that
spin alignment is slowed if not completely prevented for a specific subset of
systems. Although further modeling is necessary to understand the full
repercussions of our findings, this opens for the exciting prospect of
exploiting future LISA measurements of precessing binary BHs to infer details
on the dynamics of warped discs in gas-rich galaxies.

## Acknowledgements

The authors thank Nicola Giacobbo and the referee, Pavel Ivanov, for
discussions and comments on the manuscript. R.N. acknowledges support from
UKRI/EPSRC through a Stephen Hawking Fellowship (EP/T017287/1). E.R.
acknowledges funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No 681601 and No 864965). D.G. is supported by European Union’s
H2020 ERC Starting Grant No. 945155–GWmining, Leverhulme Trust Grant No.
RPG-2019-350, and Royal Society Grant No. RGS-R2-202004. G.R. acknowledges
support from the Netherlands Organisation for Scientific Research (NWO,
program number 016.Veni.192.233) and from an STFC Ernest Rutherford Fellowship
(grant number ST/T003855/1). This work was performed using the DiRAC Data
Intensive service at Leicester, operated by the University of Leicester IT
Services, which forms part of the STFC DiRAC HPC Facility. The equipment was
funded by BEIS capital funding via STFC capital grants ST/K000373/1 and
ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of
the National e-Infrastructure. We also acknowledge computational resources
from the University of Birmingham BlueBEAR cluster. Figures were made using
splash (Price, 2007).

## Data Availability

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

## References

  * Amaro-Seoane et al. (2017) Amaro-Seoane P., et al., 2017,  (arXiv:1702.00786)
  * Artymowicz & Lubow (1994) Artymowicz P., Lubow S. H., 1994, ApJ, 421, 651
  * Bardeen & Petterson (1975) Bardeen J. M., Petterson J. A., 1975, ApJ, 195, L65
  * Berti & Volonteri (2008) Berti E., Volonteri M., 2008, ApJ, 684, 822 (arXiv:0802.0025)
  * Bogdanović et al. (2007) Bogdanović T., Reynolds C. S., Miller M. C., 2007, ApJ, 661, L147 (arXiv:astro-ph/0703054)
  * Bonnerot et al. (2016) Bonnerot C., Rossi E. M., Lodato G., Price D. J., 2016, MNRAS, 455, 2253 (arXiv:1501.04635)
  * Bregman & Alexander (2012) Bregman M., Alexander T., 2012, ApJ, 748, 63 (arXiv:1109.5384)
  * Caproni et al. (2006) Caproni A., Abraham Z., Mosquera Cuesta H. J., 2006, ApJ, 638, 120 (arXiv:astro-ph/0510276)
  * Caproni et al. (2007) Caproni A., Abraham Z., Livio M., Mosquera Cuesta H. J., 2007, MNRAS, 379, 135 (arXiv:0705.4035)
  * Dotti et al. (2010) Dotti M., Volonteri M., Perego A., Colpi M., Ruszkowski M., Haardt F., 2010, MNRAS, 402, 682 (arXiv:0910.5729)
  * Doǧan & Nixon (2020) Doǧan S., Nixon C. J., 2020, MNRAS, 495, 1148 (arXiv:2005.00548)
  * Doǧan et al. (2015) Doǧan S., Nixon C., King A., Price D. J., 2015, MNRAS, 449, 1251 (arXiv:1502.05410)
  * Doǧan et al. (2018) Doǧan S., Nixon C. J., King A. R., Pringle J. E., 2018, MNRAS, 476, 1519 (arXiv:1801.05426)
  * Dyda & Reynolds (2020) Dyda S., Reynolds C. S., 2020,  (arXiv:2008.12381)
  * Falceta-Gonçalves et al. (2010) Falceta-Gonçalves D., Caproni A., Abraham Z., Teixeira D. M., de Gouveia Dal Pino E. M., 2010, ApJ, 713, L74 (arXiv:1003.2406)
  * Fragile et al. (2001) Fragile P. C., Mathews G. J., Wilson J. R., 2001, ApJ, 553, 955 (arXiv:astro-ph/0007478)
  * Frank et al. (2002) Frank J., King A., Raine D. J., 2002, Accretion Power in Astrophysics. Cambridge University Press
  * Gerosa et al. (2015) Gerosa D., Veronesi B., Lodato G., Rosotti G., 2015, MNRAS, 451, 3941 (arXiv:1503.06807)
  * Gerosa et al. (2020) Gerosa D., Rosotti G., Barbieri R., 2020, MNRAS, 496, 3060 (arXiv:2004.02894)
  * Haiman et al. (2009) Haiman Z., Kocsis B., Menou K., 2009, ApJ, 700, 1952 (arXiv:0904.1383)
  * Ivanov & Illarionov (1997) Ivanov P. B., Illarionov A. F., 1997, MNRAS, 285, 394
  * Kumar & Pringle (1985) Kumar S., Pringle J. E., 1985, MNRAS, 213, 435
  * Lai (1999) Lai D., 1999, ApJ, 524, 1030 (arXiv:astro-ph/9904110)
  * Lei et al. (2013) Lei W.-H., Zhang B., Gao H., 2013, ApJ, 762, 98 (arXiv:1202.4231)
  * Lense & Thirring (1918) Lense J., Thirring H., 1918, Physikalische Zeitschrift, 19, 156
  * Liska et al. (2021) Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S. B., Van Moer M., 2021, MNRAS, 507, 983 (arXiv:1904.08428)
  * Lodato & Gerosa (2013) Lodato G., Gerosa D., 2013, MNRAS, 429, L30 (arXiv:1211.0284)
  * Lodato & Price (2010) Lodato G., Price D. J., 2010, MNRAS, 405, 1212 (arXiv:1002.2973)
  * Lubow & Ogilvie (2000) Lubow S. H., Ogilvie G. I., 2000, ApJ, 538, 326 (arXiv:astro-ph/0003028)
  * Lubow et al. (2002) Lubow S. H., Ogilvie G. I., Pringle J. E., 2002, MNRAS, 337, 706 (arXiv:astro-ph/0208206)
  * Maccarone (2002) Maccarone T. J., 2002, MNRAS, 336, 1371 (arXiv:astro-ph/0209105)
  * Martin et al. (2007) Martin R. G., Pringle J. E., Tout C. A., 2007, MNRAS, 381, 1617 (arXiv:0708.2034)
  * Martin et al. (2008) Martin R. G., Reis R. C., Pringle J. E., 2008, MNRAS, 391, L15 (arXiv:0808.2139)
  * Martin et al. (2009) Martin R. G., Pringle J. E., Tout C. A., 2009, MNRAS, 400, 383 (arXiv:0907.5142)
  * Miller & Krolik (2013) Miller M. C., Krolik J. H., 2013, ApJ, 774, 43 (arXiv:1307.6569)
  * Murray (1996) Murray J. R., 1996, MNRAS, 279, 402 (arXiv:astro-ph/9511031)
  * Natarajan & Pringle (1998) Natarajan P., Pringle J. E., 1998, ApJ, 506, L97 (arXiv:astro-ph/9808187)
  * Nealon et al. (2015) Nealon R., Price D. J., Nixon C. J., 2015, MNRAS, 448, 1526 (arXiv:1501.01687)
  * Nealon et al. (2018) Nealon R., Dipierro G., Alexander R., Martin R. G., Nixon C., 2018, MNRAS, 481, 20 (arXiv:1808.05693)
  * Nelson & Papaloizou (2000) Nelson R. P., Papaloizou J. C. B., 2000, MNRAS, 315, 570 (arXiv:astro-ph/0001439)
  * Nixon et al. (2012) Nixon C., King A., Price D., Frank J., 2012, ApJ, 757, L24 (arXiv:1209.1393)
  * Nixon et al. (2013) Nixon C., King A., Price D., 2013, MNRAS, 434, 1946 (arXiv:1307.0010)
  * Ogilvie (1999) Ogilvie G. I., 1999, MNRAS, 304, 557 (arXiv:astro-ph/9812073)
  * Ogilvie & Latter (2013) Ogilvie G. I., Latter H. N., 2013, MNRAS, 433, 2403 (arXiv:1303.0263)
  * Papaloizou & Pringle (1983) Papaloizou J. C. B., Pringle J. E., 1983, MNRAS, 202, 1181
  * Price (2007) Price D. J., 2007, PASA, 24, 159 (arXiv:0709.0832)
  * Price et al. (2018) Price D. J., et al., 2018, PASA, 35, e031 (arXiv:1702.03930)
  * Pringle (1981) Pringle J. E., 1981, ARA&A, 19, 137
  * Pringle (1996) Pringle J. E., 1996, MNRAS, 281, 357
  * Raj & Nixon (2021) Raj A., Nixon C. J., 2021, ApJ, 909, 82 (arXiv:2101.05825)
  * Rees (1976) Rees M. J., 1976, in Eggleton P., Mitton S., Whelan J., eds, Vol. 73, Structure and Evolution of Close Binary Systems. p. 225
  * Sayeb et al. (2021) Sayeb M., Blecha L., Kelley L. Z., Gerosa D., Kesden M., Thomas J., 2021, MNRAS, 501, 2531 (arXiv:2006.06647)
  * Scheuer & Feiler (1996) Scheuer P. A. G., Feiler R., 1996, MNRAS, 282, 291
  * Sesana et al. (2014) Sesana A., Barausse E., Dotti M., Rossi E. M., 2014, ApJ, 794, 104 (arXiv:1402.7088)
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33
  * Siwek et al. (2020) Siwek M. S., Kelley L. Z., Hernquist L., 2020, MNRAS, 498, 537 (arXiv:2005.09010)
  * Tremaine & Davis (2014) Tremaine S., Davis S. W., 2014, MNRAS, 441, 1408 (arXiv:1308.1964)
  * Tremaine et al. (2009) Tremaine S., Touma J., Namouni F., 2009, AJ, 137, 3706 (arXiv:0809.0237)

##  Appendix A Spin-alignment implementation

Here we describe our method for calculating $d\cos\theta/dt$ as a function of
time for each disc. Importantly, although Eq. (8) could in principle be
calculated as a sum over particles, we instead calculate it over annuli as we
describe here. We begin with a standard discretisation of the disc into $N$
concentric annuli, following Lodato & Price (2010). Here particles are binned
into their respective annulus by their spherical radius because this accounts
for any warped or misaligned disc structure. We assume azimuthal symmetry and
average the properties of the particles in each annulus to recover the disc
properties as a function of radius. While we note this is not strictly
accurate as the spiral arms are not azimuthally symmetric, this best matches
the 1D approach of Gerosa et al. (2020) which we are comparing to and the
perturbations introduced by the spirals represent a small contribution to
$d\cos\theta/dt$.

With this discretisation, one has $\Sigma_{i}$ and $\bm{\hat{l}_{i}}$ for each
annulus $i$. We then write Eq. (8) in its similarly discretised form

$\displaystyle\frac{d\bm{J}}{dt}=\frac{4\pi
G}{c^{2}}\sum_{i}^{N}\frac{\bm{J}\times\bm{L}_{i}}{R_{i}^{2}}\,.$ (16)

where the angular momentum density is given by
$\bm{L}_{i}=\Sigma_{i}\sqrt{GMR_{i}}\bm{\hat{l}_{i}}$ (note this is not the
same as $m\bm{r}\times\bm{v}$ calculated on each individual particle). We then
need to calculate the time derivative of the unit spin angular momentum,

$\displaystyle\frac{d\bm{\hat{J}}}{dt}=\frac{1}{J}\left(\frac{d\textbf{{J}}}{dt}-\bm{\hat{J}}\frac{dJ}{dt}\right)\,,$
(17)

which follows from $\bm{J}=J\bm{\hat{J}}$. The last term in Eq. (17)
represents the change in the magnitude of the BH spin angular momentum.
However, by definition, the Bardeen-Petterson torque ($d\bm{J}/dt$) is always
perpendicular to $\bm{J}$. We thus use only the first term to calculate the
alignment, with

$\displaystyle\frac{\rm d\cos\theta}{\rm dt}=\frac{1}{J}\frac{\rm
d\textbf{{J}}}{\rm dt}\cdot\hat{\textbf{{L}}}_{*}\,.$ (18)

The last term is calculated from
$L_{*}=M_{\star}\bm{R}_{\star}\times\bm{V}_{\star}$, where $V_{\star}$ is the
binary’s velocity relative to the primary. We complete the above procedure for
each output of our simulation, to give $d\cos\theta/dt$ as a function of time
for each disc, as shown in Fig. 11.


$y$ in $N$. Since $W$ is dense in $V_{0}$, there is some $y^{\prime}\in V\cap
W$. By Claim B.17, $y^{\prime}$ has exactly one preimage in $U_{0}$. By Claim
B.10, $y^{\prime}$ has at least one preimage in each of $U_{1},...,U_{l}$.
Since the $U_{i}$ are disjoint, this implies $l=1$, as desired. ∎

By Claim B.18, the proof of Proposition B.1 is complete. ∎

Acknowledgments: Thank you to Assaf Hasson for many helpful discussions of the
material and for reading earlier drafts of the paper. Thanks also to Tom
Scanlon, Rahim Moosa, Ya’acov Peterzil, and Dave Marker for helpful
conversations, and to Chieu-Minh Tran for suggestions on terminology.

## References

  * [BO02] Alessandro Berarducci and Margarita Otero, _o-minimal fundamental group, homology and manifolds_ , Journal of the London Mathematical Society 65 (2002), no. 2, 257–270.
  * [Bou89] Elisabeth Bouscaren, _The group configuration—after E. Hrushovski_ , The model theory of groups (Notre Dame, IN, 1985–1987), Notre Dame Math. Lectures, vol. 11, Univ. Notre Dame Press, Notre Dame, IN, 1989, pp. 199–209.
  * [Cas23] Benjamin Castle, _Zilber’s Restricted Trichotomy in Characteristic Zero_ , Journal of the American Mathematical Society (2023).
  * [CH24] Benjamin Castle and Assaf Hasson, _Very ampleness in strongly minimal sets_ , To appear in Model Theory, 2024.
  * [CHY24] Benjamin Castle, Assaf Hasson, and Jinhe Ye, _Zilber’s trichotomy in hausdorff geometric structures_ , 2024.
  * [DG22] Alfred Dolich and John Goodrick, _Tame topology over definable uniform structures_ , Notre Dame Journal of Formal Logic 63 (2022), no. 1, 51–79.
  * [EHP19] P. Eleftheriou, A. Hasson, and Y. Peterzil, _Strongly minimal groups in o-minimal structures_ , J. Eur. Math. Soc. (JEMS) (2019), (to appear).
  * [EO04] Mário J Edmundo and Margarita Otero, _Definably compact abelian groups_ , Journal of Mathematical Logic 4 (2004), no. 02, 163–180.
  * [HK08] Assaf Hasson and Piotr Kowalski, _Strongly minimal expansions of $(\mathbb{C},+)$ definable in o-minimal fields_, Proc. Lond. Math. Soc. (3) 97 (2008), no. 1, 117–154. MR 2434093 (2009j:03059)
  * [HO10a] Assaf Hasson and Alf Onshuus, _Stable types in rosy theories_ , J. Symbolic Logic 75 (2010), no. 4, 1211–1230. MR 2767965 (2012a:03096)
  * [HO10b] by same author, _Unstable structures definable in o-minimal theories_ , Selecta Math. (N.S.) 16 (2010), no. 1, 121–143. MR 2609645 (2012g:03101)
  * [HOP10a] Assaf Hasson, Alf Onshuus, and Ya’acov Peterzil, _Definable one dimensional structures in o-minimal theories_ , Israel J. Math. 179 (2010), 297–361. MR 2735046 (2012g:03099)
  * [HOP10b] by same author, _Definable structures in o-minimal theories: one dimensional types_ , Israel J. Math. 179 (2010), 363–379. MR 2735047 (2012g:03100)
  * [HOP24] Assaf Hasson, Alf Onshuus, and Santiago Pinzon, _Strongly minimal group relics of algebraically closed valued fields_ , arXiv e-prints (2024), arXiv:2401.14618.
  * [Hru93] Ehud Hrushovski, _A new strongly minimal set_ , Ann. Pure Appl. Logic 62 (1993), no. 2, 147–166, Stability in model theory, III (Trento, 1991). MR 1226304 (94d:03064)
  * [Hru96] by same author, _The Mordell-Lang conjecture for function fields_ , J. Amer. Math. Soc. 9 (1996), no. 3, 667–690. MR 1333294 (97h:11154)
  * [HS17] A. Hasson and D. Sustretov, _Incidence systems on Cartesian powers of algebraic curves_ , ArXiv e-prints (2017).
  * [HZ96] Ehud Hrushovski and Boris Zilber, _Zariski Geometries_ , Journal of the American Mathematical Society 9 (1996), no. 1, 1–55.
  * [Joh24] Will Johnson, _Visceral theories without assumptions_ , 2024.
  * [KR14] P. Kowalski and S. Randriambololona, _Strongly minimal reducts of valued fields_ , J Symbolic logic (2014), to appear.
  * [Mar88] Gary A. Martin, _Definability in reducts of algebraically closed fields_ , J. Symbolic Logic 53 (1988), no. 1, 188–199. MR 929384 (89f:03022)
  * [Mar02] David Marker, _Model theory_ , Graduate Texts in Mathematics, vol. 217, Springer-Verlag, New York, 2002, An introduction.
  * [Mat95] Larry Mathews, _Cell Decomposition and Dimension Functions in First-Order Topological Structures_ , Proceedings of the London Mathematical Society s3-70 (1995), no. 1, 1–32.
  * [Moo04] Rahim Moosa, _A nonstandard Riemann existence theorem_ , Trans. Amer. Math. Soc. 356 (2004), no. 5, 1781–1797. MR 2031041 (2006c:03058)
  * [Moo10] by same author, _Model theory and complex geometry_ , Notices of the AMS 57 (2010), no. 2, 230–235.
  * [MP90] D. Marker and A. Pillay, _Reducts of $({\bf C},+,\cdot)$ which contain $+$_, J. Symbolic Logic 55 (1990), no. 3, 1243–1251. MR 1071326 (92d:03025)
  * [Pa01] RafałPierzchał a, _Invariance of domain in o-minimal structures_ , Ann. Polon. Math. 77 (2001), no. 3, 289–294.
  * [Pil88] Anand Pillay, _On groups and fields definable in $o$-minimal structures_, J. Pure Appl. Algebra 53 (1988), no. 3, 239–255. MR MR961362 (89i:03069)
  * [Pil96] by same author, _Geometric stability theory_ , Oxford Logic Guides, vol. 32, The Clarendon Press, Oxford University Press, New York, 1996, Oxford Science Publications. MR 1429864
  * [PS98] Ya’acov Peterzil and Sergei Starchenko, _A trichotomy theorem for o-minimal structures_ , Proc. London Math. Soc. (3) 77 (1998), no. 3, 481–523. MR 1643405 (2000b:03123)
  * [PS99] Ya’acov Peterzil and Charles Steinhorn, _Definable compactness and definable subgroups of o-minimal groups_ , J. London Math. Soc. (2) 59 (1999), no. 3, 769–786. MR MR1709079 (2000i:03055)
  * [Rab93] E. D. Rabinovich, _Definability of a field in sufficiently rich incidence systems_ , QMW Maths Notes, vol. 14, Queen Mary and Westfield College School of Mathematical Sciences, London, 1993, With an introduction by Wilfrid Hodges. MR 1213456 (94d:03065)
  * [Raz91] Vladimir Razenj, _One-dimensional groups over an $o$-minimal structure_, Ann. Pure Appl. Logic 53 (1991), no. 3, 269–277. MR MR1129780 (93b:03053a)
  * [Str94] Adam W. Strzebonski, _Euler characteristic in semialgebraic and other o-minimal groups_ , Journal of Pure and Applied Algebra 96 (1994), no. 2, 173–201.
  * [vdD98] Lou van den Dries, _Tame topology and o-minimal structures_ , London Mathematical Society Lecture Note Series, vol. 248, Cambridge University Press, Cambridge, 1998. MR 1633348
  * [Zil84] Boris Zilber, _The structure of models of uncountably categorical theories_ , Proceedings of the International Congress of Mathematicians, Vol. 1, 2 (Warsaw, 1983), PWN, Warsaw, 1984, pp. 359–368.
  * [Zil93] by same author, _Model theory and algebraic geometry_ , Proceedings of the 10th Easter Conference on Model Theory, Berlin, 1993, pp. 94720–3480.