Polluting White Dwarfs with Oort Cloud Comets

Observations point to old white dwarfs (WDs) accreting metals at a relatively constant rate over 8~Gyrs. 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 characterise the dynamics of comets around a WD with a companion having semi-major axes on the orders of 10 - 100 AU. We develop simulation techniques capable of integrating a large number ($10^8$) of objects over a 1 Gyr timescale. 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 characterise 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 timescale, (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 $\sim 10^8 \mathrm{~g~s^{-1}}$, (iii) adding a planetary-mass companion reduces the pollution rate to $\sim 10^7 \mathrm{~g~s^{-1}}$, and (iv) if the companion is stellar-mass, with $M_p \gtrsim 0.1 M_\odot$, the pollution rate reduces to $\sim 3 \times 10^5 \mathrm{~g~s^{-1}}$ 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.


INTRODUCTION
Between 25 to 50 percent of white dwarf (WD) atmospheres observed are polluted with heavy metals (e.g.Zuckerman et al. 2003Zuckerman et al. , 2010;;Koester et al. 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 destabilisation in these sources (e.g.Debes & Sigurdsson 2002;Debes et al. 2012;Mustill et al. 2014;Smallwood et al. 2018;Maldonado et al. 2020;Trierweiler et al. 2022;O'Connor et al. 2022).The bodies must be delivered to the WD Roche radius at ∼ 1 ⊙ 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 5 to 10 10 g s −1 ; the mean pollution rate is around 10 7 g s −1 with about a 1 dex spread (e.g.Blouin & Xu 2022;Johnson et al. 2022;Cunningham et al. 2022).The current minimum detection limit is about a few 10 5 g s −1 (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 timescales.
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 ★ E-mail: dang.pham@astro.utoronto.ca(DP) 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 et al. 2023).There are observations of polluted WD atmospheres with volatiles, although they are much more rare (Farihi et al. 2013;Klein et al. 2021;Doyle 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 semi-major axes via interactions with a giant planet (Vokrouhlický et al. 2019;Kaib & Volk 2022).For objects kicked into higher semi-major axes, they can be circularised by galactic tide or stellar flybys (Duncan et al. 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 11 − 10 12 objects with a total mass of ∼ 2 ⊕ with semi-major axes ranging between 3 000 AU to 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 et al. 1986;Parriott & Alcock 1998;Veras et al. 2014;Stone et al. 2015;Caiazzo & Heyl 2017).Simulations of comets into a 1 ⊙ 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 et al. (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 5 g s −1 , 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 timescale?We also investigate this question in cases where the WD has a planetary-mass or a stellarmass 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 powerlaw 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 summarise 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 summarise 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 8 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 timescale.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 summarise our findings.

Notations
In this work, we denote  * as the WD mass.We denote the orbital elements of a comet as:  for the semi-major axis,  for eccentricity,  for inclination,  for argument of pericentre, Ω for longitude of ascending node, and  for mean anomaly.Orbital elements  and  are measured relative to galactic plane.Some other quantities used to describe the comet orbits are:  for the pericentre distance,  for the apocentre distance,  for the orbital period.The comet is assumed to be a mass-less test particle.Orbital elements with the subscript  denote that they are the orbital elements for the WD companion (e.g. is the companion mass).We also regularly use the following set of Delaunay action-angle variables (quantities Λ, ,   are actions in units of specific angular momentum, and , , Ω are angles): Λ is referred to as the circular angular momentum,  as the angular momentum, and   as the  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 semi-major axis edge,  1 , ends at  2 , follows a powerlaw number density profile () ∝  − , has  Oort comets and a total cloud mass  Oort .
When referring to our own Solar System's Oort cloud, we state 'Solar System Oort cloud' explicitly.In the Solar System, estimates for  Oort typically range between 10 11 − 10 12 comets (e.g.Francis 2005; Boe et al. 2019).We assume  Oort = 10 11 in this article.It is estimated that  Oort ∼ 2 ⊕ 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 powerlaw 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  1 ,  2 , ,  Oort ,  Oort exactly like our own Solar System's Oort cloud.

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  direction. is defined as the direction orthogonal to the galactic midplane 1 , with the midplane defined at  = 0.In the Solar System, the second assumption is valid since galactic tidal terms in the  and  components are about one order of magnitude lower than the  term.Since observations of polluted WDs are typically within the Solar Neighbourhood, 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  0 . 0 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 midplane.As a result,  0 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 = 0.1 ⊙ pc −3 (Holmberg & Flynn 2000;McKee et al. 2015).
The galactic potential with these assumptions can be written as 1 Because of the way  is defined in our coordinate system, the inclination  and argument of pericentre  are measured relative to the  −  plane parallel to the galactic midplane, instead of the usual measurement relative to the ecliptic (c.f.Heisler & Tremaine 1986;Tremaine 2023).All  and  used throughout this article follow this convention.(Heisler & Tremaine 1986): (2) This potential can be averaged over the orbit of the comet and then written in Delaunay elements as: from which we find the secular (orbit-averaged) equations of motion by applying Hamilton's equations: GT and  GT are coupled differential equations.The phase space evolution described by these coupled equations is studied in detail in Heisler & Tremaine (1986); Tremaine (2023).Briefly, they show that through these equations of motion, the comet can oscillate in the (, ) 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 semi-major axis is conserved under the galactic tidal effect.Similarly,   is conserved because Ω is not in ⟨Φ GT ⟩.Since   is conserved but  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;Lidov 1962;Kozai 1962).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.

Loss Cone Theory and Comet Injection Rate
First, we describe Oort cloud comets by a distribution function  , defined as the number of comets per volume of phase space: d Oort =  (Λ, ,   , , Ω, ) dΛ d d  d dΩ d. (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).
d Oort =  (Λ, )d dΛ where ∫ d Oort =  Oort with  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  =  crit = 1 ⊙ , 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,  ≲ 1, the angular momentum can be related to the pericentre distance  as: Through this, we define a critical angular momentum once a comet has achieved a certain critical pericentre distance  crit : The loss cone is defined as the phase space region where  ⩽  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 semi-major axis  (or equivalently, Λ).Intuitively the dependency on  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 Δ over multiple orbits, slowly migrating a comet inward in  over many orbits.This is the empty loss cone case.In the other case the galactic tide is sufficiently strong to induce a large Δ and the comet is capable of reaching the loss cone within one orbit.This is the filled loss cone case.Figure 1  First, we consider the empty loss cone case, which happens when Δ ⩽  crit , where Δ ∼ |d/d| ×  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, Γ  (dimension of number of comets per unit Λ per unit time) is: where  crit ≪ Λ (highly eccentric orbit).Γ  is found by using Equation 4 and calculate the rate at which comets are pushed into the loss cone boundary at  crit .
Next, we consider the filled loss cone case when Δ ⩾  crit .The injection rate for the filled loss cone, Γ  (same dimension as Γ  ) is: where Γ  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  here denotes 'full', not the distribution function  .The loss cone is empty at small  and full at large .The transition between the two cases can be found by equating the two rates, Γ  = Γ  , yielding: When  <  crit , a comets is in the empty loss cone regime.When  ⩾  crit , it is in the full loss cone regime.In Delaunay variables, these conditions are equivalent to Λ < Λ crit and Λ ⩾ Λ crit with The total injection rate (number of comets entering a certain  crit per unit time) can be found by adding up the two injection rates, integrated over Λ: where Λ 1 = (  *  1 ) 1/2 and Λ 2 = (  *  2 ) 1/2 , with  1 the inner semi-major axis edge of the Oort cloud and  2 the outer edge.Finally, since Γ  ∝  2  (), the injection rate per unit Λ increases from  1 to  crit .On the other hand, Γ  ∝  −3/2  (), the injection rate per unit Λ decreases from  crit to  2 .At  crit , Γ  = Γ  .Since the injection rate increases until  crit and then decreases, the majority of the total injection rate is contributed from the region around  crit .

Oort cloud Distribution Function
Now, we find the explicit form of the distribution function  .First, for a dynamically relaxed distribution, the distribution is 'thermal' and the distribution in  2 is uniform (Jeans 1919;Ambartsumian 1937).This 'thermal' distribution of Oort cloud comets is seen after long term simulations of Solar System Oort cloud (c.f. Figure 7 in Vokrouhlický et al. 2019).In this case,  =  (Λ) is independent of .We can further simplify by integrating over By definition, for a spherically distributed shell with width [,  + d], we also have: where, () is the number density profile of comets.Before finding  (Λ) through these equations, we first choose, ().Long term simulations of the Solar System Oort cloud find that our own Oort cloud generally follow a powerlaw density profile (e.g.Duncan et al. 1987;Higuchi & Kokubo 2015;Vokrouhlický et al. 2019): where  ∈ [ where Λ 1 and Λ 2 are the circular angular momenta associated with  1 and  2 . is the normalisation constant: Equipped with the distribution function  (Λ), it is possible to find the total comet injection rate Γ total (dimension of number of comets per unit time) by substituting  (Λ) and integrating Equation 15.With Γ total , we can also predict the pollution rate (i.e. in g s −1 ) into a WD:

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  ⩽ 10 −2  ⊙ = 10 Jup , or stellar if it has mass  ⩾ 10 −1  ⊙ .In section 3.1, we analyse the effects of companion-induced precession on a comet with semi-major axis  and pericentre  from a companion with mass   on a circular orbit at semi-major axis   .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.

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, , induced by galactic tide and companion to study when each effect is dominant.In our case, we consider the regime where  ⩾   -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) 2 : where   is the companion's mean motion and Δ 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   (the next non-zero term occurs at the hexadecapole order, see Palacián & Yanguas 2006;Vinson & Chiang 2018), which become important as  →   .However, for our analysis in this section, the quadrupolar term is sufficient to give an estimate.When we require numerical results for   (Figure 6),   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: Precession in argument of pericentre (Δ) accompanies angular momentum change (Δ).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 .
Expanding  to first order in / (the high eccentricity limit,  ∼ 1) and taking the order of magnitude terms, we find: where we defined two dimensionless quantities  and . is the ratio of densities between the galaxy and the binary system: where  reduced =  *   /( * +   ) is the reduced mass.As  becomes smaller,  0 ≪  reduced / 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: Higher  means that comet experiences less torque from the companion. is mostly dominated by /  .When a comet has a large orbit comparing to the companion ( ≫   ), 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  (or equivalently, on the eccentricity ).As a comet migrates inward due to galactic tide, the orbit becomes more eccentricity,  increases,  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  = 5  .This is chosen because here, scattering is generally weaker comparing to both precession and galactic tide (see Figure 6).A comet semi-major axis of 10 4 AU is chosen because this is typical for incoming comets.In both tables, we give sample values for a close-in companion at   = 10 AU and a wider companion at   = 100 − 200 AU.
Figure 2 shows  on a grid of  and .The dashed lines show the contours where  = 0.1, 1, 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.
For close companions at   = 10 AU, we have  > 10 11 .In this case, to have precession dominates tide by making  ≲ 1, we require  < 10 −13 , which can be supplied by a stellar-mass companion with mass   ⩾ 0.1 ⊙ as shown in Table 1.However, a planetary-mass companion does not provide sufficient torque to overcome galactic tide.For more distant companions at   = 100 − 200 AU, the same conclusions can be reached through their values of .Therefore, a stellar-mass companion (  ⩾ 0.1 ⊙ ) is required to produce a precession barrier, which reduces the efficiency of galactic tide delivering comets from an exo-Oort cloud.

Limitations
We note that there are four limitations when using : (i)   as shown in Equation 22is only valid for  ⩾   ; that is, a comet must be strictly exterior to companion.
(ii) We only use the quadrupole term for   .As  approaches   , contributions from higher order terms become important.Thus,  should not be used when  ≈   .
(iii) We ignore angular dependencies in , , Ω.However, as seen in Figure 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.
(iv) We expand  in the limit  ≪ .In other words,  is not a good approximation for very circular orbits ( ≈ ), or for very close-in comets encountering very widely-separated companions ( is small, but   is large, so  is comparable to ).For the first case, when  ≈ , 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   ≪ .The closest Oort cloud comets in our model have  ≈ 3000 AU, so at most, we should only consider   ⩽ 300 AU for  analysis.

Precession Barrier Efficiency
Figure 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: where  GT is the number of comets that can enter a certain pericentre  with only galactic tide. 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  GT and  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   ,  * , ,   , .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 (, , Ω) are drawn randomly assuming comets are isotropically distributed.
We compare the efficiency over a grid of (, ) to the analytic values of  in Figure 3. Above the  = 10 contour, companion-induced precession does not suppress any comets experiencing galactic tide.
Here, 100% of comets can come into  due to galactic tide.As  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 colour shows the efficiency of the planet's torque at inhibiting galactic tide, where  companion is the number of comets that can enter a certain  in the existence of a companion and galactic tide and  GT is the same but there is only galactic tide (no companion).Dashed lines show contours for  = 0.1, 1, 10.
Above the  = 10 dashed line, the precession barrier efficiency is 0% and comets fully experience galactic tide and can be excited to high eccentricity, migrating inward in pericentre .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%.
decreases, the precession barrier becomes stronger and eventually at the bottom left corner, all comets are suppressed from galactic tide, preventing inward migration in .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% 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%.

Scattering Timescale
Comets not only experience a torque in Δ causing precession, but also experience a kick in energy causing a change in semi-major axis.There are two regimes where comets experience semi-major axis kicks from a massive companion: strong and diffusive scattering.In the strong scattering regime,  ≳   , comets are ejected from the system within one pericentre passage.Here, the scattering timescale is approximately  scattering ≈  comet /2; comets are kicked during their incoming pericentre passages.
In the diffusive regime,  >   , comets experience small random semi-major 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 timescale on which small diffusive kicks in semi-major axis will lead to an order unity change in binding energy (or equivalently, order unity change in semi-major axis) 3 : There are three main assumptions to using the results found by Hadden & Tremaine (2024) diffusive timescale 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,   ≪  * , the companion's mass is much smaller than the central star's mass.Hence, we cannot use this timescale to estimate the scattering timescale in the stellar-mass companion case.Third, they assume the comet is on a parabolic orbit ( = 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  = 1 diffusive timescale as written here, for comets with  ≲ 1.Finally, since we assume  ≲ 1, this timescale is only applicable in cases where  ≪  and   ≪  (the companion is much closer-in than the comet).
Next, we compare analytic expectations of  scattering to numerical simulations.We simulate two cases: a stellar-mass companion with   = 0.6 ⊙ and planetary-mass companion with   = 10 −2  ⊙ .In both cases, the companions are on a circular orbit at   = 200 AU around a central WD with mass  * = 0.6 ⊙ and the companion's initial phase is randomised.At each pericentre , where  ∈ [1, 6]  , 200 comets are initialised with  = 15 000 AU (typical semi-major axis of incoming comets to interact with companions at   = 200 AU).The comets' initial position is set at apocentre.Comets' angles are initialised either coplanar and isotropically.Simulations are stopped when comets experience a kick Δ/ 0 > 0.3 relative to the initial semi-major axis  0 , or until the simulation time reaches 10 11 years and we call this time the numerical scattering timescale.The numerical condition, Δ/ 0 > 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).
Figure 4 compares numerical scattering timescales with analytic expectations.In both panels, numerical timescales for coplanar and isotropic comets are shown.The solid coloured lines are the mean scattering timescales and the shaded areas are the timescale ranges from 200 comets.The top panel shows the timescales for a stellarmass companion and the bottom panel for a planetary-mass companion.First, there are clearly two scattering regimes.Comets sufficiently far away experience long scattering timescale according to the diffusive timescale.As the pericentre decreases, the timescale eventually converges to the strong scattering regime where comets 3 Equation 25 in Hadden & Tremaine (2024) is re-written as the timescale on which the comet experiences a strong scattering event due to diffusive kicks in semi-major axis.M p = 0.6M are ejected within one pericentre passage.Second, in the bottom panel, we find the timescale 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  scattering is consistently off by a factor of a few.We attribute this to the arbitrary Δ/ 0 > 0.3 numerical scattering condition.Third, isotropic comets have a higher scattering timescale.In addition, as the pericentre increases, coplanar and isotropic timescales 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 timescale.This is because the condition   ≪  * is strongly violated.Here, we observe that coplanar comets quickly becoming strongly scattered at  ≈ 6  , 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 4  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.
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 .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.

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 summarise their theory below and illustrate the modified loss cone theory in Figure 5.We also discuss the limitations of this modified loss cone framework.A similar effect is proposed by Teboul et al. (2024) for dense star clusters with a central black hole called "loss cone shielding".
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 semi-major axis which can eventually eject the comet.The location of the ejection loss cone is defined at: nesting on top of the tidal-disruption loss cone at where  total =  * +   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,   ≫  crit , we have  ej ≫  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 Δ due to galactic tide, it drifts inward encountering these loss cones.A comet with relatively small  experiences a small Δ ≪  ej , and drifts slowly to the edge of the ejection loss cone; the ejection loss cone is empty.Vice versa, a comet with larger  experiences a much greater Δ ≳  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, Δ ≳  crit , it can jump through both loss cones and can be tidally disrupted and pollute a WD.These regimes are illustrated in Figure 5 (following Figure 12 in O' Connor et al. 2023).The comet avoids ejection since it reaches the engulfment loss cone in one orbital period.The transition semi-major axis between the two Δ regimes is: where in the last expression we used values of  crit = 10 500 AU (Equation 13),  crit = 1 ⊙ and   = 10 AU.This transition semimajor axis corresponds to a critical circular momentum Under this framework, only comets with  >  crit,ej can experience a sufficiently large Δ ≳  ej to pollute a WD.In addition, comets that can pollute WDs are in the filled loss cone regime because  crit,ej >  crit .Therefore, the pollution rate is estimated as the filled loss cone rate integrated over semi-major axes ranging between  crit,ej ⩽  ⩽  2 , or equivalently over Λ: O' Connor et al. (2023) further simplify this expression to find: Taking a fiducial  = 3.5 (for a Solar System Oort cloud), the factor on the right hand side is approximately 0.05 for  crit = 1 ⊙ and   = 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.

Limitations
First, comets with Δ ≪  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 semi-major axis by a much smaller mass planet.A planet the size of the Earth and a planet with 10 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  <  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% effective scattering barrier, reducing pollution rate as described by Equation 34.O'Connor et al.
(2023) recognised 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 quantity4 : 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 10 Jup planet at   = 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 (Vinson & Chiang 2018;Farago & Laskar 2010).
Third, as discussed earlier, another important additional dynamics is that a comet experiences small random kicks in semi-major axis, Δ, through multiple diffusive encounters.This causes the semimajor axis to change over time which also changes how a comet experiences galactic tide Δ over time.Since Δ due to galactic tide is strongly dependent on a comet's semi-major axis, these random walks in  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 semi-major 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  <  crit,ej .Through galactic tide, the comet migrates inward in pericentre distance .At some point later in time when the comet achieves a  ≳   , the comet begins to experience random kicks in semi-major axis every pericentre passages through interactions with the planet.Due to multiple small Δ kicks, the comet is migrated to  ≳  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 , where galactic tide induces a stronger Δ 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  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 , bringing comets into a stronger Δ regime capable of bypassing the ejection barrier.As shown in Figure 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 .
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.

Timescales Comparison
To analyse the importance of the companion's scattering and precession versus galactic tide, we investigate the timescales on which these effects are important.
These timescales are measured from simulations of 10 5 comets from a Solar System-like Oort cloud: () ∝  −3.5 from 3000 to 10 5 AU.Here, we directly simulate all comets with an additional galactic tidal force at every timestep -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 .In addition, we measure independently the contribution of galactic tide and companion at a given pericentre .
First, we simulate 10 5 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 .The orbital elements of incoming comets with pericentre  are then saved.Next, we measure the galactic tidal torque over one orbit at this , giving us  GT ≈ Δ/ 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   ≈ Δ/ comet .Third, these comets are again run in a third simulation with a companion and no galactic tide to measure the scattering timescale.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, Δ/ > 0.3, or until the simulation time reaches 10 12 years.Recording the times when comets experience a scattering event gives us  scattering .These simulations are run for every  between  =   to  = 5  and for companions with masses   = 10 −2 , 0.1, 0.6 ⊙ (corresponding to a WD, M dwarf, and the upper-limit of a planetary mass companion, 10 −2  ⊙ = 10 Jup ).
We compare the rates (inverse timescales) of these effects in Figure 6.The shaded area is the 3 spread, showing the range of most comets' timescales.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 −2  ⊙ planet companion, scattering begins to dominate at  ≈ 3  , whereas this is increased to  ≈ 3.5  for a 0.1 ⊙ stellar companion.For the WD-WD binary case, scattering seems to be dominated in all  considered here.Note that the scattering timescale flattens out as  →   .This is similar to the scattering behaviour we analysed earlier: comets are ejected within one orbital period in this regime,  scattering ≈  comet /2.Similar to scattering, precession's strength becomes stronger relative to galactic tide as companion mass increases.Finally, Figure 6 shows that in all cases, as  ≃   , the dominant effect is companion-induced scattering.

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  =  switch , its orbital elements are extracted and the comet is then integrated directly with REBOUND (Rein & Liu 2012) when  ⩽  switch .When comets violate their boundary conditions (section 4.2) they are removed.Figure 7 illustrates the components of our integration.

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, 5 and 6) to quickly evolve comets.We use the fourth-order Runge-Kutta to evolve these coupled differential equations with an adaptive timestep scheme as implemented in scipy (Virtanen et al. 2020).We find that the Runge-Kutta adaptive timestep scheme reproduces well the evolution in (, ) as seen in Heisler & Tremaine (1986).
If a comet can reach a certain critical  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 Figure 8 and described in details in Appendix A.
One caveat of using these orbit-averaged equations of motion is that they are not appropriate when Δ ≫  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  switch ≈ 4  , interactions with a planetary mass companion become important.We show this earlier in Figure 6 and its discussion.At  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 semi-major axis   with   = 0.
When there is no planet, we still extract comets from the secular integration when the comets reach  switch = 10 AU.This choice is somewhat arbitrary, but it ensures that  switch ≫  crit = 1 ⊙ .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 Δ is sufficiently strong to inject the comet into the loss cone in one orbital period.

Boundary Condition
During the REBOUND part of our simulation, we enforce the following boundary conditions: (i) comets with apocentre exceeding the Hill sphere of a 0.6 ⊙ star ( > 0.8 pc) and are outbound away from the WD, (ii) comets are ejected ( < 0) and are outbound away from the WD, (iii) comets are engulfed ( ⩽ 1 ⊙ , where  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 timestep.In addition, the simulation is stopped when the simulation time reaches  = 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  ⊙ WD, the Hill radius is at 0.8 pc (O'Connor et al. 2023).In addition, boundary conditions (i) and (iii) imply that there is a maximum semi-major axis: This is the outer semi-major axis edge of the Oort cloud.Comets beyond  2 cannot be excited to high eccentricity since they will be removed for exceeding the apocentre limit.The lowest pericentre comets at this semi-major axis without exceeding the WD Hill sphere is  ∼ 5 000 AU.In addition, this upper semi-major axis limit can also be found by scaling down the Solar System's Oort cloud outer semi-major axis edge at  2 = 10 5 AU, for a  * = 0.6 ⊙ central star.
Centre-of-mass (WD and companion)

Comet
(1) Galactic tide only: Integrate comet using galactic tide equations of motion up to 1 Gyr. Companion (2) Extract comet and go to step (3) if pericentre distance reaches  = 4   .
(3) Direct integration with REBOUND: Includes galactic tide, a central star and a companion.
Extract comets when comet is engulfed at  = 1  ⊙ (Roche radius) Figure 7.A diagram to illustrate our hybrid integration scheme.The first stage is using secular galactic tide equations of motion in (Equations 4, 5 and 6) to quickly integrate a comet.If a comet can be excited to  =  switch = 4  , 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 ( > 0.8 pc), when it's ejected due to interactions with a companion ( < 0) or when it is engulfed by the WD ( ⩽ 1 ⊙ , where  is the distance between the comet and the WD).Integrations are stopped when the simulation time reaches  = 1 Gyr.

Companion Comet
Fast integration with REBOUND Only 1 particle (the comet) in simulation Force terms are analytically calculated and added to the comet's acceleration at each timestep A diagram illustrating our fast, properly time-adaptive integration method in REBOUND.There is only one particle in the simulation: the comet (red).
The positions of the WD and its companions (shown in grey) can be calculated analytically because they are in a 2-body system.After each simulation timestep, the forces are analytically computed and added to the acceleration of the comet particle.IAS15 is used as an adaptive timestep and is capable of resolve timesteps as: large timestep when the comet is far away (the timestep is a fraction of the comet's orbital period), and small timestep 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.

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 U [0, 2] where U is the uniform distribution.Comets' inclination are drawn according to cos  ∼ U [0, 1].The sign of the comet's  is not important since we set the planet at   = 0. Next, the semi-major axis is drawn such that Oort cloud comets follow a powerlaw density profile () ∝  − : d () ∝  −  2 d (38) with  ∈ [ 1 ,  2 ].  1 = 3 000 AU is the inner semi-major axis edge, set based on the Solar System's Oort cloud simulations (e.g.Duncan et al. 1987;Vokrouhlický et al. 2019). 2 = 85 000 AU is the outer semi-major 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 semi-major axis distribution.This is described in Appendix D.
After having , we draw the squared eccentricity,  2 , from: 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  (or equivalently, a minimum initial pericentre  initial,min ). 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 In the case where there is only galactic tide and a central WD, it is somewhat more arbitrary: This is to ensure that  initial,min ≫  crit = 1 ⊙ .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.

Galactic Tide Only
We perform numerical simulations to explore the pollution rate of comets into a WD's Roche limit (1 ⊙ ) over various Oort cloud powerlaw structure ().Initial conditions and sampling are done according to Section 4.3, allowing us to sample a variety of powerlaw exponents .Here, we simulate with  sim comets = 4 × 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 Figure 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 powerlaw exponent value; the number density of Oort cloud objects scale as () ∝  − .The rate on the left is the number of accreted comets per Myr.On the right axis is the rate in g s −1 , assuming a fiducial Solar System Oort cloud mass and 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00  15) is in red.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  Oort = 10 11 comets and  Oort = 2 ⊕ .The blue shaded area shows the Poisson 1 confidence interval of the solid blue line (mean rate).
number of comets with  Oort = 10 11 comets and a total cloud mass of  Oort = 2 ⊕ .To get this rate, we count the total number of comets accreted,  accreted , after the warm-up phase which is about 400 Myr (subsection 5.5).Dividing  accreted by the remaining 600 Myr of simulations yields the number of comets accreted over time as shown on the left hand side of Figure 9. 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   ≈ 5 × 10 7 − 10 8 g • s −1 , depending on .Over the course of 1 Gyr, these rates correspond to the delivery ∼ 5 × 10 −4  ⊕ 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,  accreted , the Poisson process assumption allows us to estimate the standard deviation to be  =

√
accreted .Recall that  accreted is the total number of accreted comets counted over 600 Myr.Thus, √  accreted is the uncertainty of comets entering 1 ⊙ due to our limited number of comets in our simulation (10 7 comets) over a 600 Myr timescale.In addition, we simulate a sufficiently large number of comets such that  accreted is not a small integer and the 1 interval can be meaningfully interpreted.Finally, √  accreted is based on the total number of comets in our simulation (∼ 10 7 comets).If our simulation contained 10 11 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,  accreted , from our simulation containing ∼ 10 7 comets, over a 600 Myr timescale.

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.Figure 10 shows the efficiency (Equation 27) over different  with varying companion mass:   = 0.6, 10 −1 , 10 −2  ⊙ .In contrast to Figure 3 previously, we now measure over a range of .Hence, this efficiency can include the effects of both the precession and scattering barriers produced by a companion.Especially as  →   , both precession and scattering barriers can become very important in increasing the efficiency.Efficiency is 0% when galactic tide is the only dominant effect.Vice-versa, the efficiency is 100% when precession and scattering barriers induced by a companion can suppress all galactic tidal effects.First, we use the formulations of  to predict if the precession barrier is stronger than galactic tide at  = 5  .We choose to focus at  = 5  because scattering is not important there for   = 10 −2 and 0.1 ⊙ (c.f.timescales in Figure 6).This can also be done at other  provided those  are within the limitations of our formulations, and scattering is not important.The three companion cases in Figure 10 correspond to  = 10 −10 , 10 −9 , 10 −8 .From simulations, typical incoming comets into a pericentre  = 5  = 1000 AU have semimajor axes  ∼ 15, 000 AU, giving a  ≈ 5 × 10 7 .For  = 10 −8 ,  = 25 ≫ 1 so we expect galactic tide to be dominant.For  = 10 −9 ,  ≈ 1 so we expect companion-induced angular momentum change to be important and suppress some comets' galactic tidal torque.For  = 10 −10 , the precession-barrier is dominant over galactic tide.Figure 6 at  = 5  confirms these expectations.
Second, we can predict the efficiency of the precession barrier in reducing comet engulfment by using Figure 3 with values of (, ).At  ≈ 5 × 10 7 , Figure 3 predicts about a 0% efficiency for  = 10 −8 .At  = 10 −9 , it is expected to be around 40% efficient.At  = 10 −10 , we expect about a 60% efficiency.We confirm these predictions with Figure 10 at  = 5  .
Third, scattering becomes important at  = 3 − 5  as seen in the timescale analysis in Figure 6.In Figure 10, this corresponds to the fast increase in efficiency at that  range.In both cases, a stellar-mass companion increases the efficiency by almost 100% by  = 1  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  15) is shown as a solid black line 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 red dashed line.The blue shaded area shows the Poisson 1 confidence interval of the solid blue line (mean rate).
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  ∼ 1 ⊙ will be engulfed by the WD.
Finally, for a planetary-mass (  ⩽ 10 Jup ) companion, planetinduced precession does not play a strong role.At  = 5  , the efficiency is still roughly 0%.Thus, effects from the planet at the  distance is not important.There is a quick increase in efficiency beginning at  ≈ 2  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 the next subsection.

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 semi-major 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  sim comets = 10 8 comets.
In Figure 11, we present the simulated WD pollution rate at various  in the presence of a   = 10 Jup = 10 −2  ⊙ planet at   = 10 AU.These planet mass and semi-major 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  crit = 1 ⊙ by about 1 order of magnitude.Assuming a Solar System Oort cloud with  Oort = 10 11 comets and a total cloud mass of  Oort = 2 ⊕ , the WD pollution rate in this case is  (planet)  ∼ 10 7 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 .
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  crit,ej .In the  = 4 limit, comets are more centrally distributed.Thus, more comets have  <  crit,ej .Recall that these are the comets that experience small Δ <  ej due to galactic tide and migrate slowly in  until they are ejected through encounters with the planet.Therefore, since most comets are centrally distributed, they have  <  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  <  crit,ej but not as many as in the case of  = 4. Therefore, we have more comets with  >  crit,ej .These comets are capable of experiencing a strong Δ ≳  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% effective.If the barrier is 100% effective, expect all comets with  <  crit,ej to be all ejected.However, because the ejection loss cone is not 100% 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  <  crit,ej than at  = 2, the assumption of having a 100% 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 semi-major axes,  initial in Figure 12.The initial semi-major axis is shown on the x-axis because comets experience kicks in  over time.
First, in the regime where  initial <  crit,ej , the planet significantly reduces pollution rate as predicted.However, we still see comets with  initial <  crit,ej polluting the WD.We find that the ejection loss cone barrier is not 100% effective.Therefore, the assumption in Equation 33that the pollution rate does not have any contribution with comets from  <  crit,ej gives an underestimate of the pollution rate.
Second, in the regime where  initial ⩽  crit,ej , the number of comets capable of polluting the WD increases steadily.As  initial →  crit,ej , comets need a smaller kick in Δ to deliver them over into the  ⩾  crit,ej regime where galactic tide can induce a strong enough Δ to drift them through the ejection loss cone.On the other hand, in the  initial ≪  crit,ej , comets will need to experience multiple interactions with the planet that increases their semi-major axis, but not strong enough to eject them.
Third, as  increases beyond  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 , allowing stronger Δ to migrate further in, as described in Section 3.4.1.Second, these comets have lower energy (because of large ) 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  = 50 000 AU for example, will stay there and if it cannot reach  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  initial , a lot of these additional dynamics are ignored giving an incorrect rate in that regime.As we have seen, the contribution in region  >  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 semi-major axes  initial ⩽  crit,ej as predicted by O'Connor et al. (2023).However, we find that this reduction is not 100% effective, that comets experience rich dynamics, and that beyond  initial >  crit,ej the pollution rate does not simply follow the full loss cone rate Γ f .That being said, when comparing the overall rates in Figure 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  (planet)  ∼ 10 7 g s −1 , assuming a Solar System Oort cloud.

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   on circular orbits.The Safronov number can be used to estimate how probable it is to 2.00 2.25 2.50 2.75 3.00 3.25 3.50 3.75 4.00 expect collisions of incoming comets on highly eccentric orbits: where  esc is the escape velocity from the companion's surface (at distance ) and   is the circular speed at the companion's semimajor axis.We set  = 1 ⊙ for the tidal radius where comets are captured by a WD and   = 10 AU.Since Θ ≫ 1, we predict that collisions into the tidal radius are unlikely (Tremaine 2023).
In Figure 13, we show the pollution rate of comets into one WD in a WD-WD binary separated by   = 10 AU at various .In the simulation here, we directly check every timestep 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  = 6  to accurately capture all dynamics between comets and a stellar mass companions (see Figures 4, 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 5 g s −1 .Note that this is just below the detection limit at ∼ 5×10 5 g s −1 .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  ≈   are still unlikely to collide with the WD since the Safronov number is very high.
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, Θ ∝   , would increase and collisions would be even more unlikely.Thus, we expect a decrease in pollution rate as   increases.
Finally, the results here also give an order-of-magnitude estimate of Oort cloud comet pollution rate for other stellar-mass companions, assuming   ∼  * .

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 timescale.
In Figure 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 , and arrive at  crit .The timescale 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 timescale.Thus, without a companion, an Oort cloud is capable of delivering materials into the WD tidal disruption zone at a constant rate between   ≈ 5 × 10 7 g • s −1 and 10 8 g • s −1 over a 1 Gyr timescale.
Figure 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 the previous subsection.Furthermore, we again find that the accretion rate does not significantly decline over a 1 Gyr timescale.Note that in the case where  = 2.5 (bottom panel), the accretion decreases by about a factor of 3-4 over a 1 Gyr timescale.Note that in the context of WD pollution rate which ranges 5 orders of magnitude, this reduction factor of 3-4 over a Gyr timescale 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 semi-major 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 semi-major 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 subsection 5.1, the variations over every Myr bin is because of our resolution-limited simulation (∼ 10 7 comets).This is about 4 orders of magnitude smaller than the actual Oort cloud population of 10 11 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 (Figures 14,15) are within an order of magnitude of the mean.Thus, for a full Oort cloud, these variations would be within 10 −1 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  (planet)  ≈ 1 − 3 × 10 7 g • s −1 over a 1 Gyr timescale.Depending on , this rate may be reduced by a factor of 3-4 over this timescale.

DISCUSSION
Figure 4 of Blouin & Xu (2022) presents observed WD pollution rates for WDs with ages ranging from 1 Gyr to 8 Gyr.In this Figure, they show that the pollution rate ranges between 10 5 to 10 10 g s −1 , with the majority ranging between 10 6 to 10 8 g s −1 .They also find that pollution rates decrease by no more than one order of magnitude over the course 8 Gyrs.Recent observations by Mullally et al. (2024) find giant planet candidates around two polluted WDs with ages 1.5 and 5 Gyrs.The planets have masses 1 − 7 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 6 to 10 8 g s −1 , depending on the existence of a planetary companion.These rates are above the current detection limit at ∼ 5 × 10 5 g • s −1 .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 timescale.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 timescale.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   = 10 AU, the pollution rate is ≈ 3 × 10 5 g s −1 , 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   ∼  * .
We note that an advantage of the mechanism and reservoir presented here is that it is time-independent, in contrast to other works (c.f.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 timescale with a comet injection rate of ∼ 10 8 g s −1 , an Oort cloud would only deplete by roughly 5 × 10 −3  ⊕ worth of materials, which is much less than the current Solar System Oort cloud reservoir of 2 ⊕ .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.

Robustness of the Oort cloud
The results in this work can be scaled to other exo-Oort clouds with different  Oort and  Oort by: The other parameters of an Oort cloud, ,  1 (inner edge of the Oort cloud) and  2 (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 (Figures 9,11,13).Second, we fixed the inner semi-major axis edge at  1 = 3 000 AU, based on the Solar System Oort cloud.O'Connor et al. (2023) find that varying  1 changes the pollution rate by at most a factor of 4 for 10 3 AU ⩽  1 ⩽ 10 4 AU.Third, we fixed the outer semi-major axis edge at  2 = 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,  2 would not be different in another Oort cloud around a typical WD.

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 -folding decay timescale for the Oort cloud population is 4 − 18 Gyrs.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 timescale 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  ∼ 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 timescale.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).

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 −1 .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 −1 (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% 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 mainsequence evolution is reduced by about an order of magnitude for a typical kick strength of 0.75 km s −1 .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 6 to 10 7 g s −1 , depending on the existence of a planet.If the companion is stellar-mass, the pollution rate is much lower at ∼ 3 × 10 4 g s −1 .
However, the mass and structure of exo-Oort clouds in other mainsequence 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 characterisation 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 postmain-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 () ∝  − radial density profile -like the one currently existing in the Solar System as we currently understand it -can pollute a WD.

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 disk 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% of objects (up to 8 × 10 9 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;Klein et al. 2021;Doyle et al. 2021).This includes a detection of a Kuiper Beltanalogue 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 et al. (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.

Accretion Disc Delivery
We only study mechanisms to deliver comets into the Roche radius.After a comet achieves a distance  ⩽ 1 ⊙ , 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 8 g s −1 (Rafikov 2011).However, we observe pollution rate up to 10 10 g s −1 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 8 g s −1 .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.

Very Close Companions
So far, we only considered companions with semi-major axes on the order of   = 10 or 100 AU.There are observational evidence for companions at these semi-major axes (e.g.Veras et al. 2020;Blackman et al. 2021;Mullally et al. 2024).On the other hand, Vanderburg et al. (2020) present observational evidence for a   ∼ 1 Jup planet orbiting its WD at   ∼ 4 ⊙ .In addition, Gänsicke et al. (2019); Veras (2020) present evidence for a planet at   ∼ 15 ⊙ 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 3 , assuming typical incoming comets  ∼ 10 4 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 closein 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  ⊙ 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,   ⩾ 300 AU.

Widely Separated Companions
For planets or stellar-mass companions widely separated (  ⩾ 300 AU), our analysis cannot be applied, at least for Oort clouds with an inner semi-major axis at ∼ 3 000 AU like ours.As discussed earlier in various contexts, we typically assume that comets that encounter a companion have high eccentricity,  ∼ 1 (or equivalently,  ≫   ).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 et al. (2011) provide observational evidence for a planetary companion with mass   = 7 Jup at a separation of   ≈ 2 500 AU.Zuckerman (2014) presents a catalogue of 17 WDs with companions separated by more than 10 3 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 3 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   ⩾ 10 3 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.

CONCLUSION
We have studied if an exo-Oort cloud can pollute a WD with galactic tide and a companion over a 1 Gyr timescale 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 semi-major axis they experience close-encounters with the companion.We present a fast integration method that is capable of integrating 10 8 comets over a long simulation time.
The conclusions presented below assume a Solar System Oort cloud, with total cloud mass  Oort = 2 ⊕ containing  Oort = 10 11 objects.These pollution rate results are scalable to other exo-Oort clouds with different masses and number of objects.Our main conclusions are: (i) In the absence of any companions, exo-Oort clouds like our own Solar System's can pollute WDs at a rate ∼ 5 × 10 7 − 10 8 g s −1 .
(ii) We find that the dimensionless quantity is a good indicator of the relative importance of the angular momentum change induced by a companion compared to that of galactic tide (see Figure 3).When  ≲ 1, torque from a companion dominates galactic tide and reduces comet migration.When  ≫ 1, only galactic tide is important.
(iii) Stellar-mass companions with masses   ⩾ 0.1 ⊙ reduce the WD pollution rate to ∼ 3 × 10 5 g s −1 due to their strong angular momentum and scattering barriers, and a low-likelihood of direct collisions with the WD.
(iv) 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% effective ejection barrier.The modified loss cone also cannot account for migration in comets semi-major 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 7 g s −1 , assuming a Solar System Oort cloud.
(v) The powerlaw density profile structure of the Oort cloud does not significantly affect (by more than half an order of magnitude) the pollution rate.
(vi) The pollution rate without a planet stays constant over a 1 Gyr timescale.The pollution rate with a planet can decrease by a factor of ∼ 3 over a 1 Gyr timescale.
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 3 AU).
We show that exo-Oort clouds can potentially pollute old WDs at observed rates over Gyr timescales, 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 et al. 2024).With further observations of polluted WDs and characterisation of their companions (e.g. using Gaia as shown in Sanderson et al. 2022), our results can be applied to constrain sources of WD pollution.

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 4 AU, interacting with a companion at 10 1 − 10 2 AU, and engulfed by the WD Roche radius at 1 ⊙ .To accurately resolve these length scales, we would like to have the integration timestep  when the comet is far away to be large, and  to be small when the comet can interact with the companion.That is, we would like  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 timescale 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  = 3 bodies (WD-companion-comet) simulation will always cause  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 2-body system.
To resolve this, we develop a method to simulate just the comet; the REBOUND simulation will have only  = 1 particle, as illustrated in Figure 8.At each timestep, we manually set the total acceleration the comet experiences to   ,   ,   ,   denote the distance between the comet and the companion.Similarly for  * ,  * , etc. for the distance between comet and central star.Terms like , , ,  can be quickly calculated analytically because the WD-companion system is a 2-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  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 timestep method IAS15 (Rein & Spiegel 2015) as the integrator for our  = 1 simulation.The IAS15 integrator is capable of adaptively changing timesteps to resolve close encounters.In addition, we use the recent improvement on IAS15's adaptive method (Pham et al. 2024).We do not set a minimum timestep 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 timesteps 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 (Section A) and our implementation of it in REBOUND, we plot Figure B1.In this Figure, .Verification of the fast integration method by comparing   , the root-mean-square change in  = 1/, between two integration methods.'Normal' (black dashed line) means having 3 particles, central star, companion and comet.'Fast' (red solid line) is the integration method presented in this article with only 1 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 Figure 9.3 in Tremaine (2023).
we compare   , the root-mean-square change in  = 1/, between our fast integration method and the full 3-body simulation.At each , one thousand comets are initialised isotropically, with an initial semi-major axis at  = 1 000 AU and eccentricity  = 1 − /.The planet is at 5 AU with mass 10 −3  ⊙ .Comets are integrated over one orbital period to calculate   .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 Figure 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.

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 timescale is less than the timescales on which planet-planet interactions begin to become important (for example, this would be the secular timescale 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 , 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 ⟨⟩ = 10 • .We include the Sun and the outer Solar System planets.This integration  .The efficiency of performing rejection sampling on a broken powerlaw versus rejection sampling on a uniform distribution.We can see that the using the broken powerlaw distribution to create powerlaw samples is much more effective.With the broken powerlaw rejection sampling, the number of usable samples at any particular  is about half of the total simulated samples. is run for 100 Jupiter orbits.Figure C1 shows the semi-major axis distribution of these objects.We show here that the fast method is capable of statistically reproducing the Kuiper Belt peaks induced by resonance.

APPENDIX D: SEMI-MAJOR AXIS REJECTION SAMPLING
The distribution of comet semi-major axis itself follows a density profile described by a powerlaw relationship () ∝  − .We wish to study all  between 2 ⩽  ⩽ 4. We could draw samples at fixed

Figure 1 .
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  crit is the empty loss cone, on the right is the full loss cone.This figure closely follows Figure 1 in O'Connor et al. (2023).
(following Figure 1 in O'Connor et al. 2023) is a diagram illustrating these two regimes.

Figure 3 .
Figure3.Efficiency of the precession barrier induced by a companion, measured from numerical simulations.Simulation values of   ,  * , ,   ,  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 colour shows the efficiency of the planet's torque at inhibiting galactic tide, where  companion is the number of comets that can enter a certain  in the existence of a companion and galactic tide and  GT is the same but there is only galactic tide (no companion).Dashed lines show contours for  = 0.1, 1, 10. Above the  = 10 dashed line, the precession barrier efficiency is 0% and comets fully experience galactic tide and can be excited to high eccentricity, migrating inward in pericentre .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%.

Figure 4 .
Figure 4. Numerical scattering timescales for coplanar (blue) and isotropic (orange) comets interacting with a stellar-mass (top) and planetary (bottom) companion.The analytic diffusive timescales are shown in both panels, but can only be used in the bottom panel where   ≪  * .Blue and orange shaded areas are the range of measured numerical scattering timescales.In both cases, the high diffusive timescale flattens to the strong scattering regime (where comets are ejected within one pericentre passage with timescale  comet /2) as  decreases.

Figure 5 .
Figure 5. Diagram illustrating the modified loss cone theory by O'Connor et al. (2023).The engulfment loss cone is in grey, the ejection loss cone is in orange.Comets ejected are shown as plus signs and comets capable of engulfment are shown as crosses.This figure is created following closely Figure 12 in O'Connor et al. (2023).

Figure 6 .
Figure 6.Timescale comparison of galactic tide, companion-induced precession and scattering.From the bottom panel up, these show cases of increasing importance of companions.Timescales are measured from simulations of a Solar System-like Oort cloud (() ∝  −3.5 ,  ∈ [3 000, 10 5 ] AU).The shaded area shows the 3 spread of timescales for this particular distribution of comets.

Figure 9 .
Figure 9. Accretion rate of comets into  crit = 1 ⊙ over various Oort cloud structures , where () ∝  − .The mean simulation rate is shown in blue and the analytic expectation (Equation15) is in red.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  Oort = 10 11 comets and  Oort = 2 ⊕ .The blue shaded area shows the Poisson 1 confidence interval of the solid blue line (mean rate).

Figure 10 .
Figure10.Efficiency of a companion at suppressing galactic tide and preventing comets from entering .These are numerical results from simulations of incoming Solar System-like Oort cloud comets.The three cases here are the same as those in Figure6.The values of  for each case is shown on the right side with their corresponding colours.At  = 5  , a typical value for  is ∼ 5 × 10 7 .

Figure 11 .
Figure11.Accretion rate of comets into  crit = 1 ⊙ over various Oort cloud structures , where () ∝  − .In this simulation, there is a 10 Jup planetary companion at 10 AU.The mean simulation rate is shown in blue and the analytic expectation is in red.The analytic rate in the case where there is no planet (Equation15) is shown as a solid black line for comparison.The analytic prediction in the presence of a planet that can create a sufficiently strong ejection barrier (Equation33) is shown as a red dashed line.The blue shaded area shows the Poisson 1 confidence interval of the solid blue line (mean rate).

Figure 12 .
Figure12.Distribution of accreted comets' initial semi-major axis at  = 2. Histograms are counted per  initial bin.The blue histogram shows distribution when there is galactic tide only.The orange histogram shows the case in the case with galactic tide and a 10 Jup planet.Vertical lines show  crit,ej (dashed) and  crit (solid).

Figure 13 .
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 −1 is shown on the right axis, assuming a Solar System Oort cloud.

Figure 14 .
Figure14.Accretion rate over a 1 Gyr timescale from Oort clouds with  = 2.5 and 3.5.There are no companions in this case.Accretion rates from simulations (red) is binned per Myr.Analytic rates are shown as horizontal blue lines.A moving average with a sliding window of 10 Myr is shown in grey.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  Oort = 10 11 comets and total cloud mass  Oort = 2 ⊕ .

Figure 15 .
Figure15.Accretion rate over a 1 Gyr timescale from Oort clouds with  = 2.5 and 3.5.There is a planetary-mass companion with   = 10 Jup at   = 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.
Figure B1.Verification of the fast integration method by comparing   , the root-mean-square change in  = 1/, between two integration methods.'Normal' (black dashed line) means having 3 particles, central star, companion and comet.'Fast' (red solid line) is the integration method presented in this article with only 1 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 Figure9.3 inTremaine (2023).

Figure C1 .
Figure C1.The Kuiper Belt semi-major axis distribution reproduced by the fast integration method (blue histogram) compared with the normal integration method (orange histogram).The simulation includes effects from the outer Solar System planets (and the Sun).The simulation is run for 100 Jupiter orbits.

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

Table 1 .
Order-of-magnitude values of the dimensionless quantity , which describes the ratio of densities between the galaxy and the binary system (Equation25).

Table 2 .
Order-of-magnitude  values for comets with pericentre distance  = 5  .The dimensionless quantity  describes the comet's orbit relative to the companion orbit (Equation26).