The vertical structure of debris discs and the role of disc gravity: A primer using a simplified model

TL;DR

This work develops an analytical framework to study how the gravity of a debris disk back-reacts on an interior inclined planet and how this back-reaction shapes the disk's vertical structure. Using axisymmetric Laplace–Lagrange secular theory, the authors identify two regimes—disk-dominated and planet-dominated—separated by a secular-inclination resonance, with distinct warp and scale-height behaviors: a distance-independent $\, ext{H}(R) o I_p(0)$ in the planet-dominated case, and a steep $\, ext{H}(R)\propto R^{-7/2}$ decline in the disk-dominated case. The model predicts non-Gaussian vertical density profiles and bimodal inclination distributions near resonances, and provides analytic formulas to infer planetary parameters and disk masses from observed warps and aspect-ratio trends, as demonstrated for HD 110058 and $eta$ Pictoris. The results underscore the importance of including disk gravity when interpreting debris-disk structures and offer a practical framework for guiding future, more comprehensive investigations and observations.

Abstract

Debris discs provide valuable insights into the formation and evolution of exoplanetary systems. Their structures are commonly attributed to planetary perturbations, serving as probes of as-yet-undetected planets. However, most studies of planet-debris disc interactions ignore the disc's gravity, treating it as a collection of massless planetesimals. Here, using an analytical model, we investigate how the vertical structure of a back-reacting debris disc responds to secular perturbations from an inner, inclined planet. Considering the disc's axisymmetric potential, we identify two dynamical regimes: planet-dominated and disc-dominated, which may coexist, separated by a secular-inclination resonance. In the planet-dominated regime ($M_d/m_p\ll1$), we recover the classical result: a transient warp propagates outward until the disc settles into a box-like structure centered around the planetary orbit's initial inclination $I_p(0)$, with a distance-independent aspect ratio $\mathcal{H}(R)\approx I_p(0)$. In contrast, in the disc-dominated regime ($M_d/m_p\gtrsim1$), the disc exhibits dynamical rigidity, remaining thin and misaligned, with significantly suppressed inclinations and a sharply declining aspect ratio, $\mathcal{H}(R)\propto I_p(0)R^{-7/2}$. In the intermediate regime ($M_d/m_p\lesssim1$), the system exhibits a secular-inclination resonance, leading to long-lived, warp-like structures and a bimodal inclination distribution, containing both dynamically hot and cold populations. We provide analytic formulae describing these effects as a function of system parameters. We also find that the vertical density profile is intrinsically non-Gaussian and recommend fitting observations with non-zero slopes of $\mathcal{H}(R)$. Our results may be used to infer planetary parameters and debris disc masses based on observed warps and scale heights, as demonstrated for HD 110058 and $β$ Pic.

Figures (13)

• Figure 1: An illustration depicting an edge-on view of the initial setup of the model system examined in this work. The system consists of a central star of mass $M_c$ which is orbited by a planet of mass $m_p$ and an external debris disk of mass $M_d$ (with $m_p, M_d \ll M_c$). The planetary orbit, depicted by the dashed green line, is taken to be circular with a semimajor axis $a_p$ and is slightly inclined with respect to the disk plane (typically, $I_p(0) \lesssim 20^{\circ}$). The debris disk, extending from $a_{\rm in}$ to $a_{\rm out}$, is circular and initially vertically razor-thin, lying in the equatorial plane of the star. The disk's inner edge is located at a minimum distance of $\Delta a_p$ away from the planetary semimajor axis (Equation (\ref{['eq:MMR_condition_delta_ap']})), and its initial surface density distribution is given by a truncated power-law profile (Equation (\ref{['eq:Sigma_d']})). The disk is represented as a collection of $N\gg 1$ massive planetesimals, each with a fixed, logarithmically-spaced semimajor axis $a_j$. As the system evolves, the inclinations ($I_j$) and longitudes of the ascending nodes ($\Omega_j$) of each planetesimal ($j = 1, \ldots, N$) and that of the planet ($I_p$ and $\Omega_p$) may vary due to gravitational interactions. For further details on the model system and the theoretical framework describing its long-term evolution, see Sections \ref{['sec:modelsystem']} and \ref{['sec:analytical-model']}.
• Figure 2: The ratio $A_{d,p}/A_p$, i.e., the disk-induced planetary nodal recession rate ($A_{d,p}<0$; Eq. (\ref{['eq:Adp']})) relative to that of the planetesimals due to the planet ($A_p(a_d)<0$; Eq. (\ref{['eq:A_planet']})), as a function of planetesimal semimajor axis $a_d/a_{\rm in}$. Calculations are done for various disk-to-planet mass ratios, $M_d/m_p$ (shown by different colors), assuming the fiducial disk model ($p=3/2$, $a_{\rm in} = 30$ au, $\delta =5$; Section \ref{['sec:modelsystem']}) and a planet at $a_p = 20$ au. For $10^{-2} \lesssim M_d/m_p \lesssim 10$, the nodal recession rates are equal, $A_{d,p}/A_{p} = 1$, at a unique semimajor axis $a_{\rm res}$. Interior to this, $|A_{d,p}| \lesssim |A_p|$, while exterior to it, $|A_{d,p}| \gtrsim |A_p|$. For more massive disks, i.e., $M_d/m_p \gtrsim 10$, the ratio $A_{d,p}/A_p \gtrsim 1$ across the entire disk. Conversely, for low-mass disks, i.e., $M_d/m_p \lesssim 10^{-2}$, the ratio $A_{d,p}/A_p \lesssim 1$ at all semimajor axes. Note that $A_{d,p}/A_p \propto a_d^{7/2}$ for $a_d/a_{\rm in} \gg 1$ (black dotted line; Eq. \ref{['eq:Adp_over_Ap']}). See the text (Sections \ref{['sec:dist_fn_d_p']} and \ref{['sec:dyn_reg']}) for more details.
• Figure 3: Profiles of the time-averaged planetesimal inclinations $\langle I_{d}(t) \rangle_t$ as a function of their semimajor axes $a_d/a_{\rm in}$, computed using Equation (\ref{['eq:I_forced_gen']}) for different values of disk-to-planet mass ratios $M_d/m_p$. These profiles scale linearly with the initial planetary inclination, which here we have taken to be $I_p(0) =10^{\circ}$. All other system parameters are identical to those in Figure \ref{['fig:A_a']}. For reference, the black dashed line shows the results for a massless disk, i.e., $\langle I_{d}(t) \rangle_t = I_p(0)$ (see also Eq. (\ref{['eq:Im_no_disk']})). In contrast, the black dotted line illustrates the asymptotic behavior of the inclinations when suppressed by the disk gravity, i.e., $\langle I_{d}(t) \rangle_t \propto a_d^{-7/2}$ (see also Eq. (\ref{['eq:Im_disk_sup_numest']})). Note the occurrence of a secular inclination resonance for $10^{-2} \lesssim M_d/m_p \lesssim 10$, where $\langle I_{d}(t) \rangle_t \rightarrow \infty$ (see also Eq. (\ref{['eq:lambda_If_gen']})). See the text (Section \ref{['sec:analytical_analysis']}) for details.
• Figure 4: The locations of secular-inclination resonances relative to the disk inner edge, $a_{\rm res}/a_{\rm in}$, as functions of planetary semimajor axes $a_p/a_{\rm in}$ and disk-to-planet mass ratios $M_d/m_p$. Results are obtained by numerically solving the resonance condition (Equation \ref{['eq:res_condition_gen_exp']}) for a disk model with $p=3/2$ and $\delta = a_{\rm out}/a_{\rm in} = 5$ (Section \ref{['sec:modelsystem']}). The white solid curves represent constant $M_d/m_p$ contours in $\log_{10}$ units (see the color bar). The black dashed line shows the scaling of $a_{\rm res}$ with $a_p$ for fixed $M_d/m_p$ (Equation \ref{['eq:res_loc']}). Note the robustness of establishing a resonance across a wide range of $a_p/a_{\text{inn}}$ and $M_d/m_p$. When no resonances occur within the disk, planetesimal dynamics across the entire disk will be dominated either by the planet's potential (i.e., $M_d/m_p \lesssim 10^{-3}$) or by the disk's potential (i.e., $M_d/m_p \gtrsim 1$) alone. See the text (Section \ref{['sec:dyn_reg']}) for details.
• Figure 5: The radial profiles of planetesimal inclinations $I_d$ and longitudes of ascending node $\Delta\Omega = \Omega_d - \Omega_p$ (top and bottom rows in each panel, respectively) at different snapshot times (left to right). Panels (a), (b), and (c) correspond to different planet--debris disk systems (Table \ref{['table:models']}), as indicated in subcaptions. Snapshot times are provided in absolute terms $t$ and relative to relevant secular timescales (Section \ref{['subsec:orbit-evol-massive-less-disks']}). Each model includes $N = 5,000$ planetesimals initialized on coplanar orbits ($I_d(0) = 0$) and a planet with $I_p(0) = 10^{\circ}$ and $\Omega_p(0) =0$. In each panel, the red curves show the maximum inclinations, $I_m = 2 | I_{d, {\rm forced}} |$ (Eq. \ref{['eq:I_forced_gen']}). For reference, the dashed black lines in panels (b) and (c) show the maximum inclinations in the absence of disk gravity, i.e., $I_m^{{\rm n/disk}} = 2 I_p(0)$ (Eq. \ref{['eq:Im_no_disk']}). In panel (a), the vertical black line (when present) marks the warp location $a_{\rm warp}$ (Eq. \ref{['eq:awarp_explicit']}). In panel (b), the vertical green line marks the secular inclination resonance location $a_{\rm res}$ (Eq. \ref{['eq:res_loc']}). Note the different y-axis scales for $I_d$ and $\Delta\Omega$ in each of panels (a), (b), and (c). See the text (Section \ref{['sec:analytical_analysis']}) for details.