Survivability of radio-loud planetary cores orbiting white dwarfs

TL;DR

This study models the survivability of metallic planetary cores orbiting white dwarfs by coupling gravitational tides and Lorentz drift effects, identifying conditions under which cores can survive over a billion years for potential radio detection.

Contribution

It introduces a coupled model of gravitational and electromagnetic forces affecting planetary core survivability, expanding previous models by including a Maxwell rheological framework.

Findings

High viscosity cores (>~10^{24} Pa*s) are more likely to survive.

Survivability is maximized in white dwarfs with magnetic fields in the kG range.

Cores can persist over Gyr timescales under certain physical conditions.

Abstract

The discovery of the intact metallic planetary core fragment orbiting the white dwarf SDSS J1228+1040 within one Solar radius highlights the possibility of detecting larger, unfragmented conducting cores around magnetic white dwarfs through radio emission. Previous models of this decades-old idea focussed on determining survivability of the cores based on their inward Lorentz drift towards the star. However, gravitational tides may represent an equal or dominant force. Here, we couple both effects by assuming a Maxwell rheological model and performing simulations over the entire range of observable white dwarf magnetic field strengths (10^3 -- 10^9 G) and their potential atmospheric electrical conductivities (10^{-1} -- 10^4 S/m) in order to more accurately constrain survivability lifetimes. This force coupling allows us to better pinpoint the physical and orbital parameters which allow planetary cores to survive for over a Gyr, maximizing the possibility that they can be detected. The most robust survivors showcase high dynamic viscosities (>~ 10^{24} Pa*s) and orbit within kG-level magnetic fields.

Equations

$${⟨⟨\frac{da}{dt}⟩⟩}_{\mathrm{Lorentz}}=-5⟨{\gamma }_{\star }⟩B_{\star }^2\left(\frac{R^3}{R_{\star }}\right){\left(\frac{R_{\star }}{a}\right)}^{\frac{17}{2}}\left(\frac{a^2}{M}\right)$$

$$\times {\left[1-⟨⟨\frac{d{\theta }^{\star }}{dt}⟩⟩\frac{a^{3/2}}{\sqrt{G\left(M_{\star }+M\right)}}\right]}^2,$$

$$H=-\frac{G{M}_{\star }M}{4r^2}{⟨⟨\frac{da}{dt}⟩⟩}_{\mathrm{Lorentz}}$$

$$a_{\mathrm{crit}}\approx {\left[\frac{B_{\star }{R}_{\star }^{3}R\sqrt{G{M}_{\star }}}{1\mathrm{kV}}\right]}^{2/7}.$$

Full text

††pubyear: 2019††pagerange: Survivability of radio-loud planetary cores orbiting white dwarfs–References

1 Introduction

Over two decades ago, Ferrario et al. (1997) detected striking emission lines in the spectrum of the magnetic white dwarf GD 356. Prompted by this finding, Li et al. (1998) proposed that the electromagnetic interactions between a magnetic white dwarf and a conducting planetary core could create a unipolar inductor circuit similar to Jupiter and Io’s (Goldreich & Lynden-Bell, 1969) that both drags the planet towards the star (a phenomenon known as “Lorentz drift”) and generates detectable emission. Years later, Willes & Wu (2004) and Willes & Wu (2005) followed up and investigated this idea further by providing additional details on the evolution and radio detection of these planets, based on a model from Wu et al. (2002).

Although such planets have not yet been discovered, this possibility is becoming increasingly plausible due to the rapid increase in detections of various features of white dwarf planetary systems. Since 2005, both the available data and our understanding of these systems has matured and grown substantially, as we explain in the next two subsections. Hence, the context in which researchers investigate white dwarf planetary systems has changed.

Further, during this period, investigations of detectable emission from electromagnetic interactions of planets orbiting other types of stars have been numerous; some examples since the year 2015 are Vidotto et al. (2015), Fujii et al. (2016), Katarzyński et al. (2016), Nichols & Milan (2016), Strugarek et al. (2017), Weber et al. (2017), Daley-Yates, & Stevens (2018), Lynch et al. (2018), O’Gorman et al. (2018), Strugarek (2018), and Wang & Loeb (2019). Also, these studies incorporated stellar winds, mass loss and/or planetary atmospheres, none of which exist in the scenario that we are considering.

1.1 Current state of white dwarf planetary systems

Observable features in white dwarf systems result from substellar objects which avoid engulfment after their host star has left the main sequence (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) and then are gravitationally perturbed close to the resulting white dwarf (Veras, 2016a). After entering the white dwarf’s disruption, or Roche radius, an object will break up and form a series of rings or annuli (Graham et al., 1990; Jura, 2003; Debes et al., 2012; Bear & Soker, 2013; Veras et al., 2014c; Brown et al., 2017; Veras et al., 2017b) usually before accreting on to the white dwarf atmosphere with a particular size distribution (Wyatt et al., 2014; Veras et al., 2015b).

These substellar objects which are perturbed towards the white dwarf comprise planets (Debes & Sigurdsson, 2002; Veras et al., 2013; Voyatzis et al., 2013; Mustill et al., 2014; Veras & Gänsicke, 2015; Veras et al., 2016, 2017a; Stephan, Naoz & Gaudi, 2018; Veras et al., 2018), asteroids (Bonsor et al., 2011; Debes et al., 2012; Frewen & Hansen, 2014; Bonsor & Veras, 2015; Hamers & Portegies Zwart, 2016; Petrovich & Muñoz, 2017; Stephan, Naoz & Zuckerman, 2017; Mustill et al., 2018; Smallwood et al., 2018), comets (Alcock et al., 1986; Veras et al., 2014a; Stone et al., 2015; Caiazzo & Heyl, 2017) and moons (Payne et al., 2016, 2017). In our solar system, at least five of the major planets will survive until and during the solar white dwarf phase (Sackmann et al., 1993; Duncan & Lissauer, 1998; Rybicki & Denis, 2001; Schröder & Smith, 2008; Veras, 2016b). The inner solar system will likely to be strewn with asteroid belt debris (Veras et al., 2014b) while the minor planets in the outer solar system will be re-arranged (Veras et al., 2015a, 2019a).

Observations of the consequences of these dynamical interactions are abundant. Accreted planetary debris in white dwarf photospheres are detectable spectroscopically because these stars are dense enough to stratify matter by atomic weight (Schatzman, 1958); any observed metals heavier than helium then arise from planetary debris, with few exceptions (Farihi et al., 2010; Veras et al., 2019c). Between one-quarter and one-half of all Milky Way white dwarfs are now thought to contain planetary debris (Zuckerman et al., 2003, 2010; Koester et al., 2014), a fractional range which is similar to that of the planet occurrence range around main sequence stars (Cassan et al., 2012). The composition of the debris provides exclusive and detailed insight into the bulk chemical composition of exo-planetary bodies, even enabling one to link an individual white dwarf with an individual meteorite family in the solar system (Zuckerman et al., 2007; Klein et al., 2010, 2011; Gänsicke et al., 2012; Jura & Young, 2014; Wilson et al., 2015; Xu et al., 2017; Harrison et al., 2018; Hollands et al., 2018).

These metals accrete onto the white dwarf through the aforementioned discs. Their formation has not only been theorized and modelled, but also observed: in fact about 50 such discs have been detected (Zuckerman & Becklin, 1987; Farihi, 2016). Many more discs are missed simply due to misalignment between our line-of-sight and the plane of the disc (Bergfors et al., 2014; Bonsor et al., 2017). The discs range in radius from about $0.6R_{\odot}$ to $1.2R_{\odot}$ and showcase variability (Gänsicke et al., 2008; Wilson et al., 2014; Xu & Jura, 2014; Farihi et al., 2018a; Xu et al., 2018; Swan et al., 2019) and secular changes in spectroscopic features (Manser et al., 2016; Cauley et al., 2018; Dennihy et al., 2018), both of which are indicative of dynamical activity.

All of the discs contain dust and about eight of these discs contain detectable gas (Gänsicke et al., 2006, 2007, 2008; Gänsicke, 2011; Farihi et al., 2012; Melis et al., 2012; Wilson et al., 2014; Guo et al., 2015), which is generated from a combination of sublimation and collisions (Jura, 2008; Metzger et al., 2012; Kenyon & Bromley, 2017). The sublimation radius decreases with white dwarf cooling age, and at some point will cross the disruption radius.

In addition to metal pollution and debris discs, minor planets have been observed orbiting white dwarfs. The disintegrating asteroid orbiting WD 1145+017 (Vanderburg et al., 2015) resides a distance of just one Solar radius and has spawned over 20 additional papers. The intact and metallic core object orbiting SDSS J1228+1040 instead resides at a distance of just $0.73R_{\odot}$ and hence cannot be a rubble pile but rather an object of significant density and internal strength (Manser et al., 2019). These minor planets close to white dwarfs are most likely first-generation leftovers from the main sequence phase because second-generation formation in white dwarf debris discs requires massive discs (van Lieshout et al., 2018) and/or a stellar companion (Perets, 2011; Schleicher & Dreizler, 2014; Völschow et al., 2014; Hogg et al., 2018).

1.2 Planet – white dwarf interactions

Overall, the robust state of the observations and theories outlined above has prompted researchers to revisit the potential electromagnetic interactions between a magnetic white dwarf and orbiting bodies (Farihi et al., 2017, 2018b; Rappaport et al., 2018; Bromley & Kenyon, 2019), especially now that a ferrous core fragment has been discovered orbiting SDSS J1228+1040 (Manser et al., 2019). Bromley & Kenyon (2019) in particular have considered both the Lorentz drift and Ohmic heating of asteroid-sized bodies (including dwarf planets) orbiting a variety of stars, including white dwarfs. The form of the secular orbital decay rate due to Lorentz drift from their Eq. (18) is similar to those reported in Li et al. (1998), Willes & Wu (2004) and Willes & Wu (2005). Further, Bromley & Kenyon (2019) investigated in depth the thermal evolution of the asteroids due to Ohmic heating. Nevertheless, all four of these studies assume that the Lorentz drift dominates the motion of these planets close to the white dwarf. For cores which are the sizes of terrestrial planets, this assumption is not necessarily justified (Bouvier & Cébron, 2015; Strugarek et al., 2017).

Veras et al. (2019b) demonstrated that gravitational tides between a solid planet and a white dwarf play a pivotal role in determining whether a planet migrates into the Roche radius, and if not, in what direction and to what extent does the migration occur. They adopted an arbitrary frequency dependence on the quality dissipation functions and illustrated that the resulting motion may be non-monotonic, and vary the orbit’s semimajor axis, eccentricity and inclination. The non-monotonicity arises from the chaotic interplay between the nonlinear orbital and spin equations of motion (Boué & Efroimsky, 2019) within a Maxwell rheological model (Efroimsky, 2015).

1.3 Plan for paper

The goal of this paper is to model the interplay between gravitational tides and Lorentz drag in compact white dwarf planetary systems in order to more accurately determine survivability timescales. We are keen to determine the planet’s survivability as a function primarily of the star’s magnetic field strength and atmospheric conductivity, the planet’s dynamic viscosity, and the distance from the core to the white dwarf’s Roche radius. We consider only cores large enough to not be destroyed through melting, with the intention that surviving cores can represent targets for radio emission searches (Hess & Zarka, 2011; Zarka et al., 2018).

In Section 2 we briefly summarize the relevant theoretical considerations. We report on our simulation outputs in Section 3 before discussing the applicability of our results to other types of planetary systems and white dwarfs in Sections 4-5. We summarize our results in Section 6.

2 Drivers of motion

In this section we outline the different drivers of motion for a solid metallic or iron-rich object at least the size of Ceres (which we will henceforth refer to as “planet”) orbiting a white dwarf.

2.1 Gravitational tides

First consider an unmagnetized white dwarf. As the planet approaches the white dwarf from a scattering event at a distance of at least a few au, gravitational tides between the star and planet determine the latter’s subsequent motion. The density, spin and internal strength of the planet determines the location of the star’s Roche radius, within which that planet will be destroyed. Gravitational tides may push the planet into the Roche radius or push the planet outward, depending on how the orbit is modified.

If the rheologies of both the star and planet are treated within a Maxwell model, then the character of the orbital change is dictated by 14 parameters (see Section 4 of Veras et al. 2019b). These include the four orbital parameters $a$, $e$, $i$, and $i^{\prime}$ — which represent the semimajor axis, eccentricity, inclination with respect to the planet’s equator, and inclination with respect to the star’s equator — and the following five physical parameters for both the planet (with no subscript nor superscript) and star (with a $\star$ subscript or superscript): the mass $M$, radius $R$, rotation angle about the instantaneous shortest axis of the body $\theta$ (along with the spin rate $d\theta/dt$), compliance $\mathcal{J}$ and dynamic viscosity $\eta$. The initial values of all 14 parameters are fed into the coupled set of spin and orbital equations (which are not repeated here) C1-C10 of Veras et al. (2019b) in order to determine the motion.

The evolution equations are all secular, meaning that they compute the averaged motion over long timescales. Here, “long” is a quantity, such as years, which is much greater than an orbital timescale of hours or days. For the purposes of assessing detectability, computing evolution over secular timescales is desirable. Further, white dwarfs represent the endpoints of stellar evolution for the majority of stars in the Milky Way, and hence cool quiescently over $10^{10}$ yr timescales. Integrating the equations over Gyr timescales is computationally manageable because secular equations are independent of parameters, such as the mean anomaly, which pinpoint the location of the planet along the orbit.

White dwarf discs are highly unlikely to be massive enough to affect the motion of a planet at the sizes which we are considering, although the planet may perturb the disc in interesting and observationally measurable ways (Manser et al., 2016; Xu et al., 2018; Manser et al., 2019).

2.2 Lorentz drift

Now consider a magnetic white dwarf and the scenario originally envisaged by Li et al. (1998), Willes & Wu (2004), and Willes & Wu (2005). We now know that the gas required to maintain an active magnetosphere can arise from additional sources besides the interstellar medium (as explained in Section 1.1). These studies proposed that the electric currents generated by a planet travelling through a white dwarf magnetosphere will set up a unipolar inductor circuit similar to that which is seen with Jupiter and Io, where Jupiter is akin to the white dwarf and Io is akin to the planet (Goldreich & Lynden-Bell, 1969). Then the magnetic field lines close to the white dwarf reside at the locations of radio emission zones: see Fig. 1 in each of Li et al. (1998), Willes & Wu (2004), and Willes & Wu (2005). Resulting Lorentz torques generate energy, some of which is dissipated and some of which is applied as kinetic energy to the planet.

The magnitude of the resulting inward drift on the planet, the “Lorentz drift”, may be estimated. A direct comparison of the drift rates from Li et al. (1998), Willes & Wu (2004), Willes & Wu (2005) and Bromley & Kenyon (2019) reveal that no two of their rates are equivalent, even when accounting for the different electromagnetic systems of units that were adopted111Li et al. (1998), Willes & Wu (2004), and Willes & Wu (2005) all use the electrostatic ESU system, which introduces a factor of the square of the speed of light which is not present in the SI unit system..

Now we attempt to parse these differences. The results of Bromley & Kenyon (2019) differ because they consider a static dipole field and an efficiency factor based on a functional form of the current density within an asteroid from Bidinosti et al. (2007). The drift rate form of Li et al. (1998) differs from those of Willes & Wu (2004) and Willes & Wu (2005) by a numerical factor ($\pi/2$ versus 8) as well as whether to represent the magnetic field strength by the square of the stellar magnetic moment divided by the sixth power of the white dwarf radius. Willes & Wu (2004) and Willes & Wu (2005) differ because the former includes an asynchronism term which contains a comparison of the white dwarf’s angular spin speed to the planet’s orbital angular speed.

All these considerations led us to adopt the following estimate for the Lorentz drift in SI units:

[TABLE]

[TABLE]

where the magnetic field strength is given by $B_{\star}$ and the vertically-averaged electrical conductivity in the white dwarf atmosphere is given by $\left\langle\gamma_{\star}\right\rangle$. All notations are consistent with those in Veras et al. (2019b).

In equation (1), we have (i) assumed that the drift rate is secular and slow enough to affect only the long-term change in the semimajor axis, (ii) adopted a numerical coefficient which is roughly in between those of Li et al. (1998) and Willes & Wu (2004), and (iii) included the asynchronism term222Despite including the asynchronism term, Willes & Wu (2004) assumed that the stellar spin period is much greater than the planetary orbital period. Now we know that white dwarf spin periods generally range from hours to days (Hermes et al., 2017).. The asynchronism term provides an additional coupling between the orbital and spin periods of motion along with those from Veras et al. (2019b).

In order to simulate the planet’s evolution, we added equation (1) to Eq. C1 of Veras et al. (2019b). However, equation (1) requires us to establish values for two additional and important parameters: $B_{\star}$ and $\left\langle\gamma_{\star}\right\rangle$. We devote the next two subsections to obtaining reasonable ranges for these parameters.

2.2.1 White dwarf magnetic field strengths ($B_{\star}$)

Observations of isolated magnetic white dwarfs show $B_{\star}=10^{3}-10^{9}$ G and suggest $\dot{B}_{\star}\approx 0$ (Ferrario et al., 2015). We explored this entire range of magnetic field strengths here, and further assumed $\dot{B}_{\star}\approx 0$. The incidence of magnetic white dwarfs is different amongst different spectral classes (Giammichele et al., 2012; Dufour et al., 2013; Kawka & Vennes, 2014; Hollands et al., 2015, 2017; Kawka et al., 2019), engulfed planets might even trigger a magnetic field in a white dwarf (Farihi et al., 2011; Briggs et al., 2018), and a magnetic field of at least tens of kG may explain the photometric transits and accretion rates seen in WD 1145+017 (Farihi et al., 2017, 2018b). The highest observed magnetic fields ($\sim 10^{9}$ G) may result from stellar mergers (Briggs et al., 2015, 2018). We simply treated the value of $B_{\star}$ to be independent of other stellar parameters and the dynamical history of the system.

2.2.2 Conductivities of white dwarf atmospheres ($\left\langle\gamma_{\star}\right\rangle$)

Both the core of a white dwarf and the metallic core of a planet are conductive. Hence, any closed current loop connecting the magnetic white dwarf to the planet’s core must generate heat in the white dwarf atmosphere. The conductivity value of the atmosphere therefore not only affects the Lorentz drift, but also energy deposition, which could lead to potentially detectable phenomena such as emission lines.

Although most white dwarfs are of spectral type DA (hydrogen line-dominated), investigations of their atmospheric conductivity values are not widespread. Therefore, we use a sufficiently wide range of conductivities to encompass plausible values. Table 2 of Srećković et al. (2010) provides static electric conductivity values of DB (helium line-dominated) white dwarf atmospheres for a range of equilibrium temperatures between 8000 K – 75,000 K and plasma mass densities between $5\times 10^{-7}$ g/cm and $5\times 10^{-3}$ g/cm. The conductivities for all these values range over six orders of magnitude, from $\left\langle\gamma_{\ast}\right\rangle=6.88\times 10^{-2}$–$9.84\times 10^{4}$ Siemens/metre (or, equivalently, Mho/metre or 1/Ohm$\times$metre). However, given that the vast majority of white dwarfs are observed with effective temperatures under about 20,000 K, for our purposes we decreased this upper bound by an order of magnitude.

Even so, the entire range of conductivities is still substantial and should be explored. Generally these conductivities decrease with increasing white dwarf cooling age; Fig. 4 of Mazevet et al. (2007) indicates that this trend continues for (older) white dwarfs with effective temperatures under 8000 K. As previously mentioned, we do not know how the conductivities of DA (hydrogen-dominated) white dwarf atmospheres compare with the results of Srećković et al. (2010). All these considerations led us to adopt the range $\left\langle\gamma_{\ast}\right\rangle=10^{-1}$ – $10^{4}$ S/m.

By using this range, we can determine how much energy is deposited in a white dwarf atmosphere. The heating rate ($H$) near one of the magnetic poles is given by Eq. 7 of Li et al. (1998) as

[TABLE]

where $r$ is the distance between the centres of the white dwarf and planetary core.

Hence, by ignoring the asynchronism term from equation (1) and considering only snapshots in time, we can provide a representative range of heating rates with just a minimal set of reasonable physical parameters ($M_{\star}=0.6M_{\star}$, $R_{\star}=8900$ km, $\rho=8$ g/cm3, $r_{\rm Roche}=0.68R_{\odot}$; this choice of Roche radius is justified in Section 3.1). Figure 1 displays the energy deposited for planetary cores of mass $M=0.5M_{\oplus}$ (solid lines) and $M=5.0M_{\oplus}$ (dashed lines). During planetary evolution, the heating rate will evolve with time, perhaps non-monotonically. Here, we simply plot power as a function of separation. The lower (blue) and upper (green) curves may be treated as bounds on the heat generated.

3 Simulations

Now we describe the setup and results of our simulations. The primary challenge was selecting the most revealing regions of our large parameter space to explore. The numerical implementation involved coupling equation (1) to Eqs. C1-C10 of Veras et al. (2019b). We do not repeat the details of the implementation here, except to say that we used the same generous range of the index $q$ in order to compute the eccentricity functions in their Eq. (B1). We assumed that the rheologies of both the white dwarf and planet adhere to the Maxwell model.

3.1 Parameter choices

The Lorentz drift is most strongly dependent on $a$, $B_{\star}$ and $\left\langle\gamma_{\ast}\right\rangle$, the latter two parameters due to their wide ranges of plausible values. Gravitational tides between planets and white dwarfs are most strongly dependent on $a$, $M$ and $\eta$. Although all of the other parameters affect the outcomes, here we are interested only in order-of-magnitude results. Therefore, we fixed $e(0)=0.2$, $i(0)=i^{\prime}(0)=1^{\circ}$, $M_{\star}=M_{\star}(0)=0.6M_{\star}$, $R_{\star}=R_{\star}(0)=8900$ km, $\theta(0)=\theta^{\star}(0)=0^{\circ}$, $d\theta(0)/dt=d\theta^{\star}/dt(0)=360^{\circ}/(20\ {\rm h})$, $\mathcal{J}=\mathcal{J}^{\star}=0$ (1/Pa), and $\eta^{\star}=10^{7}$ Pa$\cdot$s. These values are all motivated from Veras et al. (2019b).

Regarding the physical properties of the planet, it must be a conductor in order to be detectable in radio emission. In this sense, we can better think of the planet as a planetary core. Such a core might have once hosted a gas giant, or represented most of the volume of a Super-Earth or Earth-like planet. As mentioned in Section 1.1, the scattering of large bodies like these towards a white dwarf can occur with many different potential orbital architectures (Debes & Sigurdsson, 2002; Veras et al., 2013; Voyatzis et al., 2013; Mustill et al., 2014; Veras & Gänsicke, 2015; Veras et al., 2016, 2017a; Stephan, Naoz & Gaudi, 2018; Veras et al., 2018). Because planetary cores are thought to be iron rich, Li et al. (1998) assumed that $\rho=9.4-13$ g/cm3 and Willes & Wu (2005) assumed $\rho=13.3$ g/cm3 in their Fig. 2.

However, given that the density of iron is just under $8$ g/cm3, and some cores of the solar system bodies are thought to be mixtures of iron and silicates, we adopted $\rho=8$ g/cm3. Both Li et al. (1998) and Willes & Wu (2005), papers that were written before the revelation that super-Earths are common throughout the Galaxy333See the NASA Exoplanet Archive at exoplanetarchive.ipac.caltech.edu, assumed the planet mass is 40 per cent of the mass of the Earth. Here we adopted two values for $M$: $M=\left\{0.5M_{\oplus},5.0M_{\oplus}\right\}$. Consequently, the radii of our two planet cores are 4500 km and 9600 km (the latter actually being larger than the white dwarf radius).

Unlike for the much smaller core fragment discovered orbiting a white dwarf in Manser et al. (2019), here we did not assume that the object has any internal strength. Our planets’ density, along with the mass of the white dwarf, sets the Roche radius of the white dwarf. By assuming that our planets are solid and spinning, we combined Eq. 1 from Bear & Soker (2013) with Table 1 from Veras et al. (2017b) to obtain a Roche radius (for both planets) of $r_{\rm Roche}=0.00316$ au $=53.2R_{\star}$ $=0.68R_{\odot}$. This location, which corresponds to an orbital period of about just two hours, is much closer than the standard Roche radius for rubble piles, at about $1.0R_{\odot}$ and corresponding to an orbital period of about 4.5 hours (Vanderburg et al., 2015). The smaller core fragment from Manser et al. (2019) has an orbital period of about two hours, but is intact at its location (whereas our planets would start to break up).

The dynamic viscosity ($\eta$) of planetary cores within our solar system have been estimated to range from $10^{21}-10^{28}$ Pa$\cdot$s (Hurford et al., 2018; Patthoff et al., 2019). Despite this large range, the lower bound is a useful constraint, because without it, the range of consideration would increase by at least five orders of magnitude for planetary objects which are not solid cores. Here, we assumed three values for planet viscosity that span this range: $\eta=10^{21},10^{24},10^{27}$ Pa$\cdot$s.

For the remaining three variables, as explained earlier we adopt the ranges $B_{\star}=10^{3}-10^{9}$ G and $\left\langle\gamma_{\ast}\right\rangle=10^{-1}$ – $10^{4}$ S/m, and we did so at order-of-magnitude resolution. We made no a-priori assumptions about $a$, except to sample it at multiples of $r_{\rm Roche}$. We can also consider the critical $a$ beyond which maser emission would no longer be detectable. Willes & Wu (2005) suggested that the induced potential across a detectable planet would exceed about 1 kV.444Table 2 and Figs. 3-4 of Willes & Wu (2004) suggest that much smaller induced potentials would also be detectable, but we are being conservative here. Rewriting their equation (1) in the SI unit system then gives555Equation (1) of Wu et al. (2002) also gives this induced potential, but with the asynchronism parameter (which we ignore here for our estimate).

[TABLE]

Hence, for our $0.5M_{\oplus}$ planet, $a_{\rm crit}=7.6-400r_{\rm Roche}$ for $B_{\star}=10^{3}-10^{9}$ G. For our $5.0M_{\oplus}$ planet, $a_{\rm crit}=9.5-500r_{\rm Roche}$ for $B_{\star}=10^{3}-10^{9}$ G. Fortunately, these critical semimajor axis values are high enough to not actually limit survivability; at these values of $a_{\rm crit}$ we found that the planets survive for at least 1 Gyr across the parameter space.

We adopted integration timescales of 1 Gyr because of both computational speed and due to 1 Gyr being sufficiently long compared to a typical white dwarf cooling age. Models for scattering planetary cores towards white dwarfs currently remain too unconstrained to place reliable bounds on when they arrive and how frequently they arrive, properties which would have helped us determine the ideal integration times.

3.2 Simulation results

Our primarily goal is to determine approximately how long planets which are subject to both gravitational tides and the Lorentz drift can survive orbiting a white dwarf within $a_{\rm crit}$ (equation 3), rather than provide details of their motion in individual cases. These details would likely change with more sophisticated models.

Surviving planets may alter their orbits in non-obvious ways due to tides alone (Veras et al., 2019b), and especially now when coupled with an electromagnetic force. Nevertheless, conveniently we found that the planets do not move outward enough to exceed $a_{\rm crit}$.

We assumed that a planet survives if its orbit never intersected with the sphere of radius $r_{\rm Roche}$; the planet was therefore classified as destroyed if its orbital pericentre was shifted inside of $r_{\rm Roche}$. Because we seek just order-of-magnitude timesacle estimates, we visualize some of our results through colour-contoured number grids (Figs. 2-4). The $x$-axis on each grid is $\log{B_{\star}}$ and the $y$-axis is $\log{\left\langle\gamma_{\ast}\right\rangle}$. Although the product $B_{\star}^{2}\left\langle\gamma_{\star}\right\rangle$ could have instead be presented on one axis as a single degree of freedom, we favour the clarity of our current presentation.

The numbers within the grids in Figs. 2-4 correspond to survival timescales $t_{\rm surv}$ according to

• •

“0”:       $t_{\rm surv}<10$ yr.

• •

“1-8”:    $10^{1-8}\leq t_{\rm surv}<10^{2-9}$ yr.

• •

“9”:       $t_{\rm surv}>10^{9}$ yr.

• •

“$\ast$”:        Equation set becomes stiff due to the asynchronism parameter.

The background colours for the numbers were chosen according to the planet’s level of safety or danger. Dark (0) or light red (1-3) values indicate that destruction occurs quickly enough such that these planets are unlikely to be observed, whereas light green (8) or dark green (9) indicates that the planet may survive for at least one per cent of the white dwarf cooling age.

Consider first the dependence of $t_{\rm surv}$ on dynamic viscosity. In Fig. 2, we paint a picture of the $\left\{B_{\star},\left\langle\gamma_{\star}\right\rangle,M,a_{0}\right\}$ parameter space for the limiting case of $\eta=10^{21}$ Pa$\cdot$s (the lowest supposed dynamic viscosity for a planetary core). This low value allows tides to dominate when the magnitudes of $B_{\star}$ and $\left\langle\gamma_{\star}\right\rangle$ are low (at the bottom lower left corner of each grid). When $\eta$ is then increased by three orders of magnitude to $10^{24}$ Pa$\cdot$s (Fig. 3), tides no longer dominate, even for the lowest values of $B_{\star}$ and $\left\langle\gamma_{\star}\right\rangle$. A direct comparison of Figs. 2 and 3 demonstrates a 0-3 order-of-magnitude difference in survival timescale when the core’s viscosity is increased. We also ran simulations for the $\eta=10^{27}$ Pa$\cdot$s case, but the results are indistinguishable from those in Fig. 3.

Even when tides do not dominate the motion – at the highest values of $B_{\star}$ and $\left\langle\gamma_{\star}\right\rangle$ (upper right corner of each grid) – tides still play a role in the motion. Comparison of the bottom two panels within each of Figs. 2 and 3 reveals that the survival timescales for the highest values of $B_{\star}$ are not equivalent. Although tides do not dominate the motion here, their coupling with the Lorentz drag affects the entire dynamical system.

The next dependence to note in Figs. 2 and 3 is on $a_{0}$. Both tides and the Lorentz force are steep functions of distance to the white dwarf. For low-viscosity planets (Fig. 2), we find that none of them survive for at least 1 Gyr unless $a_{0}\gtrsim 5r_{\rm Roche}$. The survival times of these planets at the highest values of $B_{\star}$ and $\left\langle\gamma_{\star}\right\rangle$ increase by at least four orders-of-magnitude just by incrementing $a_{0}$ from $4r_{\rm Roche}$ to $5r_{\rm Roche}$. A similar trend is seen for more viscous cores (Fig. 3). In general, planets survive for longer the further away they start out from the star. Because a distance of $5r_{\rm Roche}$ is under $a_{\rm crit}$ for all values of $B_{\star}$, this interval ($5r_{\rm Roche},a_{\rm crit}$) contains potentially detectable planets, and for over a large region of the remaining parameter space.

Although $t_{\rm surv}$ appears to increase monotonically with $a_{0}$, Fig. 4 illustrates that this trend is only secular. In particular, the upper panel of the figure exhibits a few cases where increasing $B_{\star}$ actually increases the survival time, and how $t_{\rm surv}$ actually decreases from the top to the middle panel for the most magnetically active system. This behaviour is due entirely to the coupling of gravitational tides with the Lorentz drift: we re-ran the same $a_{0}=8r_{\rm Roche}$ grid simulations without tides, and the result was a grid entirely filled with “9”s. In other words, even for the most magnetically active systems, the Lorentz drift on its own would not drag the planet into the white dwarf within 1 Gyr for $a_{0}\geq 8r_{\rm Roche}$.

The final correlation to consider is with planet mass, which strongly affects gravitational tides by themselves. Here, however, comparison of the right and left columns of each of Figs. 2 and 3 show differences on the order of just one order-of-magnitude for $t_{\rm surv}$ as the planet mass is varied by one order-of-magnitude.

4 Applicability to other planetary systems

The set of equations we used to obtain our results may also be applied to other planetary systems where a unipolar inductor circuit has been established between a star and a solid planet that can both be characterised with Maxwell rheologies. Our model can thus help characterise planetary systems which do not fall under the remit of the few other investigations which have so far considered the interplay between gravitational tides and magnetic torques (Bouvier & Cébron, 2015; Strugarek et al., 2017). These studies focused on protostars and hot Jupiters, T-Tauri stars and hot Jupiters, K stars and hot Jupiters, and M-dwarfs and Earths. Hot Jupiters require a fundamentally different treatment of gravitational tides than do Earths (Ogilvie, 2014; Mathis, 2018). Further, even amongst rocky bodies, there is a wide variety of potential star-planet tidal formulations; the one we have adopted from Veras et al. (2019b) contains an arbitrary frequency dependence on the quality functions, thereby avoiding pitfalls highlighted by Efroimsky & Makarov (2013) and Makarov (2015).

Although our numerical implementation would not change when applied to a different type of host star, the physical parameters of the star would. One of the most consequential changes would be for $\eta^{\star}$, whose value for white dwarfs ($10^{7}$ Pa$\cdot$s) is many orders of magnitude lower than what we might expect for other types of stars. Hence, in other systems, both stellar and planetary tidal forces may be comparable, complicating the trends observed in this study when coupled with the Lorentz drift.

5 Applicability to GD 356 and GD 394

Although $\eta^{\star}$ could change by many orders of magnitude for other types of stars, we do not expect it to vary significantly amongst different white dwarfs. Two white dwarfs which represent tantalising candidates for potentially hosting planetary cores are GD 356 (Ferrario et al., 1997; Wickramasinghe et al., 2010) and GD 394 (Dupuis et al., 2000; Wilson et al., 2019). GD 356 is known to be magnetic, whereas a magnetic field has not yet been detected in GD 394.

The energy source which heats GD 356 remains unknown. However, Wickramasinghe et al. (2010) strengthened the unipolar circuit hypothesis by using Spitzer observations to robustly constrain the mass of a planetary companion to be under 12 Jupiter masses. Also, that star harbours a magnetic field strength on the order of 10 MG. Based on our simulations, this (high) value places further constraints on a potential companion: it must reside beyond about $3-4r_{\rm Roche}$ in order to survive for longer than about $10^{7-8}$ yr. Further, the spin period of the white dwarf is about two hours, an order of magnitude shorter than the period we adopted here. Consequently, the region of parameter space where the tidal and Lorentz contributions are comparable would differ from the region identified in the simulations performed here, but most likely lie within $3-4r_{\rm Roche}$ and hence negligibly affect a detectable planet.

For GD 394, flux variations in the EUV led Dupuis et al. (2000) to propose that the white dwarf is accreting locally, in a surface accretion spot. More recent observations (Wilson et al., 2019) show no evidence for this spot, perhaps indicating that it is no longer currently active. Regardless, heating from the unipolar circuit does occur locally, at the “feet” of the flux tube. Therefore, searching for localized heating in other white dwarfs is a potential pathway for discovering magnetic planetary cores.

6 Summary

White dwarf planetary systems are ubiquitous, with an occurrence rate similar to that of main sequence planetary systems. Although observations of these evolved systems are dominated by planetary remnants in the form of metal pollution and circumstellar dust and gas, individual objects have now been seen orbiting these types of stars (Vanderburg et al., 2015; Manser et al., 2019). The most recent find is a high-density metallic core fragment orbiting within one Solar radius of the white dwarf.

Prospects for finding additional planets have received a major boost with an order of magnitude increase in the known population of white dwarfs in the year 2018 alone from Gaia Data Release #2 (Gentile Fusillo et al., 2019), as well as the successful launch of TESS and commissioning of LSST (Cortes & Kipping, 2018; Lund et al., 2018). Further, the spectroscopic technique utilized by Manser et al. (2019) to discover the core fragment may be applied to other white dwarfs which contain gaseous discs. This disc population is also expected to increase, particularly with the expected launch of Euclid.

As the field of post-main-sequence planetary system science grows, so does the need for models which can assess the likelihood of discovering planets at particular locations and with particular physical properties. In this paper, we have revisited the possibility of a radio-loud conducting planetary core orbiting a white dwarf through a closed unipolar circuit which survives long enough to be detected (Li et al., 1998; Willes & Wu, 2004, 2005). We provided a major update to this model by (i) incorporating gravitational tides from Veras et al. (2019b) into the computation in combination with Lorentz drift, and (ii) by performing simulations across the entire range of observable magnetic field strengths and expected electrical conductivities of white dwarf atmospheres. We found:

• •

Planetary survival times of $t_{\rm surv}>10^{8}$ yr are common in regions of parameter space where the induced potential across the planet exceeds 1 kV.

• •

Lorentz drift (equation 1) and gravitational tides (Eqs. C1-C10 of Veras et al. 2019b) are coupled. Although each force dominates in a different region of parameter space, in other regions they provide feedback on each other.

• •

Planet migration may be inwards or outwards, but nearly always is inwards.

• •

$t_{\rm surv}$ most strongly depends on the planet’s dynamic viscosity $\eta$, the star’s magnetic field strength $B_{\star}$ and the initial star-planet distance $a_{0}$. The next strongest dependencies are on the white dwarf’s electrical conductivity $\left\langle\gamma_{\star}\right\rangle$ and the planet mass $M$.

• •

$t_{\rm surv}$ also is dependent on, but usually weakly, the spins of both the white dwarf and planet, the planet’s density, the star’s mass and radius (which are well constrained for white dwarfs; Tremblay et al. 2016), the star’s dynamic viscosity (as long as it is not too many orders-of-magnitude higher than $10^{7}$ Pa$\cdot$s), both the planet’s and the star’s compliances (as long as they do not vary too far from zero), and the eccentricity and inclination of the initial orbit.

• •

Planetary tides dominate the motion when $B_{\star}\approx$ kG-MG, $\eta\sim 10^{21}$ Pa$\cdot$s and $a_{0}\leq 4r_{\rm Roche}$.

• •

Planetary tides never dominate the motion when $\eta\gtrsim 10^{24}$ Pa$\cdot$s.

• •

The highest observed magnetic fields ($\sim 10^{3}$ MG) rarely allow for $t_{\rm surv}$ to reach as high as $10^{8}$ yr, even at distances of tens of $r_{\rm Roche}$.

Acknowledgements

We thank the referee for helpful comments which have improved the manuscript. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958 through the Kavli Institute for Theoretical Physics programme “Better Stars, Better Planets”. DV gratefully acknowledges the support of the STFC via an Ernest Rutherford Fellowship (grant ST/P003850/1). The Centre for Exoplanets and Habitability is supported by the University of Warwick, and the Center for Exoplanets and Habitable Worlds is supported by the Pennsylvania State University and the Eberly College of Science.