# The case of HD 106906 debris disc: a binary's revenge

Mohammad A. Farhat, Antranik A. Sefilian, Jihad R. Touma

► To cite this version:

Mohammad A. Farhat, Antranik A. Sefilian, Jihad R. Touma. The case of HD 106906 debris disc: a binary's revenge. Monthly Notices of the Royal Astronomical Society, 2023, 521, pp.2067-2086. 10.1093/mnras/stad316. insu-04849172

HAL Id: insu-04849172

https://insu.hal.science/insu-04849172v1

Submitted on 19 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The case of HD 106906 debris disc: a binary’s revenge

Mohammad A. Farhat<sup>1,†</sup>, Antranik A. Sefilian<sup>2,3,4</sup> and Jihad R. Touma<sup>5,3</sup>*

Accepted 2023 January 25. Received 2023 January 25; in original form 2022 October 13

## ABSTRACT

Debris disc architecture presents [exo-]planetary scientists with precious clues for processes of planet formation and evolution, including constraints on planetary mass perturbers. 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 external planetary companion. Here, we revisit the long-term secular dynamical evolution of the HD 106906 disc allowing for the combined gravitational action of the planetary companion and the inner stellar binary which holds the system together. We argue that the binary is strong enough to impose a dynamical break at the disc’s location, resulting in distinctive observational signatures which we render via simulated surface density maps and vertical structure profiles. Within uncertainties on the planet’s orbit, we show that the disc can go from being fully dominated by the inner binary to significantly so, and is hardly ever outside its reach. The extent of binary dominance impacts the disc’s mean eccentricity, a metric which we map as a function of the planet’s semimajor axis and orbital eccentricity, with and without radiation pressure. We can thus constrain the planet’s orbit to ease the tension between evident axisymmetry in the millimeter, and apparent asymmetry in scattered light. We discuss phase space structure, then inclination distribution, arguing for the relevance of our results to a variety of hierarchical systems, as we set the stage for generalizations that allow for disc self-gravity and collisional evolution.

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

## 1 INTRODUCTION

HD 106906 features a 13 Myr-old tight spectroscopic stellar binary (Pecaut, Mamajek & Bubar 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, Rosa & Kalas 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 multiplanet system, or through direct mean-motion resonance between migrating planet and binary (e.g. Bate, Bonnell & Bromm 2002; Rodet et al. 2017; Jennings & Chiang 2021; Moore et al. 2023). 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 the sole gravitational perturber (e.g. Jílková & Zwart 2015; Nesvold et al. 2016; Nesvold, Naoz & Fitzgerald 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 (Naoz 2016; 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 Myr). 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 modelling.

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

<footnote n="1"><sup>1</sup>IMCCE, CNRS, Observatoire de Paris, PSL University, Sorbonne Université, 77 Avenue Denfert-Rochereau, 75014, Paris, France</footnote>

<footnote n="2"><sup>2</sup>Astrophysikalisches Institut und Universitätssternwarte, Friedrich-Schiller-Universität Jena, Schillergäßchen 2–3, 07745 Jena, Germany</footnote>

<footnote n="3"><sup>3</sup>Center for Advanced Mathematical Sciences, American University of Beirut, PO Box 11-0236, Riad El-Solh, Beirut 11097 2020, Lebanon</footnote>

<footnote n="4"><sup>4</sup>Departamento de Física, Universidad Técnica Federico Santa María, Avenida España 1680, Valparaíso, Chile</footnote>

<footnote n="5"><sup>5</sup>Department of Physics, American University of Beirut, PO Box 11-0236, Riad El-Solh, Beirut 11097 2020, Lebanon</footnote>

<footnote>* E-mail: jt00@aub.edu.lb</footnote>

<footnote>† Alexander von Humboldt Fellow.</footnote>
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 loses in proximity, rendering its quadrupolar action comparable to the wide planetary companion, HD 106906 b, over a substantial range of debris semimajor 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 parametrized 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, Touma & Namouni (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 modelling 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, that is, km-sized, planetesimals constituting the HD 106906 debris disc.

<table><caption>**Table 1.** System parameters.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>Value</th><th>Reference</th></tr></thead><tbody><tr><td>Age</td><td>Myr</td><td>13 ± 2</td><td>Pecaut et al. (2012)</td></tr><tr><td>m<sub>A</sub></td><td>M<sub>☉</sub></td><td>1.34</td><td>Rodet et al. (2017)</td></tr><tr><td>m<sub>B</sub></td><td>M<sub>☉</sub></td><td>1.37</td><td>Rodet et al. (2017)</td></tr><tr><td>a<sub>AB</sub></td><td>au</td><td>0.36-0.58</td><td>Rodet et al. (2017)</td></tr><tr><td>e<sub>AB</sub></td><td>—</td><td>0.66 ± 0.002</td><td>De Rosa & Kalas (2019)</td></tr><tr><td>m<sub>b</sub></td><td>M<sub>J</sub></td><td>11 ± 2</td><td>Bailey et al. (2013)</td></tr><tr><td>a<sub>b</sub></td><td>au</td><td>850<sup>+560</sup><sub>-260</sub></td><td>Nguyen et al. (2020)</td></tr><tr><td>e<sub>b</sub></td><td>—</td><td>0.44<sup>+0.28</sup><sub>-0.31</sub></td><td>Nguyen et al. (2020)</td></tr><tr><td>i<sub>b</sub></td><td>deg</td><td>39<sup>+20</sup><sub>-15</sub></td><td>Bryan et al. (2021)</td></tr></tbody></table>

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 Laplace surface 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_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's orbital plane and the plane of the outer companion's orbit. It appears as a dimensionless parameter controlling those equilibria, and in our particular setting takes the form

$$r_L^5 = a_b^3 a_{AB}^2 \frac{m_A m_B}{m_{AB}^2} \frac{m_{AB}}{m_b} [1 - e_b^{2t}]^{\frac{3}{2}}, \quad (1)$$

where $m_b$, $a_b$, and $e_b$ refer to the outer planet's mass, semimajor axis and eccentricity, and $m_A$, $m_B$, and $a_{AB}$ to the masses and semimajor axis of the inner binary, with $m_{AB} = m_A + m_B$.<fnref n="1" /> Taking masses at face value (Table 1), we locate $r_L$ over a range of planet eccentricity and semimajor axis, and for distinct inner binary semimajor axis within observational uncertainties. The outcome is shown in Fig. 1 with a grey band marking the observed radial extent of the debris disc.

What clearly stands out in this display is that, within uncertainties on companion semimajor 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 semimajor 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 a near circular planetary orbit of a wide enough semimajor 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

<footnote n="1">¹Equation (1) is obtained by dimensionally solving for the radius at which the gravitational perturbation of the inner stellar binary balances that of the outer companion (more on that in Section 3).</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Laplace radius revealed. Plotted versus planet eccentricity $e_b$ is the location of the Laplace radius (equation 1) in the HD 106906 system, computed for various values of the planetary semimajor axis as shown in the legend. Calculations are done for three different values of the stellar binary's semimajor 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 semimajor axes associated with the debris disc (the grey 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_b = 850$ au, $e_b = 0.44$) reported in Nguyen et al. (2020).</figcaption></figure>

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.

The Laplace surface in question captures an idealized configuration to which a dissipative particulate disc is expected to settle. However, the binary, and disc along with it, are a mere 13 Myr old, so it becomes important to assess relevant time-scales of 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 emphasized with spatial considerations. In so doing, it seemed reasonable to isolate the contribution of each perturber before combining them into a continuous whole:

(i) 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_{\text{nodal}} = 2\pi/\dot{\Omega}$, where the nodal precession rate, $\dot{\Omega}$, is given by (Farago & Laskar 2010)

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

where $n_{\text{AB}}$ is the mean motion of the inner binary, while $a$, $e$, and $i$ refer to the orbital semimajor axis, eccentricity, and inclination of the planetesimal in question. In Fig. 2, we display in red the range of binary-induced nodal precession time-scales for relevant combinations of debris disc orbital elements (i.e. $a \in [20, 150]$ au; $e \in [0, 0.8]$; $i \in [0, 80]^{\circ}$) and variations of $n_{\text{AB}}$. Evidently, planetesimals closer to the stellar binary precess much faster [$\sim 10^6$ yr around 30 au] than those in the outskirts of the disc [$\sim 10^8$ yr around 150 au]. For particles around the Laplace radius, binary-induced precession time-scale is interestingly comparable to the age of the system, $\simeq 10-30$ Myr.

(ii) We then model the binary as a point mass, and consider the forcing of a planetary companion which when sufficiently inclined (and eccentric) can drive planetesimals into Kozai–Lidov cycling

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Time-scales for distinct regimes of debris dynamics. In red, we show the period of binary-induced nodal precession as a function of planetesimal semimajor axis, allowing for a range of semimajor axis for the inner binary and variations of debris orbital elements (equation 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 same range of orbital parameters (equation 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 time-scales. Interestingly this period is comparable to the age of the 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.</figcaption></figure>
with an average period of (Antognini 2015):

$$P_{\text{EKL}} \approx \frac{256}{15\pi} \sqrt{\frac{10}{\epsilon_o}} t_{\text{sec}}, \qquad (3)$$

where the secular time-scale $t_{\text{sec}}$ is given by

$$t_{\text{sec}}^{-1} = \frac{m_b}{a_b^3(1 - e_b^2)^{3/2}} \sqrt{\frac{Ga^3}{m_{\text{AB}}}}, \qquad (4)$$

and the dimensionless coefficient $\epsilon_o$ is defined as

$$\epsilon_o = \frac{e_b}{1 - e_b^2} \frac{a}{a_b}. \qquad (5)$$

In Fig. 2, we display in blue the range of time-scales associated with this regime. As expected, the time-scale 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 [≈ 2 × 10⁷ yr around 150 au].

(iii) Last but not least, we consider the combined field of both inner and outer perturbers, which imposes a typical oscillation frequency on a planetesimal situated in between.<fnref n="2" /> The position-dependent period of this oscillation is marked by the black curve in Fig. 2. Evidently, this curve interpolates, as it should, between binary-induced nodal precession before the Laplace radius (shown with uncertainties in 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 behaviour of our disc around the Laplace surface.

In sum, Figs 1 and 2 show how the HD 106906 debris disc is spatiotemporally 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:

(i) $\hat{n}$ in the direction of the planetesimal's orbital angular momentum.

(ii) $\hat{u}$ in the direction of the periapse of its orbit.

(iii) $\hat{v} = \hat{n} \times \hat{u}$.

Analogous triads are defined in the reference frame of the inner stellar binary ($\hat{n}_{\text{AB}}$, $\hat{u}_{\text{AB}}$, $\hat{v}_{\text{AB}}$), and for the outer planetary companion ($\hat{n}_b$, $\hat{u}_b$, $\hat{v}_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 semisecular resonances (the ejection being a prime example of the latter in hierarchical systems; Touma & Wisdom 1998; Touma & Sridhar 2015). Such motion is cleanly parametrized with the normalized angular momentum vector, $j = \sqrt{1 - e^2\hat{n}}$, and the Lenz vector, $e = e\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 to recover the Hamiltonian governing the orbit averaged dynamics of disc particles. The procedure, which is briefly outlined in Appendix A, generalizes related work on hierarchical triples (e.g. Tremaine et al. 2009; Correia et al. 2011; Hamers 2020; Farhat & Touma 2021) by allowing for an inner binary and yields the secular Hamiltonian:

$$H_S = -\frac{Gm_{\text{AB}}}{2a} + \frac{Gm_{\text{AB}}}{a} (\Psi_{\text{AB}} + \Psi_{\text{b,quad}} + \Psi_{\text{b,oct}}). \qquad (6)$$

Here, the inner binary's gravitational field is captured to quadrupolar order through the averaged disturbing function

$$\Psi_{\text{AB}} = \frac{3}{8} \frac{\epsilon_{\text{AB}}}{(1 - e^2)^{5/2}} \left[ 5(j \cdot e_{\text{AB}})^2 - (1 - e_{\text{AB}}^2)(j \cdot n_{\text{AB}})^2 + (1 - e^2) \left( \frac{1}{3} - 2e_{\text{AB}}^2 \right) \right], \qquad (7)$$

with dimensionless parameter

$$\epsilon_{\text{AB}} = \alpha^2 \frac{m_{\text{A}}m_{\text{B}}}{m_{\text{AB}}^2}, \quad \text{where} \quad \alpha = \frac{a_{\text{AB}}}{a}. \qquad (8)$$

Note that, although the stellar binary in the HD 106906 system lives on an eccentric orbit (De Rosa & Kalas 2019), we limit $\Psi_{\text{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, Boué & Laskar 2016).

On the other end, we expand the potential of the eccentric inclined planet to octupolar order keeping in mind recent work by Farhat & Touma (2021) which gives vivid illustrations of how going beyond quadrupoles is essential to Laplace-like dynamics with eccentric outer perturbers. The planet's quadrupolar contribution is given by

$$\Psi_{\text{b,quad}} = \frac{3\epsilon_q}{8} \left[ \frac{1}{3} - 2e^2 + 5e_n^2 - j_n^2 \right], \qquad (9)$$

and its orbit averaged octupolar perturbation by

$$\Psi_{\text{b,oct}} = \frac{15}{64} \epsilon_q \epsilon_o \left[ e_u (8e^2 - 1 - 35e_n^2 + 5j_n^2) + 10j_u j_n e_n \right], \qquad (10)$$

introducing the dimensionless parameter

$$\epsilon_q = \frac{m_b a^3}{m_{\text{AB}} a_b^3 (1 - e_b^2)^{3/2}}, \qquad (11)$$

together with the previously defined [equation (5)] octupolar coefficient $\epsilon_o$. In equations (9, 10) subscripts affected to $j$ and $e$ reflect bases vectors of the planetary reference frame on which $j$ and $e$ are projected, respectively. Using the expressions above, we can readily express the Laplace radius of equation (1) in terms of the ratio of inner to outer quadrupolar dimensionless parameters via $r_L^5 = a^5 \epsilon_{AB} / \epsilon_q$.

It is perhaps worth emphasizing here that one should use a multipolar expansion for the planet's perturbation with caution. Considering the nominal orbital configuration of HD 106906 (refer to Section 3.1), it is safe to assume that the system is sufficiently hierarchical for the multipolar expansion to converge. However, this

<footnote n="2">²The frequency in question is that of small amplitude oscillations around the equilibrium orbits maintained over the range of semimajor axis of interest (see section 3 of Tremaine et al. 2009).</footnote>
hierarchy may be compromised with a closer-in and highly eccentric planet, in which case the multipolar series diverges and fails. In this case, one can either resort to exact numerical averaging of the Hamiltonian [as was done in Beust (2016), Saillenfest et al. (2017), Batygin & Morbidelli (2017) for the case of Planet 9], or to the numerical averaging of the spherical harmonics of the planet's potential when modelled as a Gaussian ring (Farhat & Touma 2021).

With $H_S$ in hand, we have a minimal model with which to follow the dynamics of HD 106906-like debris in the secular orbit-averaged regime, one where the orientation and shape of planetesimal orbits are fully described by the evolution of $e$ and $j$, all the while their semimajor axis is conserved. Equations of motion are best expressed in the vectorial form (Milankovitch 1939; Allan & Cook 1964):

$$ 
\begin{aligned}
\frac{d\mathbf{j}}{d\tau} &= -\mathbf{j} \times \nabla_j \Psi - \mathbf{e} \times \nabla_e \Psi, \\
\frac{d\mathbf{e}}{d\tau} &= -\mathbf{e} \times \nabla_j \Psi - \mathbf{j} \times \nabla_e \Psi,
\end{aligned}
\qquad (12) $$

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

$$ 
\begin{aligned}
\frac{d\mathbf{j}}{d\tau} &= \frac{3\epsilon_{AB}}{4(1-e^2)^{\frac{5}{2}}} \left[ (1-e_{AB}^2)(\mathbf{j} \cdot \hat{\mathbf{n}}_{AB})(\mathbf{j} \times \hat{\mathbf{n}}_{AB}) \right. \\
&\qquad \left. -5(\mathbf{j} \cdot \mathbf{e}_{AB})(\mathbf{j} \times \mathbf{e}_{AB}) \right] + \frac{3}{4}\epsilon_{q,jn}\mathbf{j} \times \hat{\mathbf{n}}_{b} - \frac{15}{4}\epsilon_{q}e_{n}\mathbf{e} \times \hat{\mathbf{n}}_{b} \\
&\qquad -\frac{75}{64}\epsilon_{q}\epsilon_{o}\left\{ \left[ 2(e_{u}j_{n}+e_{n}j_{u})\mathbf{j} + 2(-7e_{n}e_{u}+j_{u}j_{n})\mathbf{e} \right] \right. \\
&\qquad \left. \times \hat{\mathbf{n}}_{b} + \left[ 2e_{n}j_{n}\mathbf{j} + \left(-7e_{n}^2 + j_{n}^2 + \frac{8}{5}e^2 - \frac{1}{5}\right)\mathbf{e} \right] \times \hat{\mathbf{u}}_{b} \right\},
\end{aligned}
\qquad (13) $$

$$ 
\begin{aligned}
\frac{d\mathbf{e}}{d\tau} &= \frac{3}{4}\epsilon_{AB} \left\{ \frac{1}{(1-e^2)^{\frac{5}{2}}} \left[ (1-e_{AB}^2)(\mathbf{j} \cdot \hat{\mathbf{n}}_{AB})(\mathbf{e} \times \hat{\mathbf{n}}_{AB}) \right. \right. \\
&\qquad \left. -5(\mathbf{j} \cdot \mathbf{e}_{AB})(\mathbf{e} \times \mathbf{e}_{AB}) \right] - \frac{1}{2(1-e^2)^{\frac{5}{2}}} \left[ (1-6e_{AB}^2)(1-e^2) \right. \\
&\qquad \left. -5(1-e_{AB}^2)(\mathbf{j} \cdot \hat{\mathbf{n}}_{AB})^2 + 25(\mathbf{j} \cdot \mathbf{e}_{AB})^2 \right] \mathbf{j} \times \mathbf{e} \\
&\qquad + \frac{3}{4}\epsilon_{q,jn}\mathbf{e} \times \hat{\mathbf{n}}_{b} - \frac{15}{4}\epsilon_{q}e_{n}\mathbf{j} \times \hat{\mathbf{n}}_{b} + \frac{3}{2}\epsilon_{q}\mathbf{j} \times \mathbf{e} \\
&\qquad -\frac{75}{64}\epsilon_{q}\epsilon_{o}\left\{ \left[ 2(e_{u}j_{n}+e_{n}j_{u})\mathbf{e} + 2(-7e_{u}e_{n}+j_{u}j_{n})\mathbf{j} \right] \right. \\
&\qquad \left. \times \hat{\mathbf{n}}_{b} + \left[ 2e_{n}j_{n}\mathbf{e} + \left(-7e_{n}^2 + j_{n}^2 + \frac{8}{5}e^2 - \frac{1}{5}\right)\mathbf{j} \right] \times \hat{\mathbf{u}}_{b} \right. \\
&\qquad \left. +\frac{16}{5}e_{u}\mathbf{j} \times \mathbf{e} \right\}.
\end{aligned}
\qquad (14) $$

In this perturbative model, the planet's presence is felt through both the axisymmetric (essentially quadrupolar) terms with $\epsilon_q$, and the symmetry-breaking terms (here limited to octupolar order) with $\epsilon_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 $\epsilon_{AB}$.

## 3.1 The eccentric Laplace surface of HD 106906 system: a brief interlude

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:

(i) 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;

(ii) 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 time-scales in mind (Fig. 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 or less significant 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 $j$, $e$, $\hat{n}_b$, $\hat{u}_b$, and $\hat{n}_{AB}$ lie in the same plane. This implies that the nodes of the three players are aligned. The practical reader 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!

Imposing nodal alignment, we recover conditions for stationary debris orbits in the form of two scalar equations in two unknowns: the debris eccentricity $e$, and the debris inclination with respect to the inner plane $i$. For a mutual inclination $i_b$ between the stellar binary's plane and the planet's orbit, conditions reduce to:

$$ 
\begin{aligned}
0 &= \frac{3}{4} \frac{\epsilon_{AB}}{(1 - e^2)^{\frac{5}{2}}} \left[ \frac{1}{2}(1 - e^2) \sin 2i (4e_{AB}^2 + 1) \right] \\
&\quad - \frac{3}{8} \epsilon_q (1 + 4e^2) \sin 2\phi \mp \frac{75}{64} \epsilon_q \epsilon_o e \left[ -\frac{7}{2}(1 + e^2) \cos \phi \sin 2\phi \right. \\
&\quad \left. + (2 + 5e^2) \sin^3 \phi + \frac{1}{5}(1 - 8e^2) \sin \phi \right],
\end{aligned}
\qquad (15) $$

$$ 
\begin{aligned}
0 &= \frac{9}{8} \frac{\epsilon_{AB}e}{(1 - e^2)^2} \left[ -\sin^2 i (4e_{AB}^2 + 1) + e_{AB}^2 + 2/3 \right] \\
&\quad + \frac{3}{4} \epsilon_q e \sqrt{1 - e^2 (1 - 4 \sin^2 \phi)} \\
&\quad \pm \frac{75}{64} \epsilon_q \epsilon_o (1 - e^2)^{\frac{1}{2}} \left[ (3e^2 - 1) \cos^3 \phi \right. \\
&\quad \left. + \left( \frac{1}{5} - \frac{24}{5}e^2 \right) \cos \phi + \left( 1 + \frac{15}{2}e^2 \right) \sin \phi \sin 2\phi \right],
\end{aligned}
\qquad (16) $$

where $\phi = i_b - i$. With nodal alignment, equilibria can be either apsidally aligned or anti-aligned with the outer planet's periapse. We
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The Laplace Surface in HD 106906. Planetesimal rings with angular momentum along $\hat{n}$, are torqued to precess about the inner stellar binary with angular momentum along $\hat{n}_{AB}$, and the outer planet with angular momentum along $\hat{n}_b$. The mutual inclination between the inner and outer perturbers is defined by the angle $i_b$. This misaligned competition promotes precession about the equilibrium Laplace plane. With increasing semimajor axis, equilibrium planes span the warped surface shown the upper panel, with an alignment transition indicated at the Laplace radius $r_L$ (equation 1). The colour codes for planetesimal eccentricity which is pumped as one moves outward on the surface. In the bottom panel, we map planetesimal inclination on the Laplace surface for different values of planetary inclination. The colour renders the corresponding eccentricity as in the upper panel. Curves are obtained by solving equations (15–16) for Laplace equilibria.</figcaption></figure>

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 semimajor axis yields eccentric Laplace surfaces shown in Fig. 3. Surfaces were computed for various inclinations of the outer companion to the 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, semimajor axis, and varying inclination (Table 1). They can be further mapped for varying planet eccentricity and semimajor 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.

## 4 DYNAMICS AROUND HD 106906 LAPLACE SURFACE

Having isolated the Laplace surface controlled by inner binary and outer planet (Fig. 3), our goal now is to study the evolution of a debris disc by numerically integrating the equations of motion (13–14). We perform a suite of numerical simulations to follow the evolution of the orbital elements of a distribution of 24 100 planetesimals situated between the disc’s inner and outer edges, $a_{in} = 30$ au and $a_{out} = 150$ au, respectively. Each planetesimal is parametrized by its secularly constant semimajor axis $a$, and four variable elements: $e, i, \Omega$, and $\omega$. The equations of motion are solved using an explicit Runge–Kutta (4, 5) integrator with adaptive time steps, and constrained by $j^2 + e^2 = 1$ and $j \cdot e = 0$.

We list key ingredients in reported simulations:

(i) as the system is sufficiently hierarchical, we ignore the weak mutual perturbation and the consequently ultra-slow precession of inner stellar binary and outer planet ($\simeq 10^{10}$ yr).<fnref n="3" /> Thus ($\hat{n}_{AB}, \hat{e}_{AB}$) and ($\hat{n}_b, \hat{e}_b$) are fixed by orbital parameters as given in Table 1, with the further prescription of ($\Omega_b, \omega_b$) = (0, 90)° for the planet. The Laplace radius of this fiducial system is located at $r_L = 84$ au;

(ii) 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 finds its way to its current location. Instead, the planet is assumed in its observed position for the duration of reported simulations, and unless otherwise stated, with the mean orbital elements of Table 1;

(iii) 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;

(iv) in this first exploration of an already rich dynamical system, we further ignore the effect of debris disc self-gravity on Laplace surface dynamics, and dedicate space in the discussion to assess its potential impact on the dynamics;

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

(vi) the debris disc is assumed initially axisymmetric, rather cold, and coplanar with the binary’s orbital planar. Planetesimal orbits are initiated with uniformly distributed eccentricities and inclinations, ranging between 0 and 0.05, and 0° and 2°, respectively. We consider and report on secular dynamics over a range of semimajor axes which

<footnote n="3">³In some configurations, that of Planet 9 in particular, the secular evolution of either perturber could be fast enough to be of concern. The matter was discussed at length in Farhat & Touma (2021) who showed that, in case the secular evolution of a perturber is significantly slower than that of disc particles, the effect is rather minor, with equilibria and dynamics around them responding adiabatically to the slow quasi-periodic evolution of the perturber.</footnote>
are of relevance to observations, but nothing prevents extending that
range in either direction while respecting model assumptions; and
(vii) in dialogue with 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 semimajor 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 modelling (Section 4.2).

**4.1 Orbital evolution and surface density**

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

(i) on orbital elements

(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 modelling that ignores the binary's non-
spherically symmetric contributions as we shall discuss below.

(b) Eccentricities: the inclination break around the Laplace
radius is further manifested in the evolution of disc eccentric-
ities. 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.

(c) Nodes and Apses: by 13 Myr, one notes how nodes of
planetary 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 accom-
panied 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.

(ii) 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 Fig. 5 maps of surface density using
the orbital element distribution of planetesimals portrayed in Fig. 4.
The procedure which we follow extends what is presented in Sefilian,
Rafikov & Wyatt (2021) to account for planetesimal inclination [refer
to Appendix B for further details]. By 13 Myr, 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 which reflects eccentricity excitation together with 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 (Ω, *i*)-plane, defined
by ω = 90°, that is, for debris apsidally aligned with the planet.
We restrict plane crossings by imposing the condition of ˙ω > 0,
and we follow the orbits evolving according to equations (13–14),
starting with initial conditions prescribed by an energy hyper-surface.
Each surface of section is parametrized by the secularly conserved
semimajor axis. We present in Fig. 6 samples of sections over
a range of debris semimajor axes which is representative of the
disc studied thus far. We further colour-code for the eccentricity
at section crossings. For each semimajor axis, we construct the
section at an energy corresponding to an initially circular non-
inclined debris orbit.

In Fig. 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
<figure><img src="image_5.png" /></figure>

**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 the inside and the planetary companion on the outside (equations 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 distinct 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. Dynamical regimes translate into morphological asymmetries, warps, and spiral arms visualized in Fig. 5. The evolution is to be compared with that in Fig. 7 when we mute the effect of the stellar binary.

<figure><img src="image_6.png" /></figure>

**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 B using the distributions of orbital elements resulting from the simulations of Fig. 4. The viewing angles are fixed such that: $\theta_i = 85^{\circ}$, that is, the disc is viewed nearly edge-on, $\theta_i = 13.7^{\circ}$ as inferred from a position angle of 283.7$^{\circ}$ (Kalas et al. 2015), and $\theta_a = 0^{\circ}$ (see Appendix B for further details). Each of these panels have been normalized individually, and the green star represents the centre of mass of the inner binary. The disc around its current age (13 Myr) 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 Myr, 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 (Fig. 9), with both scenarios further resolved vertically in the histograms of Fig. 8. An animated version of this morphological evolution is provided within the supplementary material online.

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 can 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 (colour 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 (~5$^{\circ}$).

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$). The libration island by $\Omega = 0$ is now around $i \approx 10^{\circ}$, while those around $\Omega = \pm\pi$ remain at low inclinations. Consequently, debris which is initially in the plane of the inner binary can now reach 30$^{\circ}$. Increasing the semimajor axis further, we enter the outer planet's dominated regime. Phase space is now encroached upon by area-filling chaotic zones
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Poincaré sections for debris particles at different semimajor 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 (Ω, I)-plane for ω = π/2. Crossings are collected when ω > 0. Colour coding corresponds to eccentricity. Of significance to us in these maps are the separatrices 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 semimajor axis are destined to live on the separatrix, following its inclination and, to a lesser degree, eccentricity excursions. We provide in equation (21) a closed form solution to this structure, and at each semimajor axis, we plot in black the associated trajectory crossings.</figcaption></figure>

then high order resonant chains around a prominent island of quasi-periodic motion. Our initially cold disc particles (i.e. on near-circular-coplanar orbits) are seen to circulate on the phase-space trajectory which bounds this regular region, a separatrix of sorts ferrying disc particles with increasing mean inclination as we increase *a*, all the while maintaining low to moderate eccentricity.

From what preceded, it is evident that a one parameter family of phase-space trajectories is critical to the long-term dynamics of debris discs which are initially on near-circular orbits and coplanar with the inner binary. With such initial configurations and their long-term evolution in mind, it would seem reasonable to seek an approximate expression for this family. This is what we proceed to do next as we associate members of this family, at any given semimajor axis *a*, with the 'separatrix' of an intimately related integrable dynamical system.

Tremaine et al. (2009) took a first, and already challenging, step towards analytically identifying Laplace equilibria and nonlinear dynamics around them in the limit of circular orbits under quadrupolar perturbations. We can safely work within this limit given how debris orbits maintain relatively low eccentricity in the course of their evolution. As shown in Appendix B of Tremaine et al. (2009), particle dynamics (now fully described by the evolution of the angular momentum vector *j*) can be mapped into the dynamics of an integrable free-rigid body with a suitably defined inertia tensor. Adopting this methodology, and taking $\hat{n}_{AB} = (1, 0, 0)$, and $\hat{n}_b = (\cos i_b, \sin i_b, 0)$, equation (13) is rewritten as

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

where the symmetric inertia tensor *T* is given by

$$ \mathbf{T} = \begin{pmatrix} \varepsilon_q \cos^2 i_b + \varepsilon_{AB} \varepsilon_q \cos i_b \sin i_b & 0 \\ \varepsilon_q \cos i_b \sin i_b & \varepsilon_q \sin^2 i_b & 0 \\ 0 & 0 & 0 \end{pmatrix}. \qquad (18) $$

Denoting by *R* the orthogonal transformation that diagonalizes *T*, we rotate to the principal-axis coordinate system via $\mathbf{j} = \mathbf{R} \cdot \mathbf{J}$. Under this transformation, the reduced form of the Hamiltonian in

equation (6) is written as (Tremaine et al. 2009):

$$ \hat{H} = -\frac{1}{2} \mathbf{j}^{\mathrm{T}} \cdot \mathbf{T} \cdot \mathbf{j} = -\frac{1}{2} \sum_{k=1}^{3} t_k J_k^2 = -\frac{1}{2} (t_2 J_2^2 + t_3 J_3^2), \qquad (19) $$

with eigenvalues $t_1 = 0 < t_2 < t_3$ for $i_b > 0$. This Hamiltonian generates the (Lie-Poisson) dynamics of the transformed angular momentum **J** via

$$ \frac{d\mathbf{J}}{d\tau} = -\mathbf{J} \times \nabla_j \hat{H}, \qquad (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 surface, an elliptical cylinder. From equation (20), the latter is bounded by the limits $\hat{H}_{\min} = -t_3/2$ and $\hat{H}_{\max} = 0$. Between those limits, at the critical energy $\hat{H}_c = -t_2/2$, the cylinder is tangent to the sphere, and their intersection defines the sought after 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. \qquad (21) $$

We solve for orbits with angular momentum **J** on this separatrix, then we project back to orbits with **j** using the inverse transformation. Solutions yield the phase curves shown in black in Fig. 6. Though obtained under the assumption of circular orbits under dual quadrupolar forcing, those curves are in remarkable agreement with the phase-space trajectory obtained by integrating the full equations of motion. Specifically, the amplitude of inclination excursions is matched between the two methods over the semimajor axis range of interest.

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 neighbouring
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Similar to Fig. 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 Fig. 4, and are subject to periodic excursions in inclination of amplitude $2i_b$.</figcaption></figure>

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 revisit modelling assumptions and explore physical scenarios through which we put the results above in dialogue with previous work, assess implications for disc observables (asymmetries and vertical structure in particular), improve constraints on the planet's orbit, and highlight the broad applicability of Laplace surface dynamics.

### 5.1 Binary off

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

(i) 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_b|\sin(\nu t/2)|$ (Murray & Dermott 1999). Here $\nu$ is the secular oscillation frequency, is dependent on the simulated disc lifetime (Dawson, Murray-Clay & Fabrycky 2011). By 13 Myr, the age of HD 106906, all planetesimals have undergone their first oscillation cycle, before reaching a peak inclination amplitude equal to $2i_b$. By 300 Myr, all planetesimals have completed multiple 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 raise questions for earlier approaches to the modelling of this system;

(ii) 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 Fig. 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 equation (B1). The discrepancy between the two settings emerges very clearly beyond 5 Myr 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 ~500 Myr. The cold-low inclination-population persists at all in all phases of the simulation with a mean around ~5°. In contrast, the hot population has a mean which shift back and forth in time, ranging between 20° and 70°. In the absence of the inner binary's quadrupole, and as reported in the previous section, the whole disc is subject to unhampered inclination perturbations by the planet, and it is only a matter of time before the initially flat disc loses its cold population. As seen in Fig. 8, the density peak around small inclinations smoothly transitions in time towards high inclinations and persists there, leaving no remarkable features on the cold end of the distribution;

(iii) morphologically, when the inner quadrupole is switched off, the warp is much stronger and extends to the inner regions by the age of the system (Fig. 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).

### 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
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Inclination distributions of the disc particles derived from our direct orbital integrations (Fig. 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 B. 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° 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$ Myr to much higher inclinations as time evolves.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Similar to Fig. 5, but muting the effect of the stellar binary. The maps are thus constructed from the orbital elements distributions shown in Fig. 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.</figcaption></figure>

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}_d \ge 0.16$ (Crotts et al. 2021). In contrast with 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 all.<fnref n="4" />

In our study thus far, we took the best-fit planetary parameters listed in Table 1 for granted, ignoring reported uncertainties (Nguyen

et al. 2020). However, within the wide uncertainty intervals in $a_b$ and $e_b$, the Laplace radius of the system can take any value between the reported disc's inner and outer edges. In Fig. 10, we simulate the disc for various values of $r_L$ ranging between 50 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_L$ to 50 au, which would result from $(a_b, e_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 percenter closer in, the time-scale of the secular forcing by the planet also becomes shorter, thus the disc structure observed at 100 Myr in Fig. 5 can now be seen at 13 Myr. Pushing $r_L$ closer to the outer edge of the disc, which can be attained with a planet of $(a_b, e_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

<footnote n="4">⁴While 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.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).</footnote>
<figure><img src="image_11.png" /><figcaption>**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_L$ very close to the inner and outer edges of the disc (see also Fig. 1), thus placing the disc completely under the dominance of either the stellar binary or the planet. With $r_L \rightarrow a_{in}$, we end up with a planet-dominated disc which features warps and asymmetries, and for $r_L \rightarrow a_{out}$ we have a fairly rigid, binary-dominated disc which maintains its initial axisymmetric structure.</figcaption></figure>

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 favouring a rounder and more distant planet.<fnref n="5" />

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 equation (1), which in turn helped constrain the planet's orbit with a broad range of uncertainties. Evidently, we have here a rather powerful criterion which deserves 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 the 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 secular dynamical simulations exploring a range in ($a_b, e_b$). We did so in the hope of constraining the planet's orbital parameters by confronting simulated disc eccentricity with ALMA observations of disc geometry. 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).

Following the first approach, and allowing for variations of planetary semimajor 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}_d$ and contour its surface in Fig. 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}_d \le 0.05$) to null eccentricity. Meanwhile, a maximum disc eccentricity of $\bar{e}_d = 0.276$ can be attained with a planet having ($a_b, e_b$) = (600 au, 0.73). The nominal orbital parameters of the planet from Nguyen et al. (2020) are ($a_b, e_b$) = (850 au, 0.44), corresponding to $\bar{e}_d = 0.07$. Exciting the disc eccentricity to $\bar{e}_d \approx 0.1$ requires the planetary pericenter to be no farther than 415 au. For $\bar{e}_d \ge 0.2$, the study constrains the

planet to a narrow window of the parameter space with a maximum pericenter distance of 240 au. 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, Lamy & Soter 1979; Wyatt 2008). Detailed modelling 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 given by (Burns et al. 1979):

$$a' = \frac{a(1 - e^2)(1 - \beta)}{1 - e^2 - 2\beta(1 + e \cos f)}, \qquad (22)$$

$$e' = \frac{\sqrt{e^2 + 2\beta e \cos f + \beta^2}}{1 - \beta}, \qquad (23)$$

$$\omega' = \omega + \arctan\left(\frac{\beta \sin f}{e + \beta \cos f}\right), \qquad (24)$$

$$i' = i, \qquad (25)$$

$$\Omega' = \Omega, \qquad (26)$$

where $f$ is the true anomaly of the parent planetesimal. Given the orbit averaged nature of our simulations, we retrieve the distribution of planetesimal true anomalies from a uniform distribution of mean anomalies via the Kepler equation.

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, Duchêne & Matthews 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 \propto \beta^{1.19}$. This distribution is only slightly different from that of

<footnote n="5">⁵Fehr et al. (2022) went through a similar exercise and reached a similar conclusion on the planet's orbit based on their ALMA observations of the disc.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Constraining the orbital parameter space of HD 106906 b with disc eccentricity estimates. We contour the surface of disc eccentricity $\bar{e}_d$ computed for a range of values of the planetary semimajor axis $a_b$ and eccentricity $e_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, in principle, be matched with those recovered from ALMA observations.</figcaption></figure>

a collisional cascade, which is characterized by $q \approx 3.5$ (Dohnanyi 1969). Our $\beta$-distribution extends to a maximum value, $\beta_{max}$, which marks a marginally bound orbit,<fnref n="6" /> and is defined for every parent planetesimal as:

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

In Fig. 12, for each prescribed $\beta$ value or distribution, we mark the contour in $(a_b, e_b)$-space that defines the parameter space required to generate a disc eccentricity $\bar{e}_d \ge 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. For $\beta \ge 0.33$, a planet with semimajor 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 \propto \beta^{1.19})$, we are able to constrain $(a_b, e_b)$ to within almost half of the parameter space (the area above the red contour curve in Fig. 12). We find that the nominal values in Nguyen et al. (2020), $(a_b, e_b) = (850 \text{ au}, 0.44)$, are within this constrained area. We finally note that for this nominal orbit, the resulting disc eccentricity is ~0.07 when considering gravitational perturbations alone, and ~0.19 when accounting for dust particles with radiation pressure. Thus, the difference between the reported asymmetry using scattered light (Kalas et al. 2015; Crotts et al. 2021), and the axisymmetry in

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Same as Fig. 11, but allowing for radiation pressure. After spawning dust grains from planetesimal progenitors that evolved within our Laplace dynamical model, we compute the disc’s eccentricity $\bar{e}_d$ for a wide range of planetary semimajor axis and eccentricity covering the observational uncertainty intervals (Table 1). Different curves correspond to different values of $\beta$, the ratio between stellar radiation pressure and gravitational forcing. Prescribed $\beta$-values are randomly selected; while the distribution is sampled from the dust grain size distribution inferred in Crotts et al. (2021). Contours bound from below the area of the orbital parameter space in which the disc’s eccentricity is pumped above 0.16.</figcaption></figure>

ALMA imaging (Kral et al. 2020) may very well be a natural outcome of radiation pressure.

## 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 Fig. 13 the population density for the inclination distribution of the Kuiper belt and $\beta$ Pictoris. We further include 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 B (i.e. equation (B1) with $p = 0.5$). Our simulated HD 106906 disc features a similar distribution, with distinct cold and hot populations. This bimodal distribution is an inherent aspect of Laplace dynamics as elucidated thus far.

<footnote n="6">⁶To be sure, dust grains are likely to encounter HD 106906 b before hitting this limit, which would suggest a smaller upper-bound on $\beta$. This effect is worthy of consideration along with others ignored in this preliminary look at radiation pressure.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Inclination distributions as derived and modelled from observations of the Kuiper belt (Brown 2001), ALMA observations of β 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.</figcaption></figure>

The cold population is predominantly attributed to the inner part of the disc living in the binary-dominated regime; whereas the hot population mainly is mainly the planet-dominated regime. We fitted this distribution to two overlapping Rayleigh distributions with a mean and standard deviation of 1.89° and 1.24° for the cold population, and 28.08° and 18.39° for the broader hot 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, and this occurs when the planet's actual semimajor axis and eccentricity are towards the upper and lower ends of their estimated range, respectively (see Fig. 1 and last panel of Fig. 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°. In the opposite case where the planet overshadows the binary's gravitational field, the vertical distribution is expected in the longterm to show a single peak which is now centered around large values, as seen when we mute the inner binary in Fig. 8. Thus, both the number and positions of peaks in the vertical distribution of planetesimals could help constrain the planet's orbit. Their mapping will have to await future high-resolution ALMA observations, and can then indirectly complement future RV and astrometric measurements of HD 106906 b. (see also Section 5.2, where we discuss $\tilde{e}_d$).<fnref n="7" />

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.1 M_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 B), together with the mass distribution controlled by $p$ (equation B1). 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_{in} \ll a_{out}$. Finally, we emphasize that our results derive from a collisionless model, a limitation that must be assessed before making any definitive conclusions. 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, that is, those with $a \lesssim r_L$, and thus affect the peak of the cold population in Fig. 13. This is so because relative velocities of planetesimals would be higher in the inner parts than in the outer parts (recall that $v_{rel} \propto v_K \propto r^{-1/2}$, where $v_K$ is the Keplerian speed). The extent of depletion further 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

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. Fig. 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 Fig. 14, where we plot the surface density of our simulated HD 106906 disc as viewed on the sky plane for various combinations of $\theta_i$ and $\theta_a$ (see Appendix B). Note that each of the panels in this figure is normalized individually to bring out the fainter features, and $\theta_i$ is kept fixed<fnref n="8" /> at the fiducial value of 13.7°.

Looking at Fig. 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 centre of the disc's inner cavity, as can be seen most prominently for the nearly face-on discs, that is, $\theta_i < 30°$. 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 (Fig. 4, see also Nesvold et al. 2017).

Another common feature in many of the views displayed in Fig. 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°$, starts off within the disc, that is, $a_{in} \ll r \ll a_{out}$, and extends out to the outer disc parts. This can be understood by noting that by this time, that is, $t = 13$ Myr, the stellar binary has homogenized the apsidal angles of planetesimals, together with their longitudes of ascending nodes in the inner disc parts, that is, $a \lesssim r_L$, while those at $r \gtrsim r_L$ have not had

<footnote n="7">⁷The axisymmetric structure of the disc reported in Kral et al. (2020) and Fehr et al. (2022) favours a Laplace dynamical regime in which the stellar binary dominates the full disc, which would in turn favour a dominant cold population in any future sufficiently resolved characterization of the disc's vertical structure.</footnote>

<footnote n="8">⁸Recall that $\theta_i$ measures the angle in the sky-plane from the positive x-axis to the disc's ascending node on the sky (see Appendix B), and thus varying its value will not affect the disc structure.</footnote>
<figure><img src="image_15.png" /><figcaption>**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 800×800 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, that is, with $\theta_i \lesssim 90^\circ$, and a warp for nearly edge-on discs, that is, for $\theta_i \sim 90^\circ$.</figcaption></figure>

completed one complete precession cycle (see Figs 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, Fig. 4 indicates that the pattern may survive for relatively long times, that is, $t \sim 100$ Myr.<fnref n="9" /> 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 coplanar 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 for example Fig. 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 equation (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, that is, 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 semimajor axes in triples of various configurations.

<footnote n="9">⁹This can be also understood by consulting Fig. 2 which shows that the precession time-scale for $a \gtrsim r_L$ has a weak dependence on semimajor axis.</footnote>
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<fnref n="10" /> 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 modelling 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 modelling, 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 coplanar 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 semimajor 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_d = 20 M_⊕$<fnref n="11" /> yields an apse precession rate of ≈ -0.09 Myr⁻¹ (hence a period of ≈ 70 Myr) around 70 au, a rate which, as shown in Fig. 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 Section 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 sure<fnref n="12" />!

We close this section by emphasizing that some of the above variations have already been explored around Solar system 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 emphasizing 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 neighbouring 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 well defined structure gives rise to a bimodal vertical distribution that mimics those 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

<footnote n="10">¹⁰This is as good a spot as any to raise the possibility of a hitherto undetected few Jupiter-mass planet within HD 106906's inner cavity. Such a planet would sculpt the disc's inner edge, as it couples to the binary in dictating debris dynamics by shifting $\tau_L$ further out, among other potential effects. Further consideration of this possibility is best left till when future refined observations warrant it.</footnote>

<footnote n="11">¹¹For Fehr et al. (2022), the observed 1.3 mm flux indicates a ~10$M_⊕$ reservoir of planetesimals, an estimate which surely depends on the assumed size distribution and the maximum size of planetesimals (Krivov & Wyatt 2021).</footnote>

<footnote n="12">¹²The same is true for conclusions reached by Fehr et al. (2022) assuming a disc of test particles.</footnote>
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. A lingering uncertainty involves the planet's orbit, so we mapped the general dependence of the disc's asymmetry on the planet's 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 confront our conclusions with improved constraints on the orbit of HD 106906 b, along with more resolved images of the disc. We then closed the loop on our dynamical study by probing the phase space structure of the disc. We learned that initially cold discs are destined to evolve on a robust phase space structure, specifically a separatrix which carries planetesimals on excursions of inclination and eccentricity which grow in amplitude with increasing semimajor axis. We recovered a closed form solution for this separatrix, and noted how this structure and neighbouring 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 sculped 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 & Touma (in preparation) [drawing on Sefilian (2017)]. In the process we learned: 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, emphasize 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 coplanar with the binary. This is so because, as it happens, a million year gas dissipation time-scale is at worst comparable, at best faster than secular time-scales 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 help 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. JRT 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. The authors are grateful to the referee, Hervé Beust, for insightful comments which helped improve the manuscript considerably. AAS is supported by the Alexander von Humboldt Foundation through a Humboldt Research Fellowship for postdoctoral researchers. AAS 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 R., Cook G., 1964, Proc. Royal Society of London. Series A. Mathemat. Phys. Sci., 280, 97

Antognini J. M., 2015, MNRAS, 452, 3610

Bailey V., et al., 2013, ApJ, 780, L4

Bate M. R., Bonnell I. A., Bromm V., 2002, MNRAS, 332, 65

Batygin K., Morbidelli A., 2017, AJ, 154, 229

Batygin K., Adams F. C., Brown M. E., Becker J. C., 2019, Physics Reports, 805, 1

Beust H., 2016, A&A, 590, 2

Brown M. E., 2001, AJ, 121, 2804

Bryan M. L., Chiang E., Morley C. V., Mace G. N., Bowler B. P., 2021, AJ, 162, 217

Burns J. A., Lamy P. L., Soter S., 1979, Icarus, 40, 1

Carter A. L., et al., 2021, MNRAS, 501, 1999

Charnoz S., Canup R. M., Crida A., Dones L., 2018, The Origin of Planetary Ring Systems. Cambridge Univ. Press, Cambridge, p. 517

Chiang E., Kite E., Kalas P., Graham J., Clampin M., 2009, ApJ, 693, L734

Correia A. C., Laskar J., Farago F., Boué G., 2011, Celest. Mech. Dynam. Astron., 111, 105

Correia A., Boué G., Laskar J., 2016, Celest. Mech. Dynam. Astron., 126, 189

Crotts K. A., et al., 2021, ApJ, 915, L58

Ćuk M., Hamilton D. P., Lock S. J., Stewart S. T., 2016, Nature, 539, 402

Daemgen S., Todorov K., Quanz S. P., Meyer M. R., Mordasini C., Marleau G.-D., Fortney J. J., 2017, A&A, 608, 71
Dawson R. I., Murray-Clay R. A., Fabrycky D. C., 2011, ApJ, 743, L17

De Rosa R. J., Kalas P., 2019, AJ, 157, 125

Dohnanyi J., 1969, Journal of Geophysical Research, 74, 2531

Dong R., Zhu Z., Fung J., Rafikov R., Chiang E., Wagner K., 2015, ApJ, 816, L12

Elliot J., et al., 2005, AJ, 129, 1117

Farago F., Laskar J., 2010, MNRAS, 401, 1189

Farhat M. A., Touma J. R., 2021, MNRAS, 507, 6078

Fehr A. J., et al., 2022, ApJ, 939, L56

Gomes R. S., 2003, Icarus, 161, 404

Hamers A. S., 2020, MNRAS, 494, 5492

Hughes A. M., Duchêne G., Matthews B. C., 2018, ARA&A, 56, 541

Jennings R. M., Chiang E., 2021, MNRAS, 507, 5187

Jílková L., Zwart S. P., 2015, MNRAS, 451, 804

Kalas P. G., et al., 2015, ApJ, 814, L32

Kavelaars J., et al., 2009, AJ, 137, 4917

Kral Q., Matra L., Kennedy G. M., Marino S., Wyatt M. C., 2020, MNRAS, 497, 2811

Krivov A. V., Wyatt M. C., 2021, MNRAS, 500, 718

Lagrange A. M., et al., 2016, A&A, 586, 8

Laplace P.-S., 1805, Mecanique Celeste, Vol. 4. Courcier, Paris

Lee E. J., Chiang E., 2016, ApJ, 827, L125

Levison H. F., Morbidelli A., 2003, Nature, 426, 419

Matrà L., Wyatt M. C., Wilner D. J., Dent W. R., Marino S., Kennedy G. M., Milli J., 2019, AJ, 157, 135

Maury A., et al., 2019, A&A, 621, A76

Milankovitch M., 1939, Bull. Acad. Sci. Math. Nat. A, 6, 1

Moore N. W. H., Li G., Hassenzahl L., Nesvold E. R., Naoz S., Adams F. C., 2023, ApJ, 943, 6

Murray C. D., Dermott S. F., 1999, Solar System Dynamics. Cambridge Univ. press, Cambridge

Naoz S., 2016, ARA&A, 54, 441

Nealon R., Cuello N., Gonzalez J.-F., van der Plas G., Pinte C., Alexander R., Ménard F., Price D. J., 2020, MNRAS, 499, 3857

Nesvold E. R., Naoz S., Vican L., Farr W. M., 2016, ApJ, 826, L19

Nesvold E. R., Naoz S., Fitzgerald M. P., 2017, ApJ, 837, L6

Nesvorný D., 2015, AJ, 150, 68

Nguyen M. M., Rosa R. J. D., Kalas P., 2020, AJ, 161, 22

Olofsson J., et al., 2022, MNRAS, 513, 713

Pearce T. D., Wyatt M. C., 2014, MNRAS, 443, 2541

Pecaut M. J., Mamajek E. E., Bubar E. J., 2012, ApJ, 746, L154

Petit J.-M., et al., 2011, AJ, 142, 131

Rafikov R. R., Silsbee K., 2014, ApJ, 798, L70

Rodet L., Beust H., Bonnefoy M., Lagrange A. M., Galli P. A. B., Ducourant C., Teixeira R., 2017, A&A, 602, 12

Rosengren A. J., Scheeres D. J., 2014, ApJ, 786, L45

Saillenfest M., Fouchard M., Tommei G., Valsecchi G. B., 2017, Celestial Mechanics and Dynamical Astronomy, 129, 329

Saillenfest M., Lari G., Boué G., Courtot A., 2021, A&A, 647, 92

Sefilian A. A., 2017, MSc Thesis, American University Of Beirut

Sefilian A. A., Touma J. R., 2019, AJ, 157, 59

Sefilian A. A., Rafikov R. R., Wyatt M. C., 2021, ApJ, 910, L13

Spalding C., Millholland S. C., 2020, AJ, 160, 105

Speedie J., Zanazzi J., 2020, MNRAS, 497, 1870

Statler T. S., 2001, AJ, 122, 2257

Tamayo D., Burns J. A., Hamilton D. P., Nicholson P. D., 2013, AJ, 145, 54

Tian Z., Wisdom J., 2020, Proc. National Academy of Sciences, 117, 15460

Toomre A., 1983, in Symposium-International Astronomical Union: Internal Kinematics and Dynamics of Galaxies. Cambridge Univ. Press, Cambridge, p. 177

Touma J. R., Sridhar S., 2015, Nature, 524, 439

Touma J., Wisdom J., 1998, AJ, 115, 1653

Tremaine S., Touma J., Namouni F., 2009, AJ, 137, 3706

van Holstein R. G., et al., 2021, A&A, 647, A21

Van Der Plas G., et al., 2019, A&A, 624, A33

Wagner K., Apai D., Kasper M., Robberto M., 2015, ApJ, 813, L2

Ward W. R., 1981, Icarus, 46, 97

Wyatt M. C., 2008, ARA&A, 46, 339

Wyatt M., Dermott S., Telesco C., Fisher R., Grogan K., Holmes E., Pina R., 1999, ApJ, 527, L918

Zanazzi J., Lai D., 2016, MNRAS, 464, 3945

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

### SURF Den Maps EVO 500Myrs_BIN_OFF_UP.mp4

### SURF Den Maps EVO 500Myrs_BIN_ON_UP.mp4

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

## APPENDIX A: THE SECULAR HAMILTONIAN

We are interested in the motion of a mass-less planetesimal which evolves in the combined gravitational field of a central stellar binary and an outer massive planetary companion. We provide here details on the derivation of the orbit-averaged Hamiltonian of equation (6) which drives the dynamical evolution of this planetesimal via equations [13, 14]. Following notation in the body of the text, the central binary components have masses $m_A$ and $m_B$ and total mass $m_{AB} = m_A + m_B$. We denote their instantaneous barycentric position vectors by $\mathbf{r}_A$ and $\mathbf{r}_B$ with $\mathbf{r}_{AB} = \mathbf{r}_A - \mathbf{r}_B$. We further denote by $\mathbf{r}$ the instantaneous position vector of the planetesimal relative to the barycenter of the binary. Similarly, we prescribe for the outer planet the mass $m_b$ and the position vector $\mathbf{r}_b$. With these definitions, debris dynamics is governed by the Hamiltonian

$$H = \frac{1}{2} \dot{\mathbf{r}}^2 - G \left( \frac{m_A}{|\mathbf{r} - \mathbf{r}_A|} + \frac{m_B}{|\mathbf{r} - \mathbf{r}_B|} + \frac{m_b}{|\mathbf{r} - \mathbf{r}_b|} \right). \quad (A1)$$

Denoting by $\mu$ the mass ratio $m_B/m_{AB}$, the position vectors become $\mathbf{r}_A = \mu \mathbf{r}_{AB}$ and $\mathbf{r}_B = (\mu - 1) \mathbf{r}_{AB}$, allowing us to rewrite the Hamiltonian as

$$H = \frac{1}{2} \dot{\mathbf{r}}^2 - G \left( \frac{m_A}{|\mathbf{r} - \mu \mathbf{r}_{AB}|} + \frac{m_B}{|\mathbf{r} - (\mu - 1) \mathbf{r}_{AB}|} + \frac{m_b}{|\mathbf{r} - \mathbf{r}_b|} \right). \quad (A2)$$

One then expands the potential in terms of distance ratios, $r_{AB}/r$ and $r/r_b$, respectively, and the usual Legendre polynomials, $P_\ell$. For instance:

$$\frac{1}{|\mathbf{r} - \mathbf{r}_b|} = \frac{1}{r_b} \sum_{\ell=0}^{\infty} \left( \frac{r}{r_b} \right)^{\ell} P_{\ell}(\cos \gamma), \quad (A3)$$

where $\cos \gamma = (\mathbf{r} \cdot \mathbf{r}_b)/(rr_b)$. Such expansions are reliable provided the system's architecture is sufficiently hierarchical for the series to converge. Given the nominal orbital architecture of the studied system, HD 106906, described in Table 1, we truncate the expansions at the octupolar order, that is, at order 3 in $r_{AB}/r$ and $r/r_b$, respectively. The Hamiltonian then reads

$$H = H_K + \Phi_{AB, \text{quad}} + \Phi_{AB, \text{oct}} + \Phi_{b, \text{quad}} + \Phi_{b, \text{oct}}, \quad (A4)$$

where we have defined

$$H_K = \frac{1}{2} \dot{\mathbf{r}}^2 - G \frac{m_{AB}}{r}, \quad (A5)$$
$$ \Phi_{AB, \text{quad}} = -\frac{G}{2r^3} \frac{m_A m_B}{m_{AB}} \left[ \frac{3(\mathbf{r} \cdot \mathbf{r}_{AB})^2}{r^2} - r_{AB}^2 \right], \quad (A6) $$

$$ \Phi_{AB, \text{oct}} = -\frac{G}{2r^5} \frac{m_A m_B}{m_{AB}^2} (m_A - m_B) \left[ \frac{5(\mathbf{r} \cdot \mathbf{r}_{AB})^3}{r^2} - 3(\mathbf{r} \cdot \mathbf{r}_{AB}) r_{AB}^2 \right] \quad (A7) $$

$$ \Phi_{b, \text{quad}} = -\frac{G m_b}{2 r_b^3} \left[ \frac{3(\mathbf{r} \cdot \mathbf{r}_b)^2}{r_b^2} - r^2 \right], \quad (A8) $$

$$ \Phi_{b, \text{oct}} = -\frac{G m_b}{2 r_b^5} \left[ \frac{5(\mathbf{r} \cdot \mathbf{r}_b)^3}{r_b^2} - 3(\mathbf{r} \cdot \mathbf{r}_b) r^2 \right]. \quad (A9) $$

It is evident from equation (A7) that the octupolar contribution of the inner binary vanishes when the binary components have equal masses, which is almost the case of HD 106906. We therefore ignore this contribution in our study. To proceed further, we rely on the reference frames introduced in the main text and defined by the triads $(\hat{\mathbf{n}}, \hat{\mathbf{u}}, \hat{\mathbf{v}})$, $(\hat{\mathbf{n}}_{AB}, \hat{\mathbf{u}}_{AB}, \hat{\mathbf{v}}_{AB})$, and $(\hat{\mathbf{n}}_b, \hat{\mathbf{u}}_b, \hat{\mathbf{v}}_b)$. Bound planetesimal motion is then given by

$$ r = \frac{a(1 - e^2)}{1 + e \cos(f)}, \quad \vec{r} = r(\cos(f)\hat{\mathbf{u}} + \sin(f)\hat{\mathbf{v}}), \quad (A10) $$

where a and e are the planetesimal's semimajor axis and eccentricity, and v is its true anomaly relative to $\hat{\mathbf{u}}$.

Secular motion is then captured by averaging Hamiltonian contributions in equations (A5–A9), first over P, the orbital period of the planetesimal on its osculating Keplerian orbit, then over the orbital period of the outer planet. In doing so, it is useful to keep in mind the following differential relations

$$ \frac{dt}{P} = \frac{rdE}{2\pi a} = \frac{r^2 df}{2\pi ab} = \frac{dM}{2\pi}, \quad (A11) $$

where $b = a\sqrt{1 - e^2}$ is the semiminor axis, E the eccentric anomaly, and M the mean anomaly, with analogous expressions for the outer planet. In performing those averages, we rely on the evident separation of time-scales between the orbital periods of the inner stellar binary, the planetesimal, and the planetary companion, while ignoring possible, higher order, mean-motion resonances. We then recover expressions for the orbit-averaged multipoles, parametrized by the normalized angular momentum vector, $j = \sqrt{1 - e^2}\hat{\mathbf{n}}$ and the Lenz vector $e = e\hat{\mathbf{u}}$, via a straightforward though somewhat laborious exercise which is now well documented in various publications on hierarchical triples (e.g. Tremaine et al. 2009; Correia et al. 2011; Hamers 2020; Farhat & Touma 2021).

## APPENDIX B: 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 generalized to account for orbital inclinations and viewing angles.

Maps are built with $N = 2410$ planetesimals drawn from a uniform distribution in the semimajor axis $a$ from the 24 100 planetesimals constituting our simulated debris disc. We assign a mass $m_i$, $i = 1, \dots, N$, to each planetesimal in the sample in such a way that the initial disc surface density $\Sigma_d^{t=0}$ varies with radius $r$ – or equivalently with $a$, since planetesimals are initiated on circular and coplanar orbits (Section 4) – following a power-law profile given by:

$$ \Sigma_d^{t=0}(a) = \Sigma_0 \left(\frac{a_{\text{out}}}{a}\right)^p, \quad (B1) $$

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

$$ dm(a) = 2\pi a \Sigma_d^{t=0}(a) da, \quad (B2) $$

and we compute the amount of mass per unit semimajor axis for a given initial surface density (equation B1).

After evolving the planetesimal distribution in time, at a given snapshot of our simulations, we smear the mass of each planetesimal over its orbit, that is, into an eccentric ring. This is done by spawning $N_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_p$, and mean anomalies $M$ that are randomly selected from the $[0, 2\pi]$ range. This procedure guarantees that each ring is characterized 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_t = N \times N_p$, using Kepler's equation:

$$ M = E - e \sin E. \quad (B3) $$

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

$$ \begin{aligned} x &= r(\cos \Omega \cos(\omega + f) - \sin \Omega \sin(\omega + f) \cos i), \\ y &= r(\sin \Omega \cos(\omega + f) + \cos \Omega \sin(\omega + f) \cos i), \\ z &= r \sin(\omega + f) \sin i, \end{aligned} \quad (B4) $$

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_1) \cdot \mathcal{R}_x(\theta_1) \cdot \mathcal{R}_z(\theta_a) \begin{pmatrix} x \\ y \\ z \end{pmatrix}, \quad (B5) $$

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}, \quad \mathcal{R}_z(\theta) = \begin{pmatrix} \cos \theta & -\sin \theta & 0 \\ \sin \theta & \cos \theta & 0 \\ 0 & 0 & 1 \end{pmatrix}. \quad (B6) $$

In equation (B5), the angle $\theta_1$ is the angle in the sky plane, that is, 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 P.A. = 283.7° (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_1 = \text{P.A.} + \pi/2$ (i.e. anti-clockwise from North).<fnref n="13" /> The angle $\theta_1$, on the other hand, is a measure of the angle between the normals to the sky- and

<footnote n="13">¹³A similar result could be obtained using the value of PA = -74.7° of Olofsson et al. (2022), in which case $\theta_1 = -74.7 + 90 = 15.3°$.</footnote>
disc-plane, with $\theta_i = 0^{\circ}(90^{\circ})$ indicating a face-on (edge-on) disc (default =85$^{\circ}$ in text, except stated otherwise). Finally, the angle $\theta_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_t = N \times N_p$ particles in the sky plane with $[X, Y] = [-200, 200]$ au, with a given resolution of 800 × 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 2D Gaussian having a filter size of ~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 centre of mass of the stellar binary and (ii) the logarithmic colour scale represents the surface density distribution normalized by its maximum value. In most of the presented density maps, normalization was performed individually to bring out the fainter features.

This paper has been typeset from a TEX/LATEX file prepared by the author.