Necroplanetology: Simulating the Tidal Disruption of Differentiated Planetary Material Orbiting WD 1145+017

The first exoplanet discovery showed that planets can be found orbiting post-main-sequence objects like pulsars (Wolszczan & Frail 1992) and the Kepler mission showed that most main-sequence stars have at least one planet (Cassan et al. 2012; Howard et al. 2012). Combined with the observed infrared excess, detectable dusty disks, and metal-polluted spectra of ∼4%, ∼2%, and ∼30% of white dwarfs respectively (Zuckerman et al. 2003, 2010; Koester et al. 2014; Farihi 2016), we know that planetary systems survive stellar evolution in some form for a significant fraction of stars. The "polluted" spectra refer to metal lines that must be caused by recently accreted material since the high surface gravity of a white dwarf would sink metals beneath the photosphere over the course of weeks (Schatzman 1948). There are discrepancies in the observed chemical composition, spatial and kinematic distributions, and line strengths of the polluted material compared to what would be expected for accretion from the interstellar medium, identifying the presumed source of the pollution as planetary debris orbiting the white dwarf instead (Jura 2003; Kilic & Redfield 2007; Farihi et al. 2010).

Many previous studies have used high-resolution spectroscopy of these polluted white dwarfs to reconstruct the chemistry of these former exoplanetary objects (e.g., Zuckerman et al. 2007; Gänsicke et al. 2012; Jura & Young 2014; Farihi et al. 2016; Melis & Dufour 2017; Harrison et al. 2018; Hollands et al. 2018; Doyle et al. 2019; Bonsor et al. 2020), obtaining information about their bulk composition at a level unmatched by current exoplanetary characterization techniques applied to main-sequence stars. But these chemical analyses still arrive late to the scene: the planetary material has been destroyed and accreted already so that the characteristics of individual bodies cannot be determined (Veras 2016). During Campaign 1 of the K2 mission Vanderburg et al. (2015) reported the discovery of irregular transits caused by what appeared to be an asteroid or minor planet still being disrupted by its host WD 1145+017. Follow-up photometry from the ground showed that there were multiple periodic signals of varying mean anomaly, shape, asymmetry, and depth clustered around a ∼4.5 hr orbit (e.g., Vanderburg et al. 2015; Gänsicke et al. 2016; Rappaport et al. 2016; Gary et al. 2017).

Xu et al. (2016) found circumstellar absorption in the polluted metal lines which Redfield et al. (2017) observed to vary on both minute timescales corresponding to transits coincident with the spectroscopic observations and over the course of months due to a mechanism that remains uncertain. Redfield et al. (2017) proposed a warped eccentric disk interior to the transiting material that is consistent with most of the observed features while Farihi et al. (2017) proposed a magnetospheric accretion process. However, the Farihi et al. (2017) model requires magnetic fields higher than the limits imposed by the X-ray flux emanating from WD 1145+017 (Farihi et al. 2018). Cauley et al. (2018) amends the Redfield et al. (2017) model by describing a series of 14 individual narrow disks that share the same focus but have different eccentricities and apsidal angles that increase linearly with distance from the white dwarf. Vanderburg & Rappaport (2018), Xu et al. (2019), and Fortin-Archambault et al. (2020) discuss further theoretical and observational work on the gas disk and its relationship with the transiting debris.

Gurri et al. (2017) used N-body simulations to place limits on the eccentricity and mass of the transiting objects under certain assumptions of mass ratios and restricted drift in the orbital period. Complementing this study, Veras et al. (2017) used the N-body code PKDGRAV (Richardson et al. 2000; Stadel 2001) to simulate the disruption of an individual asteroid to place constraints on the bulk density and show that the transits were best explained by a differentiated body orbiting within the Roche limit for its mantle density but outside the limit for the core. However, these simulations only tested one possible differentiated internal structure to compare to a homogeneous rubble pile. This paper employs a similar methodology but expands it to explore a variety of internal structures and generate lightcurves that can be compared to data. In Section 2 we describe the parameters for this suite of simulations and in Section 3 we present the lightcurves generated to show that structural differences lead to observationally distinguishable features in the lightcurve. Section 4 reports our best fit for a specific transit lightcurve from Gary et al. (2017). We conclude in Section 5 by discussing the likely structure of the disrupting material transiting WD 1145+017, relating our work to the Cauley et al. (2018) model, and the implications of using this method to study the interior structure of other exoplanetary bodies.

We set up the disruption simulations by randomly packing N = 3000 spheres of equal radius, creating a single rubble pile with ${M}_{b}={10}^{24}\,\mathrm{kg}$ and no spin. This higher mass is chosen for computational considerations discussed later in the text. The rubble pile structure is allowed to vary between all permutations of core volume fraction f_c = 0.15, 0.25, and 0.35; mantle/core density ratio ${\rho }_{m}/{\rho }_{c}=0.25,0.40$, and 0.55; the presence or absence of a thick crust with volume fraction f_l = 0.1 and crust/core density ratio of ρ_l/ρ_c = 0.1; and bulk density ρ = 3 or 4 g cm⁻³. By fixing both the bulk mass and density, the bulk radius is determined and the radii of the individual particles are scaled appropriately. For each simulation, the rubble pile was placed in a circular e = 0, edge-on i = 90° orbit with a semimajor axis a = 0.0054 au around a white dwarf with M_⋆ = 0.6 M_⊙ to be consistent with Veras et al. (2017). At the time the simulations were performed, this was only a fiducial mass estimate for WD 1145+017, but Izquierdo et al. (2018) report the mass of WD 1145+017 as 0.63 ± 0.03M_⊙. The parameters common to all of the final simulations are listed in Table 1 while Table 2 lists the parameters that were varied between simulations. Figure 1 is a cartoon illustration demonstrating how changing each of these parameters affects the structure of the simulated rubble pile.

Zoom In Zoom Out Reset image size

Table 2.  The Parameters which Were Varied between Simulations

	fc	ρm /ρc		fl	ρl/ρc		${{ \mathcal X }}_{\nu =110}^{2}$
Structure Number	Core Volume Fraction	Mantle/Core Density Ratio	Crust?	Crust Volume Fraction	Crust/Core Density Ratio	ρb $({\rm{g}}\,{\mathrm{cm}}^{-3}$) Bulk Density	Chi-Squared with 110 degrees of freedom
1a	0.15	0.25	No	⋯	⋯	3.0	126.33
2a	0.15	0.40	No	⋯	⋯	3.0	126.73
3	0.15	0.55	No	⋯	⋯	3.0	184.73
4	0.25	0.25	No	⋯	⋯	3.0	164.39
5	0.25	0.40	No	⋯	⋯	3.0	167.16
6	0.25	0.55	No	⋯	⋯	3.0	189.46
7	0.35	0.25	No	⋯	⋯	3.0	167.93
8	0.35	0.40	No	⋯	⋯	3.0	161.74
9	0.35	0.55	No	⋯	⋯	3.0	179.70
10	0.15	0.25	Yes	0.1	0.1	3.0	146.33
11a	0.15	0.40	Yes	0.1	0.1	3.0	128.76
12b	0.15	0.55	Yes	0.1	0.1	3.0	219.39
13	0.25	0.25	Yes	0.1	0.1	3.0	152.69
14	0.25	0.40	Yes	0.1	0.1	3.0	158.65
15b	0.25	0.55	Yes	0.1	0.1	3.0	191.29
16b	0.35	0.25	Yes	0.1	0.1	3.0	192.63
17b	0.35	0.40	Yes	0.1	0.1	3.0	190.09
18	0.35	0.55	Yes	0.1	0.1	3.0	188.69
19a	0.15	0.25	No	⋯	⋯	4.0	131.42
20a	0.15	0.40	No	⋯	⋯	4.0	129.57
21	0.15	0.55	No	⋯	⋯	4.0	171.49
22	0.25	0.25	No	⋯	⋯	4.0	159.53
23	0.25	0.40	No	⋯	⋯	4.0	135.95
24	0.25	0.55	No	⋯	⋯	4.0	157.27
25	0.35	0.25	No	⋯	⋯	4.0	169.05
26	0.35	0.40	No	⋯	⋯	4.0	188.30
27	0.35	0.55	No	⋯	⋯	4.0	176.13
28	0.15	0.25	Yes	0.1	0.1	4.0	137.11
29	0.15	0.40	Yes	0.1	0.1	4.0	140.48
30	0.15	0.55	Yes	0.1	0.1	4.0	170.63
31	0.25	0.25	Yes	0.1	0.1	4.0	173.21
32	0.25	0.40	Yes	0.1	0.1	4.0	161.53
33b	0.25	0.55	Yes	0.1	0.1	4.0	193.29
34	0.35	0.25	Yes	0.1	0.1	4.0	185.12
35	0.35	0.40	Yes	0.1	0.1	4.0	189.90
36	0.35	0.55	Yes	0.1	0.1	4.0	169.22

Notes. Ten instances of Structure 1 were run to confirm that there was no strong dependence on the random seed. Five copies of Structure 1 with ${M}_{b}\in [{10}^{18},{10}^{24}]\,\mathrm{kg}$ were run to verify that the speed of disruption followed the t_fill dependence described in Veras et al. (2017). The chi-squared statistic of each structure's best fit is listed in the rightmost column. The number of free parameters is the same for all fits so lower values of ${{ \mathcal X }}_{\nu =110}^{2}$ correspond to better fits.

^aThe five best-fitting structures are bold-faced. They all have a low core fraction f_c = 0.15 and mantle/core density ratio ρ_m /ρ_c = 0.25 or 0.4. ^bThe five worst-fitting structures are italicized. All five have a crust and four of them have core fractions f_c > 0.15. The worst-fit, Structure 12, has a low core fraction f_c = 0.15, but its mantle/core density ratio is high: ρ_m /ρ_c = 0.55.

Download table as:  ASCIITypeset image

Previous studies of tidally disrupted differentiated bodies (Leinhardt et al. 2012; Veras et al. 2017) used the same structure consistently: f_c = 0.35, ρ_m/ρ_c = 0.25 and no crust. However, Carter et al. (2015) showed that when simulating planetesimal growth via core accretion, the final core fractions spanned the complete range between 0 and 1 regardless of orbital period with most planetesimals ending up with a core fraction between 0 and 0.4, motivating the range of f_c values chosen in this study. They also showed that the relative mass fractions of the core and mantle varied widely but were strongly correlated with orbital period. Since we cannot be certain where the disrupting bodies originally formed, we adopt the three values above to roughly represent a mantle that lies between an "icy" volatile-rich and a "rocky" Earth-like density with a "metallic" iron-rich core. These names are only suggestive and a single parameter for the mantle/core density ratio cannot capture the details of the possible compositions and substructures of the differentiated layers. The two choices of bulk density were motivated by the limits placed by Veras et al. (2017) when previously modeling the WD 1145+017 system. We include a crust to see whether or not multiple layers of differentiation could be detected.

Using the leapfrog integrator of the N-body simulation code REBOUND (Rein & Liu 2012), we integrate these 36 simulations for 100 orbits with P = 4.5 hr, e = 0.0, M_⋆ = 0.6 M_⊙, and an integration timestep of dt = 10 s. We record a snapshot of the simulation after the second orbit using the Simulation Archive feature of REBOUND (Rein & Tamayo 2017) then continue to evolve this snapshot for a half-orbit while saving the positions of each particle for each timestep. By projecting the particles against the face of the white dwarf divided into ∼50,000 pixels and counting the pixels obscured by a simulation particle, we generate lightcurves for each simulation at a cadence comparable to the ground-based photometry available for the system. We smooth the scatter in the pixel-counted lightcurves using a Savitzky–Golay filter (Savitzky & Golay 1964) before fitting to the photometry. We fit these template lightcurves to the photometry by scaling the depth and duration of transit with two free parameters. Scaling the depth is equivalent to inflating all particles equally, differing from inflating only the particles separated from the bulk of the asteroid as Veras et al. (2017) did.

While the bulk density determines whether or not tidal disruption occurs, Veras et al. (2014) showed that the timescale of disruption t_fill, the time it takes for a rubble pile to completely disrupt and fill out a debris ring, depends on the bulk radius R_b. For our chosen ρ_b and M_b, the snapshot lightcurve corresponds to a t/t_fill ≈ 0.06. Theoretically, these lightcurves should match any identically structured rubble pile of a different mass at roughly the same t/t_fill. We verified this for five values of ${M}_{b}\in [{10}^{18},{10}^{24}]\,\mathrm{kg}$ and 10 random seeds for Structure 1. This verification was not done for every simulation because a smaller M_b results in smaller particles, and the face of the white dwarf must be divided into more pixels to count these particles accurately, slowing down the calculation of the lightcurve. Testing a range of masses for each of the 36 structures was not computationally tractable. The random seed tests also allowed us to ensure that the observed lightcurve differences were genuinely caused by the structural differences in the rubble pile and not simply by the particular random packing. Unfortunately, the high mass we chose for computational purposes sped up the the disintegration of the rubble pile such that only two orbits after the REBOUND snapshot had meaningful transit signatures before creating a debris ring. We only used the first transit for fitting to lightcurves as described in the following section, severely limiting our ability to study the temporal evolution of the transit feature.

Canup (2010), Leinhardt et al. (2012), and Veras et al. (2017) have previously described the geometry of mantle disruption. First the entire rubble pile is stretched into a lemon shape as tidal forces act more strongly on the rarer mantle material. When the self-gravity of the mantle material is no longer sufficient to maintain cohesion, particles stream out from the end closer to the gravitational host at the L1 point and then at the farther end at the L2 point. This streaming is initially intermittent and slow, then speeds up as the number of gravitationally bound mantle particles decreases, pushing some particles into a slightly shorter orbital period and others into a slightly longer one. If the core is safely beyond its Roche limit, it remains safely intact, but otherwise breaks up and spreads out to fill out the ring. When Veras et al. (2017) simulated this particular system, they found that there was a difference in the rate of material loss between the two streams and eventually the space between them filled out to form one ring. Our simulations are consistent with Veras et al. (2017) in this regard. The time at which we capture the lightcurve corresponds to when some of the mantle has already been stripped but before it has completely dispersed into a ring.

Figure 2 shows how changing individual structural features changes the lightcurve. Increasing the mantle/core density ratio weakens the wings relative to the core since less of the mantle has been removed from the parent body. Increasing the core fraction has a similar effect but also decreases the asymmetry of the lightcurve caused by the L1/L2 mass loss difference. This trend could be caused by either the increased gravitational potential from the core or because there is less mantle material to lose and the difference becomes less observably significant. Adding a crust has a minimal effect by only adding a little more scatter to the wings, such that the presence or absence of a crust is least constrained by fitting to the photometry.

Zoom In Zoom Out Reset image size

Gary et al. (2017) presented lightcurves of WD 1145+017 from 105 observing sessions over six months starting on 2016 January 25. These observations were taken from multiple locations and with different telescope setups. We found a recurring feature in the April observations taken with the IAC80 telescope on the Canary Islands at a cadence of ∼68 s. This feature strongly resembled our simulated lightcurves with a single well-defined dip flanked by asymmetric wings. Other features exist that appear to be blended superpositions of dips or do not persist with a coherent shape over multiple orbits. Our simulations only account for the tidal disruption of a single rubble pile and do not consider detailed physics in the sublimation of particles making them ill-suited to describing these more complex transit features (Veras et al. 2015a; Xu et al. 2018, 2019). This particular transit may have been early enough in the disruption of the transiting body or symmetrical enough in its sublimation such that its orbital evolution did not diverge significantly from purely tidal disruption. Combined with its isolation from other features, this transit observation was ideal to demonstrate our proof-of-concept comparing differentiated disruption simulations to photometry.

We use smoothed lightcurves from each structure to fit to the 2016 April 26th transit of this feature in the data from Gary et al. (2017) and crudely approximate the sublimation of the particles by scaling the transit depth by one free parameter that describes an identical effective obscuration area for all particles. Since we cannot identify exactly when this particular feature's progenitor began disrupting, we similarly scale the transit duration by another free parameter which assumes that the transit shape is preserved as the ring fills out. Therefore our simulated lightcurves can only describe the early behavior of the transiting body, before sublimation affects the orbital evolution of individual particles significantly. All fits were performed with the affine-invariant MCMC algorithm proposed by Goodman & Weare (2010) and implemented in Python by the emcee package described in Foreman-Mackey et al. (2013). The four free parameters were the time at which the transit starts, the factor by which to scale the total duration of the transit, the amplitude by which to scale the depth of the transit, and the out of transit baseline flux.

Figure 3 shows the best-fitting structure which had a bulk density of 3 g cm⁻³, a core fraction f_c = 0.15, and a mantle/core density ratio ρ_m/ρ_c = 0.25. The five best-fitting structures all shared a small core fraction f_c = 0.15. The poorest fits have higher core fractions and/or a higher mantle/core density ratio because these parameters cannot reproduce the asymmetry and strong wings that were observed in this transit. The chi-squared statistic for the best-fit from each structure is listed in Table 2. We note that these fits cannot be used to infer structural parameters with much confidence because we do not model the relative sublimation of material from different layers or account for composition or orbital evolution due to sublimation. There are some trends in better fits correlating to lower core fraction, but the different mantle densities of the best fits show that we cannot uniquely specify the structure of an asteroid using our simulations. Nevertheless, these simulations show that over similar dynamical timescales, the internal structure of a disrupting rubble pile asteroid will impact the spatial distribution of particles from these layers which should have observable consequences for the transit lightcurve.

Zoom In Zoom Out Reset image size

The WD 1145+017 system offers (currently) unique insight into the fate of post-main-sequence planetary systems. The options for post-main-sequence planetary evolution seem largely limited to feeding circumstellar disks, accreting onto the host white dwarf, or being ejected into the interstellar medium: a possible origin for the recently discovered interstellar visitor 'Oumuamua (Veras et al. 2011; de la Fuente Marcos & de la Fuente Marcos 2017; Rafikov 2018). Pulsar planets and WD 1145+017 are the rare exceptions discovered so far although other stable post-main-sequence planetary configurations could exist (Veras 2016).

Our modeling of the photometry from this singular system suggests that the bodies creating the observed transit signatures are differentiated with a very small core fraction and low-density mantle, resembling an asteroid with a partially differentiated structure and volatile-rich mantle like Vesta (Veras et al. 2017). These bodies are low mass (${M}_{b}\lesssim {10}^{21}\,\mathrm{kg}$) compared to the total accreted material inferred by accretion rates, ∼10²³ kg (Redfield et al. 2017). Each of these fragments has a bulk density high enough to avoid being immediately disrupted, but low enough that they undergo mantle disruption and maintain coherent structure for weeks (Rappaport et al. 2016). Further work is required to determine the likelihood of such bodies being created during early planet formation versus being formed as a product of post-main-sequence planetary system evolution (Perets et al. 2011; Schleicher & Dreizler 2014; van Lieshout et al. 2018).

These suggested attributes of the asteroid's internal composition are consistent with the lack of small particles found by multiple papers (Alonso et al. 2016; Zhou et al. 2016; Croll et al. 2017; Xu et al. 2018). Xu et al. (2018) show that small grains would sublimate quickly while larger grains survive, with or without a gas-rich disk to protect them. An icy mantle with highly refractive minerals could disrupt such that the small grains sublimate to form the gas disk interior to the transiting objects while the larger grains create an extended dust cloud which produces the observed transits for long periods of time. Dynamic interactions between the debris clouds could explain the variations in transit morphology without being strong enough to disperse the dust clouds entirely. Determining when the asteroids have completely disrupted to form a largely homogeneous gas disk will place constraints on the total mass of asteroids disrupting in the system while separating the individual periodicities of the transiting objects will pin down the number of individual asteroids. If the multiple disk structure proposed by Cauley et al. (2018) continues to effectively predict the spectroscopic behavior of WD 1145+017 and we assume each individual ring is the byproduct of a single asteroid's destruction, the number and evolution of the rings (Veras et al. 2015b) may provide an independent constraint on the number of disrupting bodies in addition to stability constraints imposed by simulations like those explored in Veras et al. (2016). Follow-up simulations and observations should try to constrain this number by fitting to blends of transit features and accounting for sublimation by determining an appropriate particle size distribution evolving over the disruption process.

The chromatic dips of Boyajian's star (Boyajian et al. 2018) show some resemblance to our simulated lightcurves and may be the result of disruption without the rapid sublimation of small particles. More recent discoveries in the same vein include another transiting white dwarf reported by Vanderbosch et al. (2019) and the accretion of a gas giant envelope onto a white dwarf reported by Gänsicke et al. (2019). These are the first members of a larger class of dying planetary systems that must be studied by pairing spectroscopic and photometric observations with disruption simulations, either tidal as in WD 1145+017 or rotational as Veras et al. (2020) proposes for the body transiting ZTF J0139+5245 (Vanderbosch et al. 2019). This multipronged approach would use the death of these planetary systems in action to study fundamental properties of exoplanetary bodies that are otherwise inaccessible: a study in necroplanetology.

The authors thank the anonymous referee for constructive and helpful feedback. The authors also thank Wesleyan University for computer time supported by the NSF under grant number CNS-0619508 and CNS-0959856. We acknowledge support for this work from NASA Keck funds related to polluted white dwarfs associated with RSA #1536748. G.M.D. thanks Zachory Berta-Thompson for helpful comments during the preparation of this manuscript. D.V. acknowledges the support of the STFC via an Ernest Rutherford Fellowship (grant ST/P003850/1).

Software: astropy (Astropy Collaboration et al. 2018), matplotlib (Hunter 2007), NumPy (Oliphant 2006), REBOUND (Rein & Liu 2011), emcee (Foreman-Mackey et al. 2013).