Speeding past planets? Asteroids radiatively propelled by giant branch Yarkovsky effects

TL;DR

This paper investigates how intense stellar luminosity during the giant branch phase affects asteroid dynamics through Yarkovsky effects, influencing their fate and distribution around white dwarfs.

Contribution

It models the impact of giant star luminosity on exo-asteroid migration, highlighting the importance of advanced radiative models for understanding asteroid behavior near white dwarfs.

Findings

Asteroids move too rapidly to be captured by planets during giant branch phase.

Yarkovsky effects can cause orbital inclination flips without semimajor axis change.

Resonant trapping of small asteroids by planets is generally insignificant.

Abstract

Understanding the fate of planetary systems through white dwarfs which accrete debris crucially relies on tracing the orbital and physical properties of exo-asteroids during the giant branch phase of stellar evolution. Giant branch luminosities exceed the Sun's by over three orders of magnitude, leading to significantly enhanced Yarkovsky and YORP effects on minor planets. Here, we place bounds on Yarkovsky-induced differential migration between asteroids and planets during giant branch mass loss by modelling one exo-Neptune with inner and outer exo-Kuiper belts. In our bounding models, the asteroids move too quickly past the planet to be diverted from their eventual fate, which can range from: (i) populating the outer regions of systems out to 10^4-10^5 au, (ii) being engulfed within the host star, or (iii) experiencing Yarkovsky-induced orbital inclination flipping without any Yarkovsky-induced semimajor axis drift. In these violent limiting cases, temporary resonant trapping of asteroids with radii of under about 10 km by the planet is insignificant, and capture within the planet's Hill sphere requires fine-tuned dissipation. The wide variety of outcomes presented here demonstrates the need to employ sophisticated structure and radiative exo-asteroid models in future studies. Determining where metal-polluting asteroids reside around a white dwarf depends on understanding extreme Yarkovsky physics.

Equations

$$\alpha \equiv \frac{a_{\mathrm{inner}}}{a_{\mathrm{outer}}}={\left(\frac{p}{p+1}\right)}^{2/3}<1.$$

$$\frac{d{n}_{\mathrm{crit}}}{dt}\approx \frac{5\cdot 3^{5/3}}{8}{\left|p^2\left(p+1\right)\right|}^{\frac{1}{9}}{n}_{\mathrm{pl}}^2{\left|\frac{M_{\mathrm{pl}}}{M_{\star }}{f}_p\right|}^{\frac{4}{3}}.$$

$$\frac{d{a}_{\mathrm{crit}}^{\mathrm{exterior}}}{dt}\approx -\frac{5\cdot 3^{2/3}}{4}\sqrt{\frac{G{M}_{\star }}{a_{\mathrm{pl}}}}{\left(p+1\right)}^{16/9}{p}^{-13/9}{\left|\frac{M_{\mathrm{pl}}}{M_{\star }}{f}_p\right|}^{\frac{4}{3}}.$$

$$\frac{d{a}_{\mathrm{crit}}^{\mathrm{interior}}}{dt}\approx \frac{5\cdot 3^{2/3}}{4}\sqrt{\frac{G{M}_{\star }}{a_{\mathrm{pl}}}}{\left(p+1\right)}^{-14/9}{p}^{17/9}{\left|\frac{M_{\mathrm{pl}}}{M_{\star }}{f}_p\right|}^{\frac{4}{3}}.$$

$$e_{\mathrm{crit}}=\frac{2^{1/2}{3}^{1/6}}{{\left|p^2\left(p+1\right)\right|}^{\frac{2}{9}}}{\left|\frac{M_{\mathrm{pl}}}{M_{\star }}{f}_p\right|}^{\frac{1}{3}},$$

$$f_p^{\mathrm{interior}}\equiv f_{27}+f_{\mathrm{e}}$$

$$f_{27}=-\left(p+1+\frac{1}{2}\frac{d}{d\log \alpha }\right)b_{1/2}^{\left(p+1\right)}$$

$$f_p^{\mathrm{interior}}=-\left(p+1\right)b_{1/2}^{\left(p+1\right)}-\frac{\alpha }{4}\left(b_{3/2}^{\left(p\right)}-2\alpha b_{3/2}^{\left(p+1\right)}+b_{3/2}^{\left(p+2\right)}\right).$$

$$f_p^{\mathrm{exterior}}\equiv \alpha f_{31}+f_{\mathrm{i}}$$

$$\alpha f_{31}=\alpha \left(p+\frac{3}{2}+\frac{1}{2}\frac{d}{d\log \alpha }\right)b_{1/2}^{\left(p\right)}$$

$$f_p^{\mathrm{exterior}}=\alpha \left(p+\frac{3}{2}\right)b_{1/2}^{\left(p\right)}-\frac{{\delta }_{p,1}}{2}$$

$$+\frac{{\alpha }^2}{4}\left(b_{3/2}^{\left(p-1\right)}-2\alpha b_{3/2}^{\left(p\right)}+b_{3/2}^{\left(p+1\right)}\right).$$

$$b_s^{\left(m\right)}=\frac{1}{\pi }{\int }_0^{2\pi }\frac{\cos \left(mx\right)dx}{{\left(1-2\alpha \cos x+{\alpha }^2\right)}^s}$$

$$\frac{d}{d\log \alpha }{b}_s^{(j)}=s\alpha \left(b_{s+1}^{\left(j-1\right)}-2\alpha b_{s+1}^{\left(j\right)}+b_{s+1}^{\left(j+1\right)}\right).$$

$$r_{\mathrm{Hill}}=a_{\mathrm{pl}}{\left(\frac{M_{\mathrm{pl}}}{3M_{\star }}\right)}^{1/3}.$$

$$a_{\mathrm{tc}}=a_{\mathrm{pl}}{\left[\frac{2}{A}-{\left(1\pm \sqrt{3\kappa }{r}_{\mathrm{Hill}}\right)}^2\right]}^{-1}$$

$$e_{\mathrm{tc}}=1-\frac{a_{\mathrm{pl}}A}{a_{\mathrm{tc}}}\mathrm{or}{e}_{\mathrm{tc}}=\frac{a_{\mathrm{pl}}A}{a_{\mathrm{tc}}}-1.$$

$$\frac{d{a}_{\mathrm{ast}}}{dt}\propto \Lambda a_{\mathrm{ast}}^{-1/2}$$

$$\frac{d{a}_{\mathrm{ast}}^{\mathrm{exterior}}}{dt}\approx {\Lambda }_0{\left(p+1\right)}^{-1/3}{p}^{1/3}\sqrt{\frac{G{M}_{\ast }}{a_{\mathrm{pl}}}}\left[\mathrm{au}(\mathrm{yr}/2\pi {)}^{-1}\right],$$

$$\frac{d{a}_{\mathrm{ast}}^{\mathrm{interior}}}{dt}\approx {\Lambda }_0{\left(p+1\right)}^{1/3}{p}^{-1/3}\sqrt{\frac{G{M}_{\ast }}{a_{\mathrm{pl}}}}\left[\mathrm{au}(\mathrm{yr}/2\pi {)}^{-1}\right].$$

$${\Lambda }_0=\frac{5\cdot 3^{2/3}}{4}{\left(p_{\mathrm{crit}}+1\right)}^{19/9}{p}_{\mathrm{crit}}^{-16/9}{\left|\frac{M_{\mathrm{pl}}}{M_{\star }}\left(\alpha f_{31}+f_{\mathrm{i}}\right)\right|}^{\frac{4}{3}}.$$

$${\Lambda }_0=\frac{5\cdot 3^{2/3}}{4}{\left(p_{\mathrm{crit}}+1\right)}^{-17/9}{p}_{\mathrm{crit}}^{20/9}{\left|\frac{M_{\mathrm{pl}}}{M_{\star }}\left(f_{27}+f_{\mathrm{e}}\right)\right|}^{\frac{4}{3}}.$$

$${\left(\frac{d{\vec{v}}_{\mathrm{ast}}}{dt}\right)}_{\mathrm{Yar}}=\frac{A_{\mathrm{ast}}{L}_{\star }(t)}{4\pi M_{\mathrm{ast}}c{r}_{\mathrm{ast}}^2}\left[Q_{\mathrm{PR}}\mathbb{I}+k{Q}_{\mathrm{Yar}}\mathbb{Y}(t)\right]\vec{\iota }$$

$$\vec{\iota }\equiv \left(1-\frac{{\vec{v}}_{\mathrm{ast}}\cdot {\vec{r}}_{\mathrm{ast}}}{c{r}_{\mathrm{ast}}}\right)\frac{{\vec{r}}_{\mathrm{ast}}}{r_{\mathrm{ast}}}-\frac{{\vec{v}}_{\mathrm{ast}}}{c},$$

$${\left(\frac{d{\vec{v}}_{\mathrm{ast}}}{dt}\right)}_{\mathrm{Yar}}=\frac{R_{\mathrm{ast}}^2{L}_{\star }(t)}{4M_{\mathrm{ast}}c{r}_{\mathrm{ast}}^2}\left[\mathbb{I}+\frac{1}{4}\mathbb{Y}\right]\vec{\iota }.$$

$$\mathbb{Q}\equiv \mathbb{I}+\frac{1}{4}\mathbb{Y}.$$

$$⟨\frac{d{a}_{\mathrm{ast}}}{dt}⟩\propto \frac{{\mathbb{Q}}_{21}-{\mathbb{Q}}_{12}}{\sqrt{a_{\mathrm{ast}}}\left(1-e_{\mathrm{ast}}^2\right)},$$

$$⟨\frac{d{e}_{\mathrm{ast}}}{dt}⟩\to 0,$$

$$⟨\frac{d{i}_{\mathrm{ast}}}{dt}⟩\propto \frac{\sin i\left({\mathbb{Q}}_{12}{\sin }^2{\Omega }_{\mathrm{ast}}-{\mathbb{Q}}_{21}{\cos }^2{\Omega }_{\mathrm{ast}}\right)}{a_{\mathrm{ast}}^{3/2}},$$

$$\frac{d^2{\vec{r}}_{\mathrm{pl}}}{d{t}^2}=-G{M}_{\star }(t)\frac{{\vec{r}}_{\mathrm{pl}}}{{\left|{\vec{r}}_{\mathrm{pl}}\right|}^3},$$

Full text

††pubyear: 2019††pagerange: Speeding past planets? Asteroids radiatively propelled by giant branch Yarkovsky effects–Speeding past planets? Asteroids radiatively propelled by giant branch Yarkovsky effects

1 Introduction

Upon leaving the main sequence and ascending the giant branch, a star which hosts a planetary system will subject its orbiting constituents to a violent cocktail of destructive forces. These include (i) physical expansion of the star’s envelope to au-scale distances, (ii) shedding of mass through superwinds, and (iii) luminosities which reach between $10^{3}L_{\odot}$ and $10^{4}L_{\odot}$. The dynamical and physical consequences for giant planets, terrestrial planets, minor planets, boulders, pebbles and dust are only starting to be characterised (Veras, 2016a) at a sufficiently-detailed level to connect with the abundant data available in white dwarf planetary systems.

1.1 Giant branch effects

An expanding giant branch star can engulf objects which orbit too closely. Determining the critical engulfment distance has already been the subject of intense scrutiny (Villaver & Livio, 2009; Kunitomo et al., 2011; Mustill & Villaver, 2012; Adams & Bloch, 2013; Nordhaus & Spiegel, 2013; Li et al., 2014; Villaver et al., 2014; Madappatt et al., 2016; Staff et al., 2016; Gallet et al., 2017), and is strongly dependent on both the tidal prescription and stellar model that one adopts. In general, however, this critical engulfment distance resides at a few au, and is greater for gas giants than for terrestrial planets. In the solar system, Mercury and Venus will easily be engulfed, Mars will not, and the outcome for the Earth remains unclear (Schröder & Connon Smith, 2008).

Giant branch stellar mass loss will precipitate orbital changes for all types of planets. In general, even for isotropic mass loss, a planet’s eccentricity, semimajor axis and argument of pericentre will vary (Omarov, 1962; Hadjidemetriou, 1963)111Anisotropic mass loss generates changes in all orbital elements (Veras et al., 2013a; Dosopoulou & Kalogera, 2016a, b).. However, for planets and asteroids within a few hundred au of their parent star, the adiabatic approximation may be employed (Veras et al., 2011), where orbital eccentricity is effectively conserved and semimajor axis expansion scales simply with stellar mass loss. Nevertheless, despite the invariance of the semimajor axis ratios of multiple objects in the adiabatic approximation, the central mass alterations can incite instability (Debes & Sigurdsson, 2002; Bonsor et al., 2011; Debes et al., 2012; Veras et al., 2013b; Voyatzis et al., 2013; Frewen & Hansen, 2014; Mustill et al., 2014; Veras & Gänsicke, 2015; Veras et al., 2016; Veras, 2016b; Mustill et al., 2018; Smallwood et al., 2018), which is manifest both along the giant branch and white dwarf phases. In the solar system, the giant planets are separated sufficiently far from one another to prevent post-main-sequence instability, unless there exists an additional Planet Nine-like companion (Veras, 2016b).

The time-varying giant branch luminosities affect planets and asteroids differently. Planetary atmospheres are subjected to partial or complete evaporation (Livio & Soker, 1984; Nelemans & Tauris, 1998; Soker, 1998; Villaver & Livio, 2007; Wickramasinghe et al., 2010), and their geology and surface processes are likely exposed to (as-yet-unmodelled) transformative events. Asteroids, however, are small enough to be exposed to interior water depletion (Malamud & Perets, 2016, 2017a, 2017b) and being radiatively “pushed” by the Yarkovsky effect (Veras et al., 2015) and spun-up through the YORP effect (Veras et al., 2014a). These latter two phenomena have been observed in the solar system (e.g. see Fig. 1 of Polishook et al., 2017) but are enhanced during the giant branch phases. In fact, the YORP effect could easily destroy 100m - 10km asteroids within about 7 au of their parent star through rotational fission. Giant branch luminosity also creates a drag force for pebbles and boulders (Dong et al., 2010), which becomes important for objects which are too small for the Yarkovsky effect to play a role in determining the final orbital state (Veras et al., 2015).

1.2 White dwarf planetary systems

After the star has become a white dwarf, the mutual perturbations amongst the remaining planets, asteroids, pebbles and fragments conspire to thrust enough material into the white dwarf photosphere to be observable in one-quarter to one-half of all observed white dwarfs (Zuckerman et al., 2003, 2010; Koester et al., 2014). One of these exo-asteroids is currently in the process of disintegrating, a phenomenon which has been seen in real time on a nightly basis for the last several years (Vanderburg et al., 2015). Such disintegration produces debris discs (Jura, 2003, 2008; Debes et al., 2012; Veras et al., 2014b), of which about 40 are now known (Farihi, 2016). Another exo-asteroid, one with high internal strength, has been found embedded inside one of these discs (Manser et al., 2018). The importance of these minor bodies is highlighted by how their innards after break-up are regularly observed in white dwarf photospheres (Klein et al., 2010, 2011; Gänsicke et al., 2012; Jura et al., 2012; Xu et al., 2013, 2014; Jura & Young, 2014; Wilson et al., 2015, 2016; Xu et al., 2017; Harrison et al., 2018; Hollands et al., 2018). These observations provide the most direct and extensive known measurements of the bulk chemical composition of the building blocks of exoplanets.

1.3 Plan for paper

Here, we seek to achieve a better understanding of the post-main-sequence orbital distribution of exo-asteroids by (i) considering their Yarkovsky-induced interactions with a planet during the giant branch phase of evolution, and (ii) applying numerical Yarkovsky models which can place bounds on possible behaviours. Because Yarkovsky acceleration affects only exo-asteroids and not exo-planets, the resulting relative orbital evolution at first glance might appear similar to convergent or divergent migration between two bodies within a protoplanetary disc. In Section 2, we detail the limitations of this analogy, which provides context for our modelling efforts. Sections 3 and 4 then, respectively, establish our numerical simulations and reports on their outcomes. We discuss our findings in Section 5, and summarize in Section 6.

2 Analogy with protoplanetary disc migration

In this section, we explore the possibility that the giant branch-induced Yarkovsky drift can be modelled, or, at least contained, by the formalisms which have been developed for migration within protoplanetary discs.

2.1 Capture into resonance

The architectures of planetary systems on the main-sequence arise from the movement of and interactions between protoplanets within their birth disc. The prevalence of multiple planets which are currently observed to reside close to a mean motion commensurability suggest that they drifted within the disc into such configurations. This idea has been borne out by dozens of theoretical studies, which have characterised the relative drift as “convergent migration” or “divergent migration”.

2.1.1 Previous modelling efforts

In order to investigate this type of migration, several authors have solved a set of differential equations involving $N$-body point-mass interactions, plus some perturbative accelerations (e.g. equations 17-19 of Terquem & Papaloizou 2007, equation 7 of Ogihara & Ida 2012, equation 1 of Ogihara & Kobayashi 2013, and equations 54-58 of Teyssandier & Terquem 2014). These equations do not actually incorporate the potential of the disc, but rather damping timescales which mimic dissipative disc effects. These timescales range from basic linear approximations (Lee & Peale, 2002) to empirical fits based on recent hydrodynamic simulations (Cresswell & Nelson, 2008; Xu et al., 2018).

2.1.2 Types of resonances

Capture models have focused almost exclusively on eccentricity-based first and second-order mean motion resonances, assuming that at least one planet is on a near-circular orbit222Notable exceptions are Luan (2014), who quantified the probability that the $4$:$2$ inclination-based mean motion resonance between Mimas and Tethys was the consequence of convergent migration, and El Moutamid, et al. (2017), who considered capture of massless particles in the first-order corotation eccentric resonance.. Assume that a first-order resonance is described as a $p+1$:$p$ mean motion resonance, where $p$ is a positive integer. As the value of $p$ increases, and the configuration approaches a co-orbital state, the probability of capture decreases. Nevertheless, Tadeu dos Santos et al. (2015) illustrated how planets can bypass – or migrate through – the $2$:$1$ and $3$:$2$ mean motion resonances to become trapped in the $4$:$3$ mean motion resonance. Similarly, figures 9 and 7 of Papaloizou & Szuszkiewicz (2005) respectively illustrated how entrapment into the $6$:$5$ and $9$:$8$ mean motion resonances may be possible; Quillen, et al. (2013) warned, however, that capture into a mean motion resonance such as the $7$:$6$ requires fine tuning of parameters. Known moons represent good examples of bodies which may be trapped in high-$p$ mean motion resonances: Desdemona and Portia are near the $13$:$12$ mean motion resonance and Cressida and Desdemona are near the $47$:$46$ mean motion resonance (Quillen & French, 2014).

2.1.3 Conditions for capture

Only if the migration is convergent, slow enough, minimally stochastic, and relies on near-circular orbits can the planets become trapped in mean motion resonances. Each of these criteria has already been significantly vetted in the literature, and were bolstered by recent revisitations of mean motion resonance capture theory (Quillen, 2006; Mustill & Wyatt, 2011; Goldreich & Schlichting, 2014; Batygin, 2015; Deck & Batygin, 2015). Turbulence, Brownian motion and stochasticity in general may prevent resonant capture from occurring (Murray-Clay & Chiang, 2006; Adams et al., 2008; Rein & Papaloizou, 2009, 2010). Orbital eccentricities which exceed a few hundredths also can prevent capture, a result which has been demonstrated by hydrodynamical simulations (Hands & Alexander, 2018). A final barrier to capture is migration speed: generally, if a planetary body crosses the libration width of a resonance faster than the libration timescale, then capture can be avoided (e.g. Pan & Schlichting, 2017). However, capture is probabilistic, and largely depends on angular variables and their time evolution (e.g. Folonier et al., 2014).

The appendix of Petrovich et al. (2013) provides particularly useful equations in physical units for the critical migration speed and eccentricity needed to avoid certain capture for first-order $p+1$:$p$ resonances. They assume one body is a planet and another is a test particle, and use $p$ in a different way than us, allowing the integer to take on positive or negative values for resonances where their test particle (e.g. asteroid) is respectively interior or exterior to the planet. We instead use positive values for both interior and exterior resonances, such that at a resonance the ratio $\alpha$ of the inner to outer semimajor axes is

[TABLE]

Petrovich et al. (2013) illustrated in their appendix that the critical speed beyond which capture may be avoided in terms of the planet’s mean motion $n_{\rm pl}$ relative to the asteroid, or vice versa, is

[TABLE]

Correspondingly, if the asteroid is external to the planet, then its critical decay rate is

[TABLE]

Instead, if the asteroid is internal to the planet, then its critical growth rate is

[TABLE]

For either type of resonance, the critical asteroid eccentricity beyond which capture may be avoided is

[TABLE]

where we have assumed a numerical coefficient in equation (2) from Friedland (2001). In the equations, $M_{\rm pl}$ is the mass of the planet, $M_{\star}$ is the mass of the star, $a_{\rm pl}$ is the semimajor axis of the planet, and the quantity $f_{p}$ is a function extracted from Appendix B of Murray & Dermott (1999) that depends on whether the asteroid is interior or exterior to the planet.

For first-order resonances, in the interior case, where there is an asteroid internal to the planet,

[TABLE]

where $f_{\rm e}=0$ and

[TABLE]

such that

[TABLE]

In the exterior case, where there is an asteroid external to the planet,

[TABLE]

where $f_{\rm i}=-\delta_{p,1}/2$ and

[TABLE]

such that

[TABLE]

[TABLE]

In the above relations, the Laplace coefficient is given by

[TABLE]

and satisfies the relation

[TABLE]

The conditions for capture (equations 2-5) are sufficient, but not necessary. If capture is avoided, continuing migration towards the planet will allow the asteroid to encounter an infinitely increasing number of $p+1$:$p$ resonances. For a high enough value of $p$, the gravitational reaches of these individual resonances cross one another, creating “resonant overlap”.

2.2 Capture into the planet’s Hill sphere

This resonant overlap can create gravitational instability, resulting in the asteroid’s escape from the system, or collision into the planet or star. The width of the overlap region has been estimated analytically and numerically (Wisdom, 1980; Quillen & Faber, 2006; Mustill & Wyatt, 2012; Bodman & Quillen, 2014; Shannon, et al., 2015), but is not yet known to admit an exact analytical solution.

If an asteroid is migrating towards a planet, avoids entrapment into a resonance during migration, and passes through the resonant overlap region without being scattered away, then the planet may capture the asteroid, at least temporarily. The conditions for capture have been quantified by Higuchi & Ida (2016, 2017), who specify the combination of $a$ and $e$ values required for temporary capture by a planet. Denote these values as $a_{\rm tc}$ and $e_{\rm tc}$. They are given as curves on the $a-e$ plane that satisfy $q=a_{\rm tc}\left(1\pm e_{\rm tc}\right)=r_{\rm pl}\pm r_{\rm Hill}$, where $q$ is the periastron distance, $r_{\rm pl}$ is the heliocentric distance of the planet, and

[TABLE]

First $a_{\rm tc}$ is computed through

[TABLE]

where either $A=1-r_{\rm Hill}/a_{\rm pl}$, when $r_{\rm ast}<r_{\rm pl}$ at the L1 point, or where $A=1+r_{\rm Hill}/a_{\rm pl}$, when $r_{\rm ast}>r_{\rm pl}$ at the L2 point. Also, $0\leq\kappa\leq 2$, which spans the allowed capture region. Then, $e_{\rm tc}$ is computed through

[TABLE]

2.3 Link with giant branch Yarkovsky forces

Now we link the resonant formalism above with the consequences of differential migration due to Yarkovsky forces.

Drifts from Yarkovsky forces are unlikely to produce the type of smooth migration that is often envisaged in a protoplanetary disc. The shape, size and thermal properties of an asteroid determine the migration rate (e.g. Vokrouhlicky, 1998; Brož, 2006; Wang & Hou, 2017), may change with time, and are linked with YORP spin changes. Despite these complexities, a common simplification for Yarkovsky-based studies of solar system asteroids is to assume that only their semimajor axes change. This Solar-induced drift for asteroids whose sizes are within the metre-to-kilometre range is on the order of $10^{-5}-10^{-2}$ au/Myr (e.g. Bottke et al., 2000, 2006; Gallardo et al., 2011).

2.3.1 Constant drift approximation

We can now derive some results for the current solar system by using the well-established constant drift approximation. Doing so also helps lay the groundwork for later considering the giant branch phases of evolution.

Veras et al. (2015) showed in their Eq. (A2) that the averaged secular drift due to the Yarkovsky effect has the following dependency on semimajor axis

[TABLE]

where $\Lambda$ is a factor given (or assumed) that is a function of asteroid density, size and many other parameters. For an asteroid with a radius of $R_{\rm ast}=100$ m and density of $3$ g/cm3 in the current solar system (inhabited by a $1L_{\odot}$ Sun), we obtain $\Lambda\approx 66$ m3/2/s. This value can be related to a dimensionless constant $\Lambda_{0}=2.3\times 10^{-7}$ where

[TABLE]

[TABLE]

Scaling $\Lambda_{0}=2.3\times 10^{-7}$ to different values allows us to probe different magnitudes of the Yarkovsky effect easily.

The asteroid’s migration is guaranteed to stop at the first 1st-order resonance encounter after (and if) both $da/dt<da_{\rm crit}/dt$ and $e<e_{\rm crit}$ are satisfied. In order to find the value of $p_{\rm min}$ which is the minimum that satisfies both conditions, we set equations (17-18) equal to equations (3-4), and then implicitly solve for $p=p_{\rm crit}$. For external asteroids,

[TABLE]

For internal asteroids,

[TABLE]

Then $p_{\rm min}$ is the nearest integer greater than or equal to $p_{\rm crit}$ and $a_{\rm final}$ is given by equation (1) as a function of $p_{\rm min}$ and $a_{\rm pl}$. We can then obtain the critical eccentricity (equation 5) by substituting $p$ with $p_{\rm min}$.

Figures 1 and 2 quantify these analytical results. Figure 1 demonstrates that for temporary Hill sphere capture to occur around a $10^{-5}M_{\star}$ planet, the Sun’s current luminosity would need to be increased by at least a factor of 100. An asteroid with a density of 3 g/cm3 and a radius of $R_{\rm ast}=100$ m would also need to maintain an eccentricity of under about 0.06 throughout its migration. Hill sphere capture could occur from either internal or external migration. Figure 2 then quantifies this migration boundary as a function of both planetary mass and stellar luminosity, and illustrates that temporary Hill sphere capture may occur around any type of planet provided that the stellar luminosity is sufficiently high. Such capture, however, requires orders-of-magnitude higher luminosities for giant planets than for terrestrial planets.

2.3.2 Variable drift approximation

However, applying a constant drift in only the semimajor axis throughout the giant branch phases of evolution would be unrealistic. The stellar luminosity varies non-monotonically and by three orders of magnitude while the asteroid’s orbit is expanding anyway through stellar mass loss. Further, the asteroid’s eccentricity and inclination do not remain constant, the former potentially changing by up to 0.08/Myr (Veras et al., 2015). Therefore, whether asteroid eccentricities can remain under the critical eccentricity that is given by equation (5) is unclear, as is the consequences of inclination evolution. Although at the beginning of the red giant and asymptotic giant branch phases the luminosity is low enough to potentially achieve resonant capture, the chances of capture would likely decrease towards the end of each subphase.

Modelling these complexities remains a challenge. In the forthcoming section, we pursue an extension to previous work, but one simple enough to be computationally feasible.

3 Simulating Yarkovsky drift

Having now detailed the potential limitations of the analogy with planetary migration in a disc, we set about placing bounds on Yarkovsky-induced movement of an asteroid in the vicinity of a planet during the giant branch phase.

Asteroids are orbitally perturbed from stellar radiation which is (i) absorbed, (ii) immediately reflected, and (iii) reflected after a delay. As described in detail by Veras et al. (2015), the first two processes give rise to what is known as Poynting-Robertson drag and radiation pressure. The third component gives rise to the Yarkovsky effect, which “turns on” only for objects larger than about 1-10 m in size.

Crucially, Veras et al. (2015) showed that acceleration due to Yarkovsky effect is proportional to $(1/c)$, whereas acceleration due to Poynting-Robertson drag and radiation pressure is proportional to $(1/c^{2})$, where $c$ is the speed of light. Consequently, contributions from Poynting-Robertson drag and radiation pressure are negligible in this context.

3.1 The true Yarkovsky acceleration

The acceleration due to the Yarkovsky effect is (equation 27 of Veras et al. 2015)

[TABLE]

where $\vec{v}_{\rm ast}$ is the velocity of the asteroid, $k$ is a constant between 0 and 1/4 that is strongly linked to the asteroid’s rotational state, $A_{\rm ast}$ is the momentum-carrying cross-sectional area of the asteroid’s heated surface, $L_{\star}(t)$ is the stellar luminosity (with an emphasis on time dependence), $Q_{\rm PR}$ is a constant that equals the sum of the target’s absorption efficiency and reflecting efficiency, and $Q_{\rm Yar}$ is a constant between 0 and 1 that refers to the difference between the asteroid’s absorption efficiency and albedo. The matrix $\mathbb{I}$ is the identity matrix, and $\mathbb{Y}$ is the Yarkovsky matrix, where the absolute value of each entry is less than or equal to unity. The vector $\iota$ is the relativistic direction correction:

[TABLE]

The greatest complexity in equation (21) arises from the position-, velocity- and time-dependent 3-by-3 matrix $\mathbb{Y}$, which is a function of the asteroid’s spin axis, specific angular momentum, emissivity, specific heat capacity and thermal conductivity. The particular functional forms and resulting matrices are in turn dependent on the adopted heat conduction model. The spin evolution is further dictated by the YORP effect, which is a function of the asteroid’s shape.

3.2 A simplified Yarkovsky acceleration

Modelling all of these complications is well beyond the scope of this study. Instead, we pursue a simpler approach, but one which is still more detailed than any previous giant branch Yarkovsky study (Veras et al., 2015): here we place limits on the resultant motion by adopting different constant entries for $\mathbb{Y}$. This matrix contains all of the Yarkovsky physics, and by setting its entries to their extreme values, we can bound possible asteroid motions.

Therefore, we set $Q_{\rm PR}=Q_{\rm Yar}=1$ and $k=1/4$, and assume that the entries for $\mathbb{Y}$ are constant throughout the simulation. We also assume the asteroid is a perfect sphere with radius $R_{\rm ast}$ (thereby preventing initiation of the YORP effect), such that

[TABLE]

The functional form in equation (23) highlights the strong dependence of Yarkovsky acceleration on asteroid radius.

3.2.1 Secular trends

Before choosing values of $\mathbb{Y}$ and running simulations, we can obtain some analytic results based on the simplified acceleration in equation (23). Veras et al. (2015) analytically integrated a similar acceleration, and their appendix presents the resulting secular evolution of orbital elements333In each of their Eqs. A4-A5, the first dot refers to matrix multiplication, and the second dot refers to dot product. In this context, for those equations the transpose of the $\mathbb{Q}$ should be applied instead of the traditional form of $\mathbb{Q}$.. In order to relate our results to theirs, we set

[TABLE]

Then, we can observe several important trends, which will be helpful for understanding our results:

• •

The Yarkovsky-induced semimajor axis evolution $\left\langle da_{\rm ast}/dt\right\rangle$ is completely independent of $\mathbb{Q}_{11},\mathbb{Q}_{22}$ and $\mathbb{Q}_{33}$.

• •

The Yarkovsky-induced eccentricity evolution $\left\langle de_{\rm ast}/dt\right\rangle$ does not tend towards infinity for any value of $e_{\rm ast}$.

• •

In the limit of circular, coplanar asteroid orbits,

[TABLE]

[TABLE]

• •

For non-coplanar, near-circular orbits where the $\mathbb{Q}_{21}$ and/or $\mathbb{Q}_{12}$ terms dominate, as in equation (25), then

[TABLE]

where $\Omega_{\rm ast}$ represents the longitude of ascending node of the asteroid. These trends will help set initial conditions for our simulations and explain our results as we sample limiting cases of motion.

3.2.2 Equations of motion

The Yarkovsky-induced acceleration on the asteroid must be combined with the acceleration arising from stellar mass loss and the interaction with the planet. The planet is also accelerating due to the mass loss, but is negligibly perturbed by Yarkovsky acceleration due to its size.

As explained by Hadjidemetriou (1963) and Veras et al. (2013a), this acceleration from mass loss is automatically taken into account through the equations of motion through the change of stellar mass. The final equations of motion, assuming that the asteroid does not perturb the planet, are

[TABLE]

[TABLE]

[TABLE]

We have chosen to fix the mass of the planet, which may not be a suitable approximation for planets which are primarily composed of an atmosphere and reside close enough to the giant star to undergo partial or complete evaporation (Livio & Soker, 1984; Nelemans & Tauris, 1998; Soker, 1998; Villaver & Livio, 2007; Wickramasinghe et al., 2010). Determining different evaporation prescriptions due to planet core, mantle and atmosphere compositions, as well as the type of flux emitted from the star, is beyond the scope of this study (but represents an important future consideration).

3.2.3 Integrator

In order to integrate the equations of motion, we inserted the simplified Yarkovsky acceleration (equation 23) into the combined stellar and planetary evolution code that was first used in Mustill et al. (2018). This code is a modification of the code established in Veras et al. (2013a), and is originally based on the Mercury integration suite (Chambers, 1999). The code utilises the RADAU integrator to dynamically evolve the system, and interpolates stellar mass and radius at subdivisions of each timestep. The stellar profiles of mass, radius and luminosity are obtained from the SSE stellar evolution code (Hurley et al., 2000).

There is a balance between integrator speed and accuracy. Because we modelled the possibility of entrapment into resonance, we required a high-enough accuracy to track orbital angles to within a few degrees over a few Myr. Figure A1 of Mustill et al. (2018) suggests that an accuracy of $10^{-11}$ suits our purposes, a value that we adopted.

3.3 Initial conditions

Our limited computational resources prevented us from performing a broad sweep of phase space. We instead focussed on the most relevant cases, and considered

• •

Only the asymptotic giant branch phase of an initially $2M_{\odot}$ star. This phase lasts for about 1.71 Myr, and our integrations ran for up to 2 Myr from the start of that stellar phase. At that initial time, the stellar mass had already been reduced only by about $1.2\times 10^{-4}M_{\odot}$. Throughout this stellar phase, the star lost about 68 per cent of its mass in route to becoming a white dwarf. This mass loss corresponds to an adiabatic orbital semimajor axis increase of a factor of about 3.1.

• •

A Neptune-mass planet located at 30 au from the host star at the beginning of the simulation. Its eccentricity and inclination were set initially at exactly zero to provide a reference for the asteroid evolution. We chose a Neptune-mass planet because exo-Kuiper belt objects are thought to be the primary source of white dwarf pollution. Hence, a rough solar system analogue was appropriate as a first guess. Other investigations have shown, however, that lower-mass planets are more efficient polluters (Frewen & Hansen, 2014; Mustill et al., 2018), whereas Fig. 2 of this paper reveals that higher-mass planets are more efficient resonant trappers.

• •

Asteroids in belts both internal and external to the planet. In all simulations, these asteroids were given initial semimajor axes selected from a random uniform distribution ranging from both 10-25 au and 35-50 au. The inner boundary of 10 au is far enough away from the star to ensure that at least some asteroids at those locations would survive YORP-induced rotational fission (Veras et al., 2014a), which we do not model here.

• •

200 asteroids per simulation. The asteroids were given a density of 3 g/cm3 and radii ranging from 100 m to 1000 km. These asteroids were treated as “small” particles in our integrator, a simulation designation which indicates that the asteroids do not perturb each other although they do perturb the planet (thereby introducing additional, but negligible terms in equation 28).

• •

Asteroids which were initialised at the beginning of the simulation with inclinations chosen from a random uniform distribution ranging from $0^{\circ}$ to $10^{\circ}$.

• •

Asteroids which were initialised at the beginning of the simulation with two different eccentricity distributions depending on the simulation. The first is circular asteroids, and the second are asteroids with eccentricities chosen between $0.3$ to $0.7$. We emphasise that we are not attempting to reproduce realistic exo-debris discs in these systems (which would require a more detailed Yarkovsky model), but rather are placing bounds and observing trends. The evolution of most initially eccentric asteroids that we simulated are not shown, unless they provide a revealing physical trait.

• •

Asteroids which were initialised at the beginning of the simulation with orbital angles chosen from a random uniform distribution across their entire ranges.

• •

Yarkovsky prescriptions according to the following values of $\mathbb{Q}$:

[TABLE]

[TABLE]

[TABLE]

Our choices for Models A-C arise from our intention to place bounds on the motion by considering only extreme cases. We approximate these extreme cases by focusing in on $\mathbb{Q}_{12}$ and $\mathbb{Q}_{21}$ from the analytic limits in equation (25). In terms of semimajor axis, Model A pushes asteroids outward, Model B pushes asteroids inward, and Model C does not push them outward nor inward. Alternatively, Models A and B have only a minor effect on asteroid inclination evolution, whereas Model C has a relatively large effect. We actually ran additional simulations with other combinations of $\mathbb{Q}$ entries, but found that they represent intermediate cases to Models A-C, supporting the analytics and our choices.

4 Simulation results

The asteroids in our simulations are accelerated by three major effects: (i) mass loss from the star, (ii) perturbations from the planet, and the (iii) radiation-driven Yarkovsky effect444Recall that Poynting-Robertson drag is weaker than the Yarkovsky effect by a factor of $\left(1/c\right)$ and is anyway included our formalism through $\mathbb{Q}_{11}$, $\mathbb{Q}_{22}$, and $\mathbb{Q}_{33}$.. Mass loss always pushes the asteroids away from the star (the orbital pericentre is always increasing, even in nonadiabatic motion; Veras et al. 2011), whereas the planet and the Yarkovsky effect may perturb the asteroids in either direction. Although both mass loss and planets can increase $e_{\rm ast}$ to unity (causing escape; e.g. Veras et al. 2011; Adams et al. 2013), Yarkovsky acceleration cannot. Veras et al. (2015) investigated the relative magnitudes of stellar mass loss and the Yarkovsky effect, but did not include a planet nor run numerical simulations.

The Yarkovsky drift varies in a complex manner, and is a function of all orbital elements, even in our simple approximation where the matrix elements of $\mathbb{Q}$ remain constant over time. However, equations (25)-(26) reveal that for small inclinations, the key parameters are $\mathbb{Q}_{12}$ and $\mathbb{Q}_{21}$. Although the Yarkovsky acceleration drops roughly as $1/\sqrt{a_{\rm pl}}$, the acceleration increases linearly as a function of $L_{\star}(t)$. The time dependence of $L_{\star}(t)$ is, in turn, nonlinear, and dominates changes in the Yarkovsky acceleration. Mass loss is largely adiabatic within about $10^{3}$ au (Veras et al., 2011), meaning that within this region it feeds back into the Yarkovsky acceleration only through an increase in asteroid semimajor axis.

Although each simulation contained 200 asteroids, we display the evolution of subsets of these on several of the figures for better clarity and reduced filesizes. Nevertheless, we have looked at the evolution of all 200 asteroids in each case to identify any noteworthy behaviour.

4.1 Model A

4.1.1 Semimajor axis evolution

For Model A, $\mathbb{Q}$ was chosen to maximize the asteroid’s outward drift due the Yarkovsky effect. This drift is independent of the diagonal terms (Eq. A2 of Veras et al. 2015), and is outward from equation (25). In Figs. 3-7, we illustrate the extent of this drift for asteroids ranging in radius over five orders of magnitude. In all these figures, the asteroids were initialised on circular orbits. The spread in the curves in each figure arise from the different initial values of $a_{\rm ast}$, $i_{\rm ast}$ and the orbital angles.

$R_{\rm ast}=1000$ km asteroids

Even in the extreme case of Model A, asteroids or planets with $R_{\rm ast}=1000$ km (Fig. 3) experience negligible Yarkovsky drift, and hence represent a useful standard for comparison. The evolution in Fig. 3 is what one would expect without any radiative effects: the inner disc, the planet (black curve), and the outer disc all increase their semimajor axes by a factor of about three based purely on stellar mass loss. The spread seen in the plot primarily mirrors the different semimajor axes chosen in the initial conditions.

$R_{\rm ast}=100$ km asteroids

For smaller asteroids, with $R_{\rm ast}=100$ km (Fig. 4), differences in the evolution become apparent. Both the inner and outer disc drift further outward from the $R_{\rm ast}=1000$ km case, indicating that the Yarkovsky effect is noticeable. Sharp changes in slope for inner disc asteroids close to the planet indicate that here planet-asteroid interactions play a more significant role than in Fig. 3. At kinks in the inner disc curves, resonant capture occurs, and is maintained. During this resonant capture, the asteroid’s eccentricity is increased to a greater extent than can be generated by the Yarkovsky effect. The resonant capture is maintained because these asteroids are too large to “break free” from the Yarkovsky effect.

$R_{\rm ast}=10$ km asteroids

The evolution of asteroids which are one order of magnitude smaller ($R_{\rm ast}=10$ km) exhibits planet crossing of the inner disc (Fig. 5). Even though the planet disrupts the asteroid orbits, the disruption is not great enough to qualitatively alter the final outcome. Nevertheless, the disruption in some cases highlights resonant crossings that were described in Section 2. In the lower panel of Fig. 5, we have isolated the evolution of five notable cases, which all experience temporary resonant trapping. This trapping occurs when the asteroid semimajor axis suddenly flattens out, or mirrors the planet’s evolution, seemingly unaffected by the Yarkovsky effect. These highlighted cases exhibit resonant trapping inward of the planet, outward of the planet, and coincident with the planet.

The coincident cases are of particular interest because they might showcase entrapment into the planet’s Hill sphere. Hence, we investigate further the evolution of the purple curve in Fig. 5 with Fig. 6, which highlights the timescale ($\approx 0.38-0.52$ Myr) over which the planet’s and asteroid’s semimajor axes track each other. The upper panel of Fig. 6 confirms that the asteroid is temporarily captured into a co-orbital $1$:$1$ mean-motion resonance. The resonant angle $\lambda_{\rm ast}-\lambda_{\rm pl}+\varpi_{\rm ast}-\varpi_{\rm pl}$ librates with a varying centre and amplitude, where $\lambda$ denotes mean longitude and $\varpi$ denotes longitude of pericentre. Often resonant angles librate around $0^{\circ}$ or $180^{\circ}$, but here the libration centre is here clearly under $180^{\circ}$, perhaps illustrating that a correction term due to the Yarkovsky effect would need to be included in the resonant angle. Computation of this correction term may not be trivial, even for the simplistic case of Model A, because the term could be a function of all orbital elements through Eqs. 59, A5 and A6 of Veras et al. (2015). During the capture, $e_{\rm ast}\approx 0.15$, but increases steadily, which might trigger the departure from the resonant state.

The bottom panel of Fig. 6 tracks the mutual distance between the planet and asteroid, which could help indicate if the asteroid becomes trapped within the planet’s Hill sphere. Because the Hill sphere for Neptune is under 1 au, the plot shows that the asteroid is not trapped within the Hill sphere, although does occasionally poke into it. The bottom panel’s $x$-axis is extended from the $x$-axis of the top panel to show both how the asteroid gradually drifts away from the planet over time and one particularly deep incursion into the Hill sphere at 0.54 Myr, for which there is a corresponding single libration on the top panel.

$R_{\rm ast}=1$ km and $100$ m asteroids

For even smaller asteroids in the inner disc (Fig. 7), the orbit crossing with the planet occurs too quickly for resonant capture to occur, and the semimajor axes that the asteroids attain along the white dwarf phase is orders of magnitude higher than the planet’s. Fig. 7 displays the results of two simulations on a single plot: those containing asteroids with $R_{\rm ast}=1$ km and $100$ m. In no case did the planet’s interaction with an asteroid trigger instability. Both sets of inner disc asteroids were propelled to a high enough semimajor axis for eccentricity to start playing a role in nonadiabatic evolution from mass loss (Veras et al., 2011). However, the eventual semimajor axis obtained along the white dwarf phase is primarily due to the Yarkovsky effect. In the $R_{\rm ast}=100$ m simulation, 81 of the 200 asteroids escaped the system. Here escape refers to leaving the time-dependent Hill ellipsoid of the star, as defined by Veras & Evans (2013) and Veras et al. (2014c), assuming that the star resides in the Solar Neighbourhood (at a distance of 8 kpc from the Galactic Centre).

Effect of planet

In order to better pinpoint the dependencies of the asteroid evolution on the presence of the planet and the initial eccentricities of the asteroids, we have created a figure (Fig. 8) with four plots that may be compared to one another on the same scale. The upper-left panel reproduces the $R_{\rm ast}=10$ km curves from Fig. 5 but without the planet curve. The upper-right plot shows the same simulation but without a planet. This comparison indicates that the disruption created by the planet does not qualitatively affect the final result. For the case when asteroids have initially eccentric orbits ranging from $e_{\rm ast}=0.3-0.7$ (bottom panels), then the planet has a marginal effect on the final semimajor axis distribution, despite strongly affecting a few asteroids (like the one associated with the bottom green curve).

4.1.2 Eccentricity evolution

Regarding how the eccentricity of the asteroids themselves change throughout the evolution for Yarkovsky Model A, consider Fig. 9, which includes all 200 simulated $R_{\rm ast}=10$km asteroids from the upper panels of Fig. 8. The right panel of Fig. 9, which does not include a planet, helps confirm equation (26) by illustrating that the radiation-driven Yarkovsky effect does not generate eccentricity in initially circular nearly coplanar asteroids. Nearly all of the eccentricity shown in this plot is created by stellar mass loss, at the expected level555The division between adiabatic and non-adiabatic evolution from mass loss is not sharp: non-zero eccentricity is always generated, even if relatively small, regardless of semimajor axis (Veras et al., 2011)..

The inclusion of a planet (left panel) dominates the eccentricity evolution of asteroids by orders of magnitude; the final eccentricities of most of the asteroids are greater than 0.1. Consequently, a planet exerts greater influence on the eccentricity evolution of asteroid discs than their semimajor axis evolution. This eccentricity is pumped almost immediately – within about $10^{5}$ yr – and is independent of post-main-sequence evolution. Note that the eccentricity of the planet itself (black curve) is also increased due to stellar mass loss up to about $10^{-2}$.

This asteroid eccentricity pumping may be key to understanding why Hill sphere capture does not occur in our simulations. As indicated by Fig. 1, an increase of just 0.06 may prevent capture from occurring. Without a protoplanetary disc to keep the asteroid eccentricity low, the gravitational interaction between the planet and asteroid alone may be strong enough to prevent capture from occurring. Nevertheless, not all of the asteroid eccentricities are pumped to values higher than 0.06, and the bottom panel of Fig. 5 indicates that temporary resonant capture may still occasionally occur. Discs are not actually needed to induce resonance capture (Raymond et al., 2008), but increase the chances by over an order of magnitude.

Further, the critical eccentricity condition from equation (5) holds only for small eccentricities, and along with equations (3) and (4), does not take into account the perturbations from the planet. Further, the asteroid’s speed is always changing. Capture within the Hill sphere is also more likely if there is a quick dissipation mechanism. In principle, this dissipation may occur from Yarkovsky dissipation as the star transitions into a white dwarf, but the timing of this event combined with the size of the asteroid would need to be fine-tuned.

4.1.3 Inclination evolution

In all of our simulations, our asteroids are given initial inclinations with respect to the planet’s orbital plane between $0^{\circ}$ and $10^{\circ}$. Fig. 10 illustrates how the Model A Yarkovsky effect affects the inclination evolution. The 10 km asteroids in the figure are the same as those in Fig. 9 and in the upper panel of Fig. 8.

The right panel of Fig. 10 displays a net decrease in all orbital inclinations. This behaviour is expected from the analytics (equation 27). Note additionally that mass loss has no effect on inclination evolution unless the mass loss is anisotropic (Veras et al., 2011, 2013a; Dosopoulou & Kalogera, 2016a, b; Veras et al., 2018). Therefore, the levelling-out of the inclination evolution curves entirely reflects the reduced influence of the Model A Yarkovsky effect as the asteroid semimajor axes increase in a nonlinear fashion.

The presence of a planet (left panel) spreads the curves about the $0^{\circ}$ line, although there still remains a net decrease in the inclination evolution. Like for orbital eccentricity, the initial perturbations, as seen up to about 0.6 Myr after the start of the asymptotic giant branch phase, would have realistically occurred much earlier in the system evolution and is just an artefact of our chosen starting time.

4.2 Model B

In Model B, $\mathbb{Q}$ has been chosen to maximise the asteroid’s inward drift due the Yarkovsky effect. Because of the dependence of the Yarkovsky drift on $a_{\rm ast}$, there is a positive feedback effect, accelerating the asteroid more quickly towards the central star as the semimajor axis decreases. However, stellar mass loss represents a competing effect, pushing the orbit outward.

Which effect wins depends on the size of the asteroid. Outward expansion due to mass loss is independent of the asteroid mass at the $M_{\rm ast}/M_{\star}$ level, whereas the Yarkovsky effect is proportional to the square of the asteroid radius (equation 23). This latter strong dependence dictates that asteroids with radii as large as $10$ km cannot avoid engulfment, regardless of the presence of the planet.

Figure 11 illustrates their evolution, for asteroids on initially circular orbits, and showcases the potentially destructive power of the radiation-induced Yarkovsky effect. In each case, the evolution of all 200 asteroids are shown, and none survive to reach the tip of the Asymptotic Giant Branch phase. The left panel includes the influence of the planet, and the right panel does not. Although the planet disperses incoming asteroid orbits with a flair apparent in the left panel, the presence of the planet does not actually qualitatively change the fate of these asteroids. Near-horizontal segments in the left panel indicate temporary capture within mean motion resonance, which can last long enough to delay engulfment into the star for a few $10^{5}$ yr. We do not model engulfment nor tidal effects, as none of these asteroid evolutions are necessarily realistic anyway: recall they represent just bounding cases for Yarkovsky effects along the giant branch evolution phases.

For $R_{\rm ast}=100$ km asteroids, the competing effects of Model B Yarkovsky drift and stellar mass loss are comparable (Fig. 12). Shown is the semimajor axis change from its initial value, to highlight decreases versus increases. In fact, the direction of semimajor axis drift of one of these asteroids changes depending on if it resides in the inner or outer disc. Figure 12 illustrates the evolution of both discs, with (left panel) and without (right panel) a planet. The inner disc of 95 asteroids is dragged towards the star, where 71 become engulfed. All 105 asteroids in the outer disc survive, despite initially being slightly dragged towards the star before increasing their semimajor axes. The presence of the planet slightly helps prop up the semimajor axes of the outer disc, but again does not qualitatively affect the outcome. In fact, the planet induces just one more asteroid to be engulfed than if the planet was not present.

4.3 Model C

In Model C, we quench the planar component of all Yarkovsky-based semimajor axis evolution by setting $\mathbb{Q}_{12}=\mathbb{Q}_{21}$. The model also zeroes out the other non-diagonal terms of $\mathbb{Q}$ ($\mathbb{Q}_{13}=\mathbb{Q}_{23}=\mathbb{Q}_{31}=\mathbb{Q}_{32}$) which reverts the semimajor axis drift to Poynting-Robertson drag.

Confirmation of this quenching appears in the left panel of Fig. 13, which shows the discs evolving along the giant branch due primarily to mass loss. We have also tested semimajor axis evolution for both $R_{\rm ast}=1$ km and $R_{\rm ast}=100$m asteroids in low-resolution simulations with just a few asteroids each (for computational feasibility), and the result is the same.

The right panel of the figure, however, reveals that the inclination evolution is not zeroed-out, and in fact can change significantly more than in Model A. The reason is because for circular, near-coplanar orbits in Model C,

[TABLE]

which allows the inclination to change initially more quickly than in equation (27). Another reason is that $\left\langle\frac{d\Omega_{\rm ast}}{dt}\right\rangle$ is exactly twice as fast as in Model A in this limit. When both of these differential equations are solved simultaneously, time can be eliminated, yielding

[TABLE]

[TABLE]

where $C$ is an arbitrary constant. Equations (34-35) explicitly illustrate the steeper dependence on longitude of ascending node in Model C.

These initial quick changes can have a positive feedback effect on the asteroid inclination before $a_{\rm ast}$ becomes too large from stellar mass loss. Note that the right panel of Fig. 13 shows that only some asteroids achieve inclinations of several tens of degrees, because of the spread of their initial values. Further, all of the asteroids which ended up on the most highly inclined orbits had a “quick start” (fast initial inclination change). In this sense, the duration of the giant branch phase is key to this inclination pumping. Here we sampled just one case, with an initially $2M_{\odot}$ star. More massive stars would have shorter giant branch lifetimes but higher luminosities.

Model C Yarkovsky drift illustrates how an asteroid’s orbital inclination can flip from prograde to retrograde and back again in a smooth fashion. One potential consequence of this type of motion is the formation of a near-isotropic cloud of debris, perhaps akin to so-called mini Oort-clouds (Raymond & Armitage, 2013). The geometry of these clouds could have potentially important implications for white dwarf pollution (Alcock et al., 1986; Parriott & Alcock, 1998; Veras et al., 2014d; Stone et al., 2015; Caiazzo & Heyl, 2017).

5 Discussion

None of the individual asteroid evolutions that we simulated are likely to be realistic. The reality is complex, as indicated both by observations from within the solar system and by detailed efforts to model those observations. Rather, our simulations reveal physical interplays and provide limits on the characteristics of motion that we might expect for given parameters. These approximations do, however, give some more depth to the constant semimajor axis drift approximations sometimes employed for asteroids orbiting the Sun; this generalisation is necessary in giant branch planetary systems because of the host star’s “supercharged” luminosity.

However, our models do not go nearly far enough. The entries of $\mathbb{Q}$ are functions of time, and in turn, of the spin and thermal properties of the asteroid (Veras et al., 2015). The asteroid’s spin and thermal properties may then be a strong function of its shape (Rozitis & Green, 2012, 2013; Vokrouhlický et al., 2015; Golubov et al., 2016), which, as indicated by the majority of solar system asteroids, is unlikely to be spherical.

Another consideration is the connection between the Yarkovsky and YORP effects. YORP spin-down can lead to tumbling states (the asteroid Toutatis is an example of a tumbler), and these states can “switch off” the Yarkovsky effect (Vokrouhlický et al., 2007). However, tumbling states could represent transient phenomena in systems where host stars change their luminosity on short timescales. Without dedicated investigations, this interplay remains unclear. Further, giant branch rotational fission from YORP could be ubiquitous within a particular distance (estimated to be about 7 au by Veras et al. 2014a), and the resulting fragments will be closer to the giant branch star than any considered here, and hence subject to even stronger Yarkovsky forces. Further, some asteroids at Kuiper-belt like distances, or further away, may already be spinning quickly enough to experience YORP break-up during giant branch evolution.

Our model ignores the effect of collisions, which can play a vital role in shaping the debris fields of giant branch planetary systems (Bonsor & Wyatt, 2010). These collisions can both modify the size distribution of the asteroids towards smaller sizes, and affect YORP evolution (Marzari et al., 2011) through the resetting of the spin state (with a timescale that goes as the square root of the asteroid size, see e.g. Farinella et al. 1998). A driver for collisional activity is eccentricity excitation, which we have shown is generated by the presence of a planet (Fig. 9). The consequence of significant collisional comminution will be an increase in the type of behaviour seen in Fig. 7, where smaller particles are flung into the outer reaches of the planetary system. The final orbital distribution of debris may then tend towards high values of semimajor axis when collisions are frequent and early. In this respect, by not modelling collisions, our simulations provide a type of lower bound for the extent of the eventual outward migration of debris.

A helpful result of this paper is that the giant branch Yarkovsky effect can be strong enough to render the presence of a planet as relatively unimportant. Planets certainly play a dynamical role in pumping eccentricity (Fig. 9) and changing timescales for destruction (Fig. 11). However, they struggle to retain asteroids in mean motion resonances and, in the big picture, do not qualitatively change how many asteroids are destroyed around giant branch stars or where asteroids will reside after reaching the white dwarf phase.

These statements, however, must be caveated with the fact that dissipation at different points during the evolution could result in temporary capture into mean motion resonances (Sections 2.1 and 4.1.4) or within the Hill sphere of the planet (Higuchi & Ida 2016, 2017; also see Sections 2.2-2.3).

The relevant formulae (equations 3, 4, 5, 14 and 15) are affected by post-main-sequence evolution through just $a_{\rm pl}$, $a_{\rm ast}$ and $M_{\star}$. The changes in these critical values that are induced by variations in $M_{\star}$ are no greater than a factor of a few. This factor of a few change also holds for $a_{\rm pl}$ and $a_{\rm ast}$ unless the latter value becomes so large that the planet would be too far away to produce capture anyway. Planets could play a much larger role with a time-varying $\mathbb{Q}$, which may provide dissipation at key times.

Another useful result of this paper are constraints on asteroid sizes. Even in the most extreme cases (exemplified by Models A, B, and C), the largest exo-Kuiper belt asteroids – on the order of $R_{\rm ast}=1000$ km – are not meaningfully affected by the Yarkovsky effect. On the other end of the size spectrum, almost every asteroid smaller than about 10 km is predominately perturbed by the radiation-induced Yarkovsky effect. This effect should be considered when establishing initial conditions for simulations along the white dwarf phase.

As is apparent in most of the plots in this paper, dynamical settling occurs quickly just at the onset of the white dwarf phase. Although the Yarkovsky effect does not “turn off” due to the white dwarf, the luminosity of the white dwarf becomes sub-solar within just a few Myr. Combined with the greatly expanded asteroid orbits from giant branch evolution, the white dwarf radiation-induced Yarkovsky effect is then often negligible.

In this paper, we sampled only a brief time interval within the entire evolution of a planetary system (due to computational limitations). Hence, our initial conditions were not realistic, but were chosen instead to demonstrate physical trends. Reaching the asymptotic giant branch phase would have required asteroids to first undergo all of main sequence and giant branch evolution, and maintaining $e_{\rm ast}=0$ throughout those periods may be difficult. During the white dwarf phase, instabilities amongst planets coupled with planet-asteroid interactions can lead to white dwarf pollution (Bonsor et al., 2011; Debes et al., 2012; Frewen & Hansen, 2014; Veras et al., 2016; Mustill et al., 2018; Smallwood et al., 2018). Therefore, accurately computing the relative positions of the asteroids and planets after the asymptotic giant branch phase is crucial for the estimation of white dwarf pollution rates and timescales.

Will our own Sun be polluted? The answer largely depends on the fate of the Kuiper Belt, along with whether Planet Nine exists (Veras, 2016b). One complication for the solar system that is not addressed in this paper is the Yarkovsky effect generated by a red giant branch star (as opposed to an asymptotic giant branch star). $1M_{\odot}$ main sequence stars will undergo two largely comparable periods of enhanced luminosity, at the tips of both the red giant and asymptotic giant branch phases. The flavour of Yarkovsky drift (Model A versus Model B versus Model C versus something in-between) could change both within and between these two phases.

6 Summary

Understanding how asteroids evolve during giant branch evolution crucially determines their capability to pollute the eventual white dwarf. Although the often dominant effects of Yarkvosky drift generated from post-main-sequence stellar radiation have been analytically estimated previously (Veras et al., 2015), here we provide more detail by (i) running numerical simulations with three different extreme models, and (ii) introducing a planet. The three models were chosen to place limits on the types of motion and orbital changes that we may expect, and could aid the future development of more sophisticated shape, spin and thermal inertia constructions. The range of outcomes is much greater than what is observed in the solar system, with semimajor axis changes that can vary by orders of magnitude and easily achieved orbital inclination flipping. Amidst this enhanced radiative forcing, the influence of the planet is minimised. We also analytically considered how these planets could capture asteroids within their Hill spheres (Higuchi & Ida, 2016, 2017) and into mean motion resonances (as in protoplanetary disc migration), but numerically found that these processes are infrequent and ineffectual without a fine-tuned dissipation prescription.

Acknowledgements

We are grateful to Apostolos Christou for his helpful referee report, which has improved the manuscript. DV thanks the Earth-Life Science Institute at the Tokyo Institute of Technology for their hospitality during his stay and for the initiation of this project. DV also gratefully acknowledges the support of the STFC via an Ernest Rutherford Fellowship (grant ST/P003850/1).