Chaotic rotation and evolution of asteroids and small planets in high-eccentricity orbits around white dwarfs

TL;DR

This paper investigates how chaotic rotation affects the orbital evolution and potential breakup of small bodies on highly eccentric orbits around white dwarfs, revealing a new rotational fission mechanism that influences debris distribution.

Contribution

It introduces a novel approach to model the coupled physical and orbital evolution of high-eccentricity asteroids, highlighting a self-regulating chaotic rotation process and a new rotational fission mechanism.

Findings

Chaotic rotation can lead to self-regulation, reducing orbital decay effects.

Asteroids can break apart without entering the Roche radius.

A steady stream of debris impacts white dwarfs, affecting pollution levels.

Abstract

Observed planetary debris in white dwarf atmospheres predominately originate from the destruction of small bodies on highly eccentric ($>0.99$) orbits. Despite their importance, these minor planets have coupled physical and orbital evolution which has remained largely unexplored. Here, we present a novel approach for estimating the influence of fast chaotic rotation on the orbital evolution of high-eccentricity triaxial asteroids, and formally characterize the propagation of their angular rotation velocities and orbital elements as random time processes. By employing the impulse approximation, we demonstrate that the violent gravitational interactions during periastron passages transfer energy between the orbit and asteroid's rotation. If the distribution of spin impulses were symmetric around zero, then the net result would be a secular decrease of the semimajor axis and a further increase of the eccentricity. We find evidence, however, that the chaotic rotation may be self-regulated in such a manner that these effects are reduced or nullified. We discover that asteroids on highly eccentric orbits can break themselves apart --- in a type of YORP-less rotational fission --- without actually entering the Roche radius, with potentially significant consequences for the distribution of debris and energy requirements for gravitational scattering in metal-polluted white dwarf planetary systems. This mechanism provides a steady stream of material impacting a white dwarf without rapidly depleting the number of small bodies in the stellar system.

Tables

Table 1. Table 1 Model parameters used in simulations. The semimajor axis corresponds to 1.5 au, the white dwarf mass to 0.6 ​ M ⊙ 0.6 subscript 𝑀 direct-product 0.6M_{\odot} and the asteroid mass to a value which is about 5% the mass of Ceres.

Name	Description	Units	Value
$a$	semimajor axis	m	$2.244\cdot {10}^{11}$
$e$	orbital eccentricity		variable ($>0.95$)
$M_1$	mass of primary body(white dwarf)	kg	$1.193\cdot {10}^{30}$
$M_2$	mass of asteroid	kg	$4.4\cdot {10}^{19}$
$\sigma =\left(B-A\right)/C$	triaxiality of asteroid		$0.05$
$R$	radius of asteroid	m	$2.1\cdot {10}^5$
$\xi$	coefficient of inertia		$0.35$

Equations

$${\theta }^{\prime \prime }(t)+\frac{3}{2}{n}^2\frac{\sigma \sin (2\theta (t)-2\nu (t))}{{\left(1-e\cos E(t)\right)}^3}=0,$$

$${\dot{\nu }}_{\max }/n=\frac{\sqrt{1+e}}{(1-e{)}^{\frac{3}{2}}}.$$

$$z_i=c+\varphi z_{i-1}+{ϵ}_i+{\psi }_1{ϵ}_{i-1}+{\psi }_2{ϵ}_{i-2},$$

$$\mathcal{E}[x_i^2|\{x_{i-1}\}]={\varpi }_i^2=\kappa +{\alpha }_1{x}_{i-1}^2+{\beta }_1{\varpi }_{i-1}^2,$$

$$E_{\mathrm{orb}}+E_{\mathrm{rot}}=-G\frac{m_1{m}_2}{2a}+\frac{\xi m_2}{2}{R}^2{\omega }^2=const.$$

$$\frac{d_a}{a}\simeq -\frac{\xi (2\omega d_{\omega }+d_{\omega }^2)}{n^2}{\left(\frac{R}{a}\right)}^2.$$

$$\begin{array}{r}{L}_{\mathrm{tot}}=L_{\mathrm{orb}}+L_{\mathrm{rot}}=G^{\frac{1}{2}}\frac{m_1{m}_2}{(m_1+m_2{)}^{\frac{1}{2}}}\sqrt{a(1-e^2)}\\ +\xi m_2{R}^2\omega =const.\end{array}$$

$$\frac{d_e}{e}\simeq \frac{\xi \sqrt{1-e^2}}{n{e}^2}{\left(\frac{R}{a}\right)}^2{d}_{\omega }.$$

Full text

I. Introduction

Understanding the complete history and future of a planetary system usually requires investigating the host star as it changes phases. Nearly every star in the Milky Way will or already has traversed an evolutionary sequence involving main sequence, giant branch and white dwarf phases. Snapshots of white dwarf planetary systems in particular reveal exclusive and unique insights about planetary composition (Jura & Young, 2014; Farihi, 2016; Bonsor & Xu, 2017; Harrison et al., 2018; Hollands et al., 2018; Vanderburg & Rappaport, 2018; Zuckerman & Young, 2018) and include striking examples of disintegrating (Vanderburg et al., 2015) and intact (Manser et al., 2019) orbiting minor planets.

Both major and minor planets which survive until the white dwarf phase have endured physical and orbital variations resulting from stellar evolution (Veras, 2016). Giant branch stars induce physical variations primarily through stellar mass loss, envelope expansion and increased luminosity. Gas giant planets may accrete stellar mass through its wind, altering the composition of the planetary atmospheres (Spiegel & Madhusudhan, 2012). Envelope expansion may alter planetary surfaces and interior energy budgets through tidal interactions. Enhanced stellar luminosity may shear off atmospheres (Livio & Soker, 1984; Nelemans & Tauris, 1998; Soker, 1998; Villaver & Livio, 2007; Wickramasinghe et al., 2010), melt or sublimate entire planets or some of their components, and spin up, spin down and break up minor planets through YORP-induced rotational fission (Veras et al., 2014a).

Better understood are the planetary orbital changes which are triggered by giant branch stars. Isotropic stellar mass loss expands orbits at prescribed rates (Omarov, 1962; Hadjidemetriou, 1963; Veras et al., 2011) and anisotropic mass loss (Veras et al., 2013; Dosopoulou & Kalogera, 2016a, b) is not expected to become significant for bodies within Oort-cloud distances. Orbital variations due to stellar mass loss are effectively independent of planet mass. Tidal effects between giant stars and planets crucially set the critical displacement beyond which a planet will survive engulfment (Villaver & Livio, 2009; Kunitomo et al., 2011; Mustill & Villaver, 2012; Adams & Bloch, 2013; Nordhaus & Spiegel, 2013; Valsecchi & Rasio, 2014; Villaver et al., 2014; Madappatt et al., 2016; Staff et al., 2016; Gallet et al., 2017; Rao et al., 2018; Stephan et al., 2018). Effects from giant star stellar luminosity may in fact dominate over those from gravity when altering the orbits of minor planets through the Yarkovsky effect (Veras et al., 2015a, 2019a).

The combined result of the above forces is that as a star becomes a white dwarf, it will initially clear out the inner few au of all objects, and allow only planets larger than about 10-100 km in radius to survive intact within about 10 au. Despite these large distances, observations indicate that these objects must reach and accrete onto the white dwarf (which has a typical radius of $1R_{\oplus}$ and a typical Roche radius of $1R_{\odot}$), and do so regularly: Between one-quarter and one-half of all white dwarfs contain planetary metals (Zuckerman et al., 2003, 2010; Koester et al., 2014); chemical abundances of planetary debris have been measured for white dwarfs up to 8 Gyr old (Hollands et al., 2018), and the white dwarf discs generated from minor planet break-up are recycled on short timescales of $~{}10^{4}-10^{6}$ yr (Girven et al., 2012).

Consequently, an important issue is how to perturb minor planets (asteroids, moons and comets) on eccentric-enough orbits ($e>0.99$) to reach the white dwarf. Many methods have been successfully invoked: gravitational perturbations by a single planet (Bonsor et al., 2011; Debes et al., 2012; Frewen & Hansen, 2014), multiple planets (Veras et al., 2016; Payne et al., 2016, 2017; Mustill et al., 2018; Smallwood et al., 2018), multiple stars (Bonsor & Veras, 2015; Hamers & Portegies Zwart, 2016; Petrovich & Muñoz, 2017; Stephan et al., 2017, 2018) and Galactic tides (Alcock et al., 1986; Parriott & Alcock, 1998; Veras et al., 2014b; Stone et al., 2015; Caiazzo & Heyl, 2017). The outcome of a minor planet entering the Roche radius of a white dwarf has also been explored (Debes et al., 2012; Veras et al., 2014c; Wyatt et al., 2014; Veras et al., 2015b; Brown et al., 2017), particularly with regard to the resultant disc evolution (Bochkarev & Rafikov, 2011; Rafikov, 2011a, b; Metzger et al., 2012; Rafikov & Garmilla, 2012; Kenyon & Bromley, 2017a, b; Miranda & Rafikov, 2018). Star-planet tides negligibly affect the orbital motion of minor planets (Veras et al., 2019b), such as the ones uncovered by Vanderburg et al. (2015) and Manser et al. (2019).

However, in nearly every case listed above, the physical evolution of asteroids on highly eccentric orbits was not modelled, nor was the resulting feedback on the orbital angular momentum. This paper strives to address this missing component of our understanding. In Section 2 we describe our procedure for simulating the coupled evolution of these asteroids. In Section 3 we argue that a highly-eccentric asteroid need not encounter the Roche radius before spinning itself apart. Section 4 contains a formal analysis of our simulation output in terms of random processes which describe the evolution of the physical and orbital elements. We derive functional dependencies of the orbital evolution on these random processes in Section 5 and discuss and summarize our results in Section 6. We henceforth refer to our orbiting objects as “asteroids”.

II. Chaotic rotation

Rotation of elongated bodies around their principal axes of inertia in two-body systems are mostly driven by the gradient of the force of gravitational attraction exerted by the orbiting companion. The tidal interactions are much smaller in magnitude for non-vanishing parameters of triaxiality $\sigma=(B-A)/C$, which are often found in smaller planets and asteroids. The $A$, $B$, and $C$ here and in the following designate the principal moments of inertia in increasing order. In this paper, we use the term “asteroid” in its generic meaning, which includes celestial bodies rigid enough to maintain a permanent shape, such as comets, asteroids, and minor planets. Our results are applicable to a wide range of such objects, as long as they have a significantly prolate shape and highly eccentric orbits. Rotation of triaxial asteroids is known to have complex, structured sections in the parameter space (Wisdom et al., 1984). The islands of stable equilibrium in the Poincaré sections are surrounded by bands of purely chaotic motion, which are present even for small eccentricity orbits as long as $\sigma$ is finite. In the narrow boundary zone between chaos and libration zones, small secondary resonance islands may be present (Flynn & Saha, 2005), but their physical validity remains to be confirmed. Only Hyperion and, possibly, Prometheus and Pandora, represent satellites in chaotic rotation in the Solar system (Kouprianov & Shevchenko, 2005; Melnikov & Shevchenko, 2008). But chaotic satellites are likely common around minor planets in the outer Solar System. Many TNOs have eccentric satellites, and a particularly instructive example is the satellite Thorondor around the large TNO (385446) Manwë (Grundy et al., 2014), which has an eccentricity of 0.56. The prevalence of synchronous rotation can probably be explained by the regularization effect of tidal dissipation in dynamic systems of small and moderate eccentricity. The synchronously rotating satellites of significantly elongated shape, such as Phobos and Epimetheus, could not enter their current spin-orbit state without having to cross the chaotic zone (Wisdom, 1987).

The case of high eccentricity ($e>0.95$) and significant triaxiality is different, because chaotic behavior is the only possible state for such asteroids. The islands of stable equilibrium vanish, and the tidal dissipation can not regularize the motion even in the long run. We investigate this case by means of numerical simulations for the simplest setup, in which the orbit and the asteroid’s equator are coplanar (i.e., the obliquity of the orbit on the equator is zero) and the problem becomes one-dimensional. The tumbling rotation of prolate asteroids is certainly three-dimensional in reality, and this approximation is likely to slightly overestimate the associated acceleration. Tidal forces are ignored. Our model object of choice is similar to Proteus in size, mass, and shape, but we put it into a long-period, high-eccentricity orbit around a white dwarf primary (see Table 1 for our adopted physical and orbital elements). The semimajor axis of 1.5 au was chosen for computational feasibility (see Section 4), and we expect the results to be qualitatively similar for higher values of $a$. The orbital parameters are constant, an assumption which is justified for a relatively short integration spanning less than $10\,000$ orbits in view of the relatively slow rate of orbital evolution, as described below.

The integrated ordinary differential equation (ODE) is

[TABLE]

where $t$ is time, $n$ is the mean orbital motion, $\nu(t)$ is the true anomaly in radians, $E(t)$ is the eccentric anomaly in radians, and $\theta$ denotes the orientation angle of the asteroid counted from its longest axis in the inertial space (Danby, 1962). The functions of the time-variable anomalies $\nu$ and $E$ are often replaced with well-known series in powers of eccentricity, relating them to the linearly changing mean anomaly. The main challenge of this problem is that such series become impractical at large $e$ because an exponentially growing number of terms has to be used, and the associated Kaula’s functions of eccentricity become numerically large. We have to revert to the traditional computation of $E$ using Kepler’s equation. We made use of the technique of inverse interpolation proposed by Tommasini & Olivieri (2018) to compute $E$ as a function of mean anomaly. We also employed an integration scheme with self-adaptive time steps to capture the all-important sharp variations of $\theta$ that take place during the periastron passage. The time step has to be much shorter than the typical duration of this passage, while a coarser time step is prudent to use outside this orbit phase interval, i.e., for 99.9% of the orbit. The ODE (1) can be integrated with two additional boundary conditions for $\theta(0)$ and $\theta^{\prime}(0)$.

Fig. 1 shows a typical result of numerical simulation of 1700 orbits for $e=0.99$. It shows the normalized apoastron rotation velocity $\dot{\theta}/n\equiv\omega/n$ versus orbit number, for initial conditions $\theta(0)=-0.04$ rad, $\theta^{\prime}(0)=801\;n$. The behavior is chaotic with a short Lyapunov exponential; hence, even a marginally small perturbation of the initial values brings about a completely different rotation curve. However, important observations can be made based on this limited example. The spin rate varies in a very broad range, but it seems to be confined to a certain interval of prograde rotation. It bounces back when it comes close to zero. The asteroid’s rotation is fast most of the time, but can it become arbitrarily fast? The first idea is that the velocity may be limited to the maximum orbital motion at the moment of closest approach, which is computed as

[TABLE]

A few characteristic values are 124.9 at $e=0.95$, 1410.7 at $e=0.99$, 15803.5 at $e=0.998$. The results in Fig. 1 for $e=0.99$ indicate that the spin rate can become higher than $\dot{\nu}_{\rm max}$ even on a relatively short time scale. Indeed, integrations with other initial conditions revealed that the rate of rotation may become as high as $3000\,n$. We give more consideration to the issue of the highest rate of rotation in § IV.

III. Rotational fission versus tidal break-up

Much attention has been paid in the literature to the issue of rotational fissure of Solar system asteroids (Walsh et al., 2008; Ortiz et al., 2012; Vokrouhlický et al., 2015), which may account for a significant population of binary asteroids of small mass ratio. Observations show that most minor planets with diameters smaller than 10 km have a lower bound of $\sim 2$ hr for their rotation periods (Polishook et al., 2017). Within the framework of a “rubble pile” model, where the bodies are held together predominantly by self-gravitation, the critical spin rate of fission depends on the diameter and other, less observable parameters. Larger asteroids with diameters greater than 10 km rotate slower than $\sim 4$ hr (Warner et al., 2009). It is therefore of interest to know if impulse-like excitation of rotation at perihelion can drive asteroids to critical spin rates, and how large the eccentricity should be ($\equiv e_{\rm critical}$) for that to happen.

Fig. 2 shows the values of critical eccentricity computed with these simplifying assumptions: 1) tidal fission takes place when the periastron distance $a\,(1-e)$ becomes smaller than $1\,R_{\sun}$, which is approximately the Roche radius for white dwarfs; 2) the highest rate of rotation is close to the maximum periastron angular velocity, Eq. 2. Rotational fission of larger asteroids, such as our Proteus analog, is achieved at lower eccentricity than the tidal break-up. The difference in eccentricity may seem small, but the difference in the periastron distance is significant. This result suggests that smaller parts and debris can agglomerate around white dwarfs well outside the Roche radius, to be delivered to the stellar surface by other means.

We already noted in § II that the actual spin rate can become much higher than $\dot{\nu}_{\rm max}$, as follows from our numerical simulations. A single impulse that drove the asteroid to spin-up above the critical value is enough to break it. Therefore, fission by centrifugal force is likely to be more efficient, especially for larger asteroids, because it can happen at a lower eccentricity. More discussion of the upper boundary of rotational velocity is given in § IV.

IV. Rotation velocity as a random process

Figure 1 shows a rotation velocity curve obtained by numerical integration. The large changes of velocity without any pattern or periodicity suggest that this may be a random process. This figure shows only the values of velocity computed at apoastrons where the triaxial torque is at a minimum (Eq. 1). Fig. 3, left panel, depicts how the simulated velocity changes within each orbit. For most of the time, the variation is relatively slow, but during a periastron encounter, the acceleration increases by a factor $\sim 1$ million (at $e=0.99$), leading to an almost step-wise change in velocity (Fig. 3, right panel). In view of these findings, the impulse approximation is the most suitable model for long-period asteroids with $e>0.95$.

We assume in the impulse approximation that all the dynamical interaction takes place momentarily during the closest approach to the star at the periastron. The result of this interaction is an abrupt change in rotation velocity $\omega\equiv\dot{\theta}$. It is sufficient then to measure the velocity at the apoastron where $\nu=E=\pi$, and to consider the difference between two consecutive apoastron velocities $d_{\omega,i}=\omega_{i+1}-\omega_{i}$ as the impulse magnitude. The discrete sets of tuples $\{t_{i},\omega_{i}\}$ and $\{t_{i},d_{\omega,i}\}$ are considered to be random processes in time.

We performed several numerical integrations of ODE 1 with varying initial conditions for fast prograde rotation velocities, $e=0.99$, and 9000 orbits. The solutions were sampled at apoastron times and analyzed as random processes using the Wolfram Mathematica TimeSeriesModelFit function111https://reference.wolfram.com/language/ref/TimeSeriesModelFit.html. This function determines which kind of random process best describes the given sequence of data and fits the corresponding model parameters. For the normalized rotational velocity $\omega/n$, the most frequent output is an Autoregressive Moving Average process ARMA$(1,2)$. A time series of values $z_{t}$ can be represented in this model as

[TABLE]

where $c$ is a constant, $\epsilon_{i}$ is a sequence of independent random numbers with a zero mean and a variance $\varpi^{2}$, and $c$, $\phi$, $\psi_{1}$, and $\psi_{2}$ are the fitting model parameters. The presence of a nonzero $\phi$ demonstrates the auto-regressive part of the model, in that each outcome depends on the previous value. Generally, the process is stationary if $|\phi|<1$. The one or two nonzero parameters $\xi$ indicate that the random impulses are correlated with one or two preceding impulses, which represents the moving average property. It may seem puzzling why the impulse magnitude, which depends only on the apoastron velocity and orientation angle, is correlated with the preceding impulses, but the reason will be revealed in the following analysis of $d_{\omega}$ properties.

One particular integration and model fitting for 9000 orbits produced these results: $c=4.5$, $\varpi^{2}=3010$, $\phi=0.998$, $\psi_{1}=0.009$, $\psi_{2}=-0.034$. Although the values vary between simulations, some common features are present. The $\phi$ value is close to unity, which means the process is weakly stationary. Starting the process with “unreasonable” initial conditions results in an extended transitional phase before the behavior settles down and stabilizes, which explains why long integrations covering many orbits are required to produce a reliable output. Because the root of the autoregressive polynomial is greater than 1, the process is causal. The values of $\psi_{1}$ and $\psi_{2}$ are small, thus, the impulses are weakly correlated. One of the roots of the moving average polynomial is negative, so the latter is not invertible. Finally, the positive value of $c$ indicates that the rotation velocity has a tendency to infinitely grow. However, this result is doubtful in view of complications, which we will now discuss.

Generally, we do not find the fitted ARMA(1,2) process to adequately represent the simulated velocity curve. Some important features seem to be missing in this model. For example, nothing prevents an ARMA(1,2) process with the estimated parameters to descend to negative values, whereas our simulations are always confined to a range of prograde rotation. Although a stationary process is intuitively expected, the moving average part appears to be a distorted reflection of more intricate properties of the system, which are not captured by this simple model. We believe the main reason for these shortcomings is the peculiar distribution of velocity impulses, which is neither Gaussian, nor identical for the process instances.

More progress can be made if we consider the time series of velocity impulses $d_{\omega}$ as a random process. The model fitting procedure invariably produced a Generalized Autoregressive Conditionally Heteroscedastic model GARCH(1,1). A discrete-time GARCH(1,1) process $x(t)$ is a series of statistically independent random numbers with a zero conditional mean and a conditional variance

[TABLE]

where $\kappa$, $\alpha_{1}$, and $\beta_{1}$ are positive-definite model parameters. This fit captures one additional curious property of the process under investigation, viz., the variance of a velocity pulse depends on both the variance and the realization of the previous pulse. The unconditional variance of $x(t)$ is $\kappa/(1-\alpha_{1}-\beta_{1})$, so that for a weakly stationary process $\alpha_{1}$ and $\beta_{1}$ should be between 0 and 1.

A particular numerical integration produced these model fit parameters for $d_{\omega}/n$: $\kappa=28.2$, $\alpha_{1}=0.119$, $\beta_{1}=0.872$, which represent a typical outcome of this simulation. We note the strong dependence of the current variance on the preceding variance. There is no requirement that the underlying distribution is Gaussian. If the process wanders into a domain of small impulse values, the dispersion of impulse values is likely to remain small for a period of time. In other words, the velocity is likely to be less variable in the future if the previous variations have been small. This is similar to the behavior of the financial market price volatility, for example, where extended periods of lull are separated by periods of elevated variability.

The origin of these peculiarities becomes more understandable when we map the apoastron velocity pulses versus velocity values, analogous to the Poincaré section mapping. Fig. 4, left, shows a realization obtained by numerical integration over 9000 orbits with $e=0.99$, $\theta(0)=-0.02$, $\dot{\theta}(0)=800\,n$, and other parameters from Table 1. The map confirms that the velocity is limited by $\sim 100\,n$ at the minimum, but there seems to be no upper bound. More surprisingly, the domain of possible rotation rates and impulses is confined in the latter on both sides. For any specific apoastron velocity, the distribution of periastron impulses is finite. The chaotic nature of the process is betrayed by the random position of states within the domain. When the velocity stochastically reaches the smallest possible value, the process is “cornered” in such a way that only vanishingly small velocity updates are possible, and the process lingers in this low-volatility regime for an extended period of time. The distribution of velocity impulses is not symmetric around zero. The largest positive (prograde) updates happen at a lower velocity than the largest negative (retrograde) updates. Therefore, when the velocity is smaller than $\sim 800\,n$, greater prograde impulses are possible and the velocity is likely to wander towards faster values. The situation reverses when it becomes greater than $\sim 1000\,n$, where retrograde impulses of greater magnitude become possible, which will likely succeed in pushing the velocity back. This peculiar distribution allows the process to remain stationary and in the long run, vary around some median velocity range.

We performed targeted numerical simulations of the periastron passage (instead of a continuous integration for 9000 orbits) to find out the properties of the $d_{\omega}$ distribution. For a fixed input velocity $\omega$, $20\,000$ Monte-Carlo integrations are performed with random input orientation angle $\theta$ covering a short interval around the periastron time. The resulting sampled distributions for $\omega=200\,n$ and $1000\,n$ are shown in Fig. 4, right. The width of the distributions strongly depends on the input velocity. There is also a pileup at both ends of the interval of allowable perturbations, so that impulses of greatest magnitude are more frequent than weak impulses.

These results explain why an ARMA model does not capture the properties of the $\omega(t)$ process well enough – the distribution of “noise” is drastically non-Gaussian. The GARCH model may be more successful in representing the $d_{\omega}$ process, because it does not require the distribution to be Gaussian. However, the variance as a model parameter is perhaps of limited use for the observed concave distribution of velocity updates.

V. Orbital evolution of high-eccentricity asteroids

Each periastron passage of a high-eccentricity asteroid results in an abrupt change of its rotation velocity. If we ignore the dissipative tidal forces in the asteroid and the star, the total energy of the closed system should be constant:

[TABLE]

A change in $\omega$ results in a change of the orbit semimajor axis $a$. Taking into account that the asteroid’s radius $R<<a$ and its mass $m_{2}<<m_{1}$, the change in $a$ is readily derived as

[TABLE]

The quadratic term $d_{\omega}^{2}$ may indicate that the orbit should be secularly shrinking. However, this is only true if the distribution of velocity updates $d_{\omega}$ were symmetric around zero, which is not the case. As Fig. 4 (left) indicates, the distribution may be biased toward negative values at high rotation velocities, compensating for the positive secular term or even causing the orbit to expand. Our limited numerical experiments are nor sufficient to verify this. It is safe to conclude that the orbit updates also constitute a random time process. The fractional updates of the orbit are found to be small for the model parameters in Table 1 (of the order of $10^{-9}$), but they can be much larger for super-earth exoplanets.

The total angular momentum should also be conserved in this two-body closed system:

[TABLE]

Updating the $a$, $e$, and $\omega$ by their increments, equating the change in the angular momentum to zero, and taking into account that $m_{2}<<m_{1}$, $d_{e}<<1$, $d_{a}<<1$, the following approximate equation obtains

[TABLE]

The fractional update of eccentricity $d_{e}/e$ is approximately linear in the update of rotation velocity $d_{\omega}$. The outcome of this stochastic variation is, again, uncertain, because it depends on the distribution of $d_{\omega}$. Because of the asymmetry observed in Fig. 4, left, the eccentricity is expected to grow when the asteroid rotates relatively slowly, but is likely to decline when the asteroid rapidly spins. The overall long-term evolution is unclear from our limited numerical integration.

VI. Discussion and conclusions

Exchange of angular momentum and energy between the orbit and the spin during close periastron passages of high-eccentricity asteroids causes the angular rotation velocity and the orbital elements to chaotically change. These processes can be formally modeled as random time processes. The long-term outcome for the orbit is unclear from our limited integration runs, but the indication is that the processes are weakly stationary. Any consideration of the orbital evolution is incomplete without the dissipative tidal force: the tides in the asteroid are probably negligible because of the explicit dependence of orbital action on $(R/a)^{5}$ (Makarov et al., 2018). Observational statistics for tumbling Solar system asteroids suggest that tidal dissipation inside asteroids provides long-term damping of the tumbling (chaotic, three-dimensional) motion, especially for high tensile rigidity objects (Pravec et al., 2005). However, this effect may be greatly diminished by the inverse proportionality of the tidal quality function to tidal frequency (Efroimsky, 2012) for rapidly rotating high-eccentricity asteroids around white dwarfs. The tides in the star may be more efficient on the timescale of stellar ages despite their expected low values of the tidal quality factor. Most stars rotate slower than the maximum angular velocity of the asteroid at the periastron, where practically all the action takes place. The tidal bulge raised on the star lags the direction to the asteroid, and the resulting torque gradually circularizes and shrinks the orbit.

We chose Proteus as the model object for our simulations (Table 1), which is larger, more massive, and more spherically symmetric than a typical Solar system asteroid, but the results are relevant to a broad range of objects from small comets to minor planets. The equation of motion (1) is independent of mass, but the orbital evolution is sensitive to the secondary mass. It is of interest to probe much greater triaxiality values, such as Hyperion at $\sigma\simeq 0.1$ and super-Earth exoplanets with low $\sigma$ ($\sim 10^{-5}$). Our preliminary computation for a Hyperion analog indicate that the behavior described in this paper and represented by Figs. 1, 3, and 4 is also found in smaller and more prolate bodies, but the velocity impulses are proportionally greater and their distribution around zero is more asymmetric. Such objects are likely to stochastically achieve break-up spin rates faster than the larger bodies. Conversely, they may become a viable source of white dwarf pollution at lower orbital eccentricity.

Although the behavior of rotation velocity is obviously chaotic and stochastic, modeling it with standard types of random processes is only partially successful. The reason for this is the peculiar distribution of velocity impulses, which is finite at any given velocity above the minimum prograde value. The probability of a velocity impulse sharply increases toward the boundaries of the distribution, which are not symmetric around zero. As a result, the process is weakly stationary and somewhat self-regulating. If the spin rate becomes low, the range of allowable impulses tends to zero, and the asteroid may linger in this quiescent state for an extended period of time. The high-velocity end of the distribution, on the other hand, does not seem to be bounded. The spin rates may become very high, but the probability of that is small, and it may take a long time for the asteroid to wander into this regime.

After the asteroid reaches a particular spin threshold, it will be disrupted by the centrifugal force. Our calculations illustrate that rotational fission may be more efficient than tidal disruption, especially for the largest asteroids and minor planets, because fission is triggered at lower values of eccentricity. The consequences for polluted white dwarfs are two-fold: (1) debris produced from break-up would populate a radial region which extends beyond the Roche radius, and (2) the orbital eccentricity which metal-polluting asteroids need to acheive is slightly lower than previously thought. We caution that we are unable to obtain a likelihood of rotational fission (versus tidal disruption) nor predict a specific asteroid’s long-term evolution, due to both computational limitations and the stochasticity of the systems. Nevertheless, this study represents a first step in the exploration of the coupled spin and orbital evolution of the dominant source of the observational signatures of the vast majority of known white dwarf planetary systems.

Acknowledgments

DV gratefully acknowledges the support of the STFC via an Ernest Rutherford Fellowship (grant ST/P003850/1).