Polluting white dwarfs with Oort cloud comets

ABSTRACT

Observations point to old white dwarfs (WDs) accreting metals at a relatively constant rate over 8 Gyr. Exo-Oort clouds around WDs have been proposed as potential reservoirs of materials, with galactic tide as a mechanism to deliver distant comets to the WD’s Roche limit. In this work, we characterize the dynamics of comets around a WD with a companion having semimajor axes on the orders of 10–100 au. We develop simulation techniques capable of integrating a large number (10⁸) of objects over a 1 Gyr time-scale. Our simulations include galactic tide and are capable of resolving close interactions with a massive companion. Through simulations, we study the accretion rate of exo-Oort cloud comets into a WD’s Roche limit. We also characterize the dynamics of precession and scattering induced on a comet by a massive companion. We find that (i) WD pollution by an exo-Oort cloud can be sustained over a Gyr time-scale, (ii) an exo-Oort cloud with structure like our own Solar system’s is capable of delivering materials into an isolated WD with pollution rate ∼10⁸ g s⁻¹, (iii) adding a planetary-mass companion reduces the pollution rate to ∼10⁷ g s⁻¹, and (iv) if the companion is stellar mass, with M_p ≳ 0.1 M_⊙, the pollution rate reduces to ∼3 × 10⁵ g s⁻¹ due to a combination of precession induced on a comet by the companion, a strong scattering barrier, and low likelihood of direct collisions of comets with the companion.

comets: general, Oort Cloud, white dwarfs, planets and satellites: dynamical evolution and stability

1 INTRODUCTION

Between 25 per cent to 50 per cent of white dwarf (WD) atmospheres observed are polluted with heavy metals (e.g. Zuckerman et al. 2003, 2010; Koester, Gänsicke & Farihi 2014; Wilson et al. 2019). It is generally believed that the heavy metals of polluted WDs come from its evolved planetary system, such as exomoons, exoplanets or exo-asteroid belts, through a variety of mechanisms which induce destabilization in these sources (e.g. Debes & Sigurdsson 2002; Debes, Walsh & Stark 2012; Mustill, Veras & Villaver 2014; Smallwood et al. 2018; Maldonado et al. 2020; O’Connor, Teyssandier & Lai 2022; Trierweiler et al. 2022). The bodies must be delivered to the WD Roche radius at ∼1 R_⊙ to be tidally disrupted for pollution (e.g. Veras & Heng 2020).

Observations of old polluted WDs (WDs with cooling age older than 1 Gyr) found accretion rates ranging five orders of magnitude between 5 × 10⁵ to 10¹⁰ g s⁻¹; the mean pollution rate is around 10⁷ g s⁻¹ with about a 1 dex spread (e.g. Blouin & Xu 2022; Cunningham et al. 2022; Johnson et al. 2022). The current minimum detection limit is about a few 10⁵ g s⁻¹ (Koester et al. 2014). Blouin & Xu (2022) also find that the WD accretion rate decreases by no more than one order of magnitude over 8 Gyr. Thus, their observational findings require a reservoir and mechanism that can deliver materials in the WD’s Roche radius over such long time-scales.

The chemical composition of accreted materials can be analysed to study the original reservoir of the accreted bodies. Until now, a few dozen WDs have been followed up spectroscopically to measure their element abundances (Jura & Young 2014). The majority of polluted WDs exhibit compositions resembling the bulk Earth. This suggests that the sources of materials polluted WDs must be rocky in nature (Jura 2006; Jura & Xu 2012; Xu et al. 2017; Doyle et al. 2019; Trierweiler, Doyle & Young 2023). There are observations of polluted WD atmospheres with volatiles, although they are much more rare (Farihi, Gänsicke & Koester 2013; Doyle, Desch & Young 2021; Klein et al. 2021). At the present, there is one observation of a polluted WD with an icy body composition (Xu et al. 2017). Because the majority of observations point to rocky materials as a source, the Oort cloud has often been ruled out as a reservoir for WD pollution.

The Oort cloud is a byproduct of Solar system formation, where planetesimals are either ejected or kicked into high semimajor axes via interactions with a giant planet (Vokrouhlický, Nesvorný & Dones 2019; Kaib & Volk 2022). For objects kicked into higher semimajor axes, they can be circularized by galactic tide or stellar flybys (Duncan, Quinn & Tremaine 1987; Hahn & Malhotra 1999; Higuchi & Kokubo 2015; Vokrouhlický et al. 2019). These objects form the Solar system’s Oort cloud, which is a structure containing 10¹¹–10¹² objects with a total mass of ∼2 M_⊕ with semimajor axes ranging between 3000–100 000 au (e.g. Weissman 1983; Boe et al. 2019).

Since the Solar system’s Oort cloud is a byproduct of interactions of planetesimals with planets, it is not unreasonable to expect Oort clouds to exist around other main-sequence and even WD planetary systems. As a result, several studies have investigated if the Oort cloud is a suitable reservoir for pollution. They often employ a mix of numerical and analytic methods to study the pollution rate of Oort cloud comets. Mechanisms considered consist of galactic tidal effects, stellar flybys, WD natal kicks and stellar mass loss during post-main-sequence evolution (Alcock, Fristrom & Siegelman 1986; Parriott & Alcock 1998; Veras, Shannon & Gänsicke 2014; Stone, Metzger & Loeb 2015; Caiazzo & Heyl 2017). Simulations of comets into a 1 R_⊙ tidal radius in existing literature we are aware of do not resolve pollution rate over time due to being resolution limited in the number of simulated comets. Most recently, O’Connor, Lai & Seligman (2023) perform an extensive analytic and numerical studies on the effects of galactic tide, planetary perturbations, WD natal kicks, and stellar mass loss on the pollution rate of WDs with exo-Oort clouds as a source of materials. They find that combining these effects together point to an accretion rate of a few 10⁵ g s⁻¹, which is just at or slightly below the detection limit. They argue that this is potentially why we do not observe many volatiles-rich polluted WD atmospheres.

In this article, we contribute to this existing line of investigation by answering the following question through numerical simulations: can an Oort cloud similar to the one we have in the present-day Solar system pollute a WD over a Gyr time-scale? We also investigate this question in cases where the WD has a planetary-mass or a stellar-mass companion. Since we currently have no knowledge of extrasolar Oort clouds, we simply assume throughout this article that other Oort clouds have the total number of objects and total mass like the current Solar system Oort cloud. We also study how various Oort cloud power-law density profiles affect pollution. All of our results can be easily scaled to another Oort cloud with different mass, number of objects, and density profile.

We start by presenting the analytic theory of WD pollution through galactic tide and WD companion (planet or star). Specifically, in Section 2, we summarize the analytic vertical galactic tide model (Heisler & Tremaine 1986), the loss cone theory, leading to the expected analytic pollution rate as predicted by O’Connor et al. (2023). In Section 3, we study additional dynamics comets experience when there is a companion. Specifically, we analyse the dynamics of galactic tide together with precession and scattering induced by a companion. We also summarize the loss cone shielding model proposed by O’Connor et al. (2023) leading to a prediction of the WD pollution rate in the presence of a companion in that framework.

In Section 4, we describe our simulation methodology which is capable of integrating a large number of comets (10⁸ comets) over a long time (1 Gyr), initial conditions, and boundary conditions.

In Section 5, we present pollution rate over various Oort cloud structures, in the presence of galactic tide only. Then, we present the pollution rate in the presence of a stellar companion, and in the presence of a planetary companion. Finally, we show the pollution rate over a 1 Gyr time-scale. In Section 5, we also compare and discuss our results with analytic expectations from O’Connor et al. (2023).

In Section 6, we discuss advantages and major concerns of using an Oort cloud as a potential reservoir for WD pollution, such as if an Oort cloud can survive post-main-sequence evolution and that the majority of observed WDs are volatiles poor. We also discuss our results in contexts of observations of close-in binaries, wide binaries. In Section 7, we summarize our findings.

1.1 Notations

In this work, we denote M_* as the WD mass. We denote the orbital elements of a comet as: a for the semimajor axis, e for eccentricity, I for inclination, ω for argument of pericentre, Ω for longitude of ascending node, and l for mean anomaly. Orbital elements ω and I are measured relative to galactic plane. Some other quantities used to describe the comet orbits are: q for the pericentre distance, Q for the apocentre distance, and P for the orbital period. The comet is assumed to be a mass-less test particle. Orbital elements with the subscript p denote that they are the orbital elements for the WD companion (e.g. M_p is the companion mass). We also regularly use the following set of Delaunay action-angle variables (quantities Λ, L, L_z are actions in units of specific angular momentum, and l, ω, Ω are angles):

$$\begin{eqnarray} \Lambda &=& \sqrt{G M_* a}\qquad l \\ L &=& \Lambda \sqrt{1-e^2}\qquad \omega \\ L_z &=& L \cos I\qquad \Omega \end{eqnarray}$$

(1)

Λ is referred to as the circular angular momentum, L as the angular momentum, and L_z as the z-component of the angular momentum.

We use the terms ‘exo-Oort cloud’ and ‘Oort cloud’ interchangeably. An exo-Oort cloud is presumed to start at the inner semimajor axis edge, a₁, ends at a₂, follows a power-law number density profile n(a) ∝ a^−γ, has N_Oort comets and a total cloud mass M_Oort.

When referring to our own Solar system’s Oort cloud, we state ‘Solar system Oort cloud’ explicitly. In the Solar system, estimates for N_Oort typically range between 10¹¹–10¹² comets (e.g. Francis 2005; Boe et al. 2019). We assume N_Oort = 10¹¹ in this article. It is estimated that M_Oort ∼ 2 M_⊕ for the Solar system Oort cloud, which we will use as our fiducial value. In addition, numerical simulations show that the Solar system Oort cloud has a power-law exponent γ = 3.5 (e.g. Duncan et al. 1987; Higuchi & Kokubo 2015; Vokrouhlický et al. 2019). Note that the total mass of the Solar system Oort cloud is also quite uncertain from simulations and observations of incoming long-period comets (e.g. Weissman 1983; Boe et al. 2019).

A ‘Solar system-like Oort cloud’ is an exo-Oort cloud with a₁, a₂, γ, N_Oort, and M_Oort exactly like our own Solar system’s Oort cloud.

2 ANALYTIC THEORY: GALACTIC TIDE

2.1 Vertical tide model

Heisler & Tremaine (1986) study how the galactic tide affects Oort cloud comets through an analytic model. There are three main assumptions used in their model. First, the star–comet system moves in a circular orbit around the galaxy. Second, the most important component of the galactic tidal force is in the z-direction. z is defined as the direction orthogonal to the galactic mid-plane,¹ with the mid-plane defined at z = 0. In the Solar system, the second assumption is valid since galactic tidal terms in the x- and y-components are about one order of magnitude lower than the z term. Since observations of polluted WDs are typically within the Solar nseighbourhood, the second assumption is likewise not unreasonable to apply in nearby WD planetary systems. Third, the galactic tidal potential experienced by the comet is approximated as the potential inside a homogeneous slab with constant density ρ₀. ρ₀ is the averaged background density of gas and stars in the galaxy. As a star system orbits around the galaxy, it oscillates up and down the galactic mid-plane. As a result, ρ₀ also varies over time. Following previous work (e.g. Heisler & Tremaine 1986; O’Connor et al. 2023; Tremaine 2023), we adopt an averaged fiducial value of ρ₀ = 0.1M_⊙ pc⁻³ (Holmberg & Flynn 2000; McKee, Parravano & Hollenbach 2015).

The galactic potential with these assumptions can be written as (Heisler & Tremaine 1986):

$$\begin{eqnarray} \Phi _\mathrm{GT} = 2 \pi G \rho _0 z^2. \end{eqnarray}$$

(2)

This potential can be averaged over the orbit of the comet and then written in Delaunay elements as:

$$\begin{eqnarray} \langle \Phi _\mathrm{GT} \rangle = \frac{\pi \rho }{G M_*^2} \left(\frac{\Lambda }{L}\right)^2 (L^2-L_z^2) (L^2 + 5(\Lambda ^2-L^2)\sin ^2 \omega), \end{eqnarray}$$

(3)

from which we find the secular (orbit-averaged) equations of motion by applying Hamilton’s equations:

$$\begin{eqnarray} \left\langle \frac{{\rm {d}} L}{{\rm {d}} t} \right\rangle _\mathrm{GT} = -\frac{5 \pi G \rho }{(G M_*)^2} \left(\frac{\Lambda }{L}\right)^2(L^2-L_z^2) (\Lambda ^2-L^2) \sin (2 \omega), \end{eqnarray}$$

(4)

$$\begin{eqnarray} \left\langle \frac{{\rm{d}}\omega }{{\rm{d}} t}\right\rangle _\mathrm{GT} &=& \frac{2\pi G \rho }{(G M_*)^2}\left(\frac{\Lambda }{L}\right)^2 \left(\frac{L_z^2}{L} (L^2 + 5 (\Lambda ^2-L^2) \sin ^2 \omega \right. \\ &&\left.+\, (L^2-L_z^2) (L-5 L \sin ^2 \omega) \right) , \end{eqnarray}$$

(5)

$$\begin{eqnarray} \left\langle \frac{{\rm {d}} \Omega }{{\rm {d}} t}\right\rangle _\mathrm{GT} = \frac{2\pi G \rho }{(G M_*)^2} \left(\frac{\Lambda }{L}\right)^2 L_z ~\left(L^2 (5 \sin ^2\omega -1)-5 \Lambda ^2 \sin ^2\omega \right).\\ \end{eqnarray}$$

(6)

|$\dot{L}_\mathrm{GT}$| and |$\dot{\omega }_\mathrm{GT}$| are coupled differential equations. The phase space evolution described by these coupled equations is studied in detail in Heisler & Tremaine (1986) and Tremaine (2023). Briefly, they show that through these equations of motion, the comet can oscillate in the (L, ω) phase space due to galactic tide. There are two conserved Delaunay elements. Λ is conserved because the mean anomaly is not in 〈Φ_GT〉. Thus, the comet’s semimajor axis is conserved under the galactic tidal effect. Similarly, L_z is conserved because Ω is not in 〈Φ_GT〉. Since L_z is conserved but L is not conserved, this implies that galactic tide can excite the comet to very high eccentricity by exchanging eccentricity with inclination. This inclination-eccentricity exchange is periodic and is analogous to the von Zeipel–Lidov–Kozai mechanism (von Zeipel 1910; Kozai 1962; Lidov 1962).

2.2 Loss cone theory and comet injection rate

Heisler & Tremaine (1986) study the injection rate of Oort cloud comets into the Solar system. To do so, they use the galactic tide equations of motion, together with the loss cone theory framework as proposed by Lightman & Shapiro (1977). Without perturbations from a planetary companion, this can be used to estimate the rate of Oort cloud comets capable of being excited to any arbitrarily small pericentre distance.

First, we describe Oort cloud comets by a distribution function f, defined as the number of comets per volume of phase space:

$$\begin{eqnarray} {\rm {d}} N_\mathrm{Oort} = f(\Lambda , L, L_z, \omega , \Omega , l) ~{\rm {d}} \Lambda ~{\rm {d}} L ~{\rm {d}} L_z ~{\rm {d}} \omega ~{\rm {d}} \Omega ~{\rm {d}} l. \end{eqnarray}$$

(7)

Following previous work (Heisler & Tremaine 1986; Wiegert & Tremaine 1999; O’Connor et al. 2023), the distribution function is integrated over angles, assuming a spherical distribution of comets. This is motivated by observations in long-term Solar system Oort cloud simulations that the Oort cloud is spherically symmetric (e.g. Duncan et al. 1987; Vokrouhlický et al. 2019).

$$\begin{eqnarray} {\rm {d}} N_\mathrm{Oort} &=& f(\Lambda , L) {\rm {d}} L~ {\rm {d}} \Lambda ~ \int _{-L}^{L} {\rm {d}} L_z~ \int _0^{2\pi } {\rm {d}} \omega ~ \int _0^{2\pi } {\rm {d}} \Omega ~ \int _0^{2\pi } {\rm {d}} l \\ &=& (2\pi)^3 f(\Lambda , L) \cdot 2 L ~{\rm {d}} L~ {\rm {d}} \Lambda , \end{eqnarray}$$

(8)

where |$\int {\rm {d}} N_\mathrm{Oort} = N_\mathrm{Oort}$| with N_Oort as the total number of comets in the Oort cloud. With this description of the Oort cloud, we can proceed to find the pollution rate.

In the WD pollution case, we are interested in rate of comets that can be excited to q = q_crit = 1 R_⊙, the fiducial Roche limit of a WD that we adopt. At this distance, we assume that a comet will be tidally disrupted by the WD.

For comets with high eccentricity, e ≲ 1, the angular momentum can be related to the pericentre distance q as:

$$\begin{eqnarray} L \approx (2 G M_* q)^{1/2}. \end{eqnarray}$$

(9)

Through this, we define a critical angular momentum once a comet has achieved a certain critical pericentre distance q_crit:

$$\begin{eqnarray} L_\mathrm{crit}\equiv (2 G M_* q_\mathrm{crit})^{1/2}. \end{eqnarray}$$

(10)

The loss cone is defined as the phase space region where L ≤ L_crit. This is the tidally disrupted loss cone. Once a comet enters this loss cone, we assume it is removed from the Oort cloud (because it is tidally disrupted).

Comets can be injected into the loss cone in two regimes: the filled and empty loss cones, depending on their semimajor axis a (or equivalently, Λ). Intuitively the dependency on a is because these two regimes depend on how strong the galactic tide can affect a comet. In one case the galactic tide induces small changes in angular momentum ΔL over multiple orbits, slowly migrating a comet inward in q over many orbits. This is the empty loss cone case. In the other case the galactic tide is sufficiently strong to induce a large ΔL and the comet is capable of reaching the loss cone within one orbit. This is the filled loss cone case. Fig. 1 (following fig. 1 in O’Connor et al. 2023) is a diagram illustrating these two regimes.

Figure 1.

Diagram illustrating the loss cone theory by Heisler & Tremaine (1986). The engulfment loss cone is in grey and comets engulfed are shown as crosses. The regime on the left of a_crit is the empty loss cone, on the right is the full loss cone. This figure closely follows fig. 1 in O’Connor et al. (2023).

First, we consider the empty loss cone case, which happens when ΔL ≤ L_crit, where |$\Delta L \sim |{\rm {d}} L/{\rm {d}} t| \times P$| is the angular momentum change induced by the galactic tide per orbit. Heisler & Tremaine (1986) show that the comet injection rate in the empty loss cone regime, Γ_e (dimension of number of comets per unit Λ per unit time) is:

$$\begin{eqnarray} \Gamma _e (\Lambda) ~{\rm {d}} \Lambda = \frac{160 \pi ^3 G \rho _0}{3} \frac{L_\mathrm{crit}\Lambda ^4}{(G M_*)^2} f(\Lambda , L_\mathrm{crit}) ~{\rm {d}} \Lambda , \end{eqnarray}$$

(11)

where L_crit ≪ Λ (highly eccentric orbit). Γ_e is found by using equation (4) and calculate the rate at which comets are pushed into the loss cone boundary at L_crit.

Next, we consider the filled loss cone case when ΔL ≥ L_crit. The injection rate for the filled loss cone, Γ_f (same dimension as Γ_e) is:

$$\begin{eqnarray} \Gamma _f (\Lambda) ~{\rm {d}} \Lambda = 4 \pi ^2 (G M_*)^2\frac{L_\mathrm{crit}^2}{\Lambda ^3} f(\Lambda , L_\mathrm{crit}) ~{\rm {d}} \Lambda , \end{eqnarray}$$

(12)

where Γ_f is found by dividing the number of comets inside the loss cone by the comets’ orbital period (because these comets are lost within one orbit). Note that the subscript f here denotes ‘full’, not the distribution function f.

The loss cone is empty at small a and full at large a. The transition between the two cases can be found by equating the two rates, Γ_e = Γ_f, yielding:

$$\begin{eqnarray} a_\mathrm{crit}&=& \left(\frac{3}{20\sqrt{2}\pi } \frac{M_*}{ \rho _0} q_\mathrm{crit}^{1/2}\right)^{2/7}\approx 10~500 \mathrm{~au} \times \left(\frac{M_*}{0.6 \mathrm{ M}_\odot }\right)^{2/7} \\ &&\times \,\left(\frac{\rho _0}{0.1 \mathrm{ M}_\odot ~\mathrm{pc^{-3}}}\right)^{-2/7} \left(\frac{q_\mathrm{crit}}{1 \mathrm{ R}_\odot }\right)^{1/7}. \end{eqnarray}$$

(13)

When a < a_crit, a comets is in the empty loss cone regime. When a ≥ a_crit, it is in the full loss cone regime. In Delaunay variables, these conditions are equivalent to Λ < Λ_crit and Λ ≥ Λ_crit with

$$\begin{eqnarray} \Lambda _\mathrm{crit}\equiv (G M_* a_\mathrm{crit})^{1/2}. \end{eqnarray}$$

(14)

The total injection rate (number of comets entering a certain q_crit per unit time) can be found by adding up the two injection rates, integrated over Λ:

$$\begin{eqnarray} \Gamma _\mathrm{total} = \int _{\Lambda _1}^{\Lambda _\mathrm{crit}} \Gamma _e(\Lambda) ~{\rm {d}} \Lambda ~~+~ \int _{\Lambda _\mathrm{crit}}^{\Lambda _2} \Gamma _f(\Lambda) ~{\rm {d}} \Lambda , \end{eqnarray}$$

(15)

where Λ₁ = (GM_*a₁)^1/2 and Λ₂ = (GM_*a₂)^1/2, with a₁ the inner semimajor axis edge of the Oort cloud and a₂ the outer edge.

Finally, since Γ_e ∝ a²f(a), the injection rate per unit Λ increases from a₁ to a_crit. On the other hand, Γ_f ∝ a^−3/2f(a), the injection rate per unit Λ decreases from a_crit to a₂. At a_crit, Γ_e = Γ_f. Since the injection rate increases until a_crit and then decreases, the majority of the total injection rate is contributed from the region around a_crit.

2.3 Oort cloud distribution function

Now, we find the explicit form of the distribution function f. First, for a dynamically relaxed distribution, the distribution is ‘thermal’ and the distribution in e² is uniform (Jeans 1919; Ambartsumian 1937). This ‘thermal‘distribution of Oort cloud comets is seen after long-term simulations of Solar system Oort cloud (cf. fig. 7 in Vokrouhlický et al. 2019). In this case, f = f(Λ) is independent of L. We can further simplify by integrating over |$L = \sqrt{G M_* a (1-e^2)}$| which is now uniform between [0, (GM_*a)^1/2] = [0, Λ]:

$$\begin{eqnarray} {\rm {d}} N_\mathrm{Oort} = (2\pi)^3 f(\Lambda) ~{\rm {d}} \Lambda \int _0^{\Lambda } 2 L ~{\rm {d}} L~ = (2\pi)^3 f(\Lambda) \Lambda ^2 ~{\rm {d}} \Lambda . \end{eqnarray}$$

(16)

By definition, for a spherically distributed shell with width |$[a, a + {\rm {d}} a]$|⁠, we also have:

$$\begin{eqnarray} {\rm {d}} N_\mathrm{Oort} \equiv n(a) ~{\rm {d}} V = n(a) \times 4 \pi a^2 {\rm {d}} a , \end{eqnarray}$$

(17)

where, n(a) is the number density profile of comets.

Before finding f(Λ) through these equations, we first choose, n(a). Long-term simulations of the Solar system Oort cloud find that our own Oort cloud generally follow a power-law density profile (e.g. Duncan et al. 1987; Higuchi & Kokubo 2015; Vokrouhlický et al. 2019):

$$\begin{eqnarray} n(a) \propto a^{-\gamma } , \end{eqnarray}$$

(18)

where a ∈ [a₁, a₂], with a₁ and a₂ as the inner and outer semimajor axis edges of an Oort cloud. Following O’Connor et al. (2023), we consider 2 ≤ γ ≤ 4 in this article, with γ = 3.5 as a fiducial value for a Solar system-like Oort cloud.

With a description of n(a), we can find the distribution function f(Λ) by equating the two relations for |${\rm {d}} N$| (O’Connor et al. 2023):

$$\begin{eqnarray} f = f(\Lambda) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}C N_\mathrm{Oort} (\Lambda /\Lambda _1)^{3-2\gamma } &,\Lambda _1 \le \Lambda \le \Lambda _2 \\ 0 &,\text{otherwise}, \end{array}\right . \end{eqnarray}$$

(19)

where Λ₁ and Λ₂ are the circular angular momenta associated with a₁ and a₂. C is the normalization constant:

$$\begin{eqnarray} C = \frac{8 \pi }{(2 \pi)^3} \cdot \left\lbrace \begin{array}{@{}l@{\quad }l@{}}(\gamma -3)\left(4 \pi \Lambda _1^3 (1-(\Lambda _1/\Lambda _2)^{2\gamma -6}\right)^{-1} &,\gamma \ne 3\\ \left(8\pi ^3 \Lambda _1^3 \ln (\Lambda _2/\Lambda _1)\right)^{-1} &,\gamma = 3 \end{array}\right..\\ \end{eqnarray}$$

(20)

Equipped with the distribution function f(Λ), it is possible to find the total comet injection rate Γ_total (dimension of number of comets per unit time) by substituting f(Λ) and integrating equation (15). With Γ_total, we can also predict the pollution rate (i.e. in g s⁻¹) into a WD:

$$\begin{eqnarray} \dot{M}_Z = \Gamma _\mathrm{total} \times \frac{M_\mathrm{Oort}}{N_\mathrm{Oort}}. \end{eqnarray}$$

(21)

3 ANALYTIC THEORY: COMPANION DYNAMICS

In this section, we refer to a companion in a binary system with the WD as planetary if it has mass M ≤ 10⁻² M_⊙ = 10 M_Jup, or stellar if it has mass M ≥ 10⁻¹ M_⊙.

In Section 3.1, we analyse the effects of companion-induced precession on a comet with semimajor axis a and pericentre q from a companion with mass M_p on a circular orbit at semimajor axis a_p. Companion induced precession is compared against galactic tide induced precession. In Section 3.2, we study the efficiency of precession at preventing comets’ migration due to galactic tide through simulations.

3.1 Precession

Galactic tidal effects can be suppressed by angular momentum change induced by a companion, which is also accompanied by an apsidal precession. We compare the apsidal precession rates, |$\dot{\omega }$|⁠, induced by galactic tide and companion to study when each effect is dominant. In our case, we consider the regime where q ≥ a_p – no orbit crossings occur and the companion must be interior to the comet.

With secular forcing by a companion, a comet experiences apsidal precession (at the quadrupole order) at a rate (Farago & Laskar 2010):²

$$\begin{eqnarray} \langle \dot{\omega }_p \rangle &= \frac{M_p M_*}{(M_* + M_p)^2} \times \left(\frac{a_p}{a}\right)^{7/2} \times \frac{3 n_p}{8 (1-e^2)^{2}} \times \left(5\cos ^2 \Delta I-1\right), \end{eqnarray}$$

(22)

where n_p is the companion’s mean motion and ΔI is the mutual inclination between the comet and companion. This is a secular interaction, integrated over orbital motions of both companion and comet. There are higher order terms to |$\dot{\omega }_p$| (the next non-zero term occurs at the hexadecapole order, see Palacián & Yanguas 2006; Vinson & Chiang 2018), which become important as q → a_p. However, for our analysis in this section, the quadrupolar term is sufficient to give an estimate. When we require numerical results for |$\dot{\omega }_p$| (Fig. 6), |$\dot{\omega }_p$| is measured numerically and is not limited by this approximation.

Galactic tide likewise induces apsidal precession (equation 5). Companion-induced precession begins to dominate that of galactic tide when these rates are comparable:

$$\begin{eqnarray} \zeta \equiv \left| {\frac{\langle \dot{\omega }_{\mathrm{GT}}\rangle }{\langle \dot{\omega }_p\rangle }}\right| \sim 1. \end{eqnarray}$$

(23)

Precession in argument of pericentre (Δω) accompanies angular momentum change (ΔL). That is, the companion also induces a change in angular momentum, suppressing the angular momentum change from galactic tide and inhibiting further migration in pericentre q.

Expanding ζ to first order in q/a (the high eccentricity limit, e ∼ 1) and taking the order of magnitude terms, we find:

$$\begin{eqnarray} \zeta \sim \frac{32 \pi \sqrt{2}}{3} \times \alpha \times \beta \simeq 50 \times \alpha \times \beta , \end{eqnarray}$$

(24)

where we defined two dimensionless quantities α and β.

α is the ratio of densities between the galaxy and the binary system:

$$\begin{eqnarray} \alpha &\equiv& \frac{\rho _0}{M_\mathrm{reduced}/a_p^3} \\ &\approx& 10^{-9} \times \left(\frac{\rho _0}{0.1\,\mathrm{M_\odot ~pc^{-3}}}\right)\left(\frac{M_\mathrm{reduced}}{10^{-2}\,\mathrm{ M}_\odot }\right)^{-1} \left(\frac{a_p}{100\mathrm{~au}}\right)^3 , \end{eqnarray}$$

(25)

where M_reduced = M_*M_p/(M_* + M_p) is the reduced mass. As α becomes smaller, |$\rho _0 \ll M_\mathrm{reduced}/a_p^3$|⁠, we expect the effects of the binary to be stronger than that of galactic tide.

β compares the orbit of a comet to the companion’s orbit:

$$\begin{eqnarray} \beta &\equiv& \left(\frac{q}{a_p}\right)^{3/2} \left(\frac{a}{a_p}\right)^{7/2} = \left(1-e\right)^{3/2} \times \left(\frac{a}{a_p}\right)^5 \\ &\approx& 10^{7} \times \left(\frac{q}{a_p} \right)^{3/2} \left(\frac{a}{10^4 \mathrm{~au}}\right)^{7/2} \left(\frac{a_p}{100 \mathrm{~au}}\right)^{-7/2}. \end{eqnarray}$$

(26)

Higher β means that comet experiences less torque from the companion. β is mostly dominated by a/a_p. When a comet has a large orbit comparing to the companion (a ≫ a_p), it spends most of its orbital time far from the companion and thus, receiving less torque. β is large in this case. There is also a dependence on q (or equivalently, on the eccentricity e). As a comet migrates inward due to galactic tide, the orbit becomes more eccentricity, e increases, q decreases and thus β also decreases.

Typical values of α and β are shown in Tables 1 and 2. In Table 1, typical scenarios are evaluated: WD–WD binary, WD–M dwarf, and WD–large giant planet. Table 2 shows cases where the comet pericentre is at q = 5a_p. This is chosen because here, scattering is generally weaker comparing to both precession and galactic tide (see Fig. 6). A comet semimajor axis of 10⁴ au is chosen because this is typical for incoming comets. In both tables, we give sample values for a close-in companion at a_p = 10 au and a wider companion at a_p = 100–200 au.

Table 1.

Order-of-magnitude values of the dimensionless quantity α, which describes the ratio of densities between the galaxy and the binary system (equation 25).

α	Mp [M⊙]	M* [M⊙]	Mreduced [M⊙]	ap [au]
10−14	0.6	0.6	0.3	10
10−13	10−1	0.6	10−1	10
10−12	10−2	0.6	10−2	10
10−10	0.6	0.6	0.3	200
10−9	10−1	0.6	10−1	200
10−8	10−2	0.6	10−2	200

Table 2.

Order-of-magnitude β values for comets with pericentre distance q = 5a_p. The dimensionless quantity β describes the comet’s orbit relative to the companion orbit (equation 26).

β	ap [au]	q [au]	a [au]
1011	10	50	104
108	100	500	104
107	200	1000	104

Fig. 2 shows ζ on a grid of α and β. The dashed lines show the contours where ζ = 0.1, 1, and 10. Near and below ζ = 1, we expect companion-induced precession to be stronger than galactic tide; that is, precession can overcome galactic tide and suppress inward migration. We confirm the analytic values of ζ through simulations measuring the torque a comet experiences due to only galactic tide and only companion-induced precession.

Figure 2.

Values of ζ over a grid of α and β (see Tables 1 and 2 for typical values). Dashed lines are where ζ = 0.1, 1, and 10. In the top right corner, ζ ≫ 1 and galactic tide is stronger than companion-induced precession, and vice versa in the bottom left.

For close companions at a_p = 10 au, we have β > 10¹¹. In this case, to have precession dominates tide by making ζ ≲ 1, we require α < 10⁻¹³, which can be supplied by a stellar-mass companion with mass M_p ≥ 0.1 M_⊙ as shown in Table 1. However, a planetary-mass companion does not provide sufficient torque to overcome galactic tide. For more distant companions at a_p = 100–200 au, the same conclusions can be reached through their values of ζ. Therefore, a stellar-mass companion (M_p ≥ 0.1 M_⊙) is required to produce a precession barrier, which reduces the efficiency of galactic tide delivering comets from an exo-Oort cloud.

3.1.1 Limitations

We note that there are four limitations when using ζ:

• |$\dot{\omega }_p$| as shown in equation (22) is only valid for q ≥ a_p; that is, a comet must be strictly exterior to companion.

• We only use the quadrupole term for |$\dot{\omega }_p$|⁠. As q approaches a_p, contributions from higher order terms become important. Thus, ζ should not be used when q ≈ a_p.

• We ignore angular dependencies in I, ω, and Ω. However, as seen in Fig. 3 where we averaged results over angles, ζ without angular dependencies is still able to give an order of magnitude estimate of the strength between precession and galactic tide.

• We expand ζ in the limit q ≪ a. In other words, ζ is not a good approximation for very circular orbits (q ≈ a), or for very close-in comets encountering very widely separated companions (a is small, but a_p is large, so q is comparable to a). For the first case, when q ≈ a, the comet is on a circular orbit and mostly affected by galactic tide and not by the companion. For the second case, we must limit our analyses to companions with a_p ≪ a. The closest Oort cloud comets in our model have a ≈ 3000 au, so at most, we should only consider a_p ≤ 300 au for ζ analysis.

Figure 3.

Efficiency of the precession barrier induced by a companion, measured from numerical simulations. Simulation values of M_p, M_*, a, a_p, and q are chosen from grid values of (α, β) so that scattering does not matter. In other words, this efficiency is purely from effects of the torque induced by a planetary companion reducing the effectiveness of galactic tide. The shading shows the efficiency of the planet’s torque at inhibiting galactic tide, where N_companion is the number of comets that can enter a certain q in the existence of a companion and galactic tide and N_GT is the same but there is only galactic tide (no companion). Dashed lines show contours for ζ = 0.1, 1, and 10. Above the ζ = 10 dashed line, the precession barrier efficiency is 0 per cent and comets fully experience galactic tide and can be excited to high eccentricity, migrating inward in pericentre q. Below the ζ = 10 dashed line, we begin to see the precession barrier suppressing galactic tide and the efficiency decreases. At the very bottom left point, galactic tide is completely suppressed, and efficiency is ∼100 per cent.

3.2 Precession barrier efficiency

Fig. 3 presents the efficiency of the precession barrier at overcoming galactic tidal effects over a grid of α and β. The efficiency in this figure is defined as:

$$\begin{eqnarray} \mathrm{Efficiency} = 1-(N_\mathrm{companion}/N_\mathrm{GT}) , \end{eqnarray}$$

(27)

where N_GT is the number of comets that can enter a certain pericentre q with only galactic tide. N_companion is the same number, but there are galactic tide and a companion. As defined, the efficiency measures the effectiveness of a companion in suppressing tidal effects. A companion can do this either by inducing a precession barrier or by inducing a scattering barrier. In this subsection, we only consider the efficiency of the precession barrier.

We found N_GT and N_companion numerically. For each set of (α, β), two simulations are run: with and without a companion. In both simulations, galactic tide is included. This is evaluated on a grid of α and β, which depends on M_p, M_*, a, a_p, and e. These values are chosen so that the effect of scattering is always weaker than that of galactic tide and companion-induced precession. Other Keplerian orbital elements for a comet (I, ω, Ω) are drawn randomly assuming comets are isotropically distributed.

We compare the efficiency over a grid of (α, β) to the analytic values of ζ in Fig. 3. Above the ζ = 10 contour, companion-induced precession does not suppress any comets experiencing galactic tide. Here, 100  per cent of comets can come into q due to galactic tide. As ζ decreases, the precession barrier becomes stronger and eventually at the bottom left corner, all comets are suppressed from galactic tide, preventing inward migration in q. At ζ ∼ 1, when galactic tide is on the order of companion-induced precession, about half of the comets are suppressed from entering. Below ζ = 0.1, the barrier efficiency is 100 per cent and all comets are suppressed from galactic tide.

These efficiency behaviours found from simulations match well with what we expected from ζ over many orders of magnitude: when ζ ≲ 1 companion-induced precession suppresses galactic tide. We find that ζ is a useful indicator of where companion-induced precession is important relative to galactic tide. Specifically, the contour ζ ≈ 10 is a good indicator separating regimes of where galactic tide dominates and where precession dominates. When ζ decreases, precession begins to dominate galactic tide, and becomes increasingly more effective as ζ ≪ 1 where the efficiency tends to 100 per cent.

3.3 Scattering time-scale

Comets not only experience a torque in ΔL causing precession, but also experience a kick in energy causing a change in semimajor axis. There are two regimes where comets experience semimajor axis kicks from a massive companion: strong and diffusive scattering. In the strong scattering regime, q ≳ a_p, comets are ejected from the system within one pericentre passage. Here, the scattering time-scale is approximately T_scattering ≈ P_comet/2; comets are kicked during their incoming pericentre passages.

In the diffusive regime, q > a_p, comets experience small random semimajor axis kicks during multiple pericentre passages. These kicks vary in strength due to the phase difference between the comet and the companion during the encounters. Hadden & Tremaine (2024) investigate comet–companion interactions through an analytic mapping approach. Using the results they found, we can write the characteristic time-scale on which small diffusive kicks in semimajor axis will lead to an order unity change in binding energy (or equivalently, order unity change in semimajor axis):³

$$\begin{eqnarray} T_\mathrm{scattering} &\approx & ~1.5 \times 10^4 \mathrm{yr}~\times \left(\frac{M_p}{10^{-2} \mathrm{ M}_\odot }\right)^{-2} \left(\frac{M_*}{0.6 \mathrm{ M}_\odot }\right)^{3/2}\times \\ &&\left(\frac{a}{10^4 \mathrm{~au}}\right)^{3/2}\exp \left(7.4 \times \frac{q}{a_p}\right). \end{eqnarray}$$

(28)

There are three main assumptions to using the results found by Hadden & Tremaine (2024) diffusive time-scale for Oort cloud comets. First, the incoming comet is assumed to be coplanar to the companion-WD orbital plane. This is not true for Oort cloud comets, which are isotropically distributed. Second, M_p ≪ M_*, the companion’s mass is much smaller than the central star’s mass. Hence, we cannot use this time-scale to estimate the scattering time-scale in the stellar-mass companion case. Third, they assume the comet is on a parabolic orbit (e = 1). While this is not exactly true in our case, Oort cloud comets have highly eccentric orbits when they interact with companions. Hence, we further assume that we can use the e = 1 diffusive time-scale as written here, for comets with e ≲ 1. Finally, since we assume e ≲ 1, this time-scale is only applicable in cases where q ≪ a and a_p ≪ a (the companion is much closer-in than the comet).

Next, we compare analytic expectations of T_scattering to numerical simulations. We simulate two cases: a stellar-mass companion with M_p = 0.6 M_⊙ and planetary-mass companion with M_p = 10⁻² M_⊙. In both cases, the companions are on a circular orbit at a_p = 200 au around a central WD with mass M_* = 0.6 M_⊙ and the companion’s initial phase is randomized. At each pericentre q, where q ∈ [1, 6]a_p, 200 comets are initialized with a = 15 000 au (typical semimajor axis of incoming comets to interact with companions at a_p = 200 au). The comets’ initial position is set at apocentre. Comets’ angles are initialized either coplanar and isotropically. Simulations are stopped when comets experience a kick Δa/a₀ > 0.3 relative to the initial semimajor axis a₀, or until the simulation time reaches 10¹¹ yr and we call this time the numerical scattering time-scale. The numerical condition, Δa/a₀ > 0.3, is somewhat arbitrary but we found that this is a good indicator for when comets experience order unity changes in binding energy (strong scattering).

Fig. 4 compares numerical scattering time-scales with analytic expectations. In both panels, numerical time-scales for coplanar and isotropic comets are shown. The solid coloured lines are the mean scattering time-scales and the shaded areas are the time-scale ranges from 200 comets. The top panel shows the time-scales for a stellar-mass companion and the bottom panel for a planetary-mass companion. First, there are clearly two scattering regimes. Comets sufficiently far away experience long scattering time-scale according to the diffusive time-scale. As the pericentre decreases, the time-scale eventually converges to the strong scattering regime where comets are ejected within one pericentre passage. Second, in the bottom panel, we find the time-scale for coplanar comets interacting with a planet matches well within an order of magnitude with the analytic expectations. The slopes between numerical and analytic match well in the diffusive regimes. The analytic T_scattering is consistently off by a factor of a few. We attribute this to the arbitrary Δa/a₀ > 0.3 numerical scattering condition. Third, isotropic comets have a higher scattering time-scale. In addition, as the pericentre increases, coplanar and isotropic time-scales converge. This behaviour is expected since isotropic comets experience weaker kicks than coplanar comets. But at high pericentre distances, their kick strengths are both small. Fourth, in the top panel, we find the analytic expectation no longer gives a reliable scattering time-scale. This is because the condition M_p ≪ M_* is strongly violated. Here, we observe that coplanar comets quickly becoming strongly scattered at q ≈ 6a_p, while isotropic comets transition slower. Note that in an investigation of a test particle on an initially circular orbit around a binary, Holman & Wiegert (1999) found a stability limit around 4a_p for equal mass binaries. Their result is different from ours because their test particles are on circular orbits, while ours are highly eccentric. Hence, their test particles can experience effects from mean motion resonances, as discussed in detail in Holman & Wiegert (1999). That being said, their result and ours are reminiscent of each other.

Figure 4.

Numerical scattering time-scales for coplanar and isotropic comets interacting with a stellar-mass (top) and planetary (bottom) companion. The analytic diffusive time-scales are shown in both panels, but can only be used in the bottom panel where M_p ≪ M_*. Shaded areas are the range of measured numerical scattering time-scales. In both cases, the high diffusive time-scale flattens to the strong scattering regime (where comets are ejected within one pericentre passage with time-scale P_comet/2) as q decreases.

Once a comet reaches a pericentre within the strong scattering regime, they cannot survive multiple encounters. Thus, when there is a scattering barrier strong enough to create a strong scattering regime, comets in the empty loss cone are prevented from slowly migrating inward. However, this does not mean there are no comets entering smaller q. Comets in the full loss cone regime can still migrate within one orbital period to be engulfed by the WD because they only interact with the companion once. This is the motivation for the modified loss cone theory in Section 3.4.

3.4 Modified loss cone theory

O’Connor et al. (2023) estimate the reduction of pollution rate in the presence of a planetary companion by modifying the loss cone theory. We summarize their theory below and illustrate the modified loss cone theory in Fig. 5. We also discuss the limitations of this modified loss cone framework. A similar effect is proposed by Teboul, Stone & Ostriker (2024) for dense star clusters with a central black hole called ‘loss cone shielding’.

Figure 5.

Diagram illustrating the modified loss cone theory by O’Connor et al. (2023) with two loss cones. The engulfment loss cone is at the bottom from|$L=0$| to|$L_\mathrm{crit}$|⁠. The ejection loss cone is nested on top from |$L_\mathrm{crit}$| to |$L_\mathrm{ej}$|⁠. Comets ejected are shown as plus signs and comets capable of engulfment are shown as crosses. This figure is created following closely fig. 12 in O’Connor et al. (2023).

An ‘ejection loss cone’ is defined as the region where comets are shielded from further migrating inward due to repeated encounters with the planet. During multiple strong encounters, the comet experiences changes in semimajor axis which can eventually eject the comet. The location of the ejection loss cone is defined at:

$$\begin{eqnarray} L_\mathrm{ej} = \sqrt{2 G M_\mathrm{total} a_p} \end{eqnarray}$$

(29)

nesting on top of the tidal-disruption loss cone at

$$\begin{eqnarray} L_\mathrm{crit} = \sqrt{2 G M_\mathrm{total} q_\mathrm{crit}}, \end{eqnarray}$$

(30)

where M_total = M_* + M_p is the total mass of the WD star and the planetary companion. In the case where the planet is far away from the tidal radius, a_p ≫ q_crit, we have L_ej ≫ L_crit. The two loss cones overlap, but the ejection loss cone covers a much larger region of phase space than the tidal loss cone.

As a comet experiences a change in angular momentum ΔL due to galactic tide, it drifts inward encountering these loss cones. A comet with relatively small a experiences a small ΔL ≪ L_ej, and drifts slowly to the edge of the ejection loss cone; the ejection loss cone is empty. Vice versa, a comet with larger a experiences a much greater ΔL ≳ L_ej, and jumps through the ejection loss cone in one orbital period; the ejection loss cone is filled. In addition, if a comet experiences sufficiently large change in angular momentum, ΔL ≳ L_crit, it can jump through both loss cones and can be tidally disrupted and pollute a WD. These regimes are illustrated in Fig. 5 (following fig. 12 in O’Connor et al. 2023). The comet avoids ejection since it reaches the engulfment loss cone in one orbital period. The transition semimajor axis between the two ΔL regimes is:

$$\begin{eqnarray} a_\mathrm{crit, ej} &\sim& a_\mathrm{crit}\times \left(\frac{a_p}{q_\mathrm{crit}}\right)^{1/7} \\ &\approx& 32~000 \mathrm{~au} \times \left(\frac{a_p}{10\mathrm{~au}}\right)^{1/7} \left(\frac{q_\mathrm{crit}}{1 \mathrm{ R}_\odot }\right)^{-1/7} , \end{eqnarray}$$

(31)

where in the last expression we used values of a_crit = 10 500 au (equation 13), q_crit = 1 R_⊙, and a_p = 10 au. This transition semimajor axis corresponds to a critical circular momentum

$$\begin{eqnarray} \Lambda _\mathrm{crit, ej} = \sqrt{G M_\mathrm{total} a_\mathrm{crit, ej}}. \end{eqnarray}$$

(32)

Under this framework, only comets with a > a_{crit, ej} can experience a sufficiently large ΔL ≳ L_ej to pollute a WD. In addition, comets that can pollute WDs are in the filled loss cone regime because a_{crit, ej} > a_crit. Therefore, the pollution rate is estimated as the filled loss cone rate integrated over semimajor axes ranging between a_{crit, ej} ≤ a ≤ a₂, or equivalently over Λ:

$$\begin{eqnarray} \Gamma _\mathrm{total, planet} = \int _{\Lambda _{\mathrm{crit, ej}}}^{\Lambda _2} \Gamma _\mathrm{f}(\Lambda) ~{\rm {d}} \Lambda . \end{eqnarray}$$

(33)

O’Connor et al. (2023) further simplify this expression to find:

$$\begin{eqnarray} \frac{\Gamma _\mathrm{total, planet}}{\Gamma _\mathrm{total, GT}} \sim \left(\frac{q_\mathrm{crit}}{a_p}\right)^{(2\gamma -1)/14}. \end{eqnarray}$$

(34)

Taking a fiducial γ = 3.5 (for a Solar system Oort cloud), the factor on the right-hand side is approximately 0.05 for q_crit = 1 R_⊙ and a_p = 10 au. With this, we analytically expect the existence of a planetary companion to reduce WD pollution rate by about 1.5 orders of magnitude for γ = 3.5.

3.4.1 Limitations

First, comets with ΔL ≪ L_ej drift slowly to the edge of the ejection loss cone and then are assumed to be fully ejected when they enter into the ejection loss cone. However, comet–planet interactions can be weak and not sufficiently strong to eject comets. For example, we expect a comet to experience a much weaker kick in semimajor axis by a much smaller mass planet. A planet the size of the Earth and a planet with 10M_Jup will create loss cone barriers with very different efficiency. Another example is if a comet is very inclined relative to the planet’s orbital plane, the kick will also be much weaker. Therefore, not every comet with a < a_{crit, ej} will be ejected. In the small planet mass limit, we expect Γ_{total, planet} = Γ_{total, GT} since the planetary ejection barrier is not effective at all. In the large planet mass limit, the planet creates a 100 per cent effective scattering barrier, reducing pollution rate as described by equation (34). O’Connor et al. (2023) recognized that planet mass should affect ejection efficiency, yet equation (34) does not have a mass dependence. To identify which planet could create a sufficiently strong ejection barrier, they propose a quantity:⁴

$$\begin{eqnarray} \lambda = \left(\frac{M_p}{\mathrm{ M}_\mathrm{Jup}}\right) \left(\frac{\mathrm{ M}_\odot }{M_*}\right) \left(\frac{a}{10^4\,\mathrm{au}}\right) \left(\frac{10\,\mathrm{au}}{a_p}\right). \end{eqnarray}$$

(35)

If λ ≪ 1, the ejection barrier is weak and the comet receives negligible kicks and can safely migrate inward through multiple orbits. If λ > 1, then the ejection barrier becomes important. Using λ, a 10M_Jup planet at a_p = 10 au should be able to create a sufficiently strong ejection barrier. We will test the ejection barrier strength of a planet with this configuration later.

Second, this theory cannot account for additional dynamics that can be induced by a planetary mass companion. For example, a comet can be captured into smaller orbits and experiences more complicated dynamics which can also facilitate delivery into WDs, such as von Zeipel–Kozai–Lidov or inverse Kozai (Farago & Laskar 2010; Vinson & Chiang 2018).

Third, as discussed earlier, another important additional dynamics is that a comet experiences small random kicks in semimajor axis, Δa, through multiple diffusive encounters. This causes the semimajor axis to change over time which also changes how a comet experiences galactic tide ΔL over time. Since ΔL due to galactic tide is strongly dependent on a comet’s semimajor axis, these random walks in a induced by a planet can potentially cause a comet to experience vastly different galactic tidal effects over time. Therefore, assuming a WD pollution rate through integrating the total full loss cone rate Γ_f ranging between a fixed range of semimajor axes for all comets over all time (equation 33) might not be a suitable estimate.

Consider an example where a comet begins with an initial semimajor axis a < a_{crit, ej}. Through galactic tide, the comet migrates inward in pericentre distance q. At some point later in time when the comet achieves a q ≳ a_p, the comet begins to experience random kicks in semimajor axis every pericentre passages through interactions with the planet. Due to multiple small Δa kicks, the comet is migrated to a ≳ a_{crit, ej}. In this example, we initially count this comet to be in the empty ejection loss cone and thus, cannot pollute comets. But random interactions with a planet bring the comet to a larger a, where galactic tide induces a stronger ΔL allowing the comet to bypass the ejection loss cone and pollute the WD. In this example, having a larger planet might actually increase pollution rate because larger planets can induce stronger random walks in a at a larger range pericentre range.

In the Solar system, this mechanism for Oort cloud comets to bypass the Jupiter–Saturn ejection barrier is shown in Kaib & Quinn (2009). Here, weak perturbations by Uranus and Neptune are attributed to induce small kicks in a, bringing comets into a stronger ΔL regime capable of bypassing the ejection barrier. As shown in Fig. 4, one single planet can likewise induce small perturbations on comets. Thus, the comets in our case can similarly bypass the ejection barrier through these small kicks in a.

Finally, there is a small chance of a comet becoming unbound through a strong kick, but can still pollute a WD on its last inbound passage.

3.5 Time-scales comparison

To analyse the importance of the companion’s scattering and precession versus galactic tide, we investigate the time-scales on which these effects are important.

These time-scales are measured from simulations of 10⁵ comets from a Solar system-like Oort cloud: n(a) ∝ a^−3.5 from 3000 to 10⁵ au. Here, we directly simulate all comets with an additional galactic tidal force at every time-step – no secular integration of galactic tide as done in Section 4.1. This is to ensure that we accurately capture all companion–comet interactions at all q. In addition, we measure independently the contribution of galactic tide and companion at a given pericentre q.

First, we simulate 10⁵ comets around a centre-of-mass point mass (the total mass of a WD and its companion) with galactic tide. This simulation is run until comets’ pericentres reach a certain q. The orbital elements of incoming comets with pericentre q are then saved. Next, we measure the galactic tidal torque over one orbit at this q, giving us |$\dot{\omega }_\mathrm{GT} \approx \Delta \omega /P_\mathrm{comet}$|⁠. Second, we run another simulation on these saved comets, but now there is a companion and galactic tide is turned off. In this simulation, we measure the companion-induced precession change on the comet over one orbit, giving us |$\dot{\omega }_p \approx \Delta \omega /P_\mathrm{comet}$|⁠. Third, these comets are again run in a third simulation with a companion and no galactic tide to measure the scattering time-scale. The difference between the second and third simulations is in the third simulation comets are simulated longer. In the third simulation, comets are integrated until they experience a scattering event, Δa/a > 0.3, or until the simulation time reaches 10¹² yr. Recording the times when comets experience a scattering event gives us T_scattering. These simulations are run for every q between q = a_p to 5a_p and for companions with masses M_p = 10⁻², 0.1, and 0.6 M_⊙ (corresponding to a WD, M dwarf, and the upper limit of a planetary mass companion, 10⁻² M_⊙ = 10M_Jup).

We compare the rates (inverse time-scales) of these effects in Fig. 6. The shaded area is the 3σ spread, showing the range of most comets’ time-scales. These three cases show the increasing importance of companion-induced precession and scattering. As expected, increasing the companion mass increases the strength of scattering. With a 10⁻² M_⊙ planet companion, scattering begins to dominate at q ≈ 3a_p, whereas this is increased to q ≈ 3.5a_p for a 0.1 M_⊙ stellar companion. For the WD–WD binary case, scattering seems to be dominated in all q considered here. Note that the scattering time-scale flattens out as q → a_p. This is similar to the scattering behaviour we analysed earlier: comets are ejected within one orbital period in this regime, T_scattering ≈ P_comet/2. Similar to scattering, precession’s strength becomes stronger relative to galactic tide as companion mass increases. Finally, Fig. 6 shows that in all cases, as q ≃ a_p, the dominant effect is companion-induced scattering.

Figure 6.

Time-scale comparison of galactic tide, companion-induced precession, and scattering. From the bottom panel up, these show cases of increasing importance of companions. Time-scales are measured from simulations of a Solar system-like Oort cloud (n(a) ∝ a^−3.5, a ∈ [3000, 10⁵] au). The shaded area shows the 3σ spread of time-scales for this particular distribution of comets.

4 SIMULATION METHOD

The previous sections provide us with some understanding of what dynamical effects we should expect in this kind of system. We now turn to numerical simulations to study the long-term (1 Gyr) dynamics of Oort cloud comets under the influence of galactic tide and a companions. There are two main components in our integration scheme: first is a secular integration of the galactic tide equations of motion. If a comet can be excited to a certain q = q_switch, its orbital elements are extracted and the comet is then integrated directly with rebound (Rein & Liu 2012) when q ≤ q_switch. When comets violate their boundary conditions (Section 4.2) they are removed. Fig. 7 illustrates the components of our integration.

Figure 7.

A diagram to illustrate our hybrid integration scheme. The first stage is using secular galactic tide equations of motion in (equations 4–6) to quickly integrate a comet. If a comet can be excited to q = q_switch = 4a_p, then its orbital elements are extracted to be integrated directly using rebound. The direct integration by rebound uses the full equations of motion and can include interactions with the WD and its companion. A comet is extracted when its apocentre exceeds the WD’s Hill sphere (Q > 0.8 pc), when it is ejected due to interactions with a companion (a < 0) or when it is engulfed by the WD (d ≤ 1 R_⊙, where d is the distance between the comet and the WD). Integrations are stopped when the simulation time reaches t = 1 Gyr.

4.1 Integration scheme

The majority of Oort cloud comets at most times do not have sufficiently high eccentricity to interact with the central WD or its close companion.

In this regime, the dynamics of a comet is largely governed by its orbit around the WD + companion centre of mass and tidal effects from the galaxy. In this case, we can use the orbit-averaged equations of motion (equations 4–6) to quickly evolve comets. We use the fourth-order Runge–Kutta to evolve these coupled differential equations with an adaptive time-step scheme as implemented in scipy (Virtanen et al. 2020). We find that the Runge–Kutta adaptive time-step scheme reproduces well the evolution in (L, ω) as seen in Heisler & Tremaine (1986).

If a comet can reach a certain critical q_switch, it is removed from the secular integration and integrated with rebound to allow for full interactions with the WD and its companion. When switched over to integrating with rebound, we develop a fast simulation method where only one particle is simulated to further speed up simulation. This fast rebound integration method is illustrated in Fig. 8 and described in details in Appendix  A.

Figure 8.

A diagram illustrating our fast, properly time-adaptive integration method in rebound. There is only one particle in the simulation: the comet. The positions of the WD and its companions can be calculated analytically because they are in a two-body system. After each simulation time-step, the forces are analytically computed and added to the acceleration of the comet particle. IAS15 is used as an adaptive time-step and is capable of resolve time-steps as: large time-step when the comet is far away (the time-step is a fraction of the comet’s orbital period), and small time-step when the comet is close in (fraction of the companion’s orbital period). This ensures that all close encounters are properly handled while still maintaining fast integration speed. For further descriptions of the fast rebound integration method, see Appendix  A.

One caveat of using these orbit-averaged equations of motion is that they are not appropriate when ΔL ≫ L_crit. This is in the filled loss cone regime. In this case, within one orbital period, a comet experiences a significant change in it’s angular momentum and are lost. Thus, orbit-averaging is no longer appropriate. This is also resolved by switching over to rebound where we integrate with the full (not orbit-averaged) equations of motion.

When a comet achieves a pericentre distance q_switch ≈ 4a_p, interactions with a planetary mass companion become important. We show this earlier in Fig. 6 and its discussion. At q_switch, a comet is extracted from secular equations of motion integration and simulated in rebound. In rebound, we simulate full interactions between the comet and WD-companion system, and we use the full (not orbit-averaged) equations of motion for vertical galactic tide. The companion is set on a circular orbit at a semimajor axis a_p with I_p = 0.

When there is no planet, we still extract comets from the secular integration when the comets reach q_switch = 10 au. This choice is somewhat arbitrary, but it ensures that q_switch ≫ q_crit = 1 R_⊙. This is to make sure that engulfment into the WD is integrated with the full (not orbit-averaged) equations of motion through rebound, so that all galactic tidal dynamics are properly captured. As discussed, the secular equations of motion fails in the regime where ΔL is sufficiently strong to inject the comet into the loss cone in one orbital period.

4.2 Boundary condition

During the rebound part of our simulation, we enforce the following boundary conditions:

• comets with apocentre exceeding the Hill sphere of a 0.6 M_⊙ star (Q > 0.8 pc) and are outbound away from the WD,

• comets are ejected (a < 0) and are outbound away from the WD,

• comets are engulfed (d ≤ 1 R_⊙, where d is the distance from a comet to the WD).

Only boundary condition (iii) contributes to WD pollution. Boundary conditions (i) and (ii) remove comets because these comets are ejected from the Oort cloud reservoir. These three boundary conditions are checked at every simulation time-step. In addition, the simulation is stopped when the simulation time reaches t = 1 Gyr.

The Hill sphere of a WD is scaled down from the Solar system’s Hill sphere at ∼1 pc (e.g. Higuchi & Kokubo 2015). For a 0.6 M_⊙ WD, the Hill radius is at 0.8 pc (O’Connor et al. 2023). In addition, boundary conditions (i) and (iii)

$$\begin{eqnarray} q = a(1-e) &\le& 1 \mathrm{ R}_\odot \\ Q = a(1+e) &\le& 0.8 \mathrm{~pc} \end{eqnarray}$$

(36)

imply that there is a maximum semimajor axis:

$$\begin{eqnarray} a_2 \approx 85~000\,\mathrm{au}. \end{eqnarray}$$

(37)

This is the outer semimajor axis edge of the Oort cloud. Comets beyond a₂ cannot be excited to high eccentricity since they will be removed for exceeding the apocentre limit. The lowest pericentre comets at this semimajor axis without exceeding the WD Hill sphere is q ∼ 5000 au. In addition, this upper semimajor axis limit can also be found by scaling down the Solar system’s Oort cloud outer semimajor axis edge at a₂ = 10⁵ au, for a M_* = 0.6 M_⊙ central star.

4.3 Initial condition and sampling

The initial conditions of Oort cloud comets are generated based on a spherical cloud distribution. Comets’ argument of pericentre ω and longitude of the ascending node Ω are drawn randomly from |$\mathcal {U}[0, 2 \pi ]$| where |$\mathcal {U}$| is the uniform distribution. Comets’ inclination are drawn according to |$\cos I \sim \mathcal {U}[0, 1]$|⁠. The sign of the comet’s I is not important since we set the planet at I_p = 0.

Next, the semimajor axis is drawn such that Oort cloud comets follow a power-law density profile n(a) ∝ a^−γ:

$$\begin{eqnarray} {\rm {d}} N(a) \propto a^{-\gamma } a^2 {\rm {d}} a , \end{eqnarray}$$

(38)

with a ∈ [a₁, a₂]. a₁ = 3000 au is the inner semimajor axis edge, set based on the Solar system’s Oort cloud simulations (e.g. Duncan et al. 1987; Vokrouhlický et al. 2019). a₂ = 85 000 au is the outer semimajor axis edge, as found in the previous subsection.

To efficiently simulate all γ without re-running simulations many times, we use a rejection sampling scheme to draw samples for the semimajor axis distribution. This is described in Appendix  D.

After having a, we draw the squared eccentricity, e², from:

$$\begin{eqnarray} e^2 \sim \mathcal {U}[0, 1 - 2 q_\mathrm{initial, min}/a] \end{eqnarray}$$

(39)

for a distribution uniformly filling the energy phase space (Heisler 1990), appropriate for a dynamically relaxed Oort cloud as seen after long-term simulations of Solar system Oort cloud formation (e.g. Higuchi & Kokubo 2015; Vokrouhlický et al. 2019). Note that we impose an upper e (or equivalently, a minimum initial pericentre q_{initial, min}). q_{initial, min} set sufficiently far that a comet’s interaction with a companion or the WD is negligible initially. This is to ensure that all comets–companion–WD interactions are induced by galactic tide, rather than by random initial condition. Specifically, in the case where a companion exists, we set

$$\begin{eqnarray} q_\mathrm{initial, min} = 6 a_p. \end{eqnarray}$$

(40)

In the case where there is only galactic tide and a central WD, it is somewhat more arbitrary:

$$\begin{eqnarray} q_\mathrm{initial, min} = 15 \mathrm{~au}. \end{eqnarray}$$

(41)

This is to ensure that q_{initial, min} ≫ q_crit = 1 R_⊙. This ensures all comets will be initially integrated using the secular equation of motion and then switched over to be integrated by rebound until WD engulfment, ejection, or reaching 1 Gyr.

5 NUMERICAL WHITE DWARF POLLUTION RATE

5.1 Galactic tide only

We perform numerical simulations to explore the pollution rate of comets into a WD’s Roche limit (1 R_⊙) over various Oort cloud power-law structure (γ). Initial conditions and sampling are done according to Section 4.3, allowing us to sample a variety of power-law exponents γ. Here, we simulate with |$N_\mathrm{comets}^\mathrm{sim} = 4 \times 10^7$| comets. Because there are no companions in this case, comets are allowed to freely migrate inward in pericentre distance due to galactic tide.

First, we compare simulation rate of comet accretion with analytic prediction over various Oort cloud structures. The rates are shown in Fig. 9; the solid blue line represents simulation rates and the red-dashed line shows analytic expectations (equation 15). γ on the x-axis is the Oort cloud power-law exponent value; the number density of Oort cloud objects scale as n(a) ∝ a^−γ. The rate on the left is the number of accreted comets per Myr. On the right axis is the rate in g s⁻¹, assuming a fiducial Solar system Oort cloud mass and number of comets with N_Oort = 10¹¹ comets and a total cloud mass of M_Oort = 2 M_⊕. To get this rate, we count the total number of comets accreted, N_accreted, after the warm-up phase which is about 400 Myr (Section 5.5). Dividing N_accreted by the remaining 600 Myr of simulations yields the number of comets accreted over time as shown on the left-hand side of Fig. 9.

Figure 9.

Accretion rate of comets into q_crit = 1 R_⊙ over various Oort cloud structures γ, where n(a) ∝ a^−γ. The mean simulation rate is shown as the solid line and the analytic expectation (equation 15) is the dashed line. This accretion is due solely to the effects of galactic tide (no companion in this case). The equivalent mass accretion rate on the right is found by assuming an Oort cloud with N_Oort = 10¹¹ comets and M_Oort = 2 M_⊕. The shaded area shows the Poisson 1σ confidence interval of the solid line (mean rate).

First, we find that simulated rates match well with the analytic expectations from O’Connor et al. (2023) based on the framework by Heisler & Tremaine (1986).

Second, the average pollution rate of comets into WDs due solely to galactic tide is |$\dot{M}_Z \approx 5\times 10^7-10^8~\mathrm{g\, s^{-1}}$|⁠, depending on γ. Over the course of 1 Gyr, these rates correspond to the delivery ∼5 × 10⁻⁴ M_⊕ of materials. Thus, in the case of galactic tide alone where the only comet removal mechanism is engulfment by the WD, the Oort cloud reservoir is minimally depleted.

Finally, the 1σ blue-shaded area is found by assuming that comet engulfment is a Poisson process: comet accretion into a WD is a discrete event and comets arrive independently. Since we can count the total number of comets accreted, N_accreted, the Poisson process assumption allows us to estimate the standard deviation to be |$\sigma = \sqrt{N_\mathrm{accreted}}$|⁠. Recall that N_accreted is the total number of accreted comets counted over 600 Myr. Thus, |$\sqrt{N_\mathrm{accreted}}$| is the uncertainty of comets entering 1 R_⊙ due to our limited number of comets in our simulation (10⁷ comets) over a 600 Myr time-scale. In addition, we simulate a sufficiently large number of comets such that N_accreted is not a small integer and the 1σ interval can be meaningfully interpreted. Finally, |$\sqrt{N_\mathrm{accreted}}$| is based on the total number of comets in our simulation (∼10⁷ comets). If our simulation contained 10¹¹ comets like a full Oort cloud, the blue area would be smaller by 2 orders of magnitude. In summary, the blue shaded interval is the uncertainty of the total number of accreted comets, N_accreted, from our simulation containing ∼10⁷ comets, over a 600 Myr time-scale.

5.2 Efficiency of companion-induced precession and scattering

We now add a companion to our simulation. We investigate the effects of different companion masses on a Solar system-like Oort cloud. Fig. 10 shows the efficiency (equation 27) over different q with varying companion mass: M_p = 0.6, 10⁻¹, and 10⁻²M_⊙. In contrast to Fig. 3 previously, we now measure over a range of q. Hence, this efficiency can include the effects of both the precession and scattering barriers produced by a companion. Especially as q → a_p, both precession and scattering barriers can become very important in increasing the efficiency. Efficiency is 0 per cent when galactic tide is the only dominant effect. Vice versa, the efficiency is 100 per cent when precession and scattering barriers induced by a companion can suppress all galactic tidal effects.

Figure 10.

Efficiency of a companion at suppressing galactic tide and preventing comets from entering q. These are numerical results from simulations of incoming Solar system-like Oort cloud comets. The three cases here are the same as those in Fig. 6. The values of α for each case is shown on the right-hand side with their corresponding curve. At q = 5a_p, a typical value for β is ∼5 × 10⁷.

First, we use the formulations of ζ to predict if the precession barrier is stronger than galactic tide at q = 5a_p. We choose to focus at q = 5a_p because scattering is not important there for M_p = 10⁻² and 0.1 M_⊙ (cf. time-scales in Fig. 6). This can also be done at other q provided those q are within the limitations of our formulations, and scattering is not important. The three companion cases in Fig. 10 correspond to α = 10⁻¹⁰, 10⁻⁹, and 10⁻⁸. From simulations, typical incoming comets into a pericentre q = 5a_p = 1000 au have semimajor axes a ∼ 15 000 au, giving a β ≈ 5 × 10⁷. For α = 10⁻⁸, ζ = 25 ≫ 1 so we expect galactic tide to be dominant. For α = 10⁻⁹, ζ ≈ 1 so we expect companion-induced angular momentum change to be important and suppress some comets’ galactic tidal torque. For α = 10⁻¹⁰, the precession-barrier is dominant over galactic tide. Fig. 6 at q = 5a_p confirms these expectations.

Second, we can predict the efficiency of the precession barrier in reducing comet engulfment by using Fig. 3 with values of (α, β). At β ≈ 5 × 10⁷, Fig. 3 predicts about a 0 per cent efficiency for α = 10⁻⁸. At α = 10⁻⁹, it is expected to be around 40 per cent efficient. At α = 10⁻¹⁰, we expect about a 60 per cent efficiency. We confirm these predictions with Fig. 10 at q = 5a_p.

Third, scattering becomes important at q = 3 − 5a_p as seen in the time-scale analysis in Fig. 6. In Fig. 10, this corresponds to the fast increase in efficiency at that q range. In both cases, a stellar-mass companion increases the efficiency by almost 100 per cent by q = 1a_p due to strong precession and scattering barriers.

Fourth, combining precession and scattering effects, we predict that WDs with a stellar-mass companion are unlikely to be able to be polluted by an Oort cloud exterior to the companion. The actual rate of WD pollution rate in the presence of a stellar-mass companion is discussed later in Section 5.4. This is because the WD-star binary case is special due to the centre of mass of the system being far from the WD itself. Thus, to pollute a WD in the presence of a stellar companion, we need to also consider the efficiency of direct collisions between incoming comets and the WD itself. In contrast, with a planet companion, the centre of mass is close to the central WD and all comets migrated to q ∼ 1 R_⊙ will be engulfed by the WD.

Finally, for a planetary-mass (M_p ≤ 10 M_Jup) companion, planet-induced precession does not play a strong role. At q = 5a_p, the efficiency is still roughly 0 per cent. Thus, effects from the planet at the q distance is not important. There is a quick increase in efficiency beginning at q ≈ 2a_p due to scattering. However, this is increase is not as strong as the cases with stellar-mass companions to completely prevent further pollution. We will further analyse the pollution rate in the presence of a planetary-mass companion in Section 5.3.

5.3 Planetary-mass companion

As discussed in Section 3.4, additional effects due to having a planetary companion are complicated: having a planet can decrease pollution rate due to the ejection loss cone, but can also potentially increase pollution rate since comets can diffuse to higher semimajor axes to experience stronger galactic tide. We study the effects of having a planet through numerical simulations and compare with analytic expectations from Section 3.4. The simulation here uses the same methodology as discussed in Section 4 with |$N_\mathrm{comets}^\mathrm{sim} = 10^8$| comets.

In Fig. 11, we present the simulated WD pollution rate at various γ in the presence of a M_p = 10 M_Jup = 10⁻²M_⊙ planet at a_p = 10 au. These planet mass and semimajor axis values are chosen because this configuration gives a value of λ ≈ 10 (equation 35), which is predicted to create a strong ejection barrier. In this figure, we find that having a planet decreases the pollution rate into q_crit = 1 R_⊙ by about 1 order of magnitude. Assuming a Solar system Oort cloud with N_Oort = 10¹¹ comets and a total cloud mass of M_Oort = 2 M_⊕, the WD pollution rate in this case is |$\dot{M}_Z^\mathrm{(planet)} \sim 10^7~\mathrm{g~s^{-1}}$|⁠. We further observe that the simulation rate does not match analytic expectations. As seen in the plot, simulation rate is about 2–5 times higher than predicted analytic rate. Furthermore, simulation rate matches better to analytic expectation at low γ than at high γ.

Figure 11.

Accretion rate of comets into q_crit = 1 R_⊙ over various Oort cloud structures γ, where n(a) ∝ a^−γ. In this simulation, there is a 10M_Jup planetary companion at 10 au. The mean simulation rate is shown in as a solid line with shaded area. The analytic rate in the case where there is no planet (equation 15) is shown as a solid line without shading for comparison. The analytic prediction in the presence of a planet that can create a sufficiently strong ejection barrier (equation 33) is shown as a dashed line. The shaded area around the simulation rate shows the Poisson 1σ confidence interval of the mean rate.

First, we discuss the behaviour of the analytic expectation. We notice that the analytic rate decreases as γ increases. This can be understood intuitively through analysing where comets are distributed relative to a_{crit, ej}. In the γ = 4 limit, comets are more centrally distributed. Thus, more comets have a < a_{crit, ej}. Recall that these are the comets that experience small ΔL < L_ej due to galactic tide and migrate slowly in q until they are ejected through encounters with the planet. Therefore, since most comets are centrally distributed, they have a < a_{crit, ej} and we analytically expect most comets to be ejected, reducing the pollution rate. On the other hand, in the γ = 2 limit, comets are less centrally distributed, some still have a < a_{crit, ej} but not as many as in the case of γ = 4. Therefore, we have more comets with a > a_{crit, ej}. These comets are capable of experiencing a strong ΔL ≳ L_crit, drifting through the ejection loss cone in one orbit, polluting the WD. In other words, comet ejection is more effective at γ = 4 than at γ = 2 because there are more comets available to be ejected.

Second, we discuss why simulation rate matches with theory better at γ = 2 than at γ = 4. With the same intuition where comets are distributed, we further consider that the ejection loss cone is not 100 per cent effective. If the barrier is 100 per cent effective, expect all comets with a < a_{crit, ej} to be all ejected. However, because the ejection loss cone is not 100 per cent effective, some comets are capable of drifting through the ejection loss cone and eventually pollute the WD. Since there are more comets at γ = 4 with a < a_{crit, ej} than at γ = 2, the assumption of having a 100 per cent effective ejection loss cone leads us to overestimate the reduction of pollution rate more at γ = 4 than at γ = 2.

Next, we analyse the distribution of accreted comets’ initial semimajor axes, a_initial in Fig. 12. The initial semimajor axis is shown on the x-axis because comets experience kicks in a over time.

Figure 12.

Distribution of accreted comets’ initial semimajor axis at γ = 2. Histograms are counted per a_initial bin. The two histograms show distributions when there is galactic tide only, and when there is galactic tide and a 10 M_Jup planet. Vertical lines show a_{crit, ej} (dashed) and a_crit (solid).

First, in the regime where a_initial < a_{crit, ej}, the planet significantly reduces pollution rate as predicted. However, we still see comets with a_initial < a_{crit, ej} polluting the WD. We find that the ejection loss cone barrier is not 100  per cent effective. Therefore, the assumption in equation (33) that the pollution rate does not have any contribution with comets from a < a_{crit, ej} gives an underestimate of the pollution rate.

Second, in the regime where a_initial ≤ a_{crit, ej}, the number of comets capable of polluting the WD increases steadily. As a_initial → a_{crit, ej}, comets need a smaller kick in Δa to deliver them over into the a ≥ a_{crit, ej} regime where galactic tide can induce a strong enough ΔL to drift them through the ejection loss cone. On the other hand, in the a_initial ≪ a_{crit, ej}, comets will need to experience multiple interactions with the planet that increases their semimajor axis, but not strong enough to eject them.

Third, as a increases beyond a_{crit, ej}, the pollution rate is higher than in the case of no planet. There are several mechanisms to explain this. First, these comets can be kicked into higher a, allowing stronger ΔL to migrate further in, as described in Section 3.4.1. Second, these comets have lower energy (because of large a) and can be kicked into much smaller orbits. At that point, they can be excited to high eccentricity and pollute the WD with effects like von Zeipel–Kozai–Lidov or inverse Kozai. Third, they have lower energy, and are also easier to be ejected. However, in their last inbound passage, they have pericentre distances sufficiently low to pollute the WD. None of these additional dynamics would be possible without perturbations to the comet. Without additional perturbations like a planet, a comet at a = 50 000 au for example, will stay there and if it cannot reach q_crit during a galactic tide cycle, will never be able to do so. Hence, when we assume in equation (33) that the pollution rate strictly follows Γ_f based on a fixed distribution of a_initial, a lot of these additional dynamics are ignored giving an incorrect rate in that regime. As we have seen, the contribution in region a > a_{crit, ej} is higher than expected in equation (33) leading to another source of underestimation of pollution rate.

In summary, we find that the existence of a planetary mass companion significantly reduces the pollution rate for comets with initial semimajor axes a_initial ≤ a_{crit, ej} as predicted by O’Connor et al. (2023). However, we find that this reduction is not 100 per cent effective, that comets experience rich dynamics, and that beyond a_initial > a_{crit, ej} the pollution rate does not simply follow the full loss cone rate Γ_f. That being said, when comparing the overall rates in Fig. 11, the analytic predictions by O’Connor et al. (2023) still yield a good order of magnitude estimate for the pollution rate, although it can be off by a factor of 2–5 times. We find that the WD pollution rate in the presence of a planetary companion is reduced by one magnitude to |$\dot{M}_Z^\mathrm{(planet)} \sim 10^7~\mathrm{g~s^{-1}}$|⁠, assuming a Solar system Oort cloud.

5.4 WD–WD binary

An interesting question is how a WD–WD binary would be polluted by an Oort cloud reservoir. Unlike the planetary companion case, when the companion’s mass is comparable to the WD, the WD is very far from the centre of mass. In this scenario, pollution is governed by direct collisions between incoming comets and the WDs on their orbits.

Assume we have two WDs, separated by a_p on circular orbits. The Safronov number can be used to estimate how probable it is to expect collisions of incoming comets on highly eccentric orbits:

$$\begin{eqnarray} \Theta = \frac{v_\mathrm{esc}^2}{v_c^2} = 2 \times \frac{M_p}{M_* + M_p} \times \frac{a_p}{R} \approx 2 \times 10^3 , \end{eqnarray}$$

(42)

where v_esc is the escape velocity from the companion’s surface (at distance R) and v_c is the circular speed at the companion’s semimajor axis. We set R = 1 R_⊙ for the tidal radius where comets are captured by a WD and a_p = 10 au. Since Θ ≫ 1, we predict that collisions into the tidal radius are unlikely (Tremaine 2023).

In Fig. 13, we show the pollution rate of comets into one WD in a WD–WD binary separated by a_p = 10 au at various γ. In the simulation here, we directly check every time-step if a comet’s distance is within a WD’s tidal radius. The switching point between secular integration and direct integration with rebound is increased to q = 6a_p to accurately capture all dynamics between comets and a stellar-mass companions (see Figs 4 and 6). Here, we find that the pollution rate is significantly reduced by 2.5–3 orders of magnitude compared to the galactic tide only case. Assuming a Solar system Oort cloud, the pollution rate is about 2–4 × 10⁵ g s⁻¹. Note that this is just below the detection limit at ∼5 × 10⁵ g s⁻¹. This rate is low due to a combination of effects. First, as shown earlier, ζ ∼ 0.1 and comets experience a precession torque reducing the effectiveness of galactic tide in exciting incoming comets. Second, comets in the empty loss cone cannot slowly migrate inwards as they are strongly scattered through multiple encounters with the strong scattering barrier. Third, comets that can reach sufficiently low q ≈ a_p are still unlikely to collide with the WD since the Safronov number is very high.

Figure 13.

Accretion rate of comets into one WD in a WD–WD binary over various Oort cloud structures γ. The binary is in circular orbit and separated by 10 au. The pollution rate in g s⁻¹ is shown on the right axis, assuming a Solar system Oort cloud.

Next, we consider two extreme cases: a very close WD–WD binary and a more widely separated one. If the WD–WD separation is smaller, we expect the pollution rate to increase, peaking at half of the normal galactic tide only rate. For example, take the limit where the WDs are separated by only a few solar radii, then they are both near the centre of mass and will be impacted by all incoming comets. Since there are two stars, the pollution rate (from the galactic tide only case) will be reduced by half as comets are equally likely to collide with either stars. If the WD–WD separation is larger, the Safronov number, Θ ∝ a_p, would increase and collisions would be even more unlikely. Thus, we expect a decrease in pollution rate as a_p increases.

Finally, the results here also give an order-of-magnitude estimate of Oort cloud comet pollution rate for other stellar-mass companions, assuming M_p ∼ M_*.

5.5 Pollution over time

In this subsection, we analyse the pollution rate over time to study if the Oort cloud as a reservoir can consistently maintain WD pollution rate over a Gyr time-scale.

In Fig. 14, we show the pollution rate over time for γ = 2.5 and 3.5 in the case of galactic tide only (no companion). During the first ∼ 400 Myr, we observe a ‘warm-up’ phase. Intuitively, this is because it takes time for comets to experience galactic tide, migrate in q, and arrive at q_crit. The time-scale of 400 Myr is consistent with the typical galactic tide cycle period, which is on the order of 300 Myr (estimated in equation 18 in Heisler & Tremaine 1986). After the first 400 Myr, we find that the simulation rate matches well with analytic predictions. In addition, the rate stays constant with no signs of reduction over a 1 Gyr time-scale. Thus, without a companion, an Oort cloud is capable of delivering materials into the WD tidal disruption zone at a constant rate between |$\dot{M}_Z \approx 5\times 10^7$| and 10⁸ g s⁻¹ over a 1 Gyr time-scale.

Figure 14.

Accretion rate over a 1 Gyr time-scale from Oort clouds with γ = 2.5 and 3.5. There are no companions in this case. Analytic rates are shown as horizontal dashed lines. Accretion rates from simulations is binned per Myr. A moving average with a sliding window of 10 Myr is also shown. The accretion rate can be sustained over a 1 Gyr time period. An equivalent mass accretion rate on the right is found by assuming and Oort cloud with N_Oort = 10¹¹ comets and total cloud mass M_Oort = 2 M_⊕.

Fig. 15 shows the pollution rate over time for γ = 2.5 and 3.5 in the case of galactic tide and a planetary mass companion. The simulation rate does not match as well with the analytic prediction (equation 34). Recall that there are some limitations in the analytic predictions for WD pollution rate in the presence of a companion, since this framework does not include additional complicated effects induced by scattering. This is typically off by about a factor of 2–3, consistent with what is seen in Section 5.4.

Figure 15.

Accretion rate over a 1 Gyr time-scale from Oort clouds with γ = 2.5 and 3.5. There is a planetary-mass companion with M_p = 10M_Jup at a_p = 10 au. Accretion rates from simulations (red) is binned per Myr. Analytic expectations are shown as horizontal blue lines. A moving average with a sliding window of 10 Myr is shown in grey. The accretion rate is slightly decreased over the 1 Gyr period, but no more than 0.5 dex.

Furthermore, we again find that the accretion rate does not significantly decline over a 1 Gyr time-scale. Note that in the case where γ = 2.5 (bottom panel), the accretion decreases by about a factor of 3–4 over a 1 Gyr time-scale. Note that in the context of WD pollution rate which ranges 5 orders of magnitude, this reduction factor of 3–4 over a Gyr time-scale is not significant. The reduction factor for γ = 3.5 is even less, where the accretion rate only reduces by about a factor of 2. This can be explained by the distribution of comets. In the case of γ = 2.5, there are more comets at larger at semimajor axes than γ = 3.5 and thus, are easier to be ejected. Hence, the reservoir depletes quicker at low γ than high γ. None of this is observed in the galactic tide only case because there, comets semimajor axes are conserved and the reservoir does not get depleted by ejection.

Finally, there are variations over time in the accretion rate over time for both cases with and without planet. These variations are due to the limited number of comets in our simulation. Similar to the discussion in Section 5.1, the variations over every Myr bin is because of our resolution-limited simulation (∼10⁷ comets). This is about 4 orders of magnitude smaller than the actual Oort cloud population of 10¹¹ comets. Scaling these to a full Oort cloud, the variations would be smaller by 2 orders of magnitude. We observe that the variations per Myr bin in our simulations (Figs 14 and 15) are within an order of magnitude of the mean. Thus, for a full Oort cloud, these variations would be within 10⁻¹ of the mean rate.

In summary, we find that in the presence of a planetary mass companion and galactic tide, an Oort cloud is capable of delivering materials at a relatively constant rate |$\dot{M}_Z^\mathrm{(planet)} \approx 1-3\times 10^7~\mathrm{g\times s^{-1}}$| over a 1 Gyr time-scale. Depending on γ, this rate may be reduced by a factor of 3–4 over this time-scale.

6 DISCUSSION

Fig. 4 of Blouin & Xu (2022) presents observed WD pollution rates for WDs with ages ranging from 1 to 8 Gyr. In this figure, they show that the pollution rate ranges between 10⁵ to 10¹⁰ g s⁻¹, with the majority ranging between 10⁶ to 10⁸ g s⁻¹. They also find that pollution rates decrease by no more than one order of magnitude over the course 8 Gyr. Recent observations by Mullally et al. (2024) find giant planet candidates around two polluted WDs with ages 1.5 and 5 Gyr. The planets have masses 1–7 M_Jup and have separations 10–35 au. Therefore, observational evidence suggests that giant planets do not significantly reduce WD pollution rate.

Through simulations, we find that an Oort cloud (with total number of objects and total mass like the Solar system) is capable of delivering materials into a WD Roche radius at a rate from 5 × 10⁶ to 10⁸ g s⁻¹, depending on the existence of a planetary companion. These rates are above the current detection limit at ∼5 × 10⁵g s⁻¹. Furthermore, these pollution rates can be sustained and decrease no more than a factor of 3–4 (in the existence of a planet) over a Gyr time-scale. We find that the simulated rates found here can explain a significant portion of observed pollution rates (e.g. Blouin & Xu 2022) and sustain that rate over a Gyr time-scale. Our results further show that Oort cloud comets can pollute WDs in evolved systems with giant planets as observed by Mullally et al. (2024).

In the scenario of WD–WD binary with separation a_p = 10 au, the pollution rate is ≈3 × 10⁵ g s⁻¹, below the detection limit. This rate is significantly lower than other cases due to a combination of precession, scattering, and the low chance of collisions. This rate can be used as an order-of-magnitude estimate for other stellar-mass companion cases where M_p ∼ M_*.

We note that an advantage of the mechanism and reservoir presented here is that it is time-independent, in contrast to other works (cf. Debes et al. 2012; Mustill et al. 2014; Smallwood et al. 2018). The pollution rate from an Oort cloud would only decrease if the reservoir is significantly depleted. Over an 8 Gyr time-scale with a comet injection rate of ∼10⁸ g s⁻¹, an Oort cloud would only deplete by roughly 5 × 10⁻³M_⊕ worth of materials, which is much less than the current Solar system Oort cloud reservoir of 2 M_⊕. In the presence of other perturbers like a planet or stellar flybys, the Oort cloud reservoir would be depleted faster due to ejection. As shown previously, with a planet we find that the pollution rate can decrease at most by a factor of 3–4 in the course of 1 Gyr.

We further discuss below the scalability of the Oort cloud, the potential impacts of stellar flybys, concerns of using the Oort cloud as a reservoir for WD pollution, the effects of very close companions, wide companions, and delivery from an accretion disc into a WD.

6.1 Robustness of the Oort cloud

The results in this work can be scaled to other exo-Oort clouds with different N_Oort and M_Oort by:

$$\begin{eqnarray} \Gamma _\mathrm{total}^{\mathrm{new}} = \Gamma _\mathrm{total} \times \left(\frac{N_\mathrm{Oort}^{\mathrm{new}}}{10^{11}}\right) , \end{eqnarray}$$

(43)

$$\begin{eqnarray} \dot{M}_Z^{\mathrm{new}} = \dot{M}_Z \times \left(\frac{N_\mathrm{Oort}^{\mathrm{new}}}{10^{11}}\right) \left(\frac{M_\mathrm{Oort}^{\mathrm{new}}}{2 \mathrm{ M}_{\oplus }}\right)^{-1}. \end{eqnarray}$$

(44)

The other parameters of an Oort cloud, γ, a₁ (inner edge of the Oort cloud), and a₂ (outer edge), do not significantly affect pollution rate. First, we studied various γ in this article and find that γ does not affect the pollution rate by more than a factor of 2 (Figs 9, 11, and13). Second, we fixed the inner semimajor axis edge at a₁ = 3000 au, based on the Solar system Oort cloud. O’Connor et al. (2023) find that varying a₁ changes the pollution rate by at most a factor of 4 for 10³ au ≤ a₁ ≤ 10⁴ au. Third, we fixed the outer semimajor axis edge at a₂ = 85 000 au. This is because we have shown that this outer edge is a natural consequence of enforcing the boundary condition of the WD Hill sphere at 0.8 pc. Thus, a₂ would not be different in another Oort cloud around a typical WD.

6.2 Stellar flybys

Stellar flybys are an additional mechanism which can both reduce and increase pollution rate. First, a strong flyby (slow with small impact parameter) could potentially induce strong scattering and significantly deplete an Oort cloud. We do not consider stellar flybys in this work, thus it is unclear to us how flybys would deplete the Oort cloud reservoir. However, Higuchi & Kokubo (2015) show that in simulations of Oort clouds with galactic tide and impulsive stellar flybys over different γ structures, the e-folding decay time-scale for the Oort cloud population is 4–18 Gyr. Thus, even with stellar flybys, we still do not expect a strong (more than one order of magnitude) decrease in pollution rate within a Gyr time-scale because the reservoir is not significantly depleted.

Second, distant flybys can stochastically perturb comets and cause them to diffuse into small pericentres. Hence, weak flybys can act as another mechanism, in addition to galactic tide, to deliver comets from the Oort cloud into small pericentre. This was first explored by Heisler & Tremaine (1986) and used by O’Connor et al. (2023) to estimate that the effects distant flybys contribute to the comet injection rate is on the order of the pollution rate from galactic tide alone.

Third, a strong flyby can cause comet showers (Heisler 1990). In the Solar system, these showers increase the injection rate of Oort cloud comets into q ∼ 10 au. In the Solar system these showers have increased injection rate as much as two orders of magnitude within a few Myr – much shorter than the galactic tidal time-scale. The precession and scattering barriers we discussed earlier cannot prevent these comet showers because comets are induced into low pericentre within one orbital period. In the loss cones formulation, this is equivalent to a significant number of comets are in the filled loss cone regime and are able to bypass both angular momentum and scattering barriers. In the context of exo-Oort cloud around WDs, this is perhaps another mechanism to not only deliver materials, but also potentially explain the observed spread in pollution rate (5 orders of magnitude).

6.3 Surviving stellar evolution

One major concern of having an Oort cloud as a potential material reservoir for WD pollution is that the Oort cloud might be ejected during the evolutionary process from main sequence to WDs. A comet at 10 000 au has a typical orbital speed of ∼0.3 km s⁻¹. This is lower than the typical speed the natal, anisotropic recoil kick a WD experiences during its rapid mass loss phase, which is about 0.75 km s⁻¹ (El-Badry & Rix 2018). The star also undergoes mass loss at the same time. Thus, it is a concern if objects in the Oort cloud can be kept bounded to its central star. O’Connor et al. (2023) also investigate this question to find that a post-main-sequence evolution Oort cloud retains about 10 per cent of its original objects, with a fairly complex cloud structure. Thus, an Oort cloud can remain bound, albeit with less materials, after an anisotropic mass loss during stellar evolution.

O’Connor et al. (2023) also find that the pollution rate after main-sequence evolution is reduced by about an order of magnitude for a typical kick strength of 0.75 km s⁻¹. Hence, if we assume a more complicated post-main-sequence evolution Oort cloud where the original Oort cloud is like our current Solar system Oort cloud, the pollution rates in our work is reduced by another order of magnitude. This gives the post-stellar evolution pollution rate to be in the range between 10⁶–10⁷ g s⁻¹, depending on the existence of a planet. If the companion is stellar-mass, the pollution rate is much lower at ∼3 × 10⁴ g s⁻¹.

However, the mass and structure of exo-Oort clouds in other main sequence or post-main-sequence systems are unknown. This could point to a much wider range of Oort clouds’ total mass or number of objects. In addition, Solar system Oort cloud formation and characterization remains an active area of research. For example, simulations of Solar system Oort cloud formation (e.g. Vokrouhlický et al. 2019) include a very specific migration history of the planets following the Nice model. The Solar system’s Oort cloud total mass, number of objects, and chemical composition are heavily dependent on the early configuration of the giant planets. Thus, an exo-Oort cloud might very well be different in number of objects, mass, and chemical composition from what we currently observe of our own Oort cloud. Therefore, it is unclear how well we can extrapolate the mass of the post-main-sequence Solar system Oort cloud to other WD planetary systems.

With the uncertainty of formation models and complicated post-main-sequence evolution, it is unclear how exo-Oort clouds look around WDs. Therefore, in this article, we choose to simplify by only answering the question if an exo-Oort cloud with a n(a) ∝ a^−γ radial density profile – like the one currently existing in the Solar system as we currently understand it – can pollute a WD.

6.4 Cometary composition

Another major concern for the Oort cloud as a potential material reservoir for WD pollution is that we mostly observe volatiles-poor polluted WD atmospheres (e.g. Jura 2006; Jura & Xu 2012; Doyle et al. 2019), with accreting material composition resembling of asteroids or rocky planets in the Solar system. These observed compositions are inconsistent with the volatiles-rich, icy bulk composition of typical Solar system comets. Solar system Oort cloud comets are expected to be mostly ice because they are ejected from the protoplanetary disc beyond the ice line due to interactions with the Uranus and Neptune (Vokrouhlický et al. 2019). We discuss ways to reconcile these observations with our numerical predictions that Oort clouds, at least those like our own, should be able to pollute their WD.

First, the composition of objects in the Oort cloud in our Solar system and especially in other exo-Oort clouds might not be composed of only icy comets. In our own Solar system, for example, Vida et al. (2023) recently observe a small rocky object with origins from the Oort cloud. Simulations of the Solar system Oort cloud indicate that about 4  per cent of objects (up to 8 × 10⁹ objects) are rocky asteroids (Shannon et al. 2015). In addition, our own Oort cloud is ‘icy’ and volatile-rich because of the early configuration of the giant planets. Because of this, it is also unclear regarding other details of Oort clouds in other systems. We have observations of protoplanetary discs with diverse arrangements of giant planets configurations. These giant planets do not always stay fixed beyond the ice line. Thus, it is uncertain if objects interacting with different giant planets to be injected into exo-Oort clouds are also formed within or outside of the ice line (e.g. Doner et al. 2024). Therefore, it is problematic to assume Oort clouds around WDs or even in our own Solar system to be solely composed of volatile-rich, icy comets.

To add to the complexity, we also observe some polluted WDs with volatiles in their atmospheres (e.g. Farihi et al. 2013; Doyle et al. 2021; Klein et al. 2021). This includes a detection of a Kuiper Belt-analogue composition in a polluted WD (Xu et al. 2017). In addition, Johnson et al. (2022) observe a polluted WD with compositions that is composed from both rocky and icy bodies. They conclude that the unusual composition can be explained if the WD is polluted by two parent bodies, with a mix of Mercury-like composition and an icy Kuiper Belt-analogue. In the context our discussion, if Oort clouds are composed of both rocky asteroids and icy comets as mentioned previously, they can potentially explain all these diverse composition observations.

Brouwers, Bonsor & Malamud (2023) explore in details the accretion of comets. They find that accretion can occur in two stages: the ices may sublimate and accrete first, before refractory minerals can reach the star. Thus, the composition signature on a WD’s atmosphere may vary over time in a single accretion event, potentially also explain the composition diversity in polluted WDs’ atmospheres.

6.5 Accretion disc delivery

We only study mechanisms to deliver comets into the Roche radius. After a comet achieves a distance d ≤ 1 R_⊙, we assume that it is tidally disrupted and form an accretion disc around a WD (Koester et al. 2014). Poynting–Robertson drag can deliver materials from a accretion disc into WDs at the rate 10⁸ g s⁻¹ (Rafikov 2011). However, we observe pollution rate up to 10¹⁰ g s⁻¹ on WDs. It is difficult to explain material delivery at such rate from a accretion disc. Okuya et al. (2023) show that the Poynting-Robertson delivery rate can be enhanced in the existence of some volatiles, to bring the delivery rate above 10⁸ g s⁻¹. Relating to our discussion of volatiles, it may be necessary to have some icy bodies to enhance accretion disc delivery rate to explain some observed higher pollution rates.

6.6 Very close companions

So far, we only considered companions with semimajor axes on the order of a_p = 10 or 100 au. There are observational evidence for companions at these semimajor axes (e.g. Veras, McDonald & Makarov 2020; Blackman et al. 2021; Mullally et al. 2024). On the other hand, Vanderburg et al. (2020) present observational evidence for a M_p ∼ 1 M_Jup planet orbiting its WD at a_p ∼ 4 R_⊙. In addition, Gänsicke et al. (2019) and Veras (2020) present evidence for a planet at a_p ∼ 15 R_⊙ orbiting a volatile-rich polluted WD. These observations require us to analyse the scenario of a planet on a very small orbit.

First, we analyse if effects induced by this planet is important. For the planet found by Vanderburg et al. (2020), we have ζ ≈ 5 × 10³, assuming typical incoming comets a ∼ 10⁴ au. Since ζ ≫ 1, torque induced by this planet does not overcome galactic tide. In addition, since the planet is so close to the WD, the ejection loss cone does not significantly cover phase space more than the engulfment loss cone (equation 34). Thus, the ejection loss cone should not significantly reduce pollution rate. Therefore, the existence of a close-in planet with the configurations as found by Vanderburg et al. (2020) cannot reduce pollution rate of Oort cloud comets. We expect the pollution rate to be decreased by about a factor of 2 because comets are equally likely collide into either bodies. This is because we have two close central bodies, and the Safronov number for this planet is Θ ≈ 0.01 ≪ 1. Finally, the same conclusions apply to the system found by Gänsicke et al. (2019).

We discussed the case for close stellar companions earlier in Section 5.4, where we expect the pollution rate to be reduced by half. Note that for WDs with close stellar companions (with separations on the order R_⊙ to aus), the metal pollution on those WDs can also be explained by stellar winds from their companion (Zuckerman et al. 2003; Zuckerman 2014) and not necessarily by another reservoir like exo-Oort clouds.

We show using our analysis framework that close-in companions, either stellar or planetary, cannot significantly reduce Oort cloud comet delivery rate into a WD. The analysis can be applied for other systems with varying companion masses and separations; except for when the separations are large, a_p ≥ 300 au.

6.7 Widely separated companions

For planets or stellar-mass companions widely separated (a_p ≥ 300 au), our analysis cannot be applied, at least for Oort clouds with an inner semimajor axis at ∼3000 au like ours. As discussed earlier in various contexts, we typically assume that comets that encounter a companion have high eccentricity, e ∼ 1 (or equivalently, a ≫ a_p). Thus, we regularly expand expressions in this limit for simplifications. This assumption is also used in the analytic galactic tide loss cone theory.

Observationally, there is evidence for distant companions, both planetary and stellar mass, around WDs. Luhman, Burgasser & Bochanski (2011) provide observational evidence for a planetary companion with mass M_p = 7M_Jup at a separation of a_p ≈ 2500 au. Zuckerman (2014) presents a catalogue of 17 WDs with companions separated by more than 10³ au. In addition, the WD-planet system we mentioned earlier observed by Vanderburg et al. (2020) is in fact a hierarchical triple system. There is a distant star at ∼10³ au forming a triple system with the WD-planet binary.

Relating these observations to WD pollution, Wilson et al. (2019) find that the occurrence rate of single-polluted WDs and polluted WDs with wide stellar companions are the same. The fact that there are polluted WDs in wide binaries is in contrast with our prediction in Section 5.4. We expect that wide binaries should have significantly reduced pollution rate (below detection limit) due to the difficulty of direct collisions because of a high Safronov number. An implication of this could be that Oort clouds around wide stellar binaries are not similar to ours. In addition, the existence of a distant massive object, potentially embedded within the Oort cloud itself, could change those exo-Oort cloud structures significantly from our own Solar system Oort cloud.

One dynamical study involving wide binaries around WDs is performed by Bonsor & Veras (2015). They propose that due to galactic tidal effects, a distant binary companion periodically becomes excited to high eccentricity bringing it closer to the WD. During these close approaches, the stellar companion scatters other reservoirs, like an exo-asteroid belt or an exo-Kuiper Belt, into a WD and induces pollution.

With an abundance of companions at all mass scales with orbital separations a_p ≥ 10³ au and observations of polluted WDs in wide binaries, it remains interesting if and how these distant companions affect their Oort cloud and influence the observed WD pollution rate.

7 CONCLUSION

We have studied if an exo-Oort cloud can pollute a WD with galactic tide and a companion over a 1 Gyr time-scale through numerical and analytic methods. We analysed cases when there is galactic tide only, when the companion is a star, and when the companion is a planet. We studied the dynamics of the companion, namely precession torque induced on exo-Oort cloud comets and the kick in semimajor axis they experience close encounters with the companion. We present a fast integration method that is capable of integrating 10⁸ comets over a long simulation time.

The conclusions presented below assume a Solar system Oort cloud, with total cloud mass M_Oort = 2 M_⊕ containing N_Oort = 10¹¹ objects. These pollution rate results are scalable to other exo-Oort clouds with different masses and number of objects. Our main conclusions are:

• In the absence of any companions, exo-Oort clouds like our own Solar system’s can pollute WDs at a rate ∼5 × 10⁷ − 10⁸ g s⁻¹.

• We find that the dimensionless quantity

  $$\begin{eqnarray} \zeta \sim \frac{32 \pi \sqrt{2}}{3} \times \frac{\rho _0}{M_\mathrm{reduced}/a_p^3} \times \left(\frac{q}{a_p}\right)^{3/2} \times \left(\frac{a}{a_p}\right)^{7/2} \end{eqnarray}$$

  (45)

  is a good indicator of the relative importance of the angular momentum change induced by a companion compared to that of galactic tide (see Fig. 3). When ζ ≲ 1, torque from a companion dominates galactic tide and reduces comet migration. When ζ ≫ 1, only galactic tide is important.

• Stellar-mass companions with masses M_p ≥ 0.1 M_⊙ reduce the WD pollution rate to ∼3 × 10⁵ g s⁻¹ due to their strong angular momentum and scattering barriers, and a low likelihood of direct collisions with the WD.

• Planetary-mass companions significantly reduce WD pollution rate as predicted by O’Connor et al. (2023). However, we find the simulation pollution rate to be 2–5 times higher than predicted by O’Connor et al. (2023). We attribute this discrepancy to some limitations in the modified loss cone theory. Namely, the modified loss cone theory assumes a 100 per cent effective ejection barrier. The modified loss cone also cannot account for migration in comets semimajor axes, which occurs when comets interact with the planet in the diffusive scattering regime. That said, we still find that the modified loss cone is a reasonable order of magnitude estimate. With the existence of a planet and in the absence of stellar flybys, we find that the WD pollution rate due to Oort cloud objects is reduced to about ∼10⁷ g s⁻¹, assuming a Solar system Oort cloud.

• The power-law density profile structure of the Oort cloud does not significantly affect (by more than half an order of magnitude) the pollution rate.

• The pollution rate without a planet stays constant over a 1 Gyr time-scale. The pollution rate with a planet can decrease by a factor of ∼3 over a 1 Gyr time-scale.

We discussed the advantages of the Oort cloud and the impacts of stellar flybys. We also discussed two major concerns of using the Oort cloud as potential reservoir for WD pollution: (i) the retention of Oort cloud objects after a strong anisotropic natal kick experienced by WDs and (ii) observations of volatiles-poor WDs. There are uncertainties in our current understandings of the Solar system and extrasolar Oort cloud’s structure and composition. However, these uncertainties could potentially explain the diversity in WD pollution rate and composition. Finally, we applied our analysis framework in the context of observational evidence of close companions (orbits on the order of days) to predict their minimal effects on Oort cloud pollution, and stressed how we cannot apply this framework to widely separated companions (separation ≥10³ au).

We show that exo-Oort clouds can potentially pollute old WDs at observed rates over Gyr time-scales, depending on the existence of a companion and Oort cloud structures. However, one single reservoir and mechanism may not necessarily explain all instances of WD pollution (Veras, Mustill & Bonsor 2024). With further observations of polluted WDs and characterization of their companions (e.g. using Gaia as shown in Sanderson, Bonsor & Mustill 2022), our results can be applied to constrain sources of WD pollution.

ACKNOWLEDGEMENTS

We thank the reviewer Christopher O’Connor for helpful comments. We thank Scott Tremaine for a careful reading of the manuscript and valueable feedback. We thank Norm Murray, Sam Hadden, Dimitri Veras, Yanqin Wu, Garett Brown, and Michael Poon for insightful discussions. This research has been supported by the Natural Sciences and Engineering Research Council (NSERC) Discovery Grant RGPIN-2020-04513.

DATA AVAILABILITY

Data available on request.

Footnotes

1

Because of the way z is defined in our coordinate system, the inclination I and argument of pericentre ω are measured relative to the x − y plane parallel to the galactic mid-plane, instead of the usual measurement relative to the ecliptic (cf. Heisler & Tremaine 1986; Tremaine 2023). All I and ω used throughout this article follow this convention.

2

Gauß noticed that this is equivalent to the apsidal precession rate of a test particle induced by the quadrupole moment of a homogeneous ring with mass M_*M_p/(M_* + M_p) and radius a_p (see references in Touma, Tremaine & Kazandjian 2009).

3

Equation (25) in Hadden & Tremaine (2024) is re-written as the time-scale on which the comet experiences a strong scattering event due to diffusive kicks in semimajor axis.

4

λ is called Λ in O’Connor et al. (2023). We already used Λ in this work for the circular angular momentum.

APPENDIX A: FAST INTEGRATION METHOD WITH PLANET

One problem in direct N-body integration involves lengths spanning many orders of magnitude: typical incoming comets have semimajor axes on the order of 10⁴ au, interacting with a companion at 10^1–10² au, and engulfed by the WD Roche radius at 1 R_⊙. To accurately resolve these length-scales, we would like to have the integration time-step dt when the comet is far away to be large, and dt to be small when the comet can interact with the companion. That is, we would like dt to be adaptive in the following manners: on the order of the orbital period of the comet when far away, on the order of the orbital period of the companion when close in, and on the order of the pericentre passage time-scale when close to the WD.

However, this is an issue for a numerical integrator because we always need to resolve the smallest orbit. Thus, having an N = 3 bodies (WD–companion–comet) simulation will always cause dt to be a fraction of the orbital period of the companion, not the comet. This is not ideal since a lot of time spent during the simulation will be spent on resolve the Keplerian orbit of a two-body system.

To resolve this, we develop a method to simulate just the comet; the rebound simulation will have only N = 1 particle, as illustrated in Fig. 8. At each time-step, we manually set the total acceleration the comet experiences to

$$\begin{eqnarray} \ddot{x}_\mathrm{comet} &=& -G M_* \left(\frac{x_*}{r_*^3}\right) -G M_p \left(\frac{x_p}{r_p^3}\right), \\ \ddot{y}_\mathrm{comet} &=& -G M_* \left(\frac{y_*}{r_*^3}\right) -G M_p \left(\frac{y_p}{r_p^3}\right), \\ \ddot{z}_\mathrm{comet} &=& -G M_* \left(\frac{z_*}{r_*^3}\right) -G M_p \left(\frac{z_p}{r_p^3}\right)-4 \pi G \rho _0 z. \end{eqnarray}$$

(A1)

x_p, y_p, z_p, and r_p denote the distance between the comet and the companion. Similarly for x_*, y_*, etc. for the distance between comet and central star. Terms like x, y, z, and r can be quickly calculated analytically because the WD–companion system is a two-body system (comet is massless). Therefore, at any time, we can calculate the total acceleration a comet experiences without requiring to simulate the central star or the companion.

In addition, the last term in the z-component is from vertical galactic tide. Note that this is the full vertical galactic tidal term. This does not have the limitations of the orbit-averaged equations of motion, as mentioned earlier. Hence, we use rebound with the full vertical tide formulation to handle engulfments into the WD to properly capture all galactic tide dynamics.

We use the adaptive time-step method IAS15 (Rein & Spiegel 2015) as the integrator for our N = 1 simulation. The IAS15 integrator is capable of adaptively changing time-steps to resolve close encounters. In addition, we use the recent improvement on IAS15’s adaptive method (Pham, Rein & Spiegel 2024). We do not set a minimum time-step and allow IAS15 to properly resolve all close encounters at machine precision error in energy. When used with our fast integration method, the IAS15 integrator is able to adaptively resolve time-steps as wanted.

We show some verifications of this fast direct integration method in Appendices  B and  C.

APPENDIX B: VERIFICATION OF THE FAST N-BODY METHOD

To verify the fast direct integration scheme (Appendix  A) and our implementation of it in rebound, we plot Fig. B1. In this figure, we compare σ_x, the root-mean-square change in x = 1/a, between our fast integration method and the full three-body simulation. At each q, one thousand comets are initialized isotropically, with an initial semimajor axis at a = 1000 au and eccentricity e = 1 − q/a. The planet is at 5 au with mass 10⁻³ M_⊙. Comets are integrated over one orbital period to calculate σ_x. We can see that the fast method is able to reproduce well close encounters between comets and planet. The results in this plot can be compared with isotropic curve in fig. 9.3 in Tremaine (2023) for further verification. The fast method is 40 times faster. Finally, note that this fast integration method can be applied in multiple planet systems in other applications. We show an additional application of this with Kuiper Belt objects interactions with the outer Solar system in Appendix  C.

Figure B1.

Verification of the fast integration method by comparing σ_x, the root-mean-square change in x = 1/a, between two integration methods. ‘Normal’ (dashed line) means having three particles, central star, companion, and comet. ‘Fast’ (solid line) is the integration method presented in this article with only one particle: a comet. In the production of this plot, the fast method is faster by about 40 times. This plot can be compared with isotropic curve fig. 9.3 in Tremaine (2023).

APPENDIX C: AN APPLICATION OF THE FAST N-BODY METHOD

The fast rebound integration method can be extended to multiple planets. One major caveat in this case is that we assume that massive objects (planets) do not interact with each other. In other words, we assume the massive objects to be ‘on rails’. This allows us to analytically calculate the position of bodies, which is required for our fast integration method. That being said, this method is still a good approximation when the integration time-scale is less than the time-scales on which planet–planet interactions begin to become important (e.g. this would be the secular time-scale in the Solar system). It is possible to overcome this limitation by pre-running a simulation of the massive particles, record their positions, interpolate the positions at any t, and then we can use those interpolated positions for the fast integration method.

In the example here, we will use the simple version (no position interpolation). We showcase an integration of 1000 objects with semimajor axes initially distributed between [30, 50] au and inclinations drawn randomly from a Rayleigh distribution with 〈I〉 = 10°. We include the Sun and the outer Solar system planets. This integration is run for 100 Jupiter orbits. Fig. C1 shows the semimajor axis distribution of these objects. We show here that the fast method is capable of statistically reproducing the Kuiper Belt peaks induced by resonance.

Figure C1.

The Kuiper Belt semimajor axis distribution reproduced by the fast integration method compared with the normal integration method. The simulation includes effects from the outer Solar system planets (and the Sun). The simulation is run for 100 Jupiter orbits.

APPENDIX D: SEMIMAJOR AXIS REJECTION SAMPLING

The distribution of comet semimajor axis itself follows a density profile described by a power-law relationship n(a) ∝ a^−γ. We wish to study all γ between 2 ≤ γ ≤ 4. We could draw samples at fixed γ and run simulation and repeat that over a grid γ. However, this increases the computation costs significantly. For example, to make Fig. 9 would require running 20 simulations at different γ. To resolve this, we use a rejection sampling scheme.

We initialize comets’ semimajor axes on a broken power-law distribution (cf. equation 17):

$$\begin{eqnarray} {\rm {d}} N(a) \propto \left\lbrace \begin{array}{@{}l@{\quad }l@{}}a^{-\gamma _\mathrm{max}} \cdot a^2 ~{\rm {d}} a, & a_1 \le a \le a_\mathrm{turnover}\\ a^{-\gamma _\mathrm{min}} \cdot a^2 ~{\rm {d}} a, & a_\mathrm{turnover} \le a \le a_2, \end{array}\right. \end{eqnarray}$$

(D1)

where γ_min = 2 is the minimum γ that we can sample and γ_max = 4 is the maximum. a_turnover = 17 000 au is the broken power-law turnover point. a_turnover is determined heuristically to maximize the efficiency of rejection sampling.

We simulate |$N_\mathrm{comets}^\mathrm{sim}$| comets with semimajor axes drawn from this broken power-law distribution. We record their initial and final positions. Using rejection sampling, we select comets after running the simulations to construct results for any particular γ we are interested in.

In Fig. D1, we show the efficiency of rejection sampling from the broken power law at various γ. We also compare that with rejection sampling from a uniform distribution instead. As shown, the broken power law is much efficient for rejection sampling into a power-law distribution than a uniform distribution. We find that using the broken power law as proposed gives a 50 per cent efficiency. That is, to construct results at any particular γ, the number of comets representative in that results is |$N_\mathrm{comets}^\mathrm{eff} \sim 0.5 N_\mathrm{comets}^\mathrm{sim}$|⁠. As a consequence, we need to simulate twice the number of comets we would like to see for any particular result. However, this is still less than one simulation per γ.

Figure D1.

The efficiency of performing rejection sampling on a broken power law versus rejection sampling on a uniform distribution. We can see that the using the broken power-law distribution to create power-law samples is much more effective. With the broken power-law rejection sampling, the number of usable samples at any particular γ is about half of the total simulated samples.