Orbital relaxation and excitation of planets tidally interacting with white dwarfs

TL;DR

This paper presents a new method to analyze how planets near white dwarfs are tidally affected, determining their potential inward or outward drift and destruction, with implications for understanding white dwarf planetary systems.

Contribution

It introduces a self-consistent secular approach to model tidal interactions without constant lag assumptions, considering diverse planetary rheologies and spin states.

Findings

Massive Super-Earths are more prone to destruction than minor planets.

Low-viscosity planets are destroyed more easily than high-viscosity ones.

The boundary between planetary survival and destruction is fractal and chaotic.

Abstract

Observational evidence of white dwarf planetary systems is dominated by the remains of exo-asteroids through accreted metals, debris discs, and orbiting planetesimals. However, exo-planets in these systems play crucial roles as perturbing agents, and can themselves be perturbed close to the white dwarf Roche radius. Here, we illustrate a procedure for computing the tidal interaction between a white dwarf and a near-spherical solid planet. This method determines the planet's inward and/or outward drift, and whether the planet will reach the Roche radius and be destroyed. We avoid constant tidal lag formulations and instead employ the self-consistent secular Darwin-Kaula expansions from Bou\'{e} & Efroimsky (2019), which feature an arbitrary frequency dependence on the quality functions. We adopt wide ranges of dynamic viscosities and spin rates for the planet in order to straddle many possible outcomes, and provide a foundation for the future study of individual systems with known or assumed rheologies. We find that: (i) massive Super-Earths are destroyed more readily than minor planets (such as the ones orbiting WD 1145+017 and SDSS J1228+1040), (ii) low-viscosity planets are destroyed more easily than high-viscosity planets, and (iii) the boundary between survival and destruction is likely to be fractal and chaotic.

Tables

Table 1. Table 1: Self-consistent estimations of DA white dwarf dynamic (or shear) viscosities in the convective zone (column 3), masses (column 4), and radii (column 5) for given values of effective temperatures (column 1), and stellar luminosities (column 2). A surface gravity of g ⋆ = 10 6 subscript 𝑔 ⋆ superscript 10 6 g_{\star}=10^{6} m/s 2 was assumed. Values were taken at a depth corresponding to a convection zone mass above that layer which is ten orders of magnitude smaller than the white dwarf mass.

$T_{\mathrm{eff}}^{\star }$/K	$L_{\star }$/$L_{\odot }$	${\eta }^{\star }$/(Pa$\cdot$s)	$M_{\star }$/$M_{\odot }$	$R_{\star }$/$R_{\odot }$
5000	$9\times {10}^{-4}$	$2.49\times {10}^6$	0.579	0.01260
5500	$1.3\times {10}^{-3}$	$4.02\times {10}^6$	0.583	0.01264
6000	$1.9\times {10}^{-3}$	$5.46\times {10}^6$	0.588	0.01269
6500	$2.6\times {10}^{-3}$	$6.85\times {10}^6$	0.590	0.01272
7000	$3.5\times {10}^{-3}$	$8.23\times {10}^6$	0.592	0.01273
7500	$4.7\times {10}^{-3}$	$9.38\times {10}^6$	0.593	0.01275
8000	$6.1\times {10}^{-3}$	$9.95\times {10}^6$	0.596	0.01278
8500	$7.7\times {10}^{-3}$	$1.20\times {10}^7$	0.596	0.01278
9000	$9.8\times {10}^{-3}$	$1.33\times {10}^7$	0.598	0.01280

Table 2. Table 2: Unstylised and stylised Roman variables used in this paper, with notation that assumes a stellar primary and planetary secondary. No functional dependencies are shown. Variables with a ⋆ ⋆ \star subscript or superscript refer to a stellar physical property.

Variable	Explanation	Units	Reference
$a$	Semimajor axis of the mutual orbit	length	equation (21)
$ℬ$	Auxiliary variable	1/pressure	equation (6)
${ℬ}^{\star }$	Auxiliary variable	1/pressure	equation (8)
$d_{\star }$	Typical lengthscale of stellar convective flows	distance	equation (12)
$e$	Eccentricity of the mutual orbit	dimensionless	equation (22)
$F$	Inclination function	dimensionless	equation (20)
$g$	Acceleration due to gravity on planet surface	length/time2
$g_{\star }$	Acceleration due to gravity on stellar surface	length/time2
$G$	Eccentricity function	dimensionless	equation (17)
$𝒢$	Gravitational constant	length3/(mass$\times$time2)
$i$	Inclination of the mutual orbit with respect to the planet’s equator	angle	equation (23)
$i^{\prime }$	Inclination of the mutual orbit with respect to the star’s equator	angle	equation (24)
$J$	Bessel function of the first kind	dimensionless
$𝒥$	Compliance of the planet	1/pressure
${𝒥}^{\star }$	Compliance of the star	1/pressure
$k$	Love number of the planet	angle/time
$k^{\star }$	Love number of the star	angle/time
$l$	Positive integer used as Love number index (“degree”)	dimensionless
$L_{\star }$	Luminosity of star	power	Table 1
$m$	Integer used as tidal mode index and inclination function index (“order”)	dimensionless
$M$	Mass of planet	mass
$M_{\star }$	Mass of star	mass
$n$	Anomalistic mean motion	angle/time
$p$	Integer for tidal mode index, eccentricity function index, & inclination function index	dimensionless
$q$	Integer for tidal mode index and eccentricity function index	dimensionless
$R$	Radius of planet	distance
$R_{\star }$	Radius of star	distance
$s$	Integer used as summation index for the eccentricity function	dimensionless
$s_1$	Integer used as summation index limit for the eccentricity function	dimensionless
$t$	time	time
$T_{\mathrm{eff}}^{\star }$	Effective temperature of star	temperature	Table 1
$𝒯$	Torque	mass$\times$length2/time2
$u$	Integer used as summation index for the eccentricity and inclination functions	dimensionless
$v_{\star }$	Typical speed of stellar supergranules	length/time	equation (12)

Table 3. Table 3: As in Table 1, but with Greek variables.

Variable	Explanation	Units	Reference
$\delta$	Delta function	dimensionless
$ϵ$	Phase lag of the planet	dimensionless
${ϵ}^{\star }$	Phase lag of the star	dimensionless
$\eta$	Dynamic viscosity of the planet	pressure$\times$time
${\eta }^{\star }$	Dynamic viscosity of the star	pressure$\times$time	equation (12)
$\theta$	Rotation angle about the instantaneous shortest axis of planet	angle	equation (25)
${\theta }^{\star }$	Rotation angle about the instantaneous shortest axis of star	angle
$\xi$	Coefficient for planet’s moment of inertia ($=\frac{2}{5}$ for a homogeneous sphere)	dimensionless
${\xi }^{\star }$	Coefficient for star’s moment of inertia ($=\frac{2}{5}$ for a homogeneous sphere)	dimensionless
$\rho$	Density of planet	mass/length3
${\rho }^{\star }$	Density of star	mass/length3
${\sigma }_a$	Semimajor axis decrease divided by the initial distance to the Roche radius	dimensionless	equation (14)
${\sigma }_e$	Eccentricity decrease divided by the initial eccentricity	dimensionless	equation (15)
${\sigma }_i$	Inclination decrease divided by the initial inclination	dimensionless	equation (16)
${\chi }_{lmpq}$	Positive definite physical forcing frequency in the planet	angle/time	equation (2)
${\chi }_{lmpq}^{\star }$	Positive definite physical forcing frequency in the star	angle/time	equation (4)
${\omega }_{lmpq}$	Fourier tidal mode in the planet	angle/time	equation (1)
${\omega }_{lmpq}^{\star }$	Fourier tidal mode in the star	angle/time	equation (3)
$k_l\sin {ϵ}_l$	Quality function of the planet	angle/time	equation (5)
$k_l^{\star }\sin {ϵ}_l^{\star }$	Quality function of the star	angle/time	equation (7)

Equations

$${\omega }_{lmpq}\approx \left(l-2p+q\right)n-m\frac{d\theta }{dt},$$

$${\chi }_{lmpq}=\left|{\omega }_{lmpq}\right|.$$

$${\omega }_{lmpq}^{\star }\approx \left(l-2p+q\right)n-m\frac{d{\theta }^{\star }}{dt}$$

$${\chi }_{lmpq}^{\star }=\left|{\omega }_{lmpq}^{\star }\right|$$

$$k_l({\omega }_{lmpq})\sin {ϵ}_l({\omega }_{lmpq})$$

$$=\frac{3\mathrm{Sgn}\left({\omega }_{lmpq}\right)}{2\left(l-1\right)}\frac{\frac{{\mathcal{B}}_l}{\eta {\chi }_{lmpq}}}{{\left(\mathcal{J}+{\mathcal{B}}_l\right)}^2+{\left(\eta {\chi }_{lmpq}\right)}^{-2}}$$

$${\mathcal{B}}_l=\frac{3\left(2l^2+4l+3\right)}{4l\pi \mathcal{G}{\rho }^2{R}^2}=\frac{2l^2+4l+3}{lg\rho R}.$$

$$k_l^{\star }({\omega }_{lmpq}^{\star })\sin {ϵ}_l^{\star }({\omega }_{lmpq}^{\star })$$

$$=\frac{3\mathrm{Sgn}\left({\omega }_{lmpq}^{\star }\right)}{2\left(l-1\right)}\frac{\frac{{\mathcal{B}}_l^{\star }}{{\eta }^{\star }{\chi }_{lmpq}^{\star }}}{{\left({\mathcal{J}}^{\star }+{\mathcal{B}}_l^{\star }\right)}^2+{\left({\eta }^{\star }{\chi }_{lmpq}^{\star }\right)}^{-2}}$$

$${\mathcal{B}}_l^{\star }=\frac{3\left(2l^2+4l+3\right)}{4l\pi \mathcal{G}{\rho }_{\star }^2{R}_{\star }^2}=\frac{2l^2+4l+3}{l{g}_{\star }{\rho }_{\star }{R}_{\star }},$$

$$⟨⟨\frac{da}{dt}⟩⟩\sim R^{2l},$$

$$⟨⟨\frac{de}{dt}⟩⟩\sim R^{2l},$$

$$⟨⟨\frac{di}{dt}⟩⟩\sim R^{2l},$$

$${\eta }^{\star }={\rho }_{\star }{d}_{\star }{v}_{\star }.$$

$$\frac{r_{\mathrm{Roche}}}{R_{\odot }}\approx {\left(\frac{\rho }{3{\mathrm{g}\mathrm{cm}}^{-3}}\right)}^{-1/3}.$$

$${\sigma }_a\equiv \frac{\mathrm{final}\mathrm{distance}-\mathrm{initial}\mathrm{distance}}{\mathrm{initial}\mathrm{distance}-\mathrm{Roche}\mathrm{radius}}.$$

$${\sigma }_e\equiv \frac{\mathrm{eccentricity}\mathrm{change}}{\mathrm{initial}\mathrm{eccentricity}}$$

$${\sigma }_i\equiv \frac{\mathrm{inclination}\mathrm{change}}{\mathrm{initial}\mathrm{inclination}}.$$

$$G_{lpq}(e)=\left[1+\left(\frac{e}{1+\sqrt{1-e^{2}}}\right)^{2}\right]^{l}\sum_{s=0}^{s_{1}}\sum_{u=0}^{u_{1}}\left({\begin{array}[]{cc}2p-2l\\ s\\ \end{array}}\right)\left({\begin{array}[]{cc}-2p\\ u\\ \end{array}}\right)\left(-\frac{e}{1+\sqrt{1-e^{2}}}\right)^{s+u}J_{q-s+u}\left(\left(l-2p+q\right)e\right)$$

$$s_1=2p-2l,\mathrm{if}p-l\ge 0$$

$$s_1=\infty ,\mathrm{if}p-l<0$$

$$u_1=-2p,\mathrm{if}p\le 0$$

$$u_1=\infty ,\mathrm{if}p>0.$$

$$F_{lmp}(i)=\frac{\left(l+m\right)!}{2^{l}p!\left(l-p\right)!}\sum_{u={\rm max}\left(0,\ l-2p-m\right)}^{{\rm min}\left(2l-2p,\ l-m\right)}\left(-1\right)^{u}\left({\begin{array}[]{cc}2l-2p\\ u\\ \end{array}}\right)\left({\begin{array}[]{cc}2p\\ l-m-u\\ \end{array}}\right)\cos^{3l-m-2p-2u}\left(\frac{i}{2}\right)\sin^{m-l+2p+2u}\left(\frac{i}{2}\right).$$

$$⟨⟨\frac{da}{dt}⟩⟩=-2an\sum _{l=2}^{\infty }\sum _{m=0}^l\frac{\left(l-m\right)!}{\left(l+m\right)!}\left(2-{\delta }_{0m}\right)\sum _{p=0}^l\sum _{q=-\infty }^{\infty }{G}_{lpq}^2(e)\left(l-2p+q\right)$$

$$\times \left[{\left(\frac{R}{a}\right)}^{2l+1}\frac{M_{\star }}{M}{F}_{lmp}^2(i)k_l({\omega }_{lmpq})\sin \left[{ϵ}_l({\omega }_{lmpq})\right]+{\left(\frac{R_{\star }}{a}\right)}^{2l+1}\frac{M}{M_{\star }}{F}_{lmp}^2(i^{\prime })k_l^{\star }({\omega }_{lmpq}^{\star })\sin \left[{ϵ}_l^{\star }({\omega }_{lmpq}^{\star })\right]\right].$$

$$⟨⟨\frac{de}{dt}⟩⟩=-\frac{1-e^2}{e}n\sum _{l=2}^{\infty }\sum _{m=0}^l\frac{\left(l-m\right)!}{\left(l+m\right)!}\left(2-{\delta }_{0m}\right)\sum _{p=0}^l\sum _{q=-\infty }^{\infty }{G}_{lpq}^2(e)\left(l-2p+q\right)$$

$$\times \left[{\left(\frac{R}{a}\right)}^{2l+1}\frac{M_{\star }}{M}{F}_{lmp}^2(i)k_l({\omega }_{lmpq})\sin \left[{ϵ}_l({\omega }_{lmpq})\right]+{\left(\frac{R_{\star }}{a}\right)}^{2l+1}\frac{M}{M_{\star }}{F}_{lmp}^2(i^{\prime })k_l^{\star }({\omega }_{lmpq}^{\star })\sin \left[{ϵ}_l^{\star }({\omega }_{lmpq}^{\star })\right]\right]$$

$$+\frac{\sqrt{1-e^2}}{e}n\sum _{l=2}^{\infty }\sum _{m=0}^l\frac{\left(l-m\right)!}{\left(l+m\right)!}\left(2-{\delta }_{0m}\right)\sum _{p=0}^l\sum _{q=-\infty }^{\infty }{G}_{lpq}^2(e)\left(l-2p\right)$$

$$\times \left[{\left(\frac{R}{a}\right)}^{2l+1}\frac{M_{\star }}{M}{F}_{lmp}^2(i)k_l({\omega }_{lmpq})\sin \left[{ϵ}_l({\omega }_{lmpq})\right]+{\left(\frac{R_{\star }}{a}\right)}^{2l+1}\frac{M}{M_{\star }}{F}_{lmp}^2(i^{\prime })k_l^{\star }({\omega }_{lmpq}^{\star })\sin \left[{ϵ}_l^{\star }({\omega }_{lmpq}^{\star })\right]\right].$$

Full text

††pubyear: 2019††pagerange: Orbital relaxation and excitation of planets tidally interacting with white dwarfs–C.2

1 Introduction

The discovery of bona fide asteroids surrounding two white dwarfs on low-to-moderate eccentric orbits (Vanderburg et al., 2015; Manser et al., 2019) has refocused efforts to understand the fate of planetary systems. Debris from destroyed planetesimals both surround and accrete onto the white dwarf. This destruction is preceded by the close gravitational interaction between a white dwarf and a planetary body, a process which has been investigated only sparsely. Hence, the near-complete dearth of dedicated investigations into tidal effects between white dwarfs and planets necessitate a remedy; this paper represents an initial step towards addressing this deficiency.

1.1 Planetary system evolution

When a main-sequence exo-planet host becomes a giant branch star, violent changes ensue, transforming the planetary system (Veras, 2016). The star loses most of its mass, expanding the orbits of surviving planets, moons, asteroids and comets (Omarov, 1962; Hadjidemetriou, 1963; Veras et al., 2011; Adams et al., 2013; Veras et al., 2013a). This expansion may trigger gravitational instabilities (Debes & Sigurdsson, 2002) – even without the presence of any stellar companions – with important consequences for destruction, escape and, in general, orbital rearrangement (Bonsor et al., 2011; Debes et al., 2012; Portegies Zwart, 2013; Mustill et al., 2013; Veras et al., 2013b; Voyatzis et al., 2013; Frewen & Hansen, 2014; Mustill et al., 2014; Veras et al., 2016; Mustill et al., 2018; Smallwood et al., 2018).

Further, giant branch stars exhibit luminosities which exceed the Sun’s by several orders of magnitude. Hence, pebbles and dust are propelled and dragged at different rates and can also cross stability boundaries with planets (Dong et al., 2010; Veras et al., 2015b). The consequences for asteroids are perhaps more dramatic: they can be spun up to the point of rotational fission (Veras et al., 2014c) through the YORP (Yarkvosky-O’Keefe-Radzievskii-Paddack) effect. The survivors can be thrust inward or outward, speeding past planets, due to a supercharged Yarkovsky effect (Veras et al., 2015b, 2019). The end result is a widely dispersed debris field ranging from pebbles to planets at distances from a few au to hundreds or thousands of au.

The inner boundary of this debris is set by the physical expansion of the giant branch star itself, which can reach several au in size. Tidal effects could even draw in and engulf planets which are further away than the stellar radius. This critical tidal engulfment distance has been the subject of much interest and investigation (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). The actual destructive spiral-in process and timescale for planets within the stellar envelope has also been estimated (Jia & Spruit, 2018; MacLeod, Cantiello & Soares-Furtado, 2018), as well as the consequences for the giant star (Nelemans & Tauris, 1998; Massarotti, 2008; Carlberg et al., 2012, 2013).

The star eventually transitions into a white dwarf, which is comparable in size to the Earth but has a Roche radius which extends outward to about one Solar radius for a wide range of orbiting secondary structures (Veras et al., 2017). Within this Roche radius observations are abundant. Metallic debris is seen in the white dwarf photosphere (e.g. Zuckerman et al., 2007; Gänsicke et al., 2012; Jura & Young, 2014; Koester et al., 2014; Harrison et al., 2018; Hollands et al., 2018) and just outside is the presence of Solar radii-scale debris discs (Farihi, 2016).

Perhaps most eye-catching are the discoveries of minor planets orbiting white dwarfs. At least one planetesimal is currently disintegrating around WD 1145+017 with an orbital period of about 4.5 hours (Vanderburg et al., 2015) and one other planetesimal is intact and embedded inside the disc surrounding SDSS J1228+1040 with an orbital period of about 2.06 hours (Manser et al., 2019). The orbits of the disintegrating planetesimals are thought to be nearly circular (Gurri et al., 2017; Veras et al., 2017) whereas the embedded planetesimal’s eccentricity is likely to be many tenths (potentially 0.53), and accompanies a disc-based intensity pattern with a likely eccentricity of 0.4 (Manser et al., 2019).

How objects can be emplaced at one Solar radius from several au – and then circularised – is still subject to debate. Although dust particles can be drawn into the white dwarf through Poynting-Robertson drag, larger objects need to be gravitationally perturbed onto highly eccentric ($>0.99$) orbits. In order to reach the white dwarf’s Roche radius, Graham et al. (1990) and Jura (2003) proposed that these objects are perturbed there (from distances of a least a few au) and then break up into discs, an idea which has been realised through numerical simulations (Debes et al., 2012; Veras et al., 2014a, 2015a; Malamud et al., 2019) and purely analytic formulations (Wyatt et al., 2014; Brown et al., 2017; Kenyon & Bromley, 2017a, b). Veras et al. (2018b) argued extensively that the presence of a planet is a near-necessity to perturb asteroids close to the white dwarf in a single-star system.

Nevertheless, the major planets themselves can be perturbed close to and onto the white dwarf, particularly when the planets are comparable in mass to one another (Veras & Gänsicke, 2015; Veras et al., 2016); a deficiency with those two studies is the absence of tidal modelling, which can introduce uncertainties in the orbital behaviour of a planet passing close to the white dwarf. These planets would play major roles in perturbing smaller objects and dust onto the star, with implications for fragmentation, sublimation and the resulting debris size distributions (Wyatt et al., 2014; Brown et al., 2017).

Additional motivations for studying tides between a white dwarf and planet include the prospect of the formation of second- or third-generation planets, which may even arise from second- or third-generation debris discs (Perets, 2011; Schleicher & Dreizler, 2014; Völschow et al., 2014; Hogg et al., 2018; van Lieshout et al., 2018). Such planets would reside on near-circular orbits, and may be out of reach of detection. Further, the destruction of a planet, either through collisions with other planets or by intersecting the white dwarf Roche radius, could result in easily detectable events (Bear & Soker, 2013; Vanderburg et al., 2015) and smaller planets orbiting near the white dwarf. Finally, moons of planets can be stripped, perturbed towards the white dwarf, and as a result interact with the white dwarf in close proximity (Payne et al., 2016, 2017).

In the year 2018, the number of known white dwarfs increased by an order of magnitude (Gentile Fusillo et al., 2019). An accompanying increase in the number of previously known white dwarf planetary systems (Kleinman et al., 2013; Kepler et al., 2015, 2016; Hollands et al., 2017) is likely to follow with the next generation of spectroscopic initiatives and other missions. This current study may be particularly impactful if planets are discovered orbiting white dwarfs either spectroscopically from ground-based facilities, or through TESS, LSST (Cortes & Kipping, 2018; Lund et al., 2018), Gaia (Perryman et al., 2014) and LISA (Steffen, Wu & Larson, 2018; Tamanini & Danielski, 2018). Whether tidal effects had or will have a role in these systems can then be estimated through orbital proxies and stability analyses (Veras et al., 2018a).

1.2 Tidal formulations

Despite these numerous motivations, tidal investigations of white dwarfs have been limited to star-star interactions (Fuller & Lai, 2011, 2012, 2013, 2014; Valsecchi et al., 2012; Sravan et al., 2014; Vick et al., 2017; McNeill et al., 2019). In contrast, abundant studies of the main-sequence and giant-branch phases of evolution have analysed the star-planet interaction, and have adopted a variety of approaches with both equilibrium and dynamical tides.

Two commonly-used star-planet tidal models are the CPL (constant phase lag) model pioneered by MacDonald (1964) and Goldreich (1966), and the CTL (constant time lag) model suggested in Hut (1981) and Eggleton et al. (1998). Introduced for the ease of analytical treatment and for illustrative purposes, these models have very limited application in quantitative analysis. For example, the CPL model is both physically and mathematically inconsistent (Efroimsky & Makarov, 2013). The CTL model may (arguably) render a tolerable approximation for stars, but its application to terrestrial bodies is possible only in the situations when these bodies are hot and plastic (Makarov, 2015). Under normal conditions, these models are applicable neither to solids, nor to solids with partial melt.

Generally, employment of an incorrect tidal model results in distortions of the true tidal evolution timescales and fails to describe correctly the process of tidal capture in spin-orbit resonances (for details, see Noyelles et al., 2014). Also, the use of simplistic models sometimes leads to qualitatively wrong conclusions, like with the emergence of the so-called pseudo-synchronous rotation state (Makarov & Efroimsky, 2013). Perhaps most importantly, the use of an incorrect tidal model prevents the correct assessment of the tidal heat produced in a body.

Although some of the hitherto published far-reaching conclusions on tidal dynamics of celestial bodies were based on simplistic ad hoc models, here we wish to avoid such flaws. Hence, we rely on a versatile formalism that can accurately reflect the increasing detail with which the community has been characterising planets and stars.

In order to accommodate such detail, consistent analytical modeling of linear bodily tides should comprise two steps: (1) First, Fourier-expansion of both the perturbing potential and the induced incremental potential of the tidally perturbed body. Each of these two expansions is an infinite sum over the tidal Fourier modes $\omega_{lmpq}$ which are numbered with aid of some integers $\,l,\,m,\,p,\,q\,$ (so that, technically, each such expansion becomes a sum over $l,\,m,\,p,\,q$). Appendix C describes the nature of these indices in more detail. (2) Second, linking each Fourier component of the incremental tidal potential to a corresponding Fourier component of the perturbing potential. This link implies establishing — for each Fourier mode $\omega_{lmpq}$ — both the phase lag $\epsilon_{l}(\omega_{lmpq})$ and the ratio of the magnitudes, termed the dynamical Love number and denoted with $k_{l}(\omega_{lmpq})$. Owing to the interplay of rheology and self-gravitation, the phase lags and Love numbers possess nontrivial frequency dependencies; see Section 2.2 below.

The construction of a consistent mathematical theory of bodily tides was begun by Darwin (1879) who wrote down partial sums of the Fourier expansions of both the perturbing potential and the incremental potential of a tidally perturbed sphere. A great development of this theory was achieved by Kaula (1964) who managed to derive complete series expansions for these potentials. The paper by Efroimsky & Makarov (2013) offers a relatively simple introduction to the Darwin-Kaula theory and also explains how tidal friction and lagging should be built into that theory.

With these tools at hand, it is possible to develop similar Fourier-type expansions for the orbital elements’ rates (Kaula, 1964) and the tidal torque (Efroimsky, 2012a). Fortuitously, Boué & Efroimsky (2019) recently executed a complete rederivation of the expressions for the Keplerian elements’ tidal rates. They pointed out a minor omission in Kaula’s old treatment, and also explained that Kaula’s expression for the inclination rate contained a more important flaw. We use the Boué & Efroimsky (2019) formulae in this paper. For readability purposes for a wide audience, these formulae and a few others have been relegated to Appendices B-C.

Computation of these equations, however, pose a challenge. As the eccentricity of the orbit approaches unity, the number of terms which should be included to obtain an accurate solution increases nonlinearly. Calculations show that at an orbital eccentricity of $0.95$, more than 700 terms of $G^{2}_{31q}$ (see Appendix A for the definitions of symbols used throughout the paper) are greater than 0.01 times the maximum value.

Therefore, because the first approach of a planet to a white dwarf must be on an orbit with an eccentricity which exceeds 0.99, we cannot yet model the initial damping phase. Also, the planetary rotation becomes chaotic and unpredictable after the first passage close to the star, which makes the tidal evolution scenarios essentially probabilistic. Although eccentricities of about 0.8 can be feasibly modelled for individual systems, here we wish to perform an ensemble of simulations over a broad region of parameter space and to resemble the low-to-moderate eccentricities of the known planetesimals orbiting WD 1145+017 (Vanderburg et al., 2015) and SDSS J1228+1040 (Manser et al., 2019). Therefore, we consider orbital eccentricities in the range of $10^{-4}$-$0.4$. In contrast to high eccentricities, high inclinations do not impose computational restrictions.

1.3 Plan for paper

In Section 2 and Appendices A–C, we set up and describe our planet-star setup, the equations of motion, and the meaning of the various components of those equations. Appendix A in particular is composed of two tables which describe all of the variables used in this paper. In Section 3, we identify results which can be gleaned through an inspection of the equations alone. Section 4 is devoted to elucidating our initial conditions and parameters for the simulations. In Section 5, we run the simulations and report the results. We conclude in Section 6.

2 Setup

In order to carry out our orbital simulations, we rely on the formalism presented in Boué & Efroimsky (2019). That investigation corrects earlier work and provides equations for the evolution of semimajor axis, eccentricity and inclination of the mutual orbit of two bodies due to tides raised on both. The formalism in that paper is general enough so that (i) either or both bodies may be stars or planets, and (ii) the expansions do not diverge for high values of eccentricity nor inclination, unlike for some classical expansions of planetary disturbing functions for the three-body problem (Murray & Dermott, 1999).

Our consideration will be limited to terrestrial planets (like exo-asteroids and exo-Earths) and ice giants (like exo-Neptunes) as long as they do not have any surface continents nor surface oceans. Internal oceans are allowed because tidal dissipation in internal oceans is much lower than in mantles (Table 3 of Chen et al., 2014). Surface oceans and continents are not allowed because tidal dissipation in such planets is dominated by friction in shallow seas and cannot be described by rheological models borrowed from continuum mechanics. Our method is not applicable to gas giants (like exo-Jupiters); for descriptions of tidal friction on such objects, see the review by Mathis (2018) and references therein.

2.1 Orbital evolution owing to tides

We consider a planet orbiting a white dwarf. Physical variables with a $\star$ subscript or superscript refer to a stellar physical property, and those without a superscript or subscript refer to a planetary physical property. As usual, the notations $a$ and $e$ refer to the semimajor axis and eccentricity of the mutual orbit. For inclination, $i$ refers to the orbit’s inclination with respect to the planetary equator, while $i^{\prime}$ denotes the orbit’s inclination with respect to the stellar equator. With this notation, the rates of $a$, $e$, $i$, and $i^{\,\prime}$ are given by equations (21-24) in Appendix C. We do not consider the time evolution of the arguments of pericentre, longitudes of ascending node, and times of pericentre passages; the Lagrange-type orbital equations describing their rates are written down in Boué & Efroimsky (2019).

The four orbital equations (21–24) are secular: they are not explicit functions of the mean anomaly, true anomaly, nor mean longitude. This secular property allows us to model the system for longer times than if we considered the oscillatory and libratory behaviour on a per-orbit basis. We have used the bra-ket delimiters to indicate the secular nature of these equations. Also, these four evolution rates all receive double sets of delimiters because averaging is carried out not only over one orbit but also over one cycle of the apsidal precession.

The expansions for each of the four orbital rates contain the inclination functions $\,F_{lmp}(i)\,$ and the eccentricity functions $\,G_{lpq}(e)\,$, which are given in equations (17) and (20). These functions are not dependent on any physical parameters, and contain only a single orbital parameter each. Therefore, we precompute these functions before running any simulations. The values of $\,F_{lmp}(i)\,$ are straightforwardly computed for any inclination: the finite limits of the summation in equation (20) allow us to produce compact explicit formulae.

However, convergence of the $\,G_{lpq}(e)\,$ values becomes computationally onerous for eccentricities near unity, because of the often infinite upper limit in one of the summations in equation (17). To speed up computations, Noyelles et al. (2014) tabulated values of $\,G_{lpq}(e)\,$ for different eccentricities and values of $l$, $p$, and $q$, and used the resulting look-up table in their simulations. Here, instead, for each given set of $l$, $p$, and $q$, we fit $G_{lpq}(e)$ to an explicit function of $e$. We then insert these explicit functions into the simulation arrays so that no lookup as a function of eccentricity is necessary.

2.2 Link to physical parameters

In Appendix C, we write down the expressions (21–24) for the tidal evolution rates of four orbital parameters: the mutual orbit’s semimajor axis $a$, eccentricity $e$, inclination $i$ with respect to the planetary equator, and inclination $i^{\,\prime}$ with respect to the stellar equator. Each of those expressions is an infinite sum — over some integers $\,l,\,m,\,p,\,q\,$ — of terms proportional to the so-called “quality functions”. A quality function of degree $l$ equals the product of the degree-$l$ dynamical Love number and the sine of the degree-$l$ phase lag, both understood as functions of an $lmpq$ tidal Fourier mode.

Specifically, the quality functions of the planet, $k_{l}(\omega_{lmpq})\sin\epsilon_{l}(\omega_{lmpq})$, are functions of the tidal Fourier modes111 Be mindful that the functional form of a quality function is parametrised with the degree $\,l\,$ only, while the dependence on the other three integers comes through the arguments $\omega_{lmpq}$ or $\omega_{lmpq}^{\star}$. This simplification is available only if we approximate the body with a homogeneous sphere. The general case is more complex. For example, if we take into account the permanent oblateness of a body, its quality functions will be parametrised not only with the degree $l$, but also with the order $m$, and will read as $k_{lm}(\omega_{lmpq})\sin\epsilon_{lm}(\omega_{lmpq})\,$ and $k_{lm}^{\star}(\omega^{\star}_{lmpq})\sin\epsilon^{\star}_{lm}(\omega^{\star}_{lmpq})$. In this paper, however, we will not delve into this level of complexity.

[TABLE]

which are exerted on the planet by the star. Their absolute values are the actual forcing frequencies of deformation in the planet: 222 The physical meaning of the modes $\,\omega_{lmpq}\,$ and frequencies $\,\chi_{lmpq}\,$ can be easily understood through Eq. 15 in Efroimsky & Makarov (2013).

[TABLE]

Similarly, the quality functions of the star, $k_{l}^{\star}(\omega_{lmpq}^{\star})\sin\epsilon_{l}^{\star}(\omega_{lmpq}^{\star})$, are functions of the modes

[TABLE]

which are excited inside the star by the planet. The absolute values of these modes,

[TABLE]

are the physical frequencies of deformation inside the star.

For both the planet and the star, the shapes of the quality functions are determined by interplay of self-gravitation and rheology in the appropriate body. Therefore, for each of the two bodies, the shape of an $l^{\rm th}$ quality function depends on this body’s key kinematic and physical parameters, such as the spin rate, density, radius, surface gravity, shear viscosity, and shear compliance (which is the reciprocal of shear rigidity, or shear elastic modulus).

Assume for simplicity that the tidally-perturbed celestial body is homogeneous and near-spherical, and also that its rheology is linear (which is usually the case). A general expression for a degree-$l\,$ quality function of such a body was derived in Efroimsky (2012a, Eq. 169) and Efroimsky (2015, Eq. 40a). That expression is valid for an essentially arbitrary linear rheology. Mathematically, that expression is a simple algebraic functional of two functions representing the real and imaginary parts of the complex compliance of the material of which the celestial body is composed. Recall that these parts of the complex compliance contain the entire information about the (linear) rheology.

Under the additional assumption that the planet’s material adheres to the Maxwell model,333 This assumption restricts our planets to those where most of the tidal dissipation is taking place in the rocky material. The tidal response of terrestrial mantles is viscoelastic (Maxwell) at sufficiently low frequencies (for the Earth, up to about the Chandler frequency). At higher frequencies, the behaviour of these mantles is more adequately described by the Andrade model (Efroimsky, 2012a, b) or, even better, by its generalisation named the Sundberg-Cooper model (Renaud & Henning, 2018). It is very important that both these models are firmly rooted in physics, and reflect actual physical mechanisms of friction emerging at seismic frequencies. It is also convenient that, mathematically, they both are extensions of the Maxwell model. Because the current paper is our pilot tidal publication on the topic of planets around white dwarfs, we take the liberty of omitting some technicalities with the understanding that these may be built in the theory later. Among these omitted items is the switch to more complex rheologies which is to be left for future work.  the quality functions of the planet become Eqs. 31 and 50a from Efroimsky (2015):

[TABLE]

[TABLE]

where

[TABLE]

Here $\cal G\,$ is Newton’s gravitational constant, while $\rho$, $R$, $g$, $\eta$, and ${\mathcal{J}}$ are the planet’s mean density, radius, surface gravity, shear viscosity, and shear compliance, correspondingly.

Under the assumption that the star, too, has a Maxwell rheology, 444 Stars are commonly assumed to be viscous (Stix, 2002). Simultaneously, there exist theoretical indications that they may possess magnetic rigidity (Williams, 2004, 2005, 2006; Ogilvie, 2008; Garaud et al., 2010). We are hence motivated to treat a stellar material, at large, as a Maxwell body, though with a very small shear compliance $\cal{J}^{\star}$. Not knowing how small the compliance value may be, we assume that $\cal{J}^{\star}$ is much smaller than $1/(\eta^{\star}\chi^{\star})$, where $\eta^{\star}$ and $\chi^{\star}$ are the mean viscosity of the star and a typical tidal frequency, correspondingly. This approximation simplifies the expression for the stellar quality functions (Efroimsky, 2015, Eq. 60).  the quality functions of a star approximated with a homogeneous Maxwell sphere are

[TABLE]

[TABLE]

where

[TABLE]

with $\rho^{\star}$, $R^{\star}$, $g^{\star}$, $\eta^{\star}$ and ${\mathcal{J}}^{\star}$ being the stellar mean density, radius, surface gravity, shear viscosity, and shear compliance, accordingly.

For a comprehensive list of variable definitions, see Tables A1–A2.

2.3 Spin equations of motion

The orbital evolution equations of motion (21–24) contain explicit dependences on the spin evolution of the planet and star ($d\theta/dt$ and $d\theta^{\star}/dt$) through equations (1–3). Therefore spin evolution must be solved self-consistently with orbital evolution. We provide the general spin evolution equations for both the star and planet in equations (25–26).

These equations represent the sum of the contributions from both the triaxial torque and tidal torque. The secular triaxial torque (as opposed to the non-secular version) is conventionally accepted to equal zero (equations 29–30) because the averaging (i) would otherwise be certain to yield a non-zero secular term only in physically questionable circumstances (synchronicity and nonzero net tilt) and (ii) might yield a non-zero secular term in other largely unexplored circumstances (e.g. due to the variation of the initial mean anomaly, argument of pericentre, and $\dot{\theta}$ from oblateness terms). The secular tidal torque is given by equations (27–28).

3 Properties of the equations of motion

Some results may be gleaned just through inspection of equations (21)–(26), without needing to run any numerical simluations.

3.1 Synchronisation

The spin and orbital element evolution dependencies on multiple variables suggest that a simple characterisation of spin-orbit resonances is difficult to attain. A direct comparison of the spin periods of the white dwarf and planet, along with their mutual separation, would alone be insufficient to claim synchronicity or not. Resonant trapping of an individual system is best devoted to a dedicated study (e.g. Noyelles et al. 2014 for Mercury or Makarov et al. 2012 for GJ 581d). Further, the secular (averaged) nature of the equations applied here might miss important librational behaviours which occur on orbital timescales, and can play a role in resonant capture.

3.2 Circularisation

Strictly, a perfectly circular orbit will remain circular in the midst of tidal interactions, as can be demonstrated formally by taking the limit of equation (22) as $e\rightarrow 0$. However, realistically, this situation does not occur: even a minuscule occasional interaction capable of generating a nonzero value of $e$, will, under some circumstances, initiate a secular growth in $e$. Therefore, equation (22) should be integrated simultaneously with the other orbital and spin equations regardless of how small $e$ is thought to become.

3.3 Planarisation

When the equators of both the planet and the white dwarf are perfectly aligned with the orbit, then this situation will remain unchanged in the midst of tidal interactions. Importantly, however, in this case (when $i=i^{\prime}=0$), the inclination $F_{lmp}$ functions do not vanish for all $l,m$ and $p$. Nevertheless, the coplanar case does help reduce computational cost because in that case, only those sets of $l,m$ and $p$ for which $F_{lmp}$ does not vanish need to be included in the summations in equations (21–22).

3.4 Contributions from stellar quality functions

The semimajor axis and eccentricity rates (21–22) feature distinct terms associated with each quality function. Orbital evolution is then said to be driven by planetary tides when the planetary quality function dominates, and driven by stellar tides when the stellar quality function dominates. When both tidally interacting bodies are of similar types (such as in a binary star system or a binary asteroid system), then a degree-$l\,$ term owing to the tides in one body and an analogous degree-$l\,$ terms owing to the tides in another body may make comparable contributions to the overall evolution.

However, for star-planet systems, the stellar quality function coefficients $M/M_{\star}$ and $(R_{\star}/a)^{2l+1}$ are small. The second term becomes even smaller for white dwarfs, because their radii are comparable to the Earth’s radius rather than the Sun’s radius. The smallness of these coefficients often — but not always — makes negligible the stellar quality function terms. Exceptions include when the stellar quality function itself is large, which occurs when the stellar viscosity is large and/or when the star spins quickly. One potential manifestation of this exception is an eccentricity increase when the stellar spin is faster than $(3/2)n$ (Boué & Efroimsky, 2019).

For the numerical investigations in this paper, the stellar quality function remains small — primarily because of our adopted white dwarf viscosity — greatly helping to facilitate the reduction of the phase space that we explore. Hence, because the equations for spin evolution (25–26) are each a function of just one of the quality functions, we find that the white dwarf spin rate changes negligibly due to tidal effects with a planet, and can be considered fixed. Similarly, each term in the evolution of $i^{\,\prime}$ (equation 24) is a function of the white dwarf quality function, but not the planetary quality function. Hence $i^{\,\prime}$ may also be considered fixed. The integrations of the full equations of motion in Section 5 confirm the resistance of $d\theta^{\star}/dt$ and $i^{\,\prime}$ to change.

3.5 Asteroidal tides

The expressions for the orbital elements’ rates from Boué & Efroimsky (2019) apply to any two near-spherical bodies, asteroids included. As mentioned in Section 1, asteroids play a key role in white dwarf planetary systems. The asteroid (or asteroids) discovered disintegrating around WD 1145+017 (Vanderburg et al., 2015) has now been the subject of over 20 dedicated refereed papers (Alonso et al., 2016; Gänsicke et al., 2016; Rappaport et al., 2016; Xu et al., 2016; Zhou et al., 2016; Croll et al., 2017; Farihi et al., 2017; Gary et al., 2017; Gurri et al., 2017; Hallakoun et al., 2017; Kjurkchieva et al., 2017; Veras et al., 2017; Cauley et al., 2018; Farihi et al., 2018; Izquierdo et al., 2018; Rappaport et al., 2018; Xu et al., 2018; Duvvuri et al., 2019; Karjalainen et al., 2019). Further, the asteroid embedded well-within the Roche radius (for rocky compositions) of the white dwarf SDSS J1228+1040 (Manser et al., 2019) resides in one of the most dynamically active white dwarf planetary systems.

Nevertheless, an outstanding question remains about the origin and dynamical pathways of these asteroids. If an asteroid were originally spherical, could it have been tidally torqued into the white dwarf Roche radius? The equations here cast doubt on this scenario. The effect of tides become weaker as the planetary radius decreases: for a constant density, the planetary quality function $\,\sim 1/\mathcal{B}_{l}\sim g\rho R\sim M/R\,$, and its coefficients are $\,\sim R^{2l+1}/M\,$, giving

[TABLE]

[TABLE]

[TABLE]

where $l\geq 2$. The steep dependence of these rates on the planetary radius illustrates that the tides in near-spherical asteroids are negligible and could not push the asteroids into the white dwarf Roche radius. Other mechanisms must be invoked (see Section 1.1).

4 Initial conditions and parameters for simulations

An important feature of this paper is not only the establishment of the tidal equations, but also the determination what physical situations to integrate. Our simulations require 16 initial conditions and parameters to be established:

• •

4 for the orbit: $a(0),e(0),i(0),i^{\prime}(0)$

• •

6 for the planet: $\theta(0)$, $\frac{\textstyle d\theta}{\textstyle dt}(0)$, $M$, $R$, $\mathcal{J},\eta$

• •

6 for the star: $\theta^{\star}(0)$, $\frac{\textstyle d\theta^{\star}}{\textstyle dt\,}(0)$, $M_{\star}$, $R_{\star}$, $\mathcal{J}^{\star},\eta^{\star}$

The masses, radii, shear viscosities, and shear compliances 555 We remind the reader that the shear compliance is the inverse of the shear elastic modulus. (equalling eight variables) are all assumed to remain fixed in time throughout each individual simulation. Hence we ignore the potential change in planetary viscosity as the planet approaches the white dwarf, a phenomenon which has been speculated to occur in the TRAPPIST-1 system (Makarov et al., 2018). We also ignore the plausible possibility that a planet’s mass and radius will change due to, for example, sublimation. Such sublimation, however, has been shown to have a negligible effect on the orbital pericentre (Veras et al., 2015c) but could have a larger effect through the reduction of planetary mass. Stellar evolution is irrelevant for white dwarfs on the maximum timescale of our simulations (100 Myr).

Both the star and planet are assumed to be represented fully by the 12 physical parameters above. For example, a single value of $\eta$ must be applied for the entire planet (as opposed to, for example, separate values for its mantle and core).

We cannot, in this paper, cover the entire phase space encompassed by these parameters. Therefore, we must carefully choose which ensembles of variables to vary across our simulations, and do so with computational limitations in mind (see Appendix C).

4.1 Fixed initial conditions and parameters across simulations

In total, we choose to fix 10 variables ($\mathcal{J}$, $\mathcal{J}^{\star}$, $i(0),i^{\prime}(0)$, $\theta(0),\theta^{\star}(0),(d\theta^{\star}/dt)(0)$, $M_{\star},R_{\star},\eta^{\star}$) of the 16, for the following reasons:

• •

Compliances. As explained in Footnote 4, a star may possess some effective magnetic rigidity. We, however, assume that the inverse of rigidity – the compliance – is small ($\mathcal{J}^{\star}\ll 1/(\eta^{\star}\chi^{\star})$), such that the elastic reaction of the stellar material to the tidal stress is much less important than its viscous reaction (see equation 7). Therefore, we set $\mathcal{J}^{\star}=0$ (1/Pa).

The real compliance $\mathcal{J}$ of planets, asteroids, and comets is a better-known quantity: for the solid Earth, $1/{\mathcal{J}}~{}\approx~{}0.8~{}\times~{}10^{11}$ Pa; for ices, $1/{\mathcal{J}}~{}\approx~{}4~{}\times~{}10^{9}$ Pa; and for snow, $1/{\mathcal{J}}~{}\approx~{}10^{6}$ Pa. The corresponding values for ${\mathcal{B}}_{l}$ are such that, for small bodies and for terrestrial planets smaller than Earth, ${\mathcal{J}}~{}\ll~{}{\mathcal{B}}_{l}$, so that ${\mathcal{J}}$ plays no role in equation (5) (provided that the rheology is Maxwell). For overheated super-Earths (which are expected to obey the Maxwell model, see Makarov 2015), we still can neglect ${\mathcal{J}}$, because both ${\mathcal{J}}$ and ${\mathcal{B}}_{l}$ are much smaller than $1/(\eta\chi_{lmpq})$ in the denominator of the right-hand side of equation (5). However, in the case of colder super-Earths, the term ${\mathcal{J}}$ in the denominator of that expression cannot be neglected.

• •

Initial spin orientations. Our numerical integration results are insensitive to $\theta(0)$ and $\theta^{\star}(0)$, at least for the outputs of interest. Therefore, we set $\theta(0)=\theta^{\star}(0)=0^{\circ}$.

• •

Initial spin rate of white dwarf. Because equation (25) is a second-order differential equation, we must set the initial spin rate ($d\theta^{\star}/dt$)(0) of the white dwarf in addition to initial orientation $\theta^{\star}(0)$. As argued in Section 3.4, in most cases the white dwarf’s spin rate will change negligibly due to tidal interactions with a planet, and will negligibly affect other values of interest through the white dwarf quality function. Nevertheless, we still must set a value.

White dwarf spin rates range from $0-1$ revolution/hour, where the upper bound corresponds to the shortest rotational periods found from the observed population of single pulsating white dwarfs (Fig. 4 of Hermes et al., 2017)666These values give an equatorial rotational speed range of about 0-15 km/s.. However, a more typical upper bound is likely to be 1 revolution/day (corresponding to an equatorial rotational speed of about 0.6 km/s), as 1 revolution/hour represents extreme cases where the white dwarf’s spin was likely kicked due to a former binary companion. Hence, we adopt $(d\theta^{\star}/dt)(0)=$ 1 revolution per $20$ hours.

• •

Initial orbital inclinations. As argued in Section 3.4, the evolution of $i^{\prime}$ negligibly affects $a$ and $e$ and does not affect $i$ nor $d\theta/dt$ at all. Hence, here we arbitrarily set $i^{\prime}(0)={10^{-4}}^{\circ}$.

In our preliminary simulations, we found that the orbital inclination $i$ usually either (i) decreases to the coplanar limit sharply relative to the semimajor axis and eccentricity evolution, or (ii) increases slightly. This distinction depends on the relative strength of the two terms in equation (23), and was consistent across all of our preliminary simulations. Further, the magnitude of the inclination, even if above $90^{\circ}$, does not affect this bimodal outcome, and does not substantially affect the semimajor axis nor eccentricity evolution. Therefore, we do not learn enough by varying $i(0)$ to justify doing so. Hence, purely for demonstration purposes, we set $i(0)=140^{\circ}$.777Although this retrograde value might seem high, planets have been shown to reach near the white dwarf Roche radius on arbitrarily highly inclined orbits in both single or binary systems (Veras & Gänsicke, 2015; Hamers & Portegies Zwart, 2016; Petrovich & Muñoz, 2017; Stephan, Naoz & Zuckerman, 2017; Veras et al., 2018a).

• •

White dwarf mass. Conveniently, the mass distribution of the population of white dwarfs is highly peaked around $0.60-0.65M_{\odot}$ (Tremblay et al., 2016), and here we adopt $M_{\star}=0.60M_{\odot}$.

• •

White dwarf radius White dwarf radii are closely linked with their masses (see e.g. equations 27-28 of Nauenberg 1972, equation 15 of Verbunt & Rappaport 1988, and Boshkayev & Quevedo 2018). Consequently, given our choice of $M_{\star}=0.60M_{\odot}$, we choose $R_{\star}=0.01280R_{\odot}=8900$ km.

• •

White dwarf viscosity The star’s dynamic shear viscosity is a parameter for which observational constraints have proven largely elusive. We employ a simple approximation and treat the entire star as a viscous sphere (with a constant viscosity). Even this simple approximation, however, cannot hide the gaping uncertainty in our knowledge of white dwarfs’ viscosity. Dall’Osso & Rossi (2014) suggest that $\eta^{\star}\sim 10^{7}-10^{17}$ Pa$\cdot$s, depending on whether the main driver of the viscosity is radiation, plasma, magnetism, turbulence or a combination thereof. Despite this large range, we fix a specific value through the following procedure.

A method to reduce this uncertainty is to treat the shear viscosity as the product of the star’s average density ($\rho_{\star}$), the typical scale of convective flows ($d_{\star}$, which is comparable to the typical size of supergranules), and the typical speed in the supergranules ($v_{\star}$). An extension of equation 7.32 of Stix (2002) gives

[TABLE]

For the Sun, $v_{\odot}\approx 300$ m/s, $d_{\odot}=1.6\times 10^{7}$ m and $\rho_{\odot}~{}=~{}200$ kg/m3, yielding $\eta_{\odot}=1.0\times 10^{12}$ Pa$\cdot$s. 888 Although Stix (2002) states on page 267 that a typical cell diameter is $\,1.6\times 10^{4}$ km, we assume that he meant radius because his ensuing estimate for the mean spacing between cell centres is $\,3\times 10^{4}$ km.

We perform a similar computation for a variety of potential hydrogen-atmosphere based white dwarf hosts (known as “DA white dwarfs”) by using the stellar models from Koester (2009) and Koester (2010). DA white dwarfs are the most common types of white dwarfs. These models output convective speeds ($v_{\star}$) and the density in the envelope, which is more appropriate to use than the overall white dwarf density in order to model the convection zone. The models also output pressure, which can be combined with a given surface gravity to obtain the pressure scale height (Eq. 1 of Tremblay et al. 2013), which in turn can be used to represent $d_{\star}$.

In Table 1, we list dynamic viscosities for DA white dwarfs at a location in the convection zone where the total white dwarf mass is a factor of 10 orders of magnitude more than the convective zone mass above that layer. The effective temperature range in the table from 5000 K to 9000 K corresponds to cooling ages of about 1 Gyr – 5 Gyr (“cooling age” is the time since the white dwarf was born). The table illustrates overall that the dynamic viscosity of these white dwarfs are near the lower end of the range proposed by Dall’Osso & Rossi (2014), but are nevertheless relatively well-confined. Hence, we fix $\eta^{\star}=10^{7}$ Pa$\cdot$s.

4.2 Variable initial conditions and parameters across simulations

We vary only 6 variables ($a(0),e(0),M,R,(d\theta/dt)(0),\eta$) amongst the different simulations.

• •

Initial semimajor axis and eccentricity. Planets are assumed to approach the white dwarf in the first instance on a highly eccentric ($>0.99$) orbit. Although the Darwin-Kaula expansion can model any eccentricity under unity, practically the number of terms which would need to be retained to compute an accurate solution for $e=0.99$ renders such simulations computationally infeasible. Details of the initial damping period hence remain unknown for now, and would be enlightened in a future study that could perhaps reformulate the Darwin-Kaula formalism by building an expansion over $e$ not about the value $e=0$, but about $e=1$.

Therefore, here we consider the planet at a later stage, after its eccentricity and semimajor axis have already been moderately damped. Of interest to us is the rate of the subsequent inward or outward drift and whether or not the planet will reach the white dwarf’s Roche radius $r_{\rm Roche}$. If we assume the planet is solid and spinning, then we could use the relevant Roche radius coefficient from Veras et al. (2017), and with $M_{\star}=0.60M_{\odot}$, we then obtain

[TABLE]

By using the Roche radius as a scaling, we adopt $a(0)$ values ranging from $2r_{\rm Roche}$ to $30r_{\rm Roche}$, with $10^{-4}~{}\leq~{}e(0)~{}\leq~{}0.4$. As explained in Section 3.2, the eccentricity is never realistically exactly zero.

• •

Planet mass. One advantage of utilizing the Darwin-Kaula expansion of Boué & Efroimsky (2019) is that the primary and secondary could be any objects, including ones of asteroid-size (which is particularly relevant for white dwarf planetary systems). Regardless, as explained in Section 3.5, orbital evolution of spherical asteroids is unaffected by tidal torques. Hence, we explore a selection of higher rocky masses. We vary the planet mass ($M$) from $10^{1}M_{\oplus}$ (representative of a “Super-Earth”) to $10^{-3}M_{\oplus}$ (about the mass of Haumea or Titania). The trends in the results become obvious enough that reducing the lower limit further is not worth spending the computational resources.

• •

Planet radius. We then obtain $R$ by assuming an appropriate value of $\rho$. Here, we simply adopt Earth’s density for all values of $M$.

• •

Initial spin rate of planet. As evidenced by the major planets in our solar system, planetary spin periods can range from about 10 hours to $-$243 days (with the minus indicating retrograde motion); extrasolar planets may have a larger range, but their spins are currently not as well constrained. The direction of planetary spin may alter the results, and hence we take that case into account. Here, we adopt a slightly smaller range to what is seen in the solar system, with $(d\theta/dt)(0)=$ 1 revolution per $10-10^{3}$ hours, in both the prograde and retrograde senses with respect to the orbital motion around the star.

• •

Viscosity of planet. The final parameter to specify, $\eta$, can vary significantly depending on the planet type. One recent estimate of stagnant lid planets (Thiriet et al., 2019) sets $\eta~{}=~{}10^{20}~{}-~{}10^{21}$ Pa$\cdot$s for dry mantle rheologies. On Earth, for the Tibetan crust and lithosphere, $\eta~{}=~{}10^{18}~{}-~{}10^{20}$ Pa$\cdot$s (Henriquet et al., 2019), whereas for the Earth’s mantle, $\eta~{}=~{}10^{20}~{}-~{}10^{24}$ Pa$\cdot$s (Mitrovica & Forte, 2004).

Several recent papers on the icy satellites in the solar system specify a wide range of dynamic viscosities. Table 2 of Hurford et al. (2018) lists $\eta=10^{21},10^{21},10^{20}$ and $10^{14}$ Pa$\cdot$s respectively for Europa’s iron core, brittle ice layer, silicate mantle, and ductile ice layer. Cameron et al. (2019) give similar values for Ganymede, whereas Table 1 of Patthoff et al. (2019) lists $\eta=10^{28},10^{22},10^{14}$ Pa$\cdot$s respectively for the core, upper ice layer and lower ice layer of Enceladus. The latter two values actually were chosen from wide ranges previously reported in the literature of $10^{19}-10^{26}$ Pa$\cdot$s for the upper ice layer and $10^{12}-10^{17}$ Pa$\cdot$s for the lower ice layer.

Overall, the values given in the last paragraph indicate that adopting a range of $\eta=10^{16}-10^{24}$ Pa$\cdot$s is representative of most cases, and we do so.

5 Simulation results

The simulation output of primary interest is how the planet drifts relative to the Roche radius, and whether or not the planet reaches the Roche radius. Therefore, we report the semimajor axis drift in terms of the parameter

[TABLE]

The sign of $\sigma_{a}$ indicates if the planet has drifted outward (positive) or inward (negative); $\sigma_{a}=-100\%$ indicates entering the Roche radius (with presumed destruction subsequently). Further, we define

[TABLE]

and

[TABLE]

Both $\sigma_{e}$ and $\sigma_{i}$ may be positive or negative.

Before reporting on our ensemble results, we first illustrate a couple examples of time evolution in Figure 1. The figure helps exhibit that the orbital evolution is not necessarily monotonic with time, and how the final outcomes reported can be strongly dependent on the stopping time of the simulation. Therefore, characterising the entire phase space within the space of one paper is challenging. The planet can change direction due to tides, with the semimajor axis and eccentricity evolving in a non-obvious manner. The figure also illustrates how qualitatively different behaviour can result just by decreasing the planet’s viscosity by an order of magnitude.

That figure alone flags the danger of attempting to characterise tidal effects in an individual white dwarf planetary system by appealing to simplified comparisons (Section 3.1) or by identifying a particular point on a phase portrait. Nevertheless, here we do now construct phase portraits. These are intended only for order-of-magnitude use, and to detect general trends which remain robust amidst the complexities.

Of primary interest to us is the final value of $\sigma_{a}$. We integrated the equations of motion for 100 Myr and recorded the values of $\sigma_{a}$ after this time has elapsed. For these simulations, we adopted the fixed initial conditions and parameters reported in Section 4.1 ($\mathcal{J}^{\star}=0$ (1/Pa), $i(0)=140^{\circ}$, $i^{\prime}(0)={10^{-4}}^{\circ}$, $\theta(0)=\theta^{\star}(0)=0^{\circ}$, $(d\theta^{\star}/dt)(0)=360^{\circ}/(20\ {\rm hours})$, $M_{\star}=0.60M_{\odot}$, $R_{\star}=8900$ km, $\eta^{\star}=10^{7}$ Pa$\cdot$s) and a fiducial set of initial conditions and parameters in Section 4.2 from which we vary selected parameters. This fiducial set is ($e(0)=0.2$, $M=M_{\oplus}$, $R=R_{\oplus}$, $(d\theta/dt)(0)=360^{\circ}/(20\ {\rm hours})$, and $\eta=10^{20}$ Pa$\cdot$s).

We report our results in a series of pictographs (Figs. 2-5). In every figure, we sampled 12 initial values of $a(0)$ where relevant behaviour occurs. Red crosses indicate that $a\leq r_{\rm Roche}$ at some point during the simulation, suggesting complete destruction of the planet999We do not report when the orbital pericentre $a(1-e)$ first intersects the Roche radius, because (i) the resulting destruction could be intermittent, and (ii) the eccentricity is small enough within $2r_{\rm Roche}$ such that uncertainty in the Roche radius coefficient (Veras et al., 2017) would dominate the difference between orbital pericentre and semimajor axis.. The other figure entries feature a red “V” or green caret; the former indicates the semimajor axis has decreased ($\sigma_{a}<0$) after 100 Myr, and the latter indicates that the semimajor axis has increased ($\sigma_{a}>0$) after 100 Myr. The number given inside the symbols is the order of magnitude of the change: $-1$ refers to $10\%\leq\left|\sigma_{a}\right|<100\%$; $-3$ refers to $0.1\%\leq\left|\sigma_{a}\right|<1\%$, and so forth. In one case, a [math] indicates $100\%\leq\sigma_{a}<1000\%$

The figures are ordered according to the variables which most strongly dictate the final outcome. The first figure, Fig. 2, illustrates the importance of knowing, estimating or being able to guess the planet’s dynamic viscosity. Recall that all viscosity values in the figure are within the ranges suggested for various components of solar system bodies. The critical destruction distance can vary significantly, and planets which just avoid destruction (bordering the red crosses) may be shuffled about significantly. The plot demonstrates a clear trend against destruction as the planet viscosity increases. For the highest viscosities, even planets initially within $2r_{\rm Roche}$ can survive.

The next figure (Fig. 3) reinforces the arguments presented in Section 3.5: as planetary mass (and radius) decreases, tidal effects quickly become negligible. Similarly, Super-Earths ($10M_{\oplus}$) are particularly susceptible to white dwarf tides, with a (100 Myr) destruction radius reaching out to $4r_{\rm Roche}$ for $\eta~{}=~{}10^{20}$ Pa$\cdot$s and $10r_{\rm Roche}$ for $\eta~{}=~{}10^{16}$ Pa$\cdot$s. Because asteroids are many orders of magnitude less massive than $10^{-3}M_{\oplus}$, extrapolating from the figure indicates that spherical asteroids would be virtually unaffected by white dwarf tidal torques, even when assuming our lower bound for dynamical viscosities.

Fig. 4 contains one grid for prograde planetary spins (top) and one for retrograde planetary spins (bottom): differences are highlighted by gray squares. The grids are 79 per cent coincident with each other, indicating that on an order-of-magnitude scale, the direction of spin makes little difference to the final outcome. The likely reason is that the planet’s rotation becomes quickly synchronized or pseudosynchronized. Most of the differences occur at the highest $a(0)$ values sampled and hence are small, but there are a few important exceptions close to the critical destruction distance.

The general dependence of destruction and drift with the magnitude of the spin period is not obvious, at least from the figure. Nevertheless, the top rows in each grid are well-ordered and perhaps suggests predictable asymptotic behaviour as the spin period tends towards infinity.

The final figure, Fig. 5, illustrates how the critical engulfment distance is a weak function of $e(0)$, at least for $e(0)<0.4$. For the lowest values of $e(0)$, the results are non-obvious: the boundary between being engulfed and remaining nearly stationary is sharp, and not adequately sampled with our choices of $a(0)$. These results emphasise the importance of performing dedicated studies for individual systems due to the chaotic nature of the phase space.

6 Summary

Exo-asteroids are already observed orbiting two white dwarfs in real time, and almost every known exo-planet currently orbits a star that will become a white dwarf. Planets which then survive to the white dwarf phase play a crucial role in frequently shepherding asteroids and their observable detritus onto white dwarf atmospheres, even if the planets themselves lie just outside of the narrow range of detectability. Further, planets themselves may occasionally shower a white dwarf with metal constituents through post-impact crater ejecta and when the planetary orbit grazes the star’s Roche radius (Veras & Gänsicke, 2015; Brown et al., 2017).

Here, we undertook one of the first dedicated studies of tides in a two-body system comprising a solid planet and a white dwarf. We adopted the Fourier-mode tidal formalism of Boué & Efroimsky (2019), which provides complete equations (21–24) for the secular evolution of the orbital semimajor axis, eccentricity and inclination with frequency dependencies on the quality functions. By combining these equations with secular parts of the tidal torque from Efroimsky (2012a) and Makarov (2012), and typical physical parameters for the two bodies (Section 4.1), we generated a computational framework for future, more detailed consideration of individual systems.

A broad sweep of phase space has revealed the following trends:

(i) Massive Super-Earths are more easily destroyed than low-mass planets.

(ii) Planetary survival is boosted with higher viscosities.

(iii) Orbital evolution is in general non-monotonic, and cannot usually be described by a straightforward comparison of spin period versus orbital period.

(iv) The boundary between destruction and survival appears to be chaotic.

(v) Because the magnitude of the stellar tides scales as the mass of the perturber, the orbital dynamics of the asteroids in the WD 1145+017 and SDSS J1228+1040 systems are unaffected by stellar tides.

(vi) The relatively small range of white dwarf physical parameters, as compared to those of main-sequence stars, helps constrain variable explorations and ensure that semimajor axis, planetary mass, planetary viscosity and planetary spin rates are the four key variables which influence the evolution.

(vii) Conservatively, even for low values of the planetary viscosity, planets which do not achieve a semimajor axis of under about $10r_{\rm Roche}$ will survive for at least 100 Myr of white-dwarf cooling.

(viii) The planet’s evolution is largely independent of the direction but not magnitude of the spin of the planet.

(ix) Despite the varied ways in which the orbit can be stretched (Figure 1), the critical engulfment distance is largely independent of the orbital eccentricity when it takes on low and moderate values.

Acknowledgements

All authors appreciate the helpful comments from the MNRAS referee and the four internal United States Naval Observatory (USNO) referees Paul Barrett, Christopher Dieck, George Kaplan and Robert Zavala. DV further thanks the USNO for its repeated hospitality for over a decade, and the Landessternwarte, Zentrum für Astronomie der Universität Heidelberg for a productive visit, ultimately kickstarting this work. DV also gratefully acknowledges the support of the STFC via an Ernest Rutherford Fellowship (grant ST/P003850/1). ME would like to thank Alexander Getling for a valuable consultation on the stellar viscosity. VW and SR further acknowledge support by the DFG Priority Program SPP 1992 Exploring the Diversity of Extrasolar Planets (RE 2694/5-1). PET has received funding from the European Research Council under the European Union’s Horizon 2020 research and innovation programme n. 677706 (WD3D).

Appendix A Parameter Tables

Appendix B Eccentricity and inclination functions

The eccentricity and inclination functions are utilised in the equations of motion, and may be precomputed regardless of the physical setup. Doing so helps speed up the integrations.

B.1 Eccentricity function

We utilise the useful form of the eccentricity function $G_{lpq}(e)$ which is given in the Appendix of Celletti et al. (2017). Their formulae are based on those from Giacaglia (1976). Transforming their notation into our formalism yields the following results for the eccentricity function

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

Note that because $0\leq p\leq l$, always here $s_{1}=2p-2l$ and $u_{1}=\infty$. This last relation creates an additional parameter in numerical computations which must be explored in order to achieve the desired accuracy.

B.2 Inclination function

We use the inclination function given in equation 4 of Gooding & Wagner (2008), which corrects some sign errors from the exposition in Kaula (1962) and conveniently does not include integer parts. Translating the Gooding & Wagner (2008) variables into ours yields

[TABLE]

Appendix C Equations of motion

The following equations of motion (equations 21-28) include the integers $\,l,\,m,\,p$, and $q$ explicitly through summation limits; note that $n$ here refers to mean motion (specifically, the anomalistic mean motion) and is not an index.

Known in tidal investigations as the “degree”, the positive integer $\,l\,$ in particular sets the scene: (a) it represents the index for the first summation which appears in these equations of motion, and (b) it establishes the maximum possible values of $m$ and $p$.

For point (a), the truncation of the tidal expansions by setting $l=2$ (the quadrupole approximation) represents a very common approximation for the evolution of the orbital elements and the spin. The suitability of the quadrupole approximation is discussed in Boué & Efroimsky (2019). Although for highly eccentric orbits around white dwarfs this approximation might not be adequate, it remains sufficient for the eccentricities considered here ($<$0.40).

For point (b), because $l=2$, then $p=0,1,2$ and $m=0,1,2$. The only index which then remains a free parameter is $q$. Consequently, the range of $q$ must be chosen to encompass all of the eccentricity terms of a given order. Although a typical approximation is to assume that $\left|q\right|$ is equal to the highest power of the eccentricity sampled, actually the distribution of $q$ which must be considered to achieve a given accuracy is strongly asymmetric (Noyelles et al. 2014 use $-1\leq q\leq 7$). We did not employ the typical approximation for $q$, nor consider the concept of order, which can introduce pitfalls in series expansions (Veras, 2007). We created fits to $G_{lpq}(e)$ for $e<0.4$ by choosing, quite conservatively, $-28\leq q\leq 43$. However, preliminary numerical testing indicated that we need only choose $-7\leq q\leq 7$ for our primary integrations; this range of values allowed us computationally to sample several regions of phase space.

All of the equations in this section include tidal contributions from both the star and planet, and are not truncated in any way.

C.1 Equations of orbital motion

The secular semimajor axis evolution of the mutual orbit is given by Boué & Efroimsky (2019) as

[TABLE]

[TABLE]

where the double averaging indicates averaging over both one orbital cycle and one cycle of the apsidal precession.

The secular eccentricity evolution of the mutual orbit is given by Boué & Efroimsky (2019) as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The secular inclination evolution of the mutual orbit with respect to the planetary equator is given by Boué & Efroimsky (2019) as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

With respect to the stellar equator, we instead have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

C.2 Equations of spin motion

The secular spin evolutions of the planet and star are given by Eq. 1 of Makarov (2012) as

[TABLE]

[TABLE]

and the secular part of the polar component of the tidal torque is given by equation (109) of Efroimsky (2012a) as

[TABLE]

[TABLE]

A detailed derivation of the polar component of the permanent-triaxiality-generated torque is provided in Frouard & Efroimsky (2017, Appendix D). Equation (116e) of that paper reveals that  under steady rotation  and, importantly, outside of a spin-orbit resonance , the secular part of this torque vanishes. In realistic situations, the equation further reveals that when the apsidal precession is much slower than the orbital motion, then orbital averaging is sufficient to nullify the torque. For these reasons, we equip these secular parts with only one pair of angular brackets:

[TABLE]

[TABLE]

Both these equalities no longer imply that a rotator is captured in a spin-orbit resonance. In such a resonance, the triaxiality-generated torque is leading for non-liquid bodies (much larger than the tidal torque) and becomes the driver of longitudinal libration.