Tidal Disruption of Planetesimals from an Eccentric Debris Disk Following a White Dwarf Natal Kick

The surfaces of many white dwarfs are polluted by metals, implying a recent accretion event. The tidal disruption of planetesimals is a viable source of white dwarf pollution and offers a unique window into the composition of exoplanet systems. The question of how planetary material enters the tidal disruption radius of the white dwarf is currently unresolved. Using a series of N-body simulations, we explore the response of the surrounding planetesimal debris disk as the white dwarf receives a natal kick caused by anisotropic mass loss on the asymptotic giant branch. We find that the kick can form an apse-aligned, eccentric debris disk in the range 30–240 au, which corresponds to the orbits of Neptune, the Kuiper Belt, and the scattered disk in our solar system. In addition, many planetesimals beyond 240 au flip to counterrotating orbits. Assuming an isotropic distribution of kicks, we predict that approximately 80% of white dwarf debris disks should exhibit significant apsidal alignment and a fraction of counterrotating orbits. The eccentric disk is able to efficiently and continuously torque planetesimals onto radial, star-grazing orbits. We show that the kick causes both an initial burst in tidal disruption events as well as an extended period of 100 Myr where tidal disruption rates are consistent with observed mass accretion rates on polluted white dwarfs.


INTRODUCTION
White dwarfs are the remnants of stars with main-sequence masses below 8 M ⊙ , which constitute an estimated 97% of stars in our Galaxy -including our Sun (Fontaine & Wesemael 2000).White dwarfs are extremely dense bodies with masses comparable to the Sun, despite tatsuya.akiba@colorado.edutheir sizes being closer to that of the Earth (Schatzman 1958).Most white dwarfs should have a core of carbon and oxygen surrounded by a thin, outer layer of hydrogen and helium (Burbidge & Burbidge 1954).However, an estimated 25 to 50% of observed white dwarfs show signs of metals such as calcium, magnesium, iron, and silicon in their spectra (e.g., Zuckerman et al. 2003Zuckerman et al. , 2010;;Koester et al. 2014).Metal polluted white dwarfs are found with a wide range of effective temperatures 3000 to 25 000 K (see Farihi 2016) with estimated cooling ages as old as 10 Gyr (Elms et al. 2022).
The presence of metals in white dwarf spectra is commonly attributed to the active accretion of planetary material following a tidal disruption event (e.g., Debes & Sigurdsson 2002;Jura 2003;Wang et al. 2019;Brouwers et al. 2022), given that any surface metals should sink to the core on a timescale much shorter than the white dwarf age through gravitational settling (see Jura & Young 2014).The radius at which this tidal disruption occurs is the white dwarf's Roche radius ≈ 1 R ⊙ (e.g, Li et al. 1998;Davidsson 1999;Bear & Soker 2013;Veras et al. 2014;Barber et al. 2016;Veras 2021).Once a planetesimal is tidally disrupted, the bound planetary debris forms a circumstellar disk which produces an observable infrared excess (e.g., Jura 2003;Jura et al. 2007), and many white dwarfs have been found with surrounding planetary material in this manner (Farihi 2016).The composition of tidally disrupted planetesimals can be inferred from metal abundances on polluted white dwarfs which provides a unique avenue for the study of exoplanet compositions (e.g., Zuckerman et al. 2007Zuckerman et al. , 2010;;Koester et al. 2014).
As a main-sequence star with mass below 8 M ⊙ runs out of hydrogen in its core, it turns into a red giant star to undergo subsequent fusion of heavier elements (Iben 1967).On the asymptotic giant branch, the outer layers of the star become unbound and about half of the stellar mass is lost before the core collapses into a white dwarf (Auer & Woolf 1965;Fusi-Pecci & Renzini 1976).When this mass ejection occurs anisotropically, a natal kick is imparted on the white dwarf upon its formation (Fellhauer et al. 2003).The kick magnitude is expected to be 1 to 3 km s −1 (e.g., Fregeau et al. 2009;El-Badry & Rix 2018;Hamers & Thompson 2019); the direction of the kick with respect to the planetesimal disk plane is unknown.Stone et al. (2015) considered the effect of this natal kick on exo-Oort clouds.The kick maps many comets onto radial, plunging orbits which produces a temporary burst of tidal disruption events.However, their Monte Carlo approach followed post-kick orbits on a short timescale without self-gravity, and hence does not apply to the cooler population of metal-polluted white dwarfs.In this letter, we show that the white dwarf natal kick results in the formation of an apse-aligned, eccentric debris disk of planetesimals which produces tidal disruption events at a rate consistent with observed mass accretion rates for 100 Myr.

ECCENTRIC DEBRIS DISK FORMATION AND STABILITY
When a gravitational recoil kick is imparted on a supermassive black hole, the surrounding stellar orbits in a nuclear star cluster can form an eccentric, apse-aligned disk (Akiba & Madigan 2021, 2023), and these eccentric disks exhibit stellar tidal disruption rates as high as 0.1 yr −1 gal −1 , three to four orders of magnitude higher than rates expected from isotropic distributions (Madigan et al. 2018a).The dynamics work on all scales in a near-Keplerian system and thus are directly applicable to planetesimals surrounding white dwarfs following the impartment of a natal kick.To quantify apsidal alignment, we make use of the eccentricity vector where ⃗ v is the velocity vector, ⃗ j is the (specific) angular momentum vector, M * is the white dwarf mass, and r is the unit position vector.
The eccentricity vector points from the apoapsis to the periapsis of a given orbit, with a magnitude equal to the scalar eccentricity.
The mean eccentricity vector is a measure of apsidal alignment defined by where ⃗ e i is the eccentricity vector of the i-th planetesimal and N is the number of planetesimals considered.When the white dwarf experiences an in-plane kick, planetesimals on initially circular orbits will align their eccentricity vectors such that where v kick is the natal kick speed and v circ is the initial circular speed of the planetesimals.Apsidal alignment is strongest when |⟨⃗ e⟩|= 1.
From Equation 3, this occurs at a characteristic radius For M * = 0.6 M ⊙ and v kick = 1 km s −1 , this radius occurs at r c = 240 AU.For v kick = 3 km s −1 , r c = 30 AU.In the solar system, an orbital distance of 30 AU corresponds to that of Neptune and the Kuiper Belt, a dynamically rich region of space that includes both kinematically cold and hot primordial planetesimal populations intermixed with those in orbital resonance with Neptune (e.g., Jewitt & Luu 1993;Malhotra 2019).A distance of 240 AU corresponds to that of the scattered disk, a population of planetesimals on eccentric orbits with periapses that bring them into contact with Neptune's orbit (Duncan & Levison 1997;Vokrouhlický & Nesvorný 2019).We note that this comparison does not take into account orbital expansion due to the star's mass loss (e.g., Veras et al. 2013).The existence of Kuiper Belt or scattered disk-like structures in exoplanet systems has been inferred from observations (e.g., Geiler et al. 2019;Wyatt 2020); disks of icy bodies in the outskirts of planetary systems should be common.
Dynamical stability in eccentric disks comes about via mutual gravitational torques between orbiters (Madigan et al. 2018a).When a planetesimal precesses ahead of the eccentric disk, its orbit is negatively torqued by the disk, and its angular momentum decreases which in turn increases its scalar eccentricity.This change in eccentricity works to slow down the orbit's precession and allows the rest of the disk to catch up to it.The opposite is true for an orbit that lags behind the disk: it feels a positive torque which circularizes the orbit and increases its precession speed.In this way, the eccentric disk maintains its apsidal alignment as individual orbits undergo oscillations in precession speeds and eccentricity.It is the latter oscillation in eccentricity that causes an enhancement in the rate of tidal disruption events as strong mutual torques throw planetesimals onto radial, stargrazing orbits.

Initial Setup
We use the open-source, N -body simulation package REBOUND (Rein & Liu 2012), with code units of G = 1, M * = 1, and r c = 1 where r c is the characteristic radius given by Equation 4. For concreteness, we translate to physical units assuming a white dwarf mass of M * = 0.6 M ⊙ and a kick velocity of v kick = 1 km s −1 throughout the letter.This sets r c = 240 AU.However, the dynamics are completely scalable; the presented results can be applied to a different M * or v kick by scaling lengths by r c ∝ M * /v 2 kick (see Equation 4) and timescales by t ∝ r 3/2 c .
To study the instantaneous post-kick distribution, we initialize N = 5 × 10 4 massless planetesimals in an axi-symmetric, thin disk spanning four orders of magnitude in semimajor axis space.We use a surface density profile Σ ∝ a −1.5 as motivated by Hayashi (1981).The number of planetesimals was chosen to have a sufficient number density out to Oort cloud distances of a ≈ 10 3 AU.Inclination is Rayleigh-distributed with scale parameter σ = 3 • while longitude of periapsis and mean anomaly are uniformly distributed in [0, 2π).We perform 2500 simulations in which the white dwarf instantaneously gains a velocity of v kick = 1 km s −1 with respect to its initial frame of reference at time t = 0, allowing orbits to impulsively respond.We define the planetesimal disk plane to be in the x-y plane with the angular momentum vector pointing in the +zdirection.We further define the +x-direction to be the direction of the in-plane component of the kick.We randomly sample the kick angle with respect to the disk plane from an isotropic distribution.The resulting distribution of kick angles is shown in Figure 1.We note that α is measured with respect to the +z-direction.Since α is uniform in cosine, in-plane kicks are more likely than out-of-plane kicks.

Structure of Eccentric Debris Disks
In Figure 2, we compare the distributions of orbital angular momenta and apsidal alignment of planetesimals as a function of semi-major axis at t = 0 post-kick for simulations with different kick angles with respect to the planetesimal disk.We opt to show kick angles in the range α = 0 • to 90 • since these results are completely symmetric about α = 90 • ; the distributions for α = 60 • and 120 • look identical, for instance.For each simulation, we show the distribution of the z-component of angular momentum in the top panel.The plot is color-coded by the y-component of the eccentricity vector which is the direction apsidal alignment is expected with a kick in the +x-direction.In the bottom panel for each simulation run, we show the distribution of the magnitude of the mean eccentricity vector, |⟨⃗ e⟩| as defined in Equation 2. |⟨⃗ e⟩|= 0 when the disk is axi-symmetric and deviates toward unity as the disk becomes apse-aligned.
For kicks that are nearly out-of-plane, planetesimals remain on prograde orbits and their eccentricities increase with semi-major axis.There are no patterns in the distribution of e y and |⟨⃗ e⟩|= 0 throughout the disk.When α ≥ 30 • , two dynamically important structures are observed: 1. strong apsidal alignment in the +y-direction, and 2. a significant retrograde population of planetesimals beyond 240 AU.The retrograde planetesimal population emerges at large semi-major axes as the initial speed of the planetesimals becomes smaller than the natal kick speed.In the reference frame of the kicked white dwarf, the planetesimal velocity can flip such that the post-kick orbit is retrograde with respect to the initial angular momentum axis of the disk.In particular, the in-plane kick case shows a prograde population of planetesimals that is entirely apse-aligned in the +y-direction, a retrograde population at a > 240 AU that is aligned in the same direction, and a retrograde population at a > 2000 AU that is aligned in the opposite direction.Apsidal alignment is statistically significant throughout the debris disk, and the retrograde fraction is significant beyond 240 AU.It should be emphasized that at least one retrograde planetesimal is found in approximately 80 % of our simulations.Apsidal alignment and retrograde planetesimals in white dwarf systems are thus natural consequences of natal kicks.

Eccentric Debris Disk Evolution and the Tidal Disruption of Planetesimals
To investigate the burst of tidal disruption events following the kick, we redistribute N = 5 × 10 4 massless planetesimals in the region a = 240 to 2400 AU where apsidal alignment is the strongest.This is done in order to avoid issues with small number statistics in the detection of tidal disruption events.For the most probable case of an in-plane kick, we randomly select a subset of N = 400 planetesimals in the range a = 240 to 480 AU and switch them from massless test particles to mas- sive, self-interacting particles to study the eccentric debris disk's long-term evolution.We use REBOUND's IAS15 integrator, a high-order, non-symplectic integrator with adaptive timestepping (Rein & Spiegel 2015).The low N and narrow range of semi-major axes explored are in part due to the computational limitations of IAS15, but high-accuracy integration is critical in studying eccentric disk evolution 1 .
We vary the total disk mass in the range M disk = 1 to 100 M to study the effect of disk mass on the tidal disruption rate.Tidal disruption events are detected by comparing the periapsis distance, r − = a(1 − e) with the Roche radius defined to be R Roche ≡ 1 R ⊙ at each 1 We also tested WHFast, a low-order, symplectic integrator (Rein & Tamayo 2015).Our results showed that WHFast runs significantly deviate from IAS15 runs of the same disk mass with the former exhibiting total energy and angular momentum errors of order unity within a Myr.The energy and angular momentum errors are < 10 −10 even for long-term integrations up to 100 Myr using IAS15.Due to the high-resolution required to accurately integrate eccentric orbits near their periapses accompanied by strong mutual torques between orbits, high-order integrators such as IAS15 are required to fully study the evolution of an eccentric disk.time step.When r − < R Roche , the planetesimal is considered tidally disrupted.Each massive simulation is stopped at t = 10 Myr.The only exception is the simulation with M disk = 20 M which is run until t = 100 Myr to show the long-term enhancement of tidal disruption rates.The simulation run is stopped at 100 Myr due to two computational limitations of our low N setup: 1. the number of tidally disrupted planetesimals becomes comparable to N , and 2. the two-body relaxation timescale (Rauch & Tremaine 1996), which is proportional to N for a given disk mass, becomes comparable to the simulation time (see Madigan et al. 2018b).
In Figure 3, we show planetesimal orbit projections at t = 0 pre-kick, t = 0 post-kick, and t = 100 Myr after the kick for the M disk = 20 M simulation.Initially, the pre-kick orbits are all circular and prograde.Post-kick, the orbits' eccentricity vectors coherently point in the +y-direction.We also see the emergence of retrograde orbits which are apsidally aligned in the same direction.For 100 Myr, the alignment of eccentricity vectors is maintained while the eccentric disk coherently precesses in the prograde direction.We perform a quick estimation of the fraction of kicked white dwarfs that tidally disrupt a planetesimal.At t = 0 post-kick, we compare the new periapsis distances of the redistributed 5 × 10 4 massless planetesimals in the range a = 240 to 2400 AU to R Roche .Of the massless simulations with isotropically-distributed kicks, we find that 10% of white dwarfs immediately tidally disrupt at least one planetesimal upon receiving a kick of 1 km s −1 .This fraction increases to 16% if a 3 km s −1 kick is assumed instead.Interestingly, Stone et al. (2015) predict a similar fraction of comets that tidally disrupt in exo-Oort clouds.However, this fraction is notably lower than the expected 25 to 50% of older white dwarfs that are metal-polluted.We show in the following section that eccentric debris disks can increase this fraction by inducing tidal disruption events at later times through strong mutual torques between planetesimals.

Long-term Enhancement of Tidal Disruption Rates in Eccentric Disks
In the top panel of Figure 4, we show the cumulative mass accreted by the white dwarf from planetesimal tidal disruption events during the M disk = 20 M simulation.We note that we do not include a treatment of the accretion process; we are inherently assuming that a planetesimal, once tidally disrupted, is entirely ac-creted by the white dwarf on orbital timescales t orb ≈ kyr.We see that the eccentric debris disk is able to produce tidal disruption events for an extended period of 100 Myr following the impartment of the kick.This particular disk mass simulation implies a very high accretion rate of 10 12 g s −1 .Except for a few periods or dormancy (e.g., t = 58 to 69 Myr), at least one tidal disruption event is observed every few Myr and hence the expected time-averaged accretion rate is at least 10 9 g s −1 throughout the simulation.In addition, due to the significant retrograde fraction induced by the kick, we find that approximately 40% of tidally disrupted planetesimals are on retrograde orbits when they become tidally disrupted.If the initial disk angular momentum axis is the same as that of the star's rotation axis, this mechanism would predict a significant fraction of white dwarf systems to have a circumstellar debris disk following a tidal disruption event that is retrograde with respect to the white dwarf spin axis.
In the center panel of Figure 4, we plot the eccentricity evolution of a small sample of planetesimals that become tidally disrupted.While the initial jump to high eccentricities at t = 0 is due to the natal kick, the kick alone is unable to tidally disrupt these planetesimals.The subsequent evolution to higher eccentricities is caused by the strong mutual torques between planetesimals within the eccentric disk (see Section 2).The bottom panel of Figure 4 shows the apsidal alignment evolution.Statisticallysignificant alignment is maintained throughout the simulation.We note that the oscillation in |⟨⃗ e⟩| is caused by the periodic alignment and misalignment between the eccentric disk and planetesimals near the outer edge of the disk that break off and precess in the retrograde direction (see Madigan et al. 2018a;Akiba & Madigan 2021).This figure shows that the oscillation is correlated with bursts of tidal disruption events; as the retrograde-precessing plan- etesimals come back into alignment with the eccentric disk, mutual gravitational torques increase the amplitude of the eccentricity oscillations which consequently increase tidal disruption rates.The white dwarf natal kick causes an initial burst of tidal disruption events followed by an extended 100 Myr period where tidal disruption rates from the eccentric debris disk are consistent with some of the highest mass accretion rates observed.We note that t = 100 Myr corresponds to polluted white dwarfs on the warmer end with temperatures around 20 000 K (see Farihi 2016); a higher N simulation is needed to extend this study to white dwarfs with cooling ages Gyr or older.

Tidal Disruption Rate Dependence on Eccentric Debris Disk Mass
The expected planetesimal tidal disruption rate has a steep dependence on the mass of the eccentric debris disk.The following anal-ysis assumes that the accretion process happens on the orbital timescale t orb ≈ kyr which is much shorter than the secular timescale defined by t sec ≡ (M * /M disk ) t orb .The eccentricity oscillations that cause tidal disruption events happen on the secular timescale (Madigan et al. 2018a), so the tidal disruption rate ṄTDE ∝ t −1 sec ∝ M disk .The mass accretion rate is estimated as Ṁ = m ṄTDE where m is the planetesimal mass.In our simulations, we fix N and vary M disk , so m ∝ M disk .Hence, we analytically predict Ṁ ∝ (M disk ) 2 .
In Figure 5, we show the estimated mass accretion rate from tidal disruption events for each of our massive disk simulations.For simulations with at least one tidal disruption event detected, the correlation is well-modeled by our analytic prescription Ṁ ∝ (M disk ) 2 .The M disk = 1 and 2 M simulations detect no tidal disruption events in the first 10 Myr, but this is likely due to the low N and short simulation time.For instance, M disk = 1 M implies a secular timescale of t sec ≈ 100 Myr.If instead, we extrapolate the Ṁ ∝ (M disk ) 2 scaling to lower debris disk masses, we predict that typical instantaneous accretion rates of 10 7 g s −1 on hydrogen white dwarfs are reproduced for any eccentric debris disk with mass above 0.03 M .The highest measured instantaneous accretion rates of 10 9 g s −1 are predicted for M disk > 0.3 M .Finally, the highest time-averaged accretion rates inferred from helium white dwarfs of 10 11 g s −1 require debris disk masses of at least a few M .

DISCUSSION
We perform numerical simulations to explore the formation of an eccentric debris disk following a white dwarf natal kick.We run 2500 massless simulations with isotropicallydistributed kick angles to estimate the fraction of white dwarf debris disks with significant apsidal alignment, retrograde fraction, and planetesimal tidal disruption rates.For the most likely case of an in-plane kick, we follow the evolution of the eccentric disk and tidal disruption rates by performing a series of massive N -body simulations.Our main findings are: 1.An eccentric debris disk and a retrograde planetesimal population are expected for approximately 80% of planetesimal debris disks.The retrograde population forms beyond 240 AU for a white dwarf natal kick of 1 km s −1 .
2. For a 1 km s −1 kick, 10% of white dwarfs are expected to immediately tidally disrupt a planetesimal.This fraction increases to 16% for a 3 km s −1 kick.This is consistent with earlier results of Stone et al. (2015).
3. The eccentric debris disk hosts planetesimal tidal disruption rates consistent with observed mass accretion rates for an extended 100 Myr period.Around 40% of the tidal disruption events are retrograde with respect to the initial disk angular momentum axis.
4. The predicted mass accretion rates have a steep dependence on the eccentric debris disk mass, Ṁ ∝ (M disk ) 2 .A debris disk of mass above 0.03 M is required to explain typical instantaneous mass accretion rates of 10 7 g s −1 .
The estimated fraction of white dwarf systems that exhibit eccentric debris disks, retrograde planetesimal populations, and tidal disruption events is dependent on the density profile of the pre-kick debris disk.Our results from Figure 5 show that our mechanism requires at least 0.03 M in the debris disk within the range 240 to 480 AU assuming a kick of 1 km s −1 .Throughout the letter, we have assumed that the natal kick is impulsive.In reality, the investigated range of semi-major axes lies in the complex transition region between planetesimals with a < 10 AU where mass loss seems adiabatic, and planetesimals with a > 1000 AU where the mass loss seems impulsive (Stone et al. 2015).Furthermore, stellar models predict that mass loss is pulsed (e.g., Reimers 1975;Bloecker 1995;Bonsor et al. 2011) which adds more layers of complication.In the future, our work will consider the treatment of mass loss and the complex orbital evolution of planetesimals in response to a more realistic natal kick.
Additionally, our current work only considers a pre-kick circular disk with equal-mass particles.Moving forward, we plan to repeat this numerical experiment with a scattered disk-like distribution of planetesimals accompanied by one or more giant planets.Preliminary results have shown that the initial burst of tidal disruption events is more significant for the scattered disk but apsidal alignment is weaker.A follow-up study with long-term simulations is necessary to determine whether the eccentric debris disk and high tidal disruption rates are maintained in this case.We note that higher N simulations are necessary to thoroughly explore debris disks with masses around 1 M or below.In future studies, we aim to extend our mechanism to polluted white dwarfs with cooling ages Gyr or older, and make our claims about the tidal disruption event bursts more robust through longer and higher N simulations.

Figure 1 .
Figure 1.A histogram of kick angles with respect to the orbital angular momentum axis.α = 0 • or 180 • indicates an out-of-plane kick where the kick is parallel or anti-parallel to the angular momentum vector, whereas α = 90 • means that the kick is in the planetesimal disk plane.The analytic expectation assuming an isotropic distribution is shown with the blue, solid line and the numerical results of sampling 2500 times are shown in purple bars.The most likely natal kick direction is in the planetesimal disk plane.

Figure 2 .
Figure 2. A comparison of the planetesimals' orbital angular momenta and apsidal alignment as a function of semi-major axis at t = 0 post-kick ranging from α = 3 • to 90 • .(Top:) The z-component of angular momentum normalized by the circular angular momentum, j circ = √ GM * a. Orbits above the dashed line are prograde and those below are retrograde.Orbits are more eccentric closer to the dashed line and are circular when |j z /j circ |= 1.The planetesimal orbits are color-coded by the y-component of the eccentricity vector.(Bottom:) The magnitude of the mean eccentricity vector, ⟨⃗ e⟩.The green, shaded region indicates the noise floor above which alignment is statistically significant at the 3σ-level.Retrograde orbits and apsidal alignment both emerge as the kick tends toward an in-plane direction.

Figure 4 .
Figure 4. Evolution of tidal disruption rates and apsidal alignment for the M disk = 20 M simulation.(Top:) The cumulative mass accreted by the white dwarf from tidal disruption events over time.The axis on the right shows the fraction of the debris disk that is tidally disrupted.(Center:) Eccentricity evolution of a small sample of planetesimals that become tidally disrupted.The quantity, 1 − e, is plotted on the y-axis using a logarithmic scale.The time at which tidal disruption is detected is denoted by an open black circle.(Bottom:) The evolution of the magnitude of the mean eccentricity vector.The green shaded region indicates the noise floor above which apsidal alignment is statistically significant at the 3σ-level.

Figure 5 .
Figure5.The average mass accretion rate, Ṁ from tidal disruption events as a function of the debris disk mass, M disk .The average is calculated from the first 10 Myr of each simulation.The axis on the right shows the corresponding fraction of the debris disk that is tidally disrupted in each case.For simulations with at least one tidal disruption event detection, the results are consistent with a power-law dependence Ṁ ∝ (M disk ) 2 .