# The case of HD 106906 debris disc: A binary’s revenge

Mohammad A. Farhat,2 Antranik A. Sefilian,3,4,5 Jihad R. Touma1,4

1Department of Physics, American University of Beirut PO Box 11-0236, Riad El-
Solh, Beirut 11097 2020, Lebanon
2IMCCE, CNRS, Observatoire de Paris, PSL University, Sorbonne Université, 77
Avenue Denfert-Rochereau, 75014, Paris, France
3Astrophysikalisches Institut und Universitätssternwarte, Friedrich-Schiller-
Universität Jena, Schillergäßchen 2–3, 07745 Jena, Germany
4Center for Advanced Mathematical Sciences, American University of Beirut, PO
Box 11-0236, Riad El-Solh, Beirut 11097 2020, Lebanon
5Departamento de Física, Universidad Técnica Federico Santa María, Av. España
1680, Valparaíso, Chile Alexander von Humboldt Fellow.Corresponding Author:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Debris disc architecture presents [exo-]planetary scientists with precious
clues for processes of planet formation and evolution, including constraints
on planetary mass perturbers (hidden or not). This is particularly true of the
disc in HD 106906 which in early HST, then follow up polarimetric
observations, presented asymmetries and needle-like features that have been
attributed to perturbations by a massive, and unusually distant planetary
companion external to the disc. Here, we revisit the long-term secular
dynamical evolution of the HD 106906 disc allowing for the combined
gravitational action of both the planetary companion, and for the first time
of the inner stellar binary which holds the system together. In doing so, we
draw on Laplace’s venerable insights as recently adapted for wide eccentric
perturbers, Planet 9 included. We show that within uncertainties on the
planet’s orbit, the disc can go from being fully dominated by the inner
binary’s gravitational field to significantly so, and is hardly ever outside
the reach of that binary which has been all but neglected in dedicated
modeling efforts to date. Our analysis predicts a rounded mildly warped (if at
all) disc architecture which appears to be consistent with fresh report on
ALMA observations, and can be further refined in dialog with those
observations. We also employ our results to provide dynamical constraints on
the orbit of the planetary companion.

###### keywords:

planet–disc interactions – planets and satellites: dynamical evolution and
stability – circumstellar matter – stars: individual: HD 106906

††pubyear: 2022††pagerange: The case of HD 106906 debris disc: A binary’s
revenge–A

##  1 Introduction

HD 106906 features a 13 Myr-old tight spectroscopic stellar binary (Pecaut et
al., 2012; De Rosa & Kalas, 2019), a directly imaged planet (Bailey et al.,
2013; Daemgen et al., 2017), and a debris disc in between (Kalas et al., 2015;
Lagrange et al., 2016). In recent years, improved constraints on the planet’s
mass and orbit were accompanied with attempts at understanding disc
architecture (asymmetries and apparent eccentricity in HST then polarimetric
observations; e.g., Crotts et al., 2021), then implications for the origin of
that planet in the first place (Nguyen et al., 2020).

To be sure, planet mass and distance present challenges for run-of-the-mill
planet formation scenarios, with theorists considering the by now familiar
pathways around similar such conundrums: in-situ formation in the course of
binary genesis, capture in the host stellar cluster, formation close-in then
ejection via instability in a multi-planet system, or through direct mean-
motion resonance between migrating planet and binary (e.g., Bate et al., 2002;
Rodet et al., 2017; Jennings & Chiang, 2021). An exploration of the latter
scenario by Rodet et al. (2017) was alone in allowing for the stellar binary’s
gravitational field, mainly on the migrating/ejected planet, but also on the
debris disc in ways that presage our present work [a matter we shall come back
to below].

Otherwise, all studies of debris dynamics in the HD 106906 disc to date have
focused on the planet as as sole perturber gravitational perturber (e.g.,
Jílková & Zwart, 2015; Nesvold et al., 2016, 2017), while ignoring the role of
the binary in shaping the structure of that ring of debris, if not controlling
it altogether. In other words, the inner binary present in HD 106906 was
assumed sufficiently tightly bound for its quadrupole to be rather benign when
compared to the planet’s perturbations. This left the disc in the throngs of
an eccentric Kozai-Lidov trap (Antognini, 2015; Nesvold et al., 2017), through
which it was natural to obtain and argue for significant excitation of
eccentricities and inclinations over the rather young age of the system (13
Myrs). This further appeared consistent with available observations of the
disc’s asymmetries (Kalas et al., 2015; Crotts et al., 2021; van Holstein et
al., 2021), which more often than not were sensitive to dust and its
multifarious perturbations, and were not particularly adept at detecting the
underlying distribution of planetesimals, themselves the subject of modeling.

Figure 1: Laplace Radius Revealed. Plotted versus planet eccentricity $e_{\rm
b}$ is the location of the Laplace radius (Eq. 1) in the HD 106906 system,
computed for various values of the planetary semi-major axis as shown in the
legend. Calculations are done for three different values of the stellar
binary’s semi-major axis as indicated on top of each panel. The variations are
consistent with current observational constraints (Rodet et al., 2017; Nguyen
et al., 2020, see also Table 1). For all of these variations, the resulting
curve is shown to fall within the range of semi-major axes associated with the
debris disc (the gray highlighted region), and to separate it into two
substantial chunks, one dominated by the inner binary and the other by the
planetary companion. Markers correspond to the nominal values of the planetary
parameters ($a_{\rm b}=850{\rm AU},e_{\rm b}=0.44)$ reported in Nguyen et al.
(2020).

We are fresh out of a study of an intimately related system, the Solar system,
where we investigated the fate of Trans-Neptunian debris in the combined field
of the outer gaseous giants inside and a putative super-earth on a wide
inclined and eccentric orbit outside (Farhat & Touma, 2021). In that setting,
the inner and outer perturbers were recognized as comparable and shown to
confine planetesimals into eccentric and inclined frozen orbits. It was thus
natural for us to start by allowing for the inner binary in the HD 106906
system and learn in the process that the binary makes up in mass what it
looses in proximity, rendering its quadrupolar action comparable to the wide
planetary companion, HD 106906 b, over a substantial range of debris semi-
major axes.

Subjecting the debris disc to the combined gravitational influence of the
inner binary and the outer planetary companion gives rise to the so-called
Laplace Surface. The latter, originally posited by Laplace (1805) in his
classic study of Jupiter’s satellites, interpolates a family of planes
parameterized by the satellite’s orbital distance. On these planes, now
referred to as the Laplace planes, competing torques from the Jovian bulge and
the solar tide annihilate on stationary circular orbits with “proper”
inclinations. Spatially, these planes coincide with the planet’s equator on
the inside and tend to the planet’s orbit when moving outwards, smoothly
forming the Laplace Surface.

With a view to greater realism, Tremaine et al. (2009) relaxed the assumption
of circular equilibria, and studied the stability of the classic Laplace
surface to perturbations in the satellite’s eccentricity and to variations in
planetary obliquity. Their results ushered a stream of solar system and exo-
planetary generalizations exploring: the early evolution of the Moon around a
fast spinning oblique Earth (Ćuk et al., 2016; Tian & Wisdom, 2020),
hemispherical asymmetries in Uranian satellites (Charnoz et al., 2018), warps
in circumplanetary discs (Zanazzi & Lai, 2016), and non-gravitational
perturbations such as radiation pressure on dust grains (Rosengren & Scheeres,
2014; Tamayo et al., 2013). More recently, Farhat & Touma (2021) generalized
the framework in Tremaine et al. (2009) by allowing for symmetry breaking
octupolar contributions from an eccentric external perturber. Such
contributions were shown to destroy the classical Laplace surface of circular
orbits, replacing it with an eccentric analogue together with bifurcating
families of spatially frozen orbits at large inclination and eccentricity.

In this work, we highlight the significance and draw the implications of the
eccentric Laplace surface of Farhat & Touma (2021) for debris disc dynamics in
HD 106906. We start with estimates arguing for the applicability of this
dynamics to this remarkable system. We then proceed to deploy the relevant
machinery on the debris disc, identifying equilibrium structures, then
following dynamical evolution starting from an initially cold disc. We subject
the evolving disc to both kinematic and photometric measures. We comment on
disc properties of HD 106906’s age, then consider long-term simulations and
associated observational signatures with similar such systems in mind.
Dynamical evolution is further situated in the appropriate phase space,
arguing that a single geometric structure captures the full dynamics, as it
prepares the stage for future collisional studies of debris discs around a
Laplace surface. The work raises serious questions for HD 106906 debris
modeling to date and should be eminently relevant for emerging ALMA
observations. While deeply rooted in HD 106906 debris disc dynamics, our study
is of broader relevance to discs around binaries with distant companions, and
may once again feed into Planet 9 fantasies closer to home (e.g., Batygin et
al., 2019).

##  2 You may wish the binary away, but the disc knows otherwise: Critical
spatio-temporal scales

We are interested in studying the long-term dynamics of large , i.e., km-
sized, planetesimals constituting the HD 106906 debris disc. Before delving
into detailed calculations, and mainly with fellow observers in mind, we start
by spatially locating the HD 106906 debris disc with respect to the competing
gravitational influences of the inner stellar binary (here captured to
quadrupolar order) and the outer super-Jupiter. What we are inquiring about is
of course germane to the LS associated with this system, and we shall
introduce and map it with sufficient detail in what follows. For now, we
simply point out that a critical role in this balancing act is played by a
distinguished radius, the so called Laplace radius, $r_{\rm L}$ for short.
This radius signals the transition from inner binary to outer-companion
dominated dynamics, a point of inflection in a family of equilibrium orbits
interpolating between the inner binary orbital plane and the plane of the
outer companion’s orbit. The Laplace radius appears as a dimensionless
parameter controlling those equilibria, and in our particular setting takes
the following form:

$r_{\rm L}^{5}=a_{\rm b}^{3}\,a_{\rm AB}^{2}\frac{m_{\rm A}m_{\rm B}}{{m_{\rm
AB}^{2}}}\frac{m_{\rm AB}}{m_{\rm b}}{[1-e_{\rm b}^{2}]}^{\frac{3}{2}},$ (1)

where $m_{\rm b}$ and $a_{\rm b}$ correspond to the outer planet’s mass and
semi-major axis, while $m_{\rm A}$, $m_{\rm B}$, and $a_{\rm AB}$ correspond
to the masses and semi-major axis of the inner stellar binary, with $m_{\rm
AB}=m_{\rm A}+m_{\rm B}$. Taking masses at face value (Table 1), we locate
$r_{\rm L}$ over a range of planet eccentricity and semi-major axis, and for
distinct inner binary semi-major axis within observational uncertainties. The
outcome is shown in Figure 1 with a grey band marking the observed radial
extent of the debris disc.

What stands out neatly in this display is that, within uncertainties on
companion semi-major axis and eccentricity, the Laplace radius for this system
is smack within the disc, moving in or out with decreasing or increasing
planet eccentricity (and the reverse with planet semi-major axis). For mean
planet orbital parameters, the location of the Laplace radius divides the disc
into two chunks, one dominated by the inner binary and the other by the outer
companion. Thus, and going by previous studies of Laplace surface dynamics
(e.g., Tremaine et al., 2009), we expect a warped structure in the outer part
[weak or strong depending on how small or large is the mutual inclination
between planet and stellar binary]. Also, with a potentially eccentric
planetary companion one can foresee in-plane eccentricity forcing in the
warped segment, and motion that largely respects initial conditions in the
inner disc, circular if that is indeed how the disc came to life. Even more
interesting, is the possibility that for near circular planetary orbits of
wide enough semi-major axis, the Laplace radius moves close to the disc’s
outer edge, leaving the inner binary in full control of the disc, dictating in
plane dynamics, all the while significantly hampering warp formation. In sum,
we can safely rule out unhampered Kozai-Lidov dynamics over the full range of
the disc, with it becoming quenched by the binary over a significant inner
range, and for certain planetary orbits, over the whole disc.

Figure 2: Timescales for distinct regimes of debris dynamics. In red, we show
the period of binary-induced nodal precession as a function of planetesimal
semi-major axis, allowing for a range of semi-major axis for the inner binary
and variations of debris orbital elements (Eq. 2). In blue, we show the period
of eccentric Kozai-Lidov oscillations driven by the outer planetary companion
when allowed to act alone, and for that a range of orbital parameters (Eq. 3).
In black we plot the period of small amplitude oscillation around stable
stationary orbits as maintained by binary and companion acting together. It
clearly interpolates between binary and companion controlled timescales.
Interestingly this period is comparable to the age of system in the outer and
inner parts of the disc, and few times longer around the Laplace transition.
The impact on the disc is ultimately a function of the initial state of the
disc, how close or far it is from the system’s Laplace surface.

The LS in question reflects the long-term evolution of the disc, its steady
state architecture. However, the binary, and disc along with it, are a mere 13
Myrs old, so it becomes important to assess relevant timescales in
planetesimal dynamics, partly to situate the disc temporally vis a vis any
steady state, but also to further capture through competing time scales what
we emphasised with spatial considerations. In so doing, it seemed reasonable
to isolate the contribution of each perturber before combining them into a
continuous whole:

  * •

We first consider the stellar binary acting alone as it causes the node of a
planetesimal to precess about the binary’s angular momentum vector with
period: $P_{\rm nodal}=\frac{2\pi}{\dot{\Omega}}$, where the nodal precession
rate, $\dot{\Omega}$, is given by (Farago & Laskar, 2010):

$\dot{\Omega}=-\frac{3}{4}n_{\rm AB}\left(\frac{a_{\rm
AB}}{a}\right)^{7/2}\frac{m_{\rm A}m_{\rm B}}{m_{\rm AB}^{2}}\frac{\cos
i}{(1-e^{2})^{2}}.$ (2)

Here $n_{\rm AB}$ is the mean motion of the inner binary, while $a,e,$ and $i$
refer to the orbital semi-major axis, eccentricity, and inclination of the
planetesimal in question. In Figure 2, we display in red the range of binary-
induced nodal precession timescales for relevant combinations of debris disc
orbital elements (i.e., $a\in[20,150]$ AU; $e\in[0,0.8]$; $i\in[0,80]$) and
variations of $n_{\rm AB}$. Evidently, planetesimals closer to the stellar
binary precess much faster [$\sim 10^{6}$ yrs around 30 AU] than those in the
outskirts of the disc [$\sim 10^{8}$ yrs around 150 AU]. For particles around
the Laplace radius, binary-induced precessional timescale is interestingly
comparable to the age of the system, $\simeq 10-30$ Myrs.

  * •

We then switch off the stellar binary, and consider the forcing of a planetary
companion which when sufficiently inclined (and eccentric) can drive
planetesimals into Kozai-Lidov cycling with an average period of (Antognini,
2015):

$P_{\rm EKL}\approx\frac{256}{15\pi}\sqrt{\frac{10}{\varepsilon_{\rm
o}}}t_{\rm sec},$ (3)

where the secular timescale $t_{\rm sec}$ is given by

$t_{\rm sec}^{-1}=\frac{m_{\rm b}}{a_{\rm b}^{3}(1-e_{\rm
b}^{2})^{3/2}}\sqrt{\frac{Ga^{3}}{m_{\rm AB}}},$ (4)

and the dimensionless coefficient $\varepsilon_{\rm o}$ is defined as

$\varepsilon_{\rm o}=\frac{e_{\rm b}}{1-e_{\rm b}^{2}}\frac{a}{a_{\rm b}}.$
(5)

In Figure 2, we display in blue the range of timescales associated with this
regime. As expected, the timescale of K-L cycling is longest where that of
binary-induced precession is shortest, and shortest in the outer parts of the
disc, where it is slightly longer than the age of the system [$\simeq 2\times
10^{7}$ yrs around 150 AU].

  * •

Last but not least, we consider the combined action of both inner and outer
perturbers, recovering the characteristic timescale of resulting which is here
naturally associated with the small amplitude oscillations around the
equilibrium orbits maintained over the range of semi-major axis of interest
(Tremaine et al., 2009; Farhat & Touma, 2021). That period is marked by a
black curve in Figure 2, which evidently interpolates, as it should, between
binary-induced nodal precession before the Laplace radius (shown with
uncertainties by the grey envelope), then Kozai-Lidov cycling at distances
significantly larger than the Laplace radius. A transitional interval
intimately combines both perturbers into dynamics that shape the behavior of
our disc around the Laplace Surface.

In sum, Figures 1 and 2 illustrate that the HD 106906 debris disc is spatio-
temporally shaped by the combined gravitational perturbation of both binary
and outer perturber, with its relatively young age suggesting a transient
architecture, which we proceed to map out below.

##  3 Towards HD 106906 Debris Dynamics: Secular Dynamical Model and
Associated Equilibria

By now, the reader should find it quite natural, if not unavoidable, to
consider the combined gravitational perturbation of binary and planet as one
explores the dynamics of the debris disc in HD 106906. The setting is not very
different from that of planetary rings, or satellites, in the field of the
oblate planet that hosts them and the Sun around which they revolve. Here, we
summarize key mathematical machinery leading on one hand to yet another
generalized Laplace surface, allowing for an inner stellar binary, with an
eccentric planetary companion, and on the other setting an economical
mathematical model with which to study dynamical evolution of debris (in HD
106906 or any similar such system) in arbitrary detail and for arbitrary
durations.

We introduce an orbital reference frame with the following triad:

  * •

${\boldsymbol{\hat{n}}}$ in the direction of the planetesimal ’s orbital
angular momentum.

  * •

${\boldsymbol{\hat{u}}}$ in the direction of the periapse of its orbit.

  * •

${\boldsymbol{\hat{v}}}={\boldsymbol{\hat{n}}}\times{\boldsymbol{\hat{u}}}$.

Analogous triads are defined in the reference frame of the inner stellar
binary $(\boldsymbol{\hat{n}}_{\rm AB},\boldsymbol{\hat{u}}_{\rm
AB},\boldsymbol{\hat{v}}_{\rm AB})$, and for the outer planetary companion
$(\boldsymbol{\hat{n}}_{\rm b},\boldsymbol{\hat{u}}_{\rm
b},\boldsymbol{\hat{v}}_{\rm b})$.

We are interested in the secular (read orbit averaged) evolution of
planetesimals in the debris disc, planetesimals which are assumed far enough
from binary and planet to be free of mean motion or semi-secular resonances
[the evection being a prime example of the later in hierarchical systems
(Touma & Wisdom, 1998; Touma & Sridhar, 2015)]. Such motion is cleanly
parameterized with the normalized angular momentum vector,
$\boldsymbol{j}=\sqrt{1-e^{2}}\boldsymbol{\hat{n}}$, and the Lenz vector,
$\boldsymbol{e}=e\boldsymbol{\hat{u}}$, with the help of which one avoids
coordinate related singularities at zero eccentricity and/or inclination. With
these vectorial elements, it is straight forward from first principles to
recover the doubly orbit-averaged Hamiltonian of the system, which is now well
documented in various publications on hierarchical triples (e.g., Tremaine et
al., 2009; Correia et al., 2011), and which we define here as:

$\displaystyle H_{S}=-\frac{Gm_{\rm AB}}{2a}+\frac{Gm_{\rm
AB}}{a}\left(\Psi_{\rm AB}+\Psi_{\rm b,quad}+\Psi_{\rm b,oct}\right).$ (6)

In this Hamiltonian, the inner binary’s gravitational field is captured to
quadrupolar order through the averaged disturbing function $\Psi_{\rm AB}$
defined as

$\displaystyle\Psi_{\rm AB}$ $\displaystyle=\frac{3}{8}\frac{\varepsilon_{\rm
AB}}{(1-e^{2})^{\frac{5}{2}}}\Bigg{[}5(\boldsymbol{j}\cdot\boldsymbol{e}_{\rm
AB})^{2}-(1-e_{\rm AB}^{2})(\boldsymbol{j}\cdot\boldsymbol{n}_{\rm AB})^{2}$
(7) $\displaystyle+\left(1-e^{2}\right)\left(\frac{1}{3}-2e_{\rm
AB}^{2}\right)\Bigg{]},$

where we have also defined

$\varepsilon_{\rm AB}=\alpha^{2}\frac{m_{\rm A}m_{\rm B}}{m_{\rm
AB}^{2}}\hskip 14.22636pt\text{with }\hskip 14.22636pt\alpha=\frac{a_{\rm
AB}}{a}.$ (8)

Here it is important to note that although the stellar binary in the HD 106906
system lives on an eccentric orbit (De Rosa & Kalas, 2019), we limit
$\Psi_{\rm AB}$ to quadrupolar order because we are considering a stellar twin
which has a negligible/null octupolar contribution, depending as it does on
the difference in mass between the stellar components (e.g., Correia et al.,
2016). On the other end, we expand the potential of the eccentric inclined
planet to octupolar order keeping in mind the recent work in Farhat & Touma
(2021), which gives vivid illustrations of how going beyond quadrupoles is
essential to Laplace-like dynamics with eccentric outer perturbers. We define:

$\Psi_{\rm b,quad}=\frac{3\varepsilon_{\rm
q}}{8}\Big{[}\frac{1}{3}-2e^{2}+5e_{n}^{2}-j_{n}^{2}\Big{]},$ (9) $\Psi_{\rm
b,oct}=\frac{15}{64}\varepsilon_{\rm q}\varepsilon_{\rm
o}\Big{[}e_{u}\big{(}8e^{2}-1-35e_{n}^{2}+5j_{n}^{2}\big{)}+10j_{u}j_{n}e_{n}\Big{]}.$
(10)

Here we have introduced the quadrupolar dimensionless coefficient

$\varepsilon_{\rm q}=\frac{m_{\rm b}a^{3}}{m_{\rm AB}a_{\rm b}^{3}(1-e_{\rm
b}^{2})^{\frac{3}{2}}},$ (11)

and as before, the octupolar coefficient $\varepsilon_{\rm o}$ as defined in
Eq. (5). In Eqs. (9,10) subscripts affected to $j$ and $e$ reflect bases
vectors of the planetary reference frame on which $\boldsymbol{j}$ and
$\boldsymbol{e}$ are projected respectively. With $H_{\rm S}$ in hand, we have
a minimal model with which to follow dynamics of HD 106906-like debris in a
secular regime. Accordingly, the orientation and shape of planetesimal orbits
are fully described by the evolution of $\boldsymbol{e}$ and $\boldsymbol{j}$,
all the while their semi-major axis is conserved.

Equations of motion are best expressed in a vectorial form generated by
$H_{\rm S}$ (Milankovitch, 1939; Allan & Cook, 1964):

$\displaystyle\frac{d\boldsymbol{j}}{d\tau}$
$\displaystyle=-\boldsymbol{j}\times\nabla_{\boldsymbol{j}}\Psi-\boldsymbol{e}\times\nabla_{\boldsymbol{e}}\Psi,$
$\displaystyle\frac{d\boldsymbol{e}}{d\tau}$
$\displaystyle=-\boldsymbol{e}\times\nabla_{\boldsymbol{j}}\Psi-\boldsymbol{j}\times\nabla_{\boldsymbol{e}}\Psi,$
(12)

with $\tau=\sqrt{{Gm_{\rm AB}}/{a^{3}}}t$, and $\Psi=\Psi_{\rm AB}+\Psi_{\rm
b,quad}+\Psi_{\rm b,oct}$. Expanding, we obtain:

$\displaystyle\frac{d\boldsymbol{j}}{d\tau}$
$\displaystyle=\frac{3\varepsilon_{\rm
AB}}{4(1-e^{2})^{\frac{5}{2}}}\Big{[}(1-e_{\rm
AB}^{2})(\boldsymbol{j}\cdot\boldsymbol{\hat{n}}_{\rm
AB})(\boldsymbol{j}\times\boldsymbol{\hat{n}}_{\rm AB})$
$\displaystyle-5(\boldsymbol{j}\cdot\boldsymbol{e}_{\rm
AB})(\boldsymbol{j}\times\boldsymbol{e}_{\rm
AB})\Big{]}+\frac{3}{4}\varepsilon_{\rm
q}j_{n}\boldsymbol{j}\times\boldsymbol{\hat{n}}_{\rm
b}-\frac{15}{4}\varepsilon_{\rm
q}e_{n}\boldsymbol{e}\times\boldsymbol{\hat{n}}_{\rm b}$
$\displaystyle-\frac{75}{64}\varepsilon_{\rm q}\varepsilon_{\rm
o}\Bigg{\\{}\Big{[}2\big{(}e_{u}j_{n}+e_{n}j_{u}\big{)}\boldsymbol{j}+2\big{(}-7e_{n}e_{u}+j_{u}j_{n}\big{)}\boldsymbol{e}\Big{]}\times\boldsymbol{\hat{n}}_{\rm
b}$
$\displaystyle+\Big{[}2e_{n}j_{n}\boldsymbol{j}+\big{(}-7e_{n}^{2}+j_{n}^{2}+\frac{8}{5}e^{2}-\frac{1}{5}\big{)}\boldsymbol{e}\Big{]}\times\boldsymbol{\hat{u}}_{\rm
b}\Bigg{\\}},$ (13) $\displaystyle\frac{d\boldsymbol{e}}{d\tau}$
$\displaystyle=\frac{3}{4}\varepsilon_{\rm
AB}\Bigg{\\{}\frac{1}{(1-e^{2})^{\frac{5}{2}}}\Big{[}(1-e_{\rm
AB}^{2})(\boldsymbol{j}\cdot\boldsymbol{\hat{n}}_{\rm
AB})(\boldsymbol{e}\times\boldsymbol{\hat{n}}_{\rm AB})$
$\displaystyle-5(\boldsymbol{j}\cdot\boldsymbol{e}_{\rm
AB})(\boldsymbol{e}\times\boldsymbol{e}_{\rm
AB})\Big{]}-\frac{1}{2(1-e^{2})^{\frac{7}{2}}}\Big{[}(1-6e_{\rm
AB}^{2})(1-e^{2})$ $\displaystyle-5(1-e_{\rm
AB}^{2})(\boldsymbol{j}\cdot\boldsymbol{\hat{n}}_{\rm
AB})^{2}+25(\boldsymbol{j}\cdot\boldsymbol{e}_{\rm
AB})^{2}\Big{]}\boldsymbol{j}\times\boldsymbol{e}\Bigg{\\}}$
$\displaystyle+\frac{3}{4}\varepsilon_{\rm
q}j_{n}\boldsymbol{e}\times\boldsymbol{\hat{n}}_{\rm
b}-\frac{15}{4}\varepsilon_{\rm
q}e_{n}\boldsymbol{j}\times\boldsymbol{\hat{n}}_{\rm
b}+\frac{3}{2}\varepsilon_{\rm q}\boldsymbol{j}\times\boldsymbol{e}$
$\displaystyle-\frac{75}{64}\varepsilon_{\rm q}\varepsilon_{\rm
o}\Bigg{\\{}\Big{[}2\big{(}e_{u}j_{n}+e_{n}j_{u}\big{)}\boldsymbol{e}+2\big{(}-7e_{u}e_{n}+j_{u}j_{n}\big{)}\boldsymbol{j}\Big{]}\times\boldsymbol{\hat{n}}_{\rm
b}$
$\displaystyle+\Big{[}2e_{n}j_{n}\boldsymbol{e}+\big{(}-7e_{n}^{2}+j_{n}^{2}+\frac{8}{5}e^{2}-\frac{1}{5}\big{)}\boldsymbol{j}\Big{]}\times\boldsymbol{\hat{u}}_{\rm
b}+\frac{16}{5}e_{u}\boldsymbol{j}\times\boldsymbol{e}\Bigg{\\}}.$ (14)

In this perturbative model, the planet’s presence is felt through both the
axisymmetric (essentially quadrupolar) terms with $\varepsilon_{\rm q}$, and
the symmetry-breaking terms (here limited to octupolar order) with
$\varepsilon_{\rm o}$. The terms in question are identical to those in Farhat
& Touma (2021), where the reader can further find copious details on their
emergence, properties, and significance. The contribution of the inner stellar
binary is captured through the terms with $\varepsilon_{\rm AB}$.

###  3.1 The Eccentric Laplace Surface of HD 106906 System: A brief interlude

Figure 3: An illustration of the Laplace Surface in HD 106906. Planetesimal
rings with angular momentum along $\boldsymbol{\hat{n}}$, are torqued to
precess about the inner stellar binary with angular momentum along
$\boldsymbol{\hat{n}}_{\rm AB}$, and the outer planet with angular momentum
along $\boldsymbol{\hat{n}}_{\rm b}$. The mutual inclination between the inner
and outer perturbers is defined by the angle $i_{\rm b}$. This misaligned
competition leads to precession about the intermediate Laplace plane. When the
latter is defined at each planetesimal position, it creates the sketched
warped surface, with an alignment transition defined by the Laplace radius
$r_{\rm L}$ (Eq. 1). Moving outwards on this surface, planetesimals are also
pumped in eccentricity as also indicated by the color coding. The bottom plot
is an accurate plot of this surface for different values of planetary
inclination. Each surface corresponds to the solution of the Laplace
equilibrium conditions of Eqs. (15-16). Color coding in this panel is
consistent with the upper panel for eccentricity.

With simultaneous perturbations from binary and distant planet, the orbit of a
planetesimal, conveniently envisioned as an eccentric and inclined ring, wants
to precess about the normal to the stellar binary’s plane and the normal to
the planetary companion’s plane simultaneously, with rates that are weighted
by the relative position of the planetesimal. The compromise is that it
precesses around the intermediate Laplace plane, the plane of the equilibrium
orbit at any given location. The envelope shaped by those planes at any
planetesimal position defines the Laplace surface. It was characterized in
full generality (instabilities included) for axisymmetric perturbations in
Tremaine et al. (2009), and later generalized for eccentric perturbations in
Farhat & Touma (2021).

Before embarking on debris dynamics which is our main concern, we report
briefly on debris statics, mapping out Laplace surfaces in representative HD
106906 systems, the fiducial one included. The purpose is two-fold:

  * •

Such orbits represent the skeleton around which debris dynamics is structured;
understanding them as a function of planetary companion orbits reflects
succinctly the location and extent of spatial warp and in plane distortion of
debris distributions;

  * •

Classically, the Laplace surface is where dust, ring systems, planetary
satellites, are expected to settle under adiabatically slow, non-conservative,
processes. With the above discussion of relevant timescales in mind (Figure
2), we expect the inner part of the debris disc to follow the Laplace surface
closely, and anything around the Laplace transition and beyond to be in a
transient non-adiabatic regime, where excitations can be either more
significant or less than indicated by equilibria, and this as a function of
disc initial conditions. In either case, the Laplace surface is again a useful
reference, with detailed observations (as and when they come in) providing an
indication of evolutionary processes that shaped the disc.

Without serious loss of generality, we limit our discussion to coplanar-
coplanar equilibria, in which all the vectors
$\boldsymbol{j},\boldsymbol{e},\boldsymbol{\hat{n}}_{\rm
b},\boldsymbol{\hat{u}}_{\rm b},$ and $\boldsymbol{\hat{n}}_{\rm AB}$ lie in
the same plane. This implies that the nodes of the three players are aligned.
The practical will agree that coplanar-coplanar equilibria are most relevant
here given the fiducial orientation of the planet’s orbit, though nothing
prevents further complexification when the situation calls for it!

Applying this condition to the equations of motion, we recover conditions for
stationary debris orbits in the form of two scalar equations with two
unknowns: the debris eccentricity $e$, and the debris inclination with respect
to the inner plane $i$. For a mutual inclination $i_{\rm b}$ between the
stellar binary’s plane and the planet’s orbit, these conditions read as:

$\displaystyle 0=$ $\displaystyle\frac{3}{4}\frac{\varepsilon_{\rm
AB}}{(1-e^{2})^{\frac{5}{2}}}\Big{[}\frac{1}{2}(1-e^{2})\sin{2i}(4e_{\rm
AB}^{2}+1)\Big{]}$ $\displaystyle-$ $\displaystyle\frac{3}{8}\varepsilon_{\rm
q}(1+4e^{2})\sin 2\phi\mp\frac{75}{64}\varepsilon_{\rm q}\varepsilon_{\rm
o}e\Bigg{[}-\frac{7}{2}(1+e^{2})\cos\phi\sin 2\phi$ $\displaystyle+$
$\displaystyle(2+5e^{2})\sin^{3}\phi+\frac{1}{5}(1-8e^{2})\sin\phi\Bigg{]},$
(15)

$\displaystyle 0=$ $\displaystyle\frac{9}{8}\frac{\varepsilon_{\rm
AB}e}{(1-e^{2})^{2}}\Big{[}-\sin^{2}{i}(4e_{\rm AB}^{2}+1)+e_{\rm
AB}^{2}+2/3\Big{]}$ $\displaystyle+$ $\displaystyle\frac{3}{4}\varepsilon_{\rm
q}e\sqrt{1-e^{2}}(1-4\sin^{2}\phi)$ $\displaystyle\pm$
$\displaystyle\frac{75}{64}\varepsilon_{\rm q}\varepsilon_{\rm
o}(1-e^{2})^{\frac{1}{2}}\Bigg{[}(3e^{2}-1)\cos^{3}\phi$ $\displaystyle+$
$\displaystyle(\frac{1}{5}-\frac{24}{5}e^{2})\cos\phi+(1+\frac{15}{2}e^{2})\sin\phi\sin
2\phi\Bigg{]},$ (16)

where $\phi=i_{\rm b}-i$. Under nodal alignment, equilibria can be either
apsidally aligned or anti-aligned with the outer planet’s periapse. We
distinguish between them by the sign of the octupole terms, with the upper
sign delivering the aligned configuration.

Solving those conditions as a function of planetesimal semi-major axis yield
eccentric Laplace surfaces shown in Figure 3. Surfaces were computed for
increasing companion inclination to binary plane. They clearly interpolate
between an inner circular disc region in near-alignment with binary plane,
with an outer eccentric and inclined region, with mean inclination and
eccentricity growing with the planet’s inclination. Careful treatment will
reveal instabilities in equilibria of a high enough inclination and
eccentricity (and bifurcations of host surfaces), the discussion of which will
take us too far afield. As indicated above, those surfaces, and the
architecture they scaffold, allow economic expression of disc warp and
elongation with distance away from the binary.

The equilibria were here recovered for fixed planet orientation, eccentricity,
semi-major axis and varying inclination (Table 1). They can be further mapped
for varying eccentricity and semi-major axis, shifting the Laplace transition
in or out, fully enclosing the disc in a binary dominated region with little
warp and eccentricity, or further exposing it to the onslaught of the planet.
It should be evident how this profile and its dynamical implications can be
used to rule out planetary configurations otherwise allowed within
uncertainties. We prefer to leave the discussion of such variations and their
signatures to secular dynamical simulations which we report on next.

Figure 4: Snapshots of the secular orbital evolution obtained via direct numerical integration of the eccentricity, inclination, nodal, and apsidal orientation of distributions of mass-less debris particles perturbed by the stellar binary on their inside and the planetary companion on their outside (Eqs. 13 and 14). The perturbers’ orbits are described in Table 1. The disc is initiated with near-circular low-inclination particles with uniform distributions of apsidal and nodal orientations. With the two competing perturbers, and as expected around a Laplace surface, the disc separates into two belts with different dynamical regimes: (i) an inner-binary-dominated belt with near-circular, near-coplanar planetesimals characterized by uniform nodal and apsidal orientations; and (ii) an outer planet-dominated belt that undergoes quasi-periodic eccentricity and inclination pumping. These regimes directly translate into morphological asymmetries, warps, and spiral arms visualized in Figure 5. The evolution is to be compared with that in Figure 7 when we mute the effect of the stellar binary.  Table 1: System orbital parameters.  Parameter | Unit | Value | Reference
---|---|---|---
Age | Myr | $13\pm 2$ | Pecaut et al. (2012)
$m_{\rm A}$ | M⊙ | 1.34 | Rodet et al. (2017)
$m_{\rm B}$ | M⊙ | 1.37 | Rodet et al. (2017)
$a_{\rm AB}$ | AU | 0.36-0.58 | Rodet et al. (2017)
$e_{\rm AB}$ | $-$ | $0.66\pm 0.002$ | De Rosa & Kalas (2019)
$m_{\rm b}$ | MJ | $11\pm 2$ | Bailey et al. (2013)
$a_{\rm b}$ | AU | $850^{+560}_{-260}$ | Nguyen et al. (2020)
$e_{\rm b}$ | $-$ | $0.44^{+0.28}_{-0.31}$ | Nguyen et al. (2020)
$i_{\rm b}$ | deg | $39^{+20}_{-15}$ | Bryan et al. (2021)

##  4 Dynamics around HD 106906 Laplace surface

Having isolated the Laplace surface controlled by inner binary and outer
planet (e.g., Figure 3), our goal now is to study the evolution of a debris
disc by numerically integrating the equations of motion. We perform a suite of
numerical simulations to follow the evolution of the orbital elements of a
distribution of planetesimals living in between the inner and outer
perturbers. Each planetesimal is parameterized by its secularly constant semi-
major axis $a$, and four variable elements: $e,i,\Omega,$ and $\omega$, with
the equations of motion constrained by $j^{2}+e^{2}=1$ and
$\boldsymbol{j}\cdot\boldsymbol{e}=0$.

We list key ingredients in reported simulations:

  * •

As the system is sufficiently hierarchical, we ignore the mutual perturbation
and the consequently ultra-slow precession of the dynamical vectors of the
stellar binary and the planet ($\simeq 10^{10}$ yrs; Antognini, 2015). Thus
$(\boldsymbol{\hat{n}}_{\rm AB},\boldsymbol{\hat{e}}_{\rm AB})$ and
$(\boldsymbol{\hat{n}}_{\rm b},\boldsymbol{\hat{e}}_{\rm b})$ are fixed by the
orbital parameters we set in the simulations and listed in Table 1, and we
further prescribe $(\Omega_{\rm b},\omega_{\rm b})=(0,90)$ degrees for the
planet. The resulting Laplace radius of our fiducial system is thus $r_{\rm
L}=84$ AU;

  * •

Simulations are agnostic of scenarios for the origin of HD 106906b, some of
which might lead to disruption of the debris disc as the planet finding its
way to its current location. Instead, the planet is assumed in its observed
position for the duration of simulations, and unless otherwise stated, with
the mean orbital elements of Table 1;

  * •

We ignore any role (gravitational or otherwise) for any gaseous component
and/or its dispersal on the dynamics of planetesimals in our simulations, and
limit its presence to an item in the discussion section below;

  * •

In this first exploration of an already rich dynamical system, we further
ignore the effect of debris disc self-gravity on Laplace surface dynamics,
relegating on assessment of its significance and potential effects to the
discussion section below;

  * •

Simulations do not allow for radiation pressure, the effect of which is
included in post-processed calculations reported below;

  * •

The debris disc is assumed initially axisymmetric, rather cold, and coplanar
with the binary’s orbital planar. Planetesimal orbits are initiated with
eccentricities ranging between 0 and 0.05 distributed uniformly, and similarly
for inclinations ranging between 0 and 2 degrees. We consider and report on
secular dynamics over a range of semi-major axes which are of relevance to
observations, but nothing prevents extending that range in either direction
while respecting model assumptions.

  * •

In dialogue with the published literature, we further perform and report on
simulations in which we turn off the quadrupolar perturbation due the inner
binary: rather trivial to do numerically, with significant dynamical
implications which we comment on in due course. Results which are pertinent to
HD 106906 are reviewed first, and are then contrasted with simulations that
only account for the binary’s monopole in the discussion section.

With the above assumptions and resulting minimal numerical model, we are able
to zero-in on key features pertaining to dynamics around the Laplace surface
and potential signatures in the disc’s morphology, current and future. Results
of simulations can be viewed through evolving distributions of orbital
elements and their histograms, then surface density/brightness maps at various
stages of evolution. Alternatively, one can view evolution within the system’s
phase space, at distinct and telling semi-major axes. We start with orbits and
their photometric signature (Section 4.1), then move to phase space structure
and the horizons it opens up for future collisional modeling (Section 4.2).

###  4.1 Orbital evolution and surface density

Considered together, results presented in Figs.[4,5] allow the following
observations:

  1. 1.

On Orbital Elements

    1. (a)

Inclinations. Currently and far into the future, the stellar binary’s
quadrupole forces the inner part of the disc to evolve in its orbital plane,
all the while making its signature felt into both the amplitude and frequency
of oscillations in the outer part. This clear break around the system’s
Laplace radius translates into a well defined and dynamically evolving warp
between an inner flat belt and an outer thicker distribution evolving around
the system’s Laplace surface. The break is unavoidable, and is certainly
missed in any modeling that ignores the binary’s non-spherically symmetric
contributions as we shall discuss below.

    2. (b)

Eccentricities. The inclination break around the Laplace radius is further
manifested in the evolution of disc eccentricities. The inner binary’s
quadrupole forces fast precession of inner debris rendering them relatively
immune to the planet’s perturbation, further maintaining the initial nearly
circular condition of the inner belt, while capping the quasi-periodic
eccentricity pumping of the outer part to relatively smaller amplitudes. This
is so over the age of the system, and essentially maintained over a few
hundred million years, except for a burst in eccentricity which reflects long-
term chaotic diffusion in a mixed phase (more on that aspect in Section 4.2).
The effect is familiar from works on the quenching of Kozai-Lidov dynamics by
strong enough quadrupolar perturbations, and is here sustained via Laplace
surface rigidity, modifying excitation amplitudes and frequencies all the
while seeding the frozen orbits of Laplace equilibria.

Figure 5: Surface density maps showing the time evolution of an HD 106906-like
debris disc. The disc is subject to the perturbations of both the inner
stellar binary and the outer planet. Maps at a given time and for a set of
viewing angles are constructed as described in Appendix A using the
distributions of orbital elements resulting from the simulations of Figure 4.
The viewing angles are fixed such that: $\theta_{\rm i}=85^{\circ}$ i.e. the
disc is viewed nearly edge-on, $\theta_{\rm l}=13.7^{\circ}$ as inferred from
a position angle of $283.7^{\circ}$ (Kalas et al., 2015), and $\theta_{\rm
a}=0$ (see Appendix A for further details). Each of these panels have been
normalised individually, and the green star represents the center of mass of
the inner binary. The disc around its current age (13 Myrs) features an inner
bar-like structure, courtesy of the inner binary dominated regime, which
persists on the long-term evolution of the disc. This structure is terminated
the Laplace radius, as the disc warps and transitions to the outer planet
dominated regime. Around the Laplace radius also, transient spiral arms are
created around the system’s current age, and they are twisted with time
following the planet’s periodic forcing. If allowed to live for 300 Myrs, the
separation between the inner and outer regimes becomes more evident: the inner
dense bar-like structure is completely separated from the surrounding
precessing tori. This figure is to be compared with the morphological
evolution of the disc when muting the inner quadrupolar forcing (Figure 9),
with both scenarios further resolved vertically in the histograms of Figure 8.
An animated version of this morphological evolution is provided within the
supplementary material online.

    3. (c)

Nodes and Apses. By 13 Myrs, one notes how nodes of planetesimal orbits
develop a gradual twist with increasing inclination: a telltale feature of the
unfolding warp of classical Laplace dynamics. By then, it is also evident that
the growth of eccentricities in the warped outer part of the disc is
accompanied by the sculpting of the longitude of the apse into one dimensional
manifolds which are suggestive of a developing two-armed spiral. This is while
a nearly axisymmetric coplanar structure is maintained within the Laplace
radius. Further evolution maintains a uniform distribution of the apses in the
inner belt, all the while the spiral of the outer part twists, winds and folds
into a complex triaxial formation.

  2. 2.

On Morphology

We now examine how the dynamics described above and which, in one form or the
other, is controlled by the Laplace surface, manifests through the morphology
of the disc in our fiducial configuration. To do this, we construct in Figure
5 maps of surface density using the orbital element distribution of
planetesimals portrayed in Figure 4. The procedure which we follow extends
what is presented in Sefilian et al. (2021) to account for planetesimal
inclination [refer to Appendix A for further details]. By 13 Myrs, the
system’s current age, the combined gravitational action of binary and outer
planet enforces a clear cut separation between an inner bar-like, binary-
dominated, structure, and an equally clear warp which is seeded around the
system’s Laplace radius. The bar like structure captures in projection the
near circular axisymmetric belt sheltered within the binary’s influence. The
warp reflects the almost linear growth in inclination, with disc particles
responding to both binary and companion torque as they drift vertically and
precess around the local Laplace plane orientation. The warp further features
a spiral arm, reflecting both eccentricity excitation and the self-
organisation of apses discussed above. The long-term dynamical evolution
further exemplifies the morphological features of Laplace dynamics. The inner
quadrupole’s presence maintains the bar-like structure, which further supports
precessing toroidal structures in the outer parts.

In sum, the binary imposes a definite break in the disc with clear
observational signatures of relevance to ongoing and future campaigns: a
distinctive structure in the vertical distribution, with a cold component on
the inside giving way to a hot, evidently warped component on the outside;
near circular orbits on the inside, with spiral arms seeded on the edge of
warps on the outside, the former projecting into a boxy surface density
extending till the Laplace radius on the sky, and the later delineating two
scythe blades pointing away and towards the observer in density maps. We have
much more to say on those observations and variations thereupon, variations
which we relegate to a hefty discussion following yet another and, to our
mind, essential perspective on the debris disc in its phase space.

###  4.2 Phase-space evolution of a cold disc: Look for the separatrix

To investigate the long-term dynamics of debris particles in phase space, we
construct Poincaré sections on the $(\Omega,i)-$plane, defined by $\omega=90$
degrees, i.e. for debris apsidally aligned with the planet. We restrict the
plane crossings by imposing the condition of $\dot{\omega}>0$, and we follow
the orbits evolving according to Eqs. (13-14), starting with initial
conditions prescribed by an energy hypersurface. Each surface of section is
parametrized by the secularly non-diffusive semi-major axis. Thus we present
in Figure 6 samples of sections in this space for different values of debris
semi-majors axes that are representative of the disc studied thus far. We
further color-code for the eccentricity at section crossings. For each semi-
major axis, we construct the section at an energy corresponding to an
initially circular non-inclined debris orbit.

In Figure 6, we start with a section at $a=50$ AU, focusing on debris which
are within the inner quadrupole dominated regime. Debris dynamics on this
section is fairly regular, with circulating trajectories covering almost all
of the section space and conserving their vertical angular momentum. A
libration island appears around low inclinations and $\Omega=0$. Noting that
the section’s energy is not the energy of the Laplace equilibrium, this
identified relative equilibrium is not exactly the Laplace equilibrium, but it
could be mapped into it adiabatically. Other relative equilibria appear at low
inclinations and $\Omega=\pm\pi$, and a clear separatrix emerges in between.
Debris initiated in the plane of the inner binary with near-circular orbits
(color coded with blue) would either librate within the low inclination
equilibria, or would live on the separatrix. Those living on the separatrix
would undergo inclination excursions of a few degrees ($\sim 5^{\circ})$.

Figure 6: Poincaré sections for debris particles at different semi-major axis
values, driven from the inside by the stellar binary and from the outside by a
planetary companion. Initial conditions are selected on an energy manifold
that corresponds to low eccentricity-low inclination orbits. Trajectories are
sectioned in $(\Omega,I)$-plane for $\omega=\pi/2$. Crossings are collected
when $\dot{\omega}>0$. Color coding corresponds to eccentricity. Of
significance to us in these maps are the separatricies joining the nodally
aligned and anti-aligned equilibria. If the disc is initiated as circular and
flat in the binary’s plane, the debris particles at any semi-major axis are
destined to live on the separatrix, following its inclination and, to a lesser
degree, eccentricity excursions. We provide in Eq. (21) a closed form solution
to this structure, and at each semi-major axis, we plot in black the
associated trajectory crossings.

Moving outwards in the disc, we display in the second panel from the left a
section at $a=76$ AU. As we are closer to the Laplace radius, the interplay
between the inner and outer perturbers is now manifest. The previously regular
space of circulating trajectories at higher inclinations is now broken up: a
chaotic zone emerges at high inclination with chains of higher order
resonances embedded within it. Notable is the persistence of libration islands
around nodal alignment and anti-alignment ($\Omega=0,\pm\pi)$, along with the
separatrix in between. However, the libration island by $\Omega=0$ is now
around $i\approx 10$ degrees, while those around $\Omega=\pm\pi$ remain at low
inclinations. Consequently, the separatrix now involves excursions of higher
amplitudes: orbits starting in the plane of the inner binary can now reach
$30$ degrees.

Increasing the semi-major axis further, we enter the outer planet’s dominated
regime. Phase space is now encroached upon by the area-filling chaotic zones,
with regular orbits becoming less prevalent. We note the smooth transport of
the low inclination libration zones from around $\Omega=\pm\pi$ towards
$\Omega=\pm\pi/2$ as we increase $a$. However, the dominant feature of the
phase space, persisting to control the dynamical fate of initially cold orbits
is the structure of the separatrix, only increasing in inclination excursion
amplitude as we increase $a$. Orbits starting in the plane of the inner binary
still live on the separatrix, which also maintains their eccentricity at low
to moderate values, even for high inclinations.

From what preceded, it is evident that structure of the separatrix is critical
to the long-term dynamics of debris discs which are initially coplanar with
the inner binary, on nearly circular orbits. With such configurations and
their long-term dynamics in mind, it would seem sufficient to recover an
approximate expression of the separatrix as a function of semi-major axis, and
this is what we attempt to do next. Tremaine et al. (2009) took a first, and
already challenging, step towards analytically identifying the equilibria and
nonlinear dynamics around them in the limit of circular orbits under
quadrupolar perturbations. We stick to this limit here because orbits on the
separatrix are of relatively low eccentricity. In this limit, particle
dynamics (now fully described by the evolution of the angular momentum vector
$\boldsymbol{j}$) can be mapped into that integrable free-rigid body dynamics
with a suitably defined inertia tensor (refer to Appendix B of Tremaine et al.
(2009) for further details). Adopting this methodology, and taking
$\boldsymbol{\hat{n}}_{\rm AB}=(1,0,0)$, and $\boldsymbol{\hat{n}}_{\rm
b}=(\cos i_{\rm b},\sin i_{\rm b},0)$, Eq. (13) is rewritten as

$\frac{d\boldsymbol{j}}{d\tau}=\frac{3}{4}\boldsymbol{j}\times\boldsymbol{T}\cdot\boldsymbol{j},$
(17)

where the symmetric inertia tensor $\boldsymbol{T}$ is given by

$\boldsymbol{T}=\begin{pmatrix}\varepsilon_{\rm q}\cos^{2}i_{\rm
b}+\varepsilon_{\rm AB}&\varepsilon_{\rm q}\cos i_{\rm b}\sin i_{\rm b}&0\\\
\varepsilon_{\rm q}\cos i_{\rm b}\sin i_{\rm b}&\varepsilon_{\rm
q}\sin^{2}i_{\rm b}&0\\\ 0&0&0\end{pmatrix},$ (18)

with eigenvalues $t_{1}=0<t_{2}<t_{3}$ for $i_{\rm b}>0.$ Denoting by
$\boldsymbol{R}$ the orthogonal transformation that diagonalizes
$\boldsymbol{T}$, we rotate to the principal-axis coordinate system via
$\boldsymbol{j}=\boldsymbol{R}\cdot\boldsymbol{J}$. Under this transformation,
the reduced form of the Hamiltonian in Eq. (6) is written as (Tremaine et al.,
2009):

$\hat{H}=-\frac{1}{2}\boldsymbol{j}^{\rm
T}\cdot\boldsymbol{T}\cdot\boldsymbol{j}=-\frac{1}{2}\sum_{k=1}^{3}t_{k}J_{k}^{2}=-\frac{1}{2}\left(t_{2}J_{2}^{2}+t_{3}J_{3}^{2}\right).$
(19)

This Hamiltonian generates the Lie-Poisson dynamics of the transformed angular
momentum $\boldsymbol{J}$ via

$\frac{d\boldsymbol{J}}{d\tau}=-\boldsymbol{J}\times\nabla_{\boldsymbol{J}}\hat{H},$
(20)

with two integrals of motion: $\hat{H}$ itself, and the magnitude of the
angular momentum $J^{2}$. Orbital trajectories and equilibria can now be
identified as the intersection between the angular momentum sphere and the
energy elliptical cylinder. From Eq. (20), the latter is bounded by the limits
$\hat{H}_{\rm min}=-t_{3}/2$ and $\hat{H}_{\rm max}=0$. Between those limits,
at the critical energy $\hat{H}_{\rm c}=-t_{2}/2$, the cylinder is tangent to
the sphere, and their intersection defines the separatrix. To obtain the
orbits on this separatrix we eliminate $J_{3}^{2}$ between the angular
momentum and energy integrals to obtain:

$(t_{3}-t_{2})J_{2}^{2}+t_{3}J_{1}^{2}=\text{const}=t_{3}-t_{2}.$ (21)

We solve for orbits with angular momentum $\boldsymbol{J}$ on this separatrix,
then we project back to orbits with $\boldsymbol{j}$ using the inverse
transformation. Solutions yield the phase curves shown in black in Figure 6.
Though have obtained under the assumption of circular orbits under dual
quadrupolar forcing, those curves are in remarkable agreement with the
separatrix structure obtained by integrating the full equations of motion.
Specifically, the amplitude of inclination excursions is matched between the
two methods at any semi-major axis.

Figure 7: Similar to Figure 4, but with the outer planet as the sole perturber
to the disc. With the muted effect of the inner stellar binary, the whole
disc, inside to outside, is dominated by the torque of the planet: debris
particles undergo amplified eccentricity oscillations compared to Figure 4,
and are subject to periodic excursions in inclination of amplitude $2i_{\rm
b}$. Figure 8: Inclination distributions of the disc particles derived from
our direct orbital integrations (Figure 4) and sliced at different times. The
two kernel-smoothed curves compare disc vertical distributions in the presence
and absence of the inner binary forcing, using the surface density
distribution described in Appendix A. The clear discrepancy between the two
distributions in each panel signifies taking into account the stellar binary
as an inner perturber rather than a single effective central star. Namely, the
disc under the two forcings exhibits a bimodal distribution in inclination,
with a robust peak around low inclinations, and a fluid peak that shifts
between 20 and 70 degrees as time evolves. This is evidently not the case when
the inner binary is switched off: the bulk of the distribution smoothly
transitions from around low inclinations at $t=0$ Myrs to much higher
inclinations as time evolves.

Phase space structure succinctly captures disc observables as it sets the
stage for long-term collisional evolution of the disc, whether for the cold
initial disc we consider here or for potentially stirred initially
distributions that one may wish to envisage. The separatrix recovered
analytically here controls collisionless cold disc dynamics and will prove
essential for the study of long-term collisional evolution of such discs. Such
dynamics will then play out over the richer neighboring structures including
the potential for chaotic diffusion to higher eccentricities and/or
inclinations. An initially hot disc will start out distributed over phase-
space with collisional dynamics reflecting kinetics over distinct dynamical
regimes, then diffusion between them.

##  5 Variations and Implications: A Discussion

Here we explore modeling assumptions and physical scenarios which put the
results above in dialogue with previous on HD 106906, allow comments on the
vertical structure of disc, explore impact of collisional evolution, and
consider dynamics over uncertainties in planet’s orbit. In other words, as the
title says, variations on the fiducial setup with implications for
observables, of the disc, but also constraints on the planet’s orbit.

Figure 9: Similar to Figure 5, but muting the effect of the stellar binary.
The maps are thus constructed from the orbital elements distributions shown in
Figure 7. Around the system’s age, the whole disc is subject to inclination
pumping, leaving little room for maintaining a cold population, and extending
the warp onto closer-in radii. Allowing the system to evolve longer, it enters
into a steady state where the disc uniformizes in orientations around clear
x-shaped and spiral overdensities. An animated version of this morphological
evolution is provided within the supplementary material online.

###  5.1 Binary Off

Neglecting the binary’s quadrupole, as was mistakenly done in earlier studies,
exposes the disc to the planets onslaught, with effects manifest over the full
extent of the disc (extent which is still hierarchical enough to fall under
secular dynamical considerations):

  1. 1.

With the binary reduced to its monopole, we end up with forced eccentricities
that are twice as large, and over the whole expanse of the disc in question.
This is evident in Fig.7 which further shows how the planet acting alone
forces planetesimals over the full range of the disc to oscillate about its
orbital plane. Inclination evolution can be well approximated with $i=2i_{\rm
b}|\sin\left(\nu t/2\right)|$ (Murray & Dermott, 1999). Here $\nu$, the
secular oscillation frequency, is dependent on the simulated disc lifetime
(Dawson et al., 2011). By 13 Myrs, the age of HD 106906, all planetesimals
have undergone their first oscillation cycle, before reaching a peak
inclination amplitude equal to $2i_{\rm b}$. By 300 Myrs, all planetesimals
are undergoing complete cycles, with those closer to the planet oscillating at
a higher frequency, creating a vertically thick inclined disc. The presence of
the inner perturbation modifies this picture in ways which will have profound
implications for observations of HD 106906, as they pose questions for earlier
approaches to the modeling of this system.

  2. 2.

To further clarify distinctive observable features which are maintained by the
tight binary in this relatively distant debris distribution, we revisit our
simulations via inclination histograms. In each panel of Figure 8, we show the
disc’s inclination distribution, in the presence and absence of the stellar
binary’s quadrupole, adopting the surface density defined in Eq. (28). The
discrepancy between the two settings emerges very clearly beyond 5 Myrs of the
disc’s evolution time. In the presence of the inner binary, the bimodal
distribution appears to be a robust feature which is sustained up to $\sim$500
Myrs. The cold – low inclination – population persists at all times of the
simulated disc’s lifetime with a mean around $\sim$5 degrees. In contrast, the
hot population smoothly shifts its mean back and forth with time, ranging
between 20 and 70 degrees. In the absence of the inner binary’s quadrupole, as
we reported in the previous section, the whole disc is subject to the untamed
periodic inclination excitation by the planet, and it is only a matter of time
before the initially flat disc loses its cold population. As seen in Figure 8,
the density peak around small inclinations smoothly transitions in time
towards high inclinations and persists there, leaving no remarkable features
in the cold end of the distribution.

  3. 3.

Morphologically, when the inner quadrupole is switched off, the warp is much
stronger extending to the inner regions by the age of the system (Figure 9).
On the long-term, a winding spiral shapes the disc’s evolution into a triaxial
distribution which leaves a distinctive x-shaped overdensity in projection.
While an inner bar is typical of Laplace surface dynamics with competing inner
and outer perturbers, the x-shape in the brightness profile of an edge-on disc
is telling of the long-term impact of an inclined massive outer planet (Pearce
& Wyatt, 2014).

Figure 10:  Surface density maps of the disc at its current age as a function
of the system’s Laplace radius. We control the latter by varying the orbital
parameters of the outer planet within the uncertainty intervals reported in
Table 1. These intervals allow for values of $r_{\rm L}$ very close to the
inner and outer edges of the disc (see also Figure 1), thus placing the disc
completely under the dominance of either the stellar binary or the planet.
With $r_{\rm L}\rightarrow a_{\rm in}$, we end up with a planet-dominated disc
which features warps and asymmetries, and for $r_{\rm L}\rightarrow a_{\rm
out}$ we have a fairly rigid, binary-dominated disc which maintains its
initial axisymmetric structure. Figure 11: Constraining the orbital parameter
space of HD 106906 b by disc eccentricity estimates. We contour the surface of
disc eccentricity $\bar{e}_{\rm d}$ computed for a range of values of the
planetary semi-major axis $a_{\rm b}$ and eccentricity $e_{\rm b}$. The latter
ranges are obtained from the uncertainty intervals in Nguyen et al. (2020).
The modelled disc is subject here to pure secular gravitational forcing from
the inner binary and the outer planet as described in our Laplace machinery.
Thus eccentricity estimates obtained here can be, in principle, matched with
those modelled from ALMA observations.

###  5.2 Disc asymmetries, radiation pressure, and constraints on HD 106906 b
orbit

The HD 106906 debris disc was observed in scattered light to feature asymmetry
in brightness, with the southeast extension appearing brighter than the
northwest extension in near-infrared intensity (Kalas et al., 2015; Lagrange
et al., 2016). Kalas et al. (2015) further suggested that the faint west side
extension creates a “needle"-like structure. Observed asymmetries suggest that
the disc is eccentric, leading to pericenter glow which explains brightness
bias, together with a faint tail extending towards the apocenter side of the
disc (Wyatt et al., 1999). More recent polarimetric data analysis over
multiple bands confirmed the asymmetry and suggested that the disc can be
fitted, in surface brightness and structure, to an eccentric ring of
eccentricity $\bar{e}_{\rm d}\geq 0.16$ (Crotts et al., 2021). In contrast to
this structure in the $\mu$m, no obvious asymmetry was detected using ALMA
imaging in the mm (Kral et al., 2020), with expectations for higher resolution
images to resolve the matter once and for all111While this work was in the
final stages of preparation, Fehr et al. (2022) reported on their ALMA
observations of the disc at a resolution of 0$\aas@@fstack{\prime\prime}$38
[39 AU]. They found that the observed structure can be fitted to an
axisymmetric, radially broad disc, without any statistical evidence of an
asymmetry or eccentricity, consistent with the analysis of Kral et al.
(2020)..

In our study thus far, we took the best-fit planetary parameters listed in
Table 1 for granted, ignoring the uncertainties reported therein (Nguyen et
al., 2020). However, within the wide uncertainty intervals in $a_{\rm b}$ and
$e_{\rm b}$, the Laplace radius of the system can take any value between the
reported disc’s inner and outer edges. In Figure 10, we simulate the disc for
various values of $r_{\rm L}$ ranging between 50 AU and 120 AU, by varying the
planetary parameters within the allowed range (Table 1), and we regenerate
surface density maps at the system’s current age. Pushing $r_{\rm L}$ to 50
AU, which would result from $(a_{\rm b},e_{\rm b})=(600$ AU, 0.72), a smaller
chunk of the disc remains controlled by the inner binary, and the larger part
becomes planet-dominated. With the planet’s percienter closer in, the
timescale of the secular forcing by the planet also becomes shorter, thus the
disc structure observed at 100 Myrs in Figure 5 can now be seen at 13 Myrs.
Pushing $r_{\rm L}$ to near the outer edge of the disc, which can be attained
with a planet of $(a_{\rm b},e_{\rm b})=(1150$ AU, 0.25), allows the binary to
overshadow the planet’s effect throughout the entire disc, forcing a near flat
and circular disc. Thus, and going by Kral et al. (2020) who see no
statistically significant indication of disc eccentricity, it would seem that
the system’s Laplace radius sits at or beyond the edge of the debris disc,
hence favoring a rounder and more distant planet 222Fehr et al. (2022) went
through this same exercise and reached a similar conclusion on the planet’s
orbit based on their ALMA observations of the disc..

In sum, we have shown how questions about the existence or not of warps and/or
asymmetries in the disc reduce to an estimate of a single powerful parameter,
the system’s Laplace radius of Eq. (1), which in turn helped constrain the
planet’s orbit with a broad range of uncertainties. A rather powerful
criterion and deserving of further refinements in dialog with emerging and
future observations of this and similar such discs.

For now, and with the indications we have available to us, we undertook a
systematic exploration of the asymmetry in our simulated disc with aim of
correlating global disc eccentricity with the parameters (mainly orbital) of
the planetary companion. We approached the question from two different
perspectives. First, we examined disc eccentricity in a suite of simulations
using our secular dynamical framework. The hope being to couple our parametric
study with ALMA characterization of the disc’s geometry, and in the process
further constrain the planet’s orbital parameters. Then, we allowed for
radiation pressure (in the addition to the purely gravitational interactions
in our model) and undertook the same parametric exploration to compare with
scattered light imaging (Kalas et al., 2015; Crotts et al., 2021).

Within the first approach, and allowing for variations of planetary semi-major
axis and eccentricity within the uncertainty intervals inferred in Nguyen et
al. (2020) and listed in Table 1, we compute the global density-weighed disc
eccentricity $\bar{e}_{\rm d}$ and contour its surface in Figure 11. As
expected, a closer and more eccentric planet generates a more eccentric disc.
For almost half of the parameter space within the uncertainty intervals, the
disc features minimal ($\bar{e}_{\rm d}\leq 0.05$) to null eccentricity.
Meanwhile, a maximum disc eccentricity of $\bar{e}_{\rm d}=0.276$ can be
attained with a planet having $(a_{\rm b},e_{\rm b})=(600$ AU, $0.73)$. The
nominal orbital parameters of the planet from Nguyen et al. (2020) are
$(a_{\rm b},e_{\rm b})=(850$ AU, $0.44)$, corresponding to $\bar{e}_{\rm
d}=0.07$. Exciting the disc eccentricity to $\bar{e}_{\rm d}\approx 0.1$
requires the planetary pericenter to be no farther than 415 AU. For
$\bar{e}_{\rm d}\geq 0.2$, the study constrains the planet to a narrow window
of the parameter space with a maximum pericenter distance of 240 AU. Thus it
would be interesting to see whether future better resolved characterization of
the disc will confirm these constraints or call for more physics into our
modelling.

On a scale smaller than that of the gravitationally interacting planetesimals,
dust particles are spawned and are propelled outwards by stellar radiation
pressure (Burns et al., 1979; Wyatt, 2008). Detailed modeling of this effect
will take us too far afield. Instead, we post process our simulations,
parametrizing the effect through the $\beta$-parameter, the ratio between
radiation pressure and the gravitational forcing, which we take to be
inversely proportional to the geometric cross section of dust grains. Thus on
the disc scale, $\beta$ features a distribution that is related to the size
distribution of dust grains, with smaller particles featuring more sensitivity
to radiation pressure. As a function of parent planetesimal orbital elements
(unprimed), the orbital elements of the dust particles (primed) are defined as
(Burns et al., 1979):

$\displaystyle a^{\prime}$
$\displaystyle=\frac{a(1-e^{2})(1-\beta)}{1-e^{2}-2\beta(1+e\cos f)},$ (22)
$\displaystyle e^{\prime}$ $\displaystyle=\frac{\sqrt{e^{2}+2\beta e\cos
f+\beta^{2}}}{1-\beta},$ (23) $\displaystyle\omega^{\prime}$
$\displaystyle=\omega+\arctan\left(\frac{\beta\sin f}{e+\beta\cos f}\right),$
(24) $\displaystyle i^{\prime}$ $\displaystyle=i,$ (25)
$\displaystyle\Omega^{\prime}$ $\displaystyle=\Omega,$ (26)

where $f$ is the true anomaly of the parent planetesimal, which is drawn in
our computations from a uniform distribution (see Lee & Chiang, 2016, for a
discussion on changing the distribution) .

Allowing for radiation pressure on dust grains, we compute the eccentricity of
the simulated disc that can be imaged in scattered starlight. We do so for
different fixed values of $\beta$ corresponding to different dust grain sizes,
which can in turn be observed at different imaging wavelengths (e.g., Hughes
et al., 2018). We also compute the disc eccentricity assuming the dust size
distribution inferred in Crotts et al. (2021). The authors obtained for the
latter a power law with index $q=3.19^{+0.11}_{-0.20}$, from which we infer
$dN/d\beta\varpropto\beta^{1.19}.$ This distribution is only slightly
different from that of a collisional cascade, which is characterized by
$q\approx 3.5$ (Dohnanyi, 1969). Our $\beta$-distribution extends to a maximum
value, $\beta_{\rm max}$, which marks a marginally bound orbit, and is defined
for every parent planetesimal as:

$\beta_{\rm max}=\frac{1-e^{2}}{2(1+e\cos f)}\approx 0.5+\mathcal{O}(e).$ (27)

In Figure 12, for each prescribed $\beta$ value or distribution, we mark the
contour in $(a_{\rm b},e_{\rm b})-$space that defines the parameter space
required to generate a disc eccentricity $\bar{e}_{\rm d}\geq 0.16$. The
latter value was obtained in Crotts et al. (2021) by fitting the disc’s spine
to an eccentric ring. Evidently, the larger the value of $\beta$ the more
propelled the dust grains are and the more eccentric the disc gets for a given
pair of planetary parameters. Consequently, for $\beta>0.4$, a planet with
semi-major axis and eccentricity of any value within the uncertainty range can
excite the disc’s eccentricity above 0.16. For the $\beta-$distribution
defined above ($dN/d\beta\varpropto\beta^{1.19}$), we are able to constrain
$(a_{\rm b},e_{\rm b})$ to within less than half of the parameter space (the
area above the red contour curve in Figure 12). We find that the nominal
values in Nguyen et al. (2020), $(a_{\rm b},e_{\rm b})=(850$ AU, $0.44)$, are
within this constrained area. We finally note that for this nominal orbit, the
produced disc eccentricity is $\sim 0.07$ when considering planetesimal
gravitational perturbations only, and $\sim 0.19$ when accounting for dust
particles under radiation pressure. Thus the inconsistency between the
reported asymmetry using scattered light (Kalas et al., 2015; Crotts et al.,
2021), and the axisymmetry in ALMA imaging (Kral et al., 2020) could be a
natural outcome of radiation pressure.

Figure 12: Same as Figure 11, but allowing for radiation pressure. After
spawning dust grains from planetesimal progenitors that evolved within our
Laplace dynamical model, we computed the disc’s eccentricity $\bar{e}_{\rm d}$
for a wide range of planetary semi-major axis and eccentricity. Different
curves correspond to different values, then a distribution, of $\beta$, the
ratio between stellar radiation pressure and gravitational forcing. Prescribed
$\beta$-values are randomly selected; the distribution is drawn from the dust
grain size distribution inferred in Crotts et al. (2021). Contours bound from
below the area of orbital parameter space that is sufficient to pump the
disc’s eccentricity above 0.16.

###  5.3 Inclinations Again: Analogues and Potential implications

The simulated inclination profile of constituent planetesimals is a direct
proxy for the vertical structure of the disc. Closer to home, the Kuiper belt
has been the focus of numerous studies attempting to model its vertical
structure (Brown, 2001; Elliot et al., 2005; Kavelaars et al., 2009). Such a
characterization is crucial to better constrain scenarios for the belt’s
formation and issuing dynamical evolution (Gomes, 2003; Levison & Morbidelli,
2003; Nesvornỳ, 2015). Through such studies, it has been established that the
Kuiper belt features two distinct populations, one with low inclination
objects (the cold population), and another with high inclination objects (the
hot population). Those populations are often modelled by the sum of two
overlapping Gaussians (Brown, 2001): a narrow Gaussian for the cold
population, and a broader Gaussian for the hot one, with possible refinements
on the former (e.g., Petit et al., 2011). A similar decomposition has been
attempted for discs beyond our solar system, particularly in the $\beta$
Pictoris system where Matrà et al. (2019) again identified hot and cold
planetesimal populations.

We plot in Figure 13 the population density for the inclination distribution
of the Kuiper belt and $\beta$ Pictoris. We accompany these by the inclination
distribution of our simulated HD 106906 system at its current age. We use the
same surface density distribution used before and described in Appendix A
(i.e., Eq. (28) with $p=0.5$). Our simulated HD 106906 disc features an
analogous density separating a cold population from a hot population. This
bimodal distribution is an inherent aspect of Laplace dynamics as elucidated
thus far. The cold population is predominantly attributed to the inner part of
the disc living in the stellar binary’s dominated regime; where as the hot
population mainly lives in the planetary dominated regime. We fitted this
distribution to two overlapping Rayleigh distributions with a mean and
standard deviation of $1.89$ and $1.24$ degrees for the cold population, and
$28.08$ and $18.39$ degrees for the broader hot population

Figure 13: Inclination distributions as derived and modelled from observations
of the Kuiper belt (Brown, 2001), ALMA observations of $\beta$ Pictoris (Matrà
et al., 2019), and our dynamical simulations of HD 106906, captured at the
system’s age (13 Myr). Each of the three distributions is fitted to two
overlapping Rayleigh distributions describing a high inclination hot
planetesimal population, and a low inclination cold population.

Accounting for uncertainties on the planetary orbit as we did in Section 5.2,
we noted significant shifts in the Laplace radius which can alter this bimodal
distribution significantly. An extreme of course is a system where the
perturbations due to the planet are negligible, which occurs when its true
semi-major axis and eccentricity are towards the upper and lower ends of their
estimated range, respectively (see Figure 1 and last panel of Figure 10). In
this case, the dominance of the inner binary over the planet leads to the
sculpting of a single, rather than two, Gaussian distribution in inclination
centered around $0^{\circ}$. In the opposite case where the planet overshadows
the binary’s gravitational field, the vertical distribution is expected to
show only a single peak, but now centered around large values, at least at the
late times, as seen when we mute the inner binary in Figure 8. Thus, both the
number and positions of peaks in the vertical distribution of planetesimals
could indirectly help constrain the outer planetary orbit, thus complementing
future RV and astrometric measurements of HD 106906 b. Yet another reason to
look forward for future high-resolution ALMA observations resolving the
vertical thickness of the HD 106906 disc (see also Section 5.2, where we
discuss $\bar{e}_{\rm d}$)333The axisymmetric structure of the disc reported
in Kral et al. (2020) and Fehr et al. (2022) favors a Laplace dynamical regime
in which the stellar binary dominates the full disc, which would in turn favor
a dominant cold population in any future sufficiently resolved
characterization of the disc’s vertical structure..

Finally, and conversely, future refinements in the mass and orbital parameters
of HD 106906 b could provide further support [or lack thereof depending on the
findings] for the dominant role of the binary in this specific system. This
represents a promising avenue, especially in light of the advent of JWST which
could directly image planets of $\gtrsim 0.1M_{\rm J}$ at separations of $\sim
100$ au and beyond (e.g., Carter et al., 2021).

Before moving on, we note that our predictions are, to some extent, model-
dependent via the prescribed disc’s innermost and outermost radii (Appendix
A), together with the mass distribution controlled by $p$ (Eq. 28). We have,
however, confirmed that this dependence is weak for all astrophysically-
motivated values of $p$ (i.e., $0<p<2$) as long as $a_{\rm in}\ll a_{\rm
out}$. Finally, we emphasise that our results derive from a collisionless
model, a limitation that must be assessed before making any definitive
conclusions. This is because in reality, collisional activity within the disc
may have a significant – and potentially dominant – influence on the
inclination distribution of the planetesimals. While this is beyond the scope
of our current work, we expect collisions to predominantly deplete the
innermost planetesimals, i.e., those with $a\lesssim r_{\rm L}$, and thus
affect the peak of the cold population in Figure 13. This is simply because
relative velocities of planetesimals would be higher in the inner parts than
in the outer parts (recall that $v_{\rm rel}\propto v_{\rm K}\propto
r^{-1/2}$, where $v_{\rm K}$ is the Keplerian speed). The amount of the
depletion also depends on the micro-physics governing the collisional activity
(Wyatt, 2008), e.g., the strength of the planetesimals, the distribution of
sizes, as well as the maximum size. There is clearly a need for a more
sophisticated analysis once and when the vertical distribution of the HD
106906 disc is sufficiently well-resolved by ALMA.

###  5.4 A variety of vantage points

Figure 14: Synthetic maps of surface density corresponding to HD 106906-like
discs at $t=13$ Myr as viewed from a variety of vantage points. Results are
shown as a function of line-of-sight inclination angle $\theta_{i}$ to the
observer and azimuthal angle $\theta_{a}$; see Section 5.4 for further
details. Each image is constructed following the procedure described in
Appendix A, has 800x800 pixels, and is scaled individually to bring out faint
details. The green star symbol in each panel represent the location of the
inner binary’s centre of mass. One can see the presence of two-armed spiral
pattern for discs which are nearly face-on, i.e., with $\theta_{i}\lesssim
90^{\circ}$, and a warp for nearly edge-on discs, i.e., for $\theta_{i}\sim
90^{\circ}$.

Motivated by the fact that HD 106906 is seen nearly-edge on, we have thus far
focused on describing the vertical structure imprinted by the Laplace surface
dynamics (e.g., Figure 5). For completeness, we now probe the appearance of
the HD 106906 disc at its current age when viewed from a variety of vantage
points. We do so to exemplify observables of Laplace surface dynamics for
similar such dynamical configurations, over a range of potential orientations.
This is illustrated in Figure 14, where we plot the surface density of our
simulated HD 106906 disc as viewed on the sky plane for various combinations
of $\theta_{\rm i}$ and $\theta_{\rm a}$ (see Appendix A). Note that each of
the panels in this figure is normalised individually to bring out the fainter
features, and $\theta_{\rm l}$ is kept fixed444Recall that $\theta_{\rm l}$
measures the angle in the sky-plane from the positive $X$-axis to the disc’s
ascending node on the sky (see Appendix A), and thus varying its value will
not affect the disc structure. at the fiducial value of $13.7^{\circ}$.

Looking at Figure 14, one can see that one of the common features evident in
each of the displayed panels is the little or no evidence of any stellocentric
offset. In other words, the stellar binary’s centre of mass, denoted by a
green star, coincides with the apparent geometric center of the disc’s inner
cavity, as can be seen most prominently for the nearly face-on discs, i.e.,
$\theta_{\rm i}<30^{\circ}$. This simply is a restatement of the fact that the
stellar binary’s perturbations tame the eccentricity excitation by the
external planet, otherwise expected to be significant in the inner parts of
the disc (Figure 4, see also Nesvold et al., 2017).

Another common feature in many of the views displayed in Figure 14 is the
presence of a two-armed spiral pattern. This pattern, as can be seen more
readily for discs viewed with $\theta_{i}\lesssim 90^{\circ}$, starts off
within the disc, i.e., $a_{\rm in}\ll r\ll a_{\rm out}$, and extends out to
the outer disc parts. This can be understood by noting that by this time,
i.e., $t=13$ Myr, the stellar binary has homogenised the apsidal angles of
planetesimals, together with their longitudes of ascending nodes in the inner
disc parts i.e., $a\lesssim r_{\rm L}$, while those at $r\gtrsim r_{\rm L}$
have not had completed one complete precession cycle (see Figures 2 and 4).
Thus, by and large, the spiral pattern demarcates the regions of dominance by
the inner binary and the exterior planet. Note that while this spiral pattern
is expected to be transient in nature, and thus should eventually wrap onto
itself given enough time, Figure 4 indicates that the pattern may survive for
relatively long times, i.e., $t\sim 100$ Myr555This can be also understood by
consulting Figure 2 which shows that the precession timescale for $a\gtrsim
r_{\rm L}$ has a weak dependence on semi-major axis.. Note that similar spiral
patterns were observed by Rodet et al. (2017) who reported on HD 106906 debris
disc dynamics, allowing for the binary and a co-planar planet in a study which
is mainly concerned with a dynamical origin for that planet.

###  5.5 Generality of Laplace surface dynamics

By applying our model to HD 106906, we have shown that Laplace surface
dynamics can play a significant role in shaping this system’s debris disc for
a wide range of system parameters – see e.g. Figure 10 and Section 4. In the
competition between an equal-mass stellar binary interior to the disc and a
super-Jupiter planetary companion exterior to it, and over a wide range of
triple parameters, the disc was shown to sustain a winding spiral within a
twisting warp both having distinctive kinematic and photometric signatures as
discussed above. Within our multiply averaged framework for hierarchical
systems, we can handle a whole slew of configurations, from unequal mass
binaries of varying separation, including a Sun-Jupiter pair, to wide
companions ranging from the planetary to the stellar.

Essentially, it all comes down to the location of the Laplace radius as
captured in Eq.1. With a solar mass external companion, the early spirals and
warps of Laplace surface dynamics are pushed to inner, typical planet forming
regions, with debris in the hundred AU range now fully exposed to the tug of
this stellar companion. On the other hand, Sun and outer solar system can work
with a Neptune-mass ninth planet to sculpt a Laplace surface in trans-
Neptunian realm, i.e. in the few hundred AU range; this range shifts inward
(outward) with decreasing (increasing) companion distance, and/or larger
(smaller) companion eccentricity. Suffice it to say that a sweep over system
parameters will readily reveal the process described here to be generic to
debris discs over a wide range of semi-major axes in triples of various
configurations.

Of interest in this context is the case of HD 100453 which features an inner
binary, a known stellar companion, and a debris disc with a warp, then a gap
and a few Jupiter mass planet hypothesized in that gap, and a substantial
gaseous component (Dong et al., 2015; Van Der Plas et al., 2019). Observations
reveal spirals and a warp with a break around 27 AU (e.g., Wagner et al.,
2015). Associating the break with a Laplace radius which is controlled by
inner planet, and outer stellar companion, we find that a 5 Jupiter mass
planet within the gap will enforce a Laplace transition around 27 AU, where
observations suggest it to be. This straightforward estimate is consistent
with the conclusion of the extensive modeling and simulations in Nealon et al.
(2020), all the while necessitating further elaboration to allow for the
gaseous component.

This state of affairs raises concerns about the dynamical modelling of debris
discs in general which often considers either an inner or an external
perturber (e.g., Wyatt et al., 1999; Chiang et al., 2009; Pearce & Wyatt,
2014; Lee & Chiang, 2016; Nesvold et al., 2017), and proceeds to fit for its
mass, location and/or eccentricity and inclination. Such systems, and
associated modeling, are worth revisiting for sure, allowing for the combined
effect of inner and outer companions (observed or presumed).

We note in passing that for the equal mass binary in HD 106906, octupolar
perturbations cancel out, leaving the included quadrupolar contribution as the
dominant term (Section 3). However, a binary composed of unequal-mass stars,
if and when on an eccentric orbit, can further disturb test-particle dynamics
at the octupolar level potentially carving the axisymmetric sub-region under
its influence. Also, and depending on the strength of the octupole, and the
spatial extent of the debris disc, it could also introduce another level of
competition between the inner and outer components of the planetary system
further shaping the Laplace surface warp. It is straightforward to generalize
the analysis above to such configurations, an extension which we defer to the
notorious category of future works.

We further note that debris disc mass and resulting self-gravity, which we
ignored in this already challenging first look at the problem, could alter the
dynamics reported herein. Disc self-gravity was already highlighted as a key
component of generalized co-planar Laplace dynamics by Sefilian & Touma (2019)
who showed that the combined effect of disc self-gravity with an inner
quadrupole can freeze orbital precession over a range of semi-major axes in a
disc of few Earth masses and moderate eccentricity. It is not difficult to
estimate (Sefilian & Touma, 2019; Sefilian et al., 2021) that an axisymmetric
disc of mass $M_{\rm d}=20M_{\earth}$666For Fehr et al. (2022), the observed
1.3 mm flux indicates a $\sim 10M_{\earth}$ reservoir of planetesimals, an
estimate which surely depends on the assumed size distribution and the maximum
size of planetesimals (Krivov & Wyatt, 2021). yields an apse precession rate
of $\approx-0.09$ ${\rm Myr}^{-1}$ (hence a period of $\approx 70$ ${\rm
Myr}$) around $70$ AU, a rate which, as shown in Figure 2, is comparable in
magnitude to that induced by the combined effect of binary and planet at that
location. The challenge here is to allow for disc self-gravity over a warping,
precessing geometry which is induced by an external eccentric and inclined
companion. We defer full consideration of this potentially important effect to
future work, and simply note for now that a massive disc, which contributes
negatively to apse precession, would usher an inward shift in the Laplace
radius, all else remaining fixed. Such a shift would leave the disc exposed to
twists and warps now unfolding in the combined field of binary, planet, and
disc! By insisting on the disc being maintained in an axisymmetric state, we
would then need to revisit Sec.5.2 constraints on the planet’s orbit to make
sure that the effective Laplace radius is once again beyond the disc’s outer
edge: a game worth playing for sure777The same is true for conclusions reached
by Fehr et al. (2022) assuming a disc of test particles.!

We close this section by emphasising that some of the above variations have
already been explored in the context of the Solar system’s planetary
satellites (Ward, 1981; Saillenfest et al., 2021), and further considered
following the work of Tremaine et al. (2009) for exo-planetary systems with
binary companions, allowing for planet migration, and disc self-gravity as and
when called for (e.g., Zanazzi & Lai, 2016; Spalding & Millholland, 2020;
Speedie & Zanazzi, 2020). Here, we are simply emphasising the same for debris
discs in configurations that allow for substantial companion eccentricity and
inclination.

##  6 Summary and outlook

Much was recognized in common between HD 106906 and our solar system. Both
have a prominent debris disc, with evident inner quadrupolar forcing (HD
106906’s stellar binary, the solar system’s planets). HD 106906 features a
directly imaged, and astrometrically constrained distant planet, surmised to
disturb a young planetesimal/debris disc, which is guessed at through the dust
cloud that shrouds it but otherwise inaccessible to our current observatories.
Our solar system features an exquisitely mapped debris disc with Trans-
Neptunian Objects showing signs of clustering dynamically attributed to
hypothesized and yet to be imaged distant planet.

We had recently generalized the Laplace surface dynamics toolbox to apprehend
and successfully recover key ingredients in support of Planet 9 shepherding of
TNOs (Farhat & Touma, 2021). With the analogies outlined here, it was not
surprising for this toolbox to be particularly handy (with minor
modifications) to address the relatively transient dynamics of the debris disc
in HD 106906, which is perturbed by the inner binary and a directly imaged
distant super-Jupiter. We had developed this machinery to map out generalized
Laplace equilibria and neighboring dynamics with a view to identifying islands
that can harbor a fairly evolved debris population in the solar system. Here,
we deployed it to map out transient dynamics shaping current observations of
HD 106906 debris disc. We further extended analysis, and simulations to allow
for a range of parameters, and system ages, permitting conclusions for a
broader ensemble of exo-planetary systems.

We revisited HD 106906 having convinced ourselves with a back-of-the-envelope
calculation that the inner binary is a significant player in the disc’s
dynamics as the outer planetary companion. For the nominal orbital parameters
of the system, the comparable effects of the inner binary and the outer planet
impose an inevitable sharp break of dynamical regimes, leaving behind an inner
flat and circular region that warps into a vertically thicker and asymmetric
belt in the outskirts, with spiral arms anchored to the edges of the warp.
This clear cut structure gives rise to a bimodal vertical distribution that
mimics that observed in the Kuiper belt and $\beta$ Pictoris. In light of
contrasting observational inferences on the disc’s asymmetry (e.g., Kral et
al., 2020; Crotts et al., 2021), we further studied the disc over the full
range of system parameters, within uncertainties, and confirmed that the disc
can range from being fully dominated by the inner binary, thus forced to live
axi-symmetrically in the binary’s plane, to being significantly so, with the
Laplace radius providing the demarcation line. The major remaining uncertainty
propagates from the planetary orbit, so we mapped the general dependence of
the disc’s asymmetry on the planetary parameters with and without radiation
pressure. Our analysis shows that, for the nominal orbital parameters of the
planet (Nguyen et al., 2020), the disc would appear eccentric in scattered
light and nearly circular in ALMA, thus explaining the discrepancy. It will be
interesting to confirm these maps with improved constrains on the orbit of HD
106906 b, along with more resolved images of the disc. We then closed back the
loop of our dynamical endeavour onto the phase space structure of the disc. We
learned that initially cold discs are destined to live around a robust phase
space structure, specifically a separatrix that carries planetesimals on
excursions of inclination and eccentricity which grow in amplitude as the
planetesimal migrates outwards. We verified the existence of a closed form
solution for this separatrix, and noted how this structure and neighboring
dynamics help capture the long-term evolution of similar such discs.

The hypothesized Planet 9 motivated the development of a toolbox which was
deployed on HD 106906 debris disc. From this work’s perspective on debris disc
dynamics, we are now tempted to revisit the Kuiper belt in the context of a
dynamically evolving solar system and in the presence of a primordial ninth
Planet. We envisage it scultped with double spirals and warps, all the while
the Laplace radius migrates with the outer planets: what role, if any, does
this migration play in the clustering of TNOs around stable high inclination
and eccentricity islands? Can the observed double peak in the inclination
distribution be partly attributed to a warping Laplace surface in the presence
of a distant outer perturber?

In this work, we assumed that HD 106906 b was around in its current position
at the get go. This is certainly the simplest choice when concerned with
secular dynamics of the debris disc. The outcome of our calculations can then
be confronted with detailed observations to decide whether disc structure is
consistent with the gentle secular forcing by a primordial distant planet
formed in-situ (Maury et al., 2019; Jennings & Chiang, 2021), or whether it
holds signs of more violent disruption that could have resulted from a planet
formed close in then disrupted in the course of migration then resonance with
the inner binary (Rodet et al., 2017).

Of course, there is the further complication of the primordial gas on one hand
and the collisional evolution of the debris on the other. The interplay
between secular debris dynamics and the gaseous component of protoplanetary
discs around binaries was clarified by Rafikov & Silsbee (2014), then further
generalized to nonlinear regimes by Sefilian and Touma (In Preparation,
drawing on Sefilian, 2017). Key in those works is the realization that gas
gravity can dramatically alter collision rates; that secular resonances can
sweep through the disc in the course of gas dissipation, altering debris disc
eccentricity (and surely inclination) distributions; and that such sweeping
together with gas drag can significantly sculpt debris distributions.

One naturally wonders how all of this carries over to HD 106906 Laplace
surface dynamics with a disc that is fairly young, apparently massive enough
to be self-gravitating, and has surely emerged from a primordial now largely
dissipated massive gaseous component. The question resonates with work on
warping self-gravitating discs by Toomre (1983), as well as the early and
insightful meditation of Ward (1981) on the role of formative gas in the
inclination of Iapetus. Both studies, with separate motivations, emphasise the
role of disc gravity in shaping the structure of the warp in the presence of
external forcing.

When it comes to HD 106906, we may very well be justified in assuming an
initial debris disc which is cold and co-planar with the binary. This is so
because, as it happens, a million year gas dissipation timescale is at worst
comparable, at best faster than secular timescales around the system’s Laplace
radius. Still, and for reasons discussed above, HD 106906 invites various
critical generalizations of Laplace surface dynamics to allow for a
dissipating gaseous component, disc self-gravity and collisional evolution.
Such generalizations will helps us recover physically secure initial
conditions for our calculations, then draw the implications of disc self-
gravity and collisional evolution for the currently observed structure. More
generally, and perhaps more importantly, they will give us the pleasure of
weaving yet another thread in the rich and layered tapestry initiated by
Laplace more than two centuries ago.

## Acknowledgements

The work resulted from a back-of-the-envelope calculation which came together
around a refined meal at Badguer, Bourj Hammoud, Lebanon. J.R.T. wishes to
express his gratitude to its hostess, Madame Arpi Mangassarian, for sustaining
an island of warmth and serenity in the midst of utter collapse. A.A.S. is
supported by the Alexander von Humboldt Foundation through a Humboldt Research
Fellowship for postdoctoral researchers. A.A.S. also acknowledges partial
financial support from ESO/Gobierno de Chile during the early stages of this
work.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Allan & Cook (1964) Allan R., Cook G., 1964, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 280, 97
  * Antognini (2015) Antognini J. M., 2015, Monthly Notices of the Royal Astronomical Society, 452, 3610
  * Bailey et al. (2013) Bailey V., et al., 2013, The Astrophysical Journal Letters, 780, L4
  * Bate et al. (2002) Bate M. R., Bonnell I. A., Bromm V., 2002, Monthly Notices of the Royal Astronomical Society, 332, L65
  * Batygin et al. (2019) Batygin K., Adams F. C., Brown M. E., Becker J. C., 2019, Physics Reports, 805, 1
  * Brown (2001) Brown M. E., 2001, The Astronomical Journal, 121, 2804
  * Bryan et al. (2021) Bryan M. L., Chiang E., Morley C. V., Mace G. N., Bowler B. P., 2021, The Astronomical Journal, 162, 217
  * Burns et al. (1979) Burns J. A., Lamy P. L., Soter S., 1979, Icarus, 40, 1
  * Carter et al. (2021) Carter A. L., et al., 2021, Monthly Notices of the Royal Astronomical Society, 501, 1999
  * Charnoz et al. (2018) Charnoz S., Canup R. M., Crida A., Dones L., 2018, The Origin of Planetary Ring Systems. Cambridge University Press, p. 517–538
  * Chiang et al. (2009) Chiang E., Kite E., Kalas P., Graham J., Clampin M., 2009, The Astrophysical Journal, 693, 734
  * Correia et al. (2011) Correia A. C., Laskar J., Farago F., Boué G., 2011, Celestial Mechanics and Dynamical Astronomy, 111, 105
  * Correia et al. (2016) Correia A., Boué G., Laskar J., 2016, Celestial Mechanics and Dynamical Astronomy, 126, 189
  * Crotts et al. (2021) Crotts K. A., et al., 2021, The Astrophysical Journal, 915, 58
  * Ćuk et al. (2016) Ćuk M., Hamilton D. P., Lock S. J., Stewart S. T., 2016, Nature, 539, 402
  * Daemgen et al. (2017) Daemgen S., Todorov K., Quanz S. P., Meyer M. R., Mordasini C., Marleau G.-D., Fortney J. J., 2017, Astronomy & Astrophysics, 608, A71
  * Dawson et al. (2011) Dawson R. I., Murray-Clay R. A., Fabrycky D. C., 2011, The Astrophysical Journal Letters, 743, L17
  * De Rosa & Kalas (2019) De Rosa R. J., Kalas P., 2019, The Astronomical Journal, 157, 125
  * Dohnanyi (1969) Dohnanyi J., 1969, Journal of Geophysical Research, 74, 2531
  * Dong et al. (2015) Dong R., Zhu Z., Fung J., Rafikov R., Chiang E., Wagner K., 2015, The Astrophysical Journal Letters, 816, L12
  * Elliot et al. (2005) Elliot J., et al., 2005, The Astronomical Journal, 129, 1117
  * Farago & Laskar (2010) Farago F., Laskar J., 2010, Monthly Notices of the Royal Astronomical Society, 401, 1189
  * Farhat & Touma (2021) Farhat M. A., Touma J. R., 2021, Monthly Notices of the Royal Astronomical Society, 507, 6078
  * Fehr et al. (2022) Fehr A., et al., 2022, arXiv e-prints, p. arXiv:2209.11775
  * Gomes (2003) Gomes R. S., 2003, Icarus, 161, 404
  * Hughes et al. (2018) Hughes A. M., Duchêne G., Matthews B. C., 2018, Annual Review of Astronomy and Astrophysics, 56, 541
  * Jennings & Chiang (2021) Jennings R. M., Chiang E., 2021, Monthly Notices of the Royal Astronomical Society, 507, 5187
  * Jílková & Zwart (2015) Jílková L., Zwart S. P., 2015, Monthly Notices of the Royal Astronomical Society, 451, 804
  * Kalas et al. (2015) Kalas P. G., et al., 2015, The Astrophysical Journal, 814, 32
  * Kavelaars et al. (2009) Kavelaars J., et al., 2009, The Astronomical Journal, 137, 4917
  * Kral et al. (2020) Kral Q., Matra L., Kennedy G. M., Marino S., Wyatt M. C., 2020, Monthly Notices of the Royal Astronomical Society, 497, 2811
  * Krivov & Wyatt (2021) Krivov A. V., Wyatt M. C., 2021, MNRAS, 500, 718
  * Lagrange et al. (2016) Lagrange A. M., et al., 2016, A&A, 586, L8
  * Laplace (1805) Laplace P.-S., 1805, Mecanique Celeste, Vol. 4. Courcier, Paris
  * Lee & Chiang (2016) Lee E. J., Chiang E., 2016, The Astrophysical Journal, 827, 125
  * Levison & Morbidelli (2003) Levison H. F., Morbidelli A., 2003, Nature, 426, 419
  * Matrà et al. (2019) Matrà L., Wyatt M. C., Wilner D. J., Dent W. R., Marino S., Kennedy G. M., Milli J., 2019, The Astronomical Journal, 157, 135
  * Maury et al. (2019) Maury A., et al., 2019, Astronomy & Astrophysics, 621, A76
  * Milankovitch (1939) Milankovitch M., 1939, Bull. Acad. Sci. Math. Nat. A, 6, 1
  * Murray & Dermott (1999) Murray C. D., Dermott S. F., 1999, Solar system dynamics. Cambridge university press
  * Nealon et al. (2020) Nealon R., Cuello N., Gonzalez J.-F., van der Plas G., Pinte C., Alexander R., Ménard F., Price D. J., 2020, Monthly Notices of the Royal Astronomical Society, 499, 3857
  * Nesvold et al. (2016) Nesvold E. R., Naoz S., Vican L., Farr W. M., 2016, The Astrophysical Journal, 826, 19
  * Nesvold et al. (2017) Nesvold E. R., Naoz S., Fitzgerald M. P., 2017, The Astrophysical Journal Letters, 837, L6
  * Nesvornỳ (2015) Nesvornỳ D., 2015, The Astronomical Journal, 150, 68
  * Nguyen et al. (2020) Nguyen M. M., Rosa R. J. D., Kalas P., 2020, The Astronomical Journal, 161, 22
  * Olofsson et al. (2022) Olofsson J., et al., 2022, MNRAS, 513, 713
  * Pearce & Wyatt (2014) Pearce T. D., Wyatt M. C., 2014, Monthly Notices of the Royal Astronomical Society, 443, 2541
  * Pecaut et al. (2012) Pecaut M. J., Mamajek E. E., Bubar E. J., 2012, The Astrophysical Journal, 746, 154
  * Petit et al. (2011) Petit J.-M., et al., 2011, The Astronomical Journal, 142, 131
  * Rafikov & Silsbee (2014) Rafikov R. R., Silsbee K., 2014, The Astrophysical Journal, 798, 70
  * Rodet et al. (2017) Rodet L., Beust H., Bonnefoy M., Lagrange A. M., Galli P. A. B., Ducourant C., Teixeira R., 2017, A&A, 602, A12
  * Rosengren & Scheeres (2014) Rosengren A. J., Scheeres D. J., 2014, The Astrophysical Journal, 786, 45
  * Saillenfest et al. (2021) Saillenfest M., Lari G., Boué G., Courtot A., 2021, Astronomy & Astrophysics, 647, A92
  * Sefilian (2017) Sefilian A. A., 2017, MSc Thesis,American University Of Beirut
  * Sefilian & Touma (2019) Sefilian A. A., Touma J. R., 2019, AJ, 157, 59
  * Sefilian et al. (2021) Sefilian A. A., Rafikov R. R., Wyatt M. C., 2021, ApJ, 910, 13
  * Spalding & Millholland (2020) Spalding C., Millholland S. C., 2020, The Astronomical Journal, 160, 105
  * Speedie & Zanazzi (2020) Speedie J., Zanazzi J., 2020, Monthly Notices of the Royal Astronomical Society, 497, 1870
  * Statler (2001) Statler T. S., 2001, AJ, 122, 2257
  * Tamayo et al. (2013) Tamayo D., Burns J. A., Hamilton D. P., Nicholson P. D., 2013, The Astronomical Journal, 145, 54
  * Tian & Wisdom (2020) Tian Z., Wisdom J., 2020, Proceedings of the National Academy of Sciences, 117, 15460
  * Toomre (1983) Toomre A., 1983, in Symposium-International Astronomical Union. pp 177–185
  * Touma & Sridhar (2015) Touma J. R., Sridhar S., 2015, Nature, 524, 439
  * Touma & Wisdom (1998) Touma J., Wisdom J., 1998, AJ, 115, 1653
  * Tremaine et al. (2009) Tremaine S., Touma J., Namouni F., 2009, The Astronomical Journal, 137, 3706
  * Van Der Plas et al. (2019) Van Der Plas G., et al., 2019, Astronomy & Astrophysics, 624, A33
  * Wagner et al. (2015) Wagner K., Apai D., Kasper M., Robberto M., 2015, The Astrophysical Journal Letters, 813, L2
  * Ward (1981) Ward W. R., 1981, Icarus, 46, 97
  * Wyatt (2008) Wyatt M. C., 2008, ARA&A, 46, 339
  * Wyatt et al. (1999) Wyatt M., Dermott S., Telesco C., Fisher R., Grogan K., Holmes E., Pina R., 1999, The Astrophysical Journal, 527, 918
  * Zanazzi & Lai (2016) Zanazzi J., Lai D., 2016, Monthly Notices of the Royal Astronomical Society, p. stw2629
  * van Holstein et al. (2021) van Holstein R. G., et al., 2021, Astronomy & Astrophysics, 647, A21

##  Appendix A Constructing Maps of Surface Density

Here, we describe our recipe to construct maps of disc surface density using
the distribution of the planetesimal orbital elements. The procedure we follow
is generally similar in spirit to that outlined in Sefilian et al. (2021) for
razor-thin discs (see Appendix C therein), but now generalised to account for
orbital inclinations and viewing angles.

We first begin by assigning a mass $m_{i}$, $i=1,..,N$, to each planetesimal
constituting our simulated debris disc. This is done in such a way that the
initial disc surface density $\Sigma_{\rm d}^{t=0}$ varies with radius $r$ –
or equivalently with semi-major axis $a$, since planetesimals are initiated on
circular and coplanar orbits (Section 4) – following a power-law profile given
by:

$\Sigma_{\rm d}^{t=0}(a)=\Sigma_{0}\left(\frac{a_{\rm out}}{a}\right)^{p},$
(28)

for $a_{\rm in}\leq a\leq a_{\rm out}$. Unless stated otherwise, we use
$p=0.5$, and we fix the innermost and outermost particles to have $a_{\rm
in}=30{\rm au}$ and $a_{\rm out}=150{\rm au}$. Given that the semi-major axis
distribution remains constant within the secular regime, we make use of the
relationship (Statler, 2001):

$dm(a)=2\pi a\Sigma_{\rm d}^{t=0}(a)da,$ (29)

and we compute the amount of mass per unit semi-major axis for a given initial
surface density (Equation 28).

After evolving the planetesimal distribution in time, at a given snapshot of
our simulations, we smear the mass of each planetesimal over its orbit, i.e.,
into an eccentric ring. This is done by spawning $N_{\rm p}=5\times 10^{4}$
new particles for each of the $N$ parent planetesimals, such that the new
particles share the same orbital elements as the parent planetesimals, have
masses $m_{i}/N_{\rm p}$, and mean anomalies $M$ that are randomly selected
from the $[0,2\pi]$ range. This procedure guarantees that each ring is
characterised with a non-uniform linear density reflecting the fact that
particles on eccentric orbits spend more time at their apocentre than at their
pericentre. Next, we compute the eccentric anomaly $E$ for each of the
considered particles, which in total sum up to $N_{\rm t}=N\times N_{\rm p}$,
using Kepler’s equation:

$M=E-e\sin E.$ (30)

This allows us to obtain the position of each particle via standard formulae:

$\displaystyle x$ $\displaystyle=$ $\displaystyle
r\left(\cos\Omega\cos(\omega+f)-\sin\Omega\sin(\omega+f)\cos i\right),$
$\displaystyle y$ $\displaystyle=$ $\displaystyle
r\left(\sin\Omega\cos(\omega+f)+\cos\Omega\sin(\omega+f)\cos i\right),$
$\displaystyle z$ $\displaystyle=$ $\displaystyle r\sin(\omega+f)\sin i,$ (31)

where, as usual, one has $r=a(1-e\cos E)$. For ease of description in what
follows, we refer to the coordinate system $(x,y,z)$ as the “disc-plane”
system.

Next, we transform from the “disc-plane” system to the “sky-plane” coordinate
system, denoted by $(X,Y,Z)$. We do this in such a way that the sky plane is
represented by the $(X,Y)-$plane, with the positive $X-$ and $Y-$axes pointing
west and north, respectively, and the positive $Z-$axis pointing along the
observer’s line of sight. This transformation can be done via:

$\begin{pmatrix}X\\\ Y\\\ Z\end{pmatrix}=\mathcal{R}_{z}(\theta_{\rm
l})\cdot\mathcal{R}_{x}(\theta_{\rm i})\cdot\mathcal{R}_{z}(\theta_{\rm
a})\begin{pmatrix}x\\\ y\\\ z\end{pmatrix},$ (32)

where $\mathcal{R}_{x}(\theta)$ and $\mathcal{R}_{z}(\theta)$ represent
counter-clockwise rotations about the $x-$ and $z-$axes such that:

$\mathcal{R}_{x}(\theta)=\begin{pmatrix}1&0&0\\\ 0&\cos\theta&-\sin\theta\\\
0&\sin\theta&\cos\theta\end{pmatrix},~{}~{}~{}\mathcal{R}_{z}(\theta)=\begin{pmatrix}\cos\theta&-\sin\theta&0\\\
\sin\theta&\cos\theta&0\\\ 0&0&1\end{pmatrix}.$ (33)

In Equation (32), the angle $\theta_{\rm l}$ is the angle in the sky plane,
i.e., in $(X,Y)$, measured from the positive $X-$axis to the disc’s apparent
ascending node on the sky. It can be determined by the apparent major axis of
the disc, which, for HD 106906, has a position angle ${\rm
P.A.}=283.7^{\circ}$ (Kalas et al., 2015). Given that the pericenter side of
the HD 106906 disc appears to be towards the East-West direction, one then has
$\theta_{\rm l}={\rm P.A.}+\pi/2$ (i.e., anti-clockwise from North)888A
similar result could be obtained using the value of $PA=-74.7^{\circ}$ of
Olofsson et al. (2022), in which case $\theta_{\rm l}=-74.7+90=15.3^{\circ}$..
The angle $\theta_{\rm i}$, on the other hand, is a measure of the angle
between the normals to the sky- and disc-plane, with $\theta_{\rm
i}=0^{\circ}(90^{\circ})$ indicating a face-on (edge-on) disc (default
$=85^{\circ}$ in text, except stated otherwise). Finally, the angle
$\theta_{\rm a}$ is measured in the disc plane from the disc’s ascending node
on the sky to the positive $x-$axis (default $=0$ in text, unless stated
otherwise).

Next, we bin the positions of all $N_{\rm t}=N\times N_{\rm p}$ particles in
the sky plane with $[X,Y]=[-200,200]$ AU, with a given resolution of
$800\times 800$ pixels. We then compute the disc’s surface density $\Sigma$ as
projected on the sky by finding the total mass in each bin and dividing it by
the bin’s area.

Finally, we note that we convolve images of surface density with a 2-D
Gaussian having a filter size of $\sim 3$ AU. In all surface density maps
presented throughout the bulk of this paper, (i) the origin is shown by a
green star representing the center of mass of the stellar binary and (ii) the
logarithmic colour scale represents the surface density distribution
normalised by its maximum value. In most of the presented density maps,
normalization was performed individually to bring out the fainter features.