The giant nature of WD 1856 b implies that transiting rocky planets are rare around white dwarfs

White dwarfs (WDs) have roughly Earth-sized radii - a fact long recognized to facilitate the potential discovery of sub-Earth sized planets via transits, as well atmospheric characterization including biosignatures. Despite this, the first (and still only) transiting planet discovered in 2020 was a roughly Jupiter-sized world, found using TESS photometry. Given the relative paucity of giant planets compared to terrestrials indicated by both exoplanet demographics and theoretical simulations (a"bottom-heavy"radius distribution), this is perhaps somewhat surprising. Here, we quantify the surprisingness of this fact accounting for geometric bias and detection bias assuming 1) a bottom-heavy Kepler derived radius distribution, and 2) a top-heavy radial velocity inspired radius distribution. Both are concerning, with the latter implying rocky planets are highly unusual and the former implying WD 1856 b would have to be highly surprising event at the<0.5% level. Using an HBM, we infer the implied power-law radius distribution conditioned upon WD 1856 b and arrive at a top-heavy distribution, such that 0.1-2 REarth planets are an order-of-magnitude less common than 2-20 REarth planets in the period range of 0.1-10 days. The implied hypothesis is that transiting WD rocky planets are rare. We discuss ways to reconcile this with other evidence for minor bodies around WDs, and ultimately argue that it should be easily testable.


INTRODUCTION
WD 1856 b (more formally WD 1856+534 b) is the first transiting exoplanet discovered around a white dwarf (WD) star; indeed it remains the only such world even three years after its announcement (Vanderburg et al. 2020).In many ways, the detection of such a planet was highly anticipated, with indirect evidence for their existence having been amassing for years prior to WD 1856 b.For example, 25-50% of WDs exhibit the spectral fingerprints of metal contamination (Koester, Gänsicke, & Farihi 2014), indicative of in-falling planetary debris (Jura & Young 2014;Bonsor et al. 2020).Such debris was even observed in transit for WD 1145+017 (Vanderburg et al. 2015).
Encouraged by these early signs, searches for WD exoplanets had been attempted in previous years (Fulton et al. 2014;van Sluijs & Van Eylen 2018).Fueling these searches was the potential opportunity to spectroscopically characterize the atmosphere of temperate, rocky planets around such stars using existing facilities to exquisite precision (Limbach et al. 2022).Indeed, Earth-like planets around WDs represent the only class of planet for which biosig-⋆ E-mail: dkipping@astro.columbia.edunatures could be detected with existing facilities, namely JWST (Kaltenegger et al. 2020).Although planets of any size around WDs provide important clues to post main sequence evolution (Veras 2021), the potential for observational evidence of life -arguably the most impactful possible scientific discovery -with existing facilities is completely unique and thus has an outsized influence in these discussions.
Despite this excitement, there are perhaps some reasons for concern hinted at in the very discovery of WD 1856 b.Specifically, it is a giant exoplanet with a radius of (10.4 ± 1.0) R ⊕ (Vanderburg et al. 2020) 1 , which is argued in what follows to be intuitionally surprising.
First, it appears unlikely that WD 1856 b is a so-called "second-generation" planet forming from a WD disk post-MS (main sequence).van Lieshout et al. (2018) consider this formation channel and find that minor planets, up to potentially Earth-mass, are plausible -but a ≳ 100 M ⊕ is certainly not anticipated.Accordingly, WD 1856 b likely formed along soon after the star itself formed, and subsequently survived the RGB and AGB phases (Vanderburg et al. 2020).In this way, the planet's current orbital period of 1.4 days is surely not the MS orbit, since this would be deep inside the star.Instead, the planet likely lay beyond ≳ AU and migrated inwards via planet-planet scattering (Nordhaus & Spiegel 2013;Mustill et al. 2018), Kozai-Lidov oscillations (Stephan, Naoz, & Zuckerman 2017) or even common-envelope evolution (Lagos et al. 2021).
With this point established, it's not constructive to compare the radius-period tuple of WD 1856 b with demographics of MS stars.However, assuming that the planet's current radius is approximately representative of its MS radius, then it is useful to compare its radius (in isolation) to the exoplanet demographics around FGK stars.A common result from exoplanet demographics studies is that Jupitersized/mass planets are less common than their sub-Jovian counterparts (e.g.Howard et al. 2012, Wittenmyer et al. 2016, He, Ford, & Ragozzine 2019), a result also supported by theoretical simulations of planet formation (Mordasini et al. 2012).Despite this, one might argue that WD 1856 b's giant radius is not surprising though because giant planets typically have much deeper transits, as well as longer durations and enhanced transit probabilities (Beatty & Gaudi 2008).However, the roughly Earth-like radius of the host star here certainly erodes away the first (and most major) of those advantages, and thus at least intuitionally there is something surprising in that the first WD transiting exoplanet is a giant planet.
This argument is qualitative and hand-wavy as presented thus far, but ultimately motivated us to investigate this question more rigorously.In this work, such a calculation is performed where the surprisingness of WD 1856 b's giant radius is quantified in Section 3. Following this, we ask what kind of radius distribution is required to dissolve this tension using a Hierarchical Bayesian Model (HBM) in Section 4.However, both of these sections requires a completeness model for TESS, which is introduce first in Section 2. Although a single data point, WD 1856 b in this case, has limited power to make claims about an entire population, we discuss the provocative implications of our work in Section 5, which fortunately will be testable in coming years.
1 Assuming what is very likely a circular orbit given the compact orbital period of 1.4 days.

Overview
A basic component in considering the likelihood of detecting a planet such as WD 1856 b is the detection probability, also known as the sensitivity, true positive rate or completeness.Vanderburg et al. (2020) report that WD 1856 b was initially detected in just a single sector of TESS data (Sector 14) and thus in this work we strictly consider what other exoplanet radii and orbits would be equally detectable in that sector.However, this requires us to determine some function which evaluates completeness given an input set of planet parameters.
The basic starting point in this discussion is signal-tonoise ratio (SNR), assumed in what follows to be the cumulative SNR of multiple transits within a given time frame.We expect the true positive rate to directly map onto SNR, but before we discuss this mapping procedure, we first need to define SNR.Given the grazing nature of WD 1856 b, a box-like transit assumption is a poor approximation and thus so too are traditional transit SNR formulae of the form δ / √ W , where δ is the transit depth and W is the full width half maximum duration.This issue is explored in detail in Kipping (2023), where it is shown that if limb darkening is non-negligible (as is true in TESS bandpass), the numerical integration is the only path to accurately computing the SNR.
Accordingly, our first task is to generate an efficient algorithm into which we can take in a set of planet input parameters and quickly return the corresponding SNR -assuming the same photometric precision and data volume as measured in TESS Sector 14 for WD 1856.

Computing SNRs
To compute the transit SNR of a given simulated transit we follow a series of steps.First, given a set of input parameters, a Mandel & Agol (2002) light curve is generated for a single transit from −0.025 days to +0.025 days (which if sufficient to span even our longest period planets on equatorial transits) using a cadence of one hundredth the integration time, which itself is set equal to that of TESS's Sector 14 data (2 minutes).Integration time effects are handled using numerical resampling following Kipping (2010) with N resamp.= 30.We use a quadratic limb darkening law, interpolating the Claret et al. (2020) WD DA Z = 0 table of coefficients (Table 100C) for log g = 7.915 and T eff = 4710 (best reported values from Vanderburg et al. 2020) to give u 1 = 0.0101 and u 2 = 0.4091.
Next, the simulated light curve is adjusted for blended light, following the prescription of Kipping & Tinetti (2010).However, the challenge here is that the blend factor was never reported in Vanderburg et al. (2020), as their final set of parameters instead using photometry from GTC and Spitzer.These were obtained with much tighter angular apertures for which contamination was simply not an issue.If we similarly assume the GTC photometry is unblended 2 , The giant nature of WD 1856 b 3 we can jointly fit the TESS photometry3 with the GTC photometry to infer the most probable blend factor of each sector, assuming the radius of the planet is approximately achromatic.
We thus obtained the corrected GTC photometry (A.Vanderburg, priv. comm.), detrended the TESS photometry ourselves using CoFiAM (Kipping et al. 2012), and then jointly fit both data sets using a Mandel & Agol (2002) light curve model and MultiNest (Feroz, Hobson, & Bridges 2009).The maximum a-posteriori blend factor for Sector 14 was determined to be (F WD 1856 + ∑ F blend,i )/F WD 1856 = 9.580.With this blend factor, we can now correct the simulated light curves to the correct depths that TESS would observe.
Next, the light curve is spline interpolated and the interpolation function is numerically integrated to compute the SNR of a single transit following Equation (15) of Kipping (2023).This formula requires inputting a noise value as computed over the same time unit basis as the photometry (in our case days), σ .Hence, σ will equal the typical measurement uncertainty of the TESS Sector 14 WD 1856 photometry divided by the square root of the number of points per day, or simply σ = σ 1 √ cadence.In our case, we set σ 1 to be the median uncertainty of the whole sector's set of reported uncertainties.
To account for multiple transits, we multiply the singletransit SNR by √ N, where N is the number of epochs occurring within Sector 14, calculated using the simulated ephemeris parameters.To account for data gaps, we only count epochs for which the expected time of transit minimum has a nearest Sector 14 data point within one cadence.The final SNR is then stored to file.
In practice, the above was too computationally expensive to practically integrate into a Bayesian inference model with large numbers of calls, as used later in Section 4. We thus decided the best approach was to calculate a grid of SNRs based on a broad range of input parameters, and then use this grid to train an interpolation model that would be far faster to evaluate than the SNR numerical integration each time.
Accordingly, the SNR calculation was repeated across a grid of input ratio-of-radii, p, impact parameters, b, and orbital periods, P. We define a log-uniform grid in p such that the minimum value corresponds to 0.1 R ⊕ and the maximum is 20 R ⊕ (assuming R ⋆ = 0.0131 R ⊙ ; Vanderburg et al. 2020).The impact parameter grid is also log-uniform from 0.01 to (1+ p max ).Finally, the orbital period grid is log-uniform from 0.1 to 10 days.We carefully selected the grid to be approximately cubical but reach ∼100, 000 grid points of non-zero SNR (which was sufficiently dense but not computationally prohibitive).We elected to use 54 grid points for p and b, and 47 for P, giving 134,514 grid positions but only 103,447 led to non-zero SNR (corresponding to non-transiting geometries).

SNR of WD 1856 b
At this point, it's worth discussing the SNR of WD 1856 b specifically.From Vanderburg et al. (2020), we know that the TESS SPOC pipeline identified WD 1856 b from Sector 14. Vanderburg et al. (2020) do not report the Sector 14 TESS SNR, but the ExoFOP DV report states a multi-event statistic (a proxy for SNR) of MES= 12.1.
Let us now compare this to value calculated using the same sector but with the Kipping (2023) SNR definition.That work presents both a discrete and continuous version of the SNR expression; the former is intended for real data and the latter for predicting the SNR of model functions.Starting with the discrete version, we can calculate the ∆χ 2 between the best-fitting transit model (taken from our earlier light curve fits in Section 2.2) and a null (flat-line) model on a point-by-point basis, which gives SNR= 18.80 for Sector 14.
Alternatively, we can interpolate the best-fitting model and then numerically integrate the function setting the noise to be equal to the sector-median TESS reported uncertainties, which gives SNR=19.51 for Sector 14.The value from the continuous method is slightly different to the discrete version, but still within 10%.This happens because the latter tracks the actual observed noise behaviour whereas the former assumes expectation value Gaussian statistics (i.e.averaged over many experiments versus a single stochastic realization).
In any case, both theoretical SNRs are ∼60% higher than that reported in the TESS Sector 14 DV report.This may be due to a possible difference in how SNR is defined.In particular, the Kipping ( 2023) is in principle the maximum possible value under all definitions and obviously post-dates the TESS pipeline development.Further, our SNR values benefit from a refined best-fit model thanks to the GTC data and 12 other TESS sectors, will inevitably leads to a differences between the best-fitting model we use to calculate SNR and that used by the TESS pipeline using Sector 14 data alone.Recall -SNR is an act of model comparison and thus always requires a formal model.
Perhaps the most likely explanation is simply that the SPOC pipeline is optimized for much longer transits than that of WD 1856 b and thus plausibly underestimates the significance of transits like this.In any case, the SNR is certainly sufficiently high that we should expect that the completeness probability was very high, as discussed further in the next subsection.

Completeness Curve
A basic question is -how should we map the true-positive probability (TPP), often called the completeness, to SNR?We should expect TPP to be a monotonically increased function with respect to SNR, with boundary conditions of TPP= 0 at SNR= 0 and TPP= 1 at SNR= ∞.This actual function that represents this curve really comes down to the detailed mechanics of how the TESS pipeline identifies, vets and flags TESS Objects of Interest (TOIs).
The precise shape of this function was a subject of intense discussion and iteration in the previous Kepler Mission, since it strongly affects occurrence rate calculations.For example, Howard et al. (2012) attempted the first such occurrence rate calculation prior to any effort to measure this completeness curve and thus made the simplifying assumption of a Heaviside Theta function with a trigger value at SNR= 10.This was later replaced with the linear ramp of Fressin et al. (2013).Later, Christiansen et al. (2016), presented a calibrated curve that resembled a logistic function and equally the cumulative density function of a gamma distribution.However, it should be noted that it was never discovered why these functions truly worked from a first principles perspective (J.Christiansen, priv. comm.), nor whether the trained input parameters would be in anyway generalizeable to future missions like TESS.
Unfortunately, to our knowledge, there is no published TESS general-purpose completeness curve.Attempting to write an TESS pipleine emulator and calibrate its completeness function is a formidable task far beyond the scope of this work.Instead, for the sake of simplicity, we will assume that the completeness function is a Heaviside Theta function with a trigger value at SNR= 10.This is the same assumption made by Howard et al. (2012) working with Kepler data prior to the any Kepler completeness maps being published.We justify this choice that the TESS pipeline uses the same MES (multiple event statistic -a proxy to SNR) threshold as Kepler of 7.1 to flag candidates (Guerrero et al. 2021), and for Kepler at least the completeness function reaches 90% beyond 10 (Christiansen et al. 2016).Whilst far from perfect, we consider this the most pragmatic approach in the current situation.

Interpolating Detection Probabilities
Equipped with our grid of 103,447 training SNRs, we are ready to define an interpolation function.We initially tried out-of-the-box interpolators and machine learning packages but found that they did not reliably reproduce the expected behaviour smoothly and so took a slightly different approach.
We begin by extracting a single slice of constant orbital period from our training data and filling out any missing grid positions to create a regular {p, b} → SNR grid mapping.We train use a Hermite cubic interpolator on this grid.Using this function, we input b as its first unique grid value and then freely vary p until we obtain SNR = 10 (our critical threshold for a detection).This p value is saved as the corresponding critical p, p crit , and then we repeat for all b grid positions.Near the edges of the grid, we switch to 2nd and 1st order Hermite interpolation as appropriate.The resulting list of {p crit , b} positions define the iso-SNR10 contour for this particular choice of P. We thus repeat for all P grid positions to create 47 such contours.
For a fixed choice of b, we expect p crit to monotonically increase with P. We largely see this behaviour but the contours get very close together so to ensure this behaviour against numerical error we regressed quadratic functions in each b-slice and use these predictions to slightly refine the contours.The final contours are then saved to file and used in what follows.

Simulating a Catalog of Detected WD Transiting Exoplanets
Armed with the completeness model from Section 2, we are now well-positioned to produce a simulated catalog of detected WD transiting exoplanets.This is accomplished with three ingredients.First, we need to simulate a population of WD exoplanets with random but representative physical/orbital parameters.Second, we need down-sample to only those which would transit WD 1856.Third, we need to down-sample further to only those which would be detectable using one sector of TESS data, as occurred for WD 1856 b (Vanderburg et al. 2020).
Step three has already been solved for in Section 2.
Step two is trivially through simple geometry.In our implementation, a simulated planet has assigned orbital period, which is converted into a/R ⋆ via Kepler's Third Law and assuming M P ≪ M ⋆ .The planet is then assigned a uniform random impact parameter between 0 and a/R ⋆ , corresponding to a uniform distribution in the cosine of inclination angle (i.e. a isotropic distribution).Only planets for which b < (1 + R P /R ⋆ ) are flagged with a successful transit flag.We prefer to use this brute force approach rather than generating a random Bernoulli variable with the transit probability because it allows us to assign and track and actual b value to each planet which is then consistently used in the step three completeness calculation.
This just leaves us with step one, generating the actual tuples of period-radius.In what follows, we only consider periods between 0.1 and 10 days, since interior to this tidal destruction is almost guarenteed and exterior the transit probability drops off precipitously.The period distribution of WD exoplanets within this range is wholly unknown at this point.However, Kepler statistics broadly resemble a log-uniform distribution (Foreman-Mackey, Hogg, & Morton 2014), especially in fairly localized windows of period space like this.Further, the period distribution is really a nuisance parameter in this work and our central interest is the marginalized result for the radius in any case.
Proceeding with this period distribution, we now just need a radius distribution.Unfortunately, there is no observed inference of the radius distribution of 1 − 10 AU exoplanets.We basically have two possible options.First, we can take the well-measured radius distribution of < AU exoplanets (specifically using Kepler data) and extrapolate it out to to longer periods.Or, alternatively, we could take the mass distribution of 1-10 AU giant planets found from radial velocity surveys and extrapolate it down to smaller masses (and also convert masses to radii).In truth, neither is likely to resemble reality but they do produce radically different radii distributions which at least serve as mirror hypotheses to compare to.A third, and arguably least worst option, is to use neither forward model and instead just simply attempt to infer the radius distribution from the single datum of WD 1856 b, which is precisely what we will do in Section 4. But we proceed with these two forward models in what follows to at least highlight the surprisingness of WD 1856 b conditioned upon them, however questionable their extrapolated models may be.

A Bottom-Heavy Radius Distribution: Extrapolating
Kepler Statistics to > AU Power-law radius distributions have been a common approach for modelling the Kepler exoplanet radius distribution since the first demographics papers (e.g.Howard et al. 2012).In recent years, these models have typically been modified to a broken power-law structure, with an inflection point, R C , occuring around the mini-Neptune regime (e.g.Burke et al. 2015), which better matches the data.In this work, we elected to use the broken power-law of He, Ford, & Ragozzine ( 2019), which has an inflection point at In that work, several models are presented and we use the "non-clustered" version -given that exoplanet multiplicity is yet to be demonstrated for WDs.This model is calibrated against the Kepler catalog using orbital periods of 3 to 300 days.The 300 day cap exemplifies why this Kepler-based model formally says nothing about > AU exoplanets, and thus we are forced to assume it holds in this longer period regime as an extrapolation.Formally, our radius distribution is where k is a normalization constant set such that the sum of probabilities equals unity.For α 1 and α 2 , each simulation draws a random realization for them based on the reported values by He, Ford, & Ragozzine (2019).Specifically, we draw from asymmetric normals using4 α 1 = 1.08 +0.31  −0.20 and α 2 = −5.30+0.53  −0.40 .In this model (and all other models), the minimum and maximum allowed radius is set to 0.1 R ⊕ and 20 R ⊕ .
For this model, and indeed the others that follows, we apply a final step of planet elimination due to tidal destruction.We remove all simulated planets that orbit within R Roche = 2.44(ρ ⋆ /ρ P ) 1/3 , where for ρ ⋆ we adopt the value of Vanderburg et al. (2020).Since we work in period space, not semi-major space, we convert between them using Kepler's Third Law.Planetary densities are estimated using forecaster (?) set in deterministic mode.Above a Saturnradius, the mass-radius relation is nearly flat (and thus degenerate) and so we simply assume all planets in this regime are approximately Saturn-mass.

A Top-Heavy Radius Distribution: Extrapolating CLS
Statistics to Sub-Neptunes Kepler inferences benefit from the enormous statistical power of the Kepler Mission -a mission designed specifically for this purpose (Borucki et al. 2003).However, like all transit surveys it suffers from a strong bias towards inner orbits and thus has poor coverage beyond an AU.Since WD 1856 b likely originated from beyond this distance, it's unclear that extrapolating the < 1 AU Kepler radius distribution is representative of where this planet spent most of its life during the MS phase.Accordingly, we need a radius distribution for planets at wider orbits.Since only transits provide radius (semi-)directly, other methods sensitive to > AU will unfortunately only provide planetary mass.Nevertheless, we can forecast the corresponding radii using a mass-radius relation (Chen & Kipping 2017).
For our mass distribution, we take the results from the California Legacy Survey (CLS; Howard et al. 2010) presented in Fulton et al. (2021).We start by digitizing their Figure 6, which shows the binned occurrence rate versus exoplanet mass in the range 1-5 AU; a highly plausible initial pre-WD position for WD 1856 b.Unfortunately, the distribution only extends down to 30 M ⊕ , since radial velocities are presently unable to detect >AU terrestrial mass planets.As a result, we will be forced to extrapolate this distribution in mass space to address the central question of this work.
To make progress, we proposed an asymmetric normal distribution as the underlying unbinned occurrence rate distribution (in log-mass space), characterized by a mean, µ, a negative standard deviation, σ − , a positive standard deviation, σ + .
The asymmetric normal was integrated over the four bins of that figure, and then we defined a likelihood function between our integrated model occurrence rate model and the reported occurrence rates in each bin.Since the reported uncertainties are asymmetric, we again use asymmetric normals for the likelihood functions.We then maximized the likelihood function with respect to µ, σ − , σ + and a yaxis scaling parameter, f , obtaining µ = 4.678, σ − = 0.697, σ + = 1.791 and f = 20.706.A comparison of our integrated best fitting model to the original Fulton et al. (2021) plot is shown in Figure 1.We did not pursue an MCMC approach here to infer posteriors for these terms, since with four data points and four unknowns we are precariously degenerate already.
Equipped with this model, we draw a random variate from the unbinned distribution, representing a fair value of log(M P ).Masses exceeding 14 M J are rejected given the upper limit constraint from Vanderburg et al. (2020).The mass is then converted into a radius using the forecaster package (Chen & Kipping 2017), again rejected any radii outside of our range of 0.1 to 20 R ⊕ for which we have computed the completeness.
We note that the resulting radius distribution is radically different from that found using Kepler statistics.If we integrate the He, Ford, & Ragozzine (2019) radius distribution, the probability of obtaining a R P < 2 R ⊕ (potentially terrestrial) exoplanet is 83%, or 75% after trimming tidally disrupted planets.In contrast, the Fulton et al. (2021) distribution gives a probability of R P < 2 R ⊕ of just 0.016%, both with and without Roche limit trimming.In other words, terrestial planets are incredibly rare in this alternative model.This is surely a consequence of the assumed asymmetric Gaussian.Since there is no data below 30 M ⊕ here, the model is simply extrapolating and fundamentally assumes a monotonically decreasing function away from the mean.In reality, the distribution could equally plateau, but terrestrial >AU exoplanets are simply beyond the sensitivity of CLS or indeed any other radial velocity survey (for now), and thus this situation is somewhat unavoidable.Regardless, the point here isn't whether the assumed distribution is truly valid for terrestrial planets, but rather it represents an aggressively "top-heavy" distribution juxtaposing the "bottomheavy" Kepler model.

Evaluating the Surprisingness
In each model, we randomly generated planets until 10 5 transiting, detected planets were found for a WD 1856-twin star.The simulated detections define a probability map of expected yield, and against that we can compare the location of WD 1856 b.The results are summarized in Figure 2.
Consider first the bottom-heavy Kepler -inspired model.In this model, any hopes we might have for Earth-sized planets around WDs are well founded.They exist in abundance, comprising 75% of the exoplanet population and 34% of the detected population.If this model holds, it would be totally reasonable to expect JWST, LSST, etc to soon reveal pale blue dots around these diminutive stars.However, in this model, the probability that the first transiting exoplanet would have a radius of greater than or equal to that of WD 1856 b is just 0.37%.It is perfectly possible that WD 1856 b is simply a rather improbable coincidence in this sense, but a much simpler explanation is that the radius distribution of close-in WD exoplanets does not resemble a bottom-heavy distribution, as assumed here.And of course, such a conclusion would have profound implications for our ongoing searches.
Switching to the top-heavy model inspired from Fulton et al. (2021), the surprisingness of detected a WD 1856 b or larger planet is a completely plausible 73% now.However, this comes at a significant cost.The probability of detecting a terrestrial planet is just 0.016%, less than 1 in 6000.So this top-heavy model is consistent with the first WD transiting exoplanet being a giant, but it projects a pessimistic forecast on the chance of ever detecting terrestrial planets around WDs.
As noted earlier, a reasonable concern here with this top-heavy model is the fact it extrapolates in mass.And indeed similarly a concern with the bottom-heavy model is that it extrapolates in semi-major axis.Regardless, both concur that terrestrial planets around WDs are not ex-pected, given WD 1856 b.As a final approach, we instead simply attempt to infer the radius distribution from the single point of WD 1856 b directly in the next section.

Setting up the HBM
In the previous section, we proposed two possible radius distributions for the WD exoplanet population and evaluated the surprisingness of the first detected transiting planet being a giant planet like WD 1856 b.However, whilst both proposed distributions have physical motivations to justify them, neither can be claimed to be representative since we have no previous population of WD exoplanets from which to calibrate our expectations.Instead of proposing a specific radius distribution and evaluating its surpringness, in this section we infer the radius distribution conditioned upon WD 1856 b.
The parameters of WD 1856 b, θ , can be inferred from transit modeling of the available photometry, as was done by Vanderburg et al. (2020).In principle, one could infer the parameters describing the radius distribution, φ , conditioned upon θ with a simple Bayesian inference model.However, those terms not only have uncertainty but also those uncertainties will be somewhat informed by whatever radius distribution is assumed.This speaks to the need for joint inference of both the underlying radius distribution and the planet parameters, i.e. a hierarchical Bayesian model (HBM).In this framework, we can write that Pr(θ , φ |D) = Pr(D|θ , φ )Pr(θ |φ )Pr(φ ) (2) The first term on the right-hand side is the likelihood function, the probability of obtaining the data given some choice of input parameters.It's easy to see that the likelihood here in our case is independent of the radius distribution parameters, φ , and thus Pr(θ , φ |D) = Pr(D|θ )Pr(θ |φ )Pr(φ ). (3) We next note that our data, D, is really three pieces of data.First, we have the fact that the planet transits the star, denoted by symbol b.Second, we have the fact that this transiting planet was detected by Sector 14 of TESS, denoted by d.And third, we have data concerning the transit shape, for which we'll retain the symbol D, thus giving us (4) We can now expand out the likelihood-like term in here, by noting that Pr For the transit data D, we could use the full TESS light curve but this would add significant computational cost to our inference and is ultimately unnecessary since in our inference these parameters are ultimately nuisance parameters we will marginalize over.Instead, we take refined transit pa-rameters from Xu et al. (2021) for this system, specifically we adopt their white light curve inferred parameters from their Table 4 of R P /R ⋆ ≡ p = (7.86 ± 0.01) and b = (7.75 ± 0.01).We assume both quantities can be described by normal tributions, such that where µ p = 7.86, σ p = 0.01, µ b = 7.75 and σ b = 0.01 (Xu et al. 2021).The term Pr( d| b, θ ) describes our completeness function, already described in detail in Section 2. The next term, Pr( b|θ ), describes the transit probability, here given by where the second-line exploits Kepler's Third Law to remove the a dependency and instead use the mean stellar density for which we fix ρ ⋆ = 3.248 × 10 8 kg m −3 (Vanderburg et al. 2020).Similarly, we fix P = 1.407939 days (Xu et al. 2021) and R ⋆ = 0.0131 R ⊙ (Vanderburg et al. 2020).The above implicitly assumes a uniform prior on impact parameter, b.
We now need to choose Pr(θ |φ ), which really means picking a parametric form for the radius distribution.Although the parameters describing that distribution will be inferred (φ ), a typical limitation of HBMs is that we must still choose the mathematical form of that distribution.Given we only have a single data point, we elected to use a power-law distribution, similar to that used by Howard et al. (2012).We did experiment with a broken power-law but found the parameters could not converge.Accordingly, we have where α is the power-law index and R min & R max are the minimum & maximum radii allowed by the distribution, for which we use the same 0.1 R ⊕ to 20 R ⊕ range over which our completeness model is defined.
Finally, for the hyper-prior, Pr(φ ) = Pr(α), we adopt an unbounded uniform prior such that in practice there's no need to actually include this term.

Results
To infer both θ and φ , which correspond to the parameters R P & b and α respectively, we turn to a Markov Chain Monte Carlo (MCMC) approach.We use a Metropolis sampler with normal proposal functions using one million accepted steps.Burn-in points were removed and mixing and convergence verified by inspection of the chains.Finally, we thinned the resulting chain to 100,000 points.A corner plot of the posteriors is presented in Figure 3.
The posterior distributions for R P and b closely match the strong prior constraint we leverage from Xu et al. (2021), as expected.After checking this basic point, we turn to α -the power-law index for the radius distribution.Since Pr(R P |α) ∝ R −α , then negative indices correspond to topheavy radius distributions (paucity of terrestrial planets), whereas positive indices are bottom-heavy (abundance of terrestrial planets).The median and ±34.15% quantile range of the marginalized a-posteriori distribution of α can be summarized as α = −1.9+1.8 −2.9 .From this value, one might naively conclude that there is a ∼1-σ preference for a top-heavy radius distribution, which implies a less significant discrepancy than that indicated by the summary statistics quoted in the inset panels of Figure 2. The confusion here stems from the fact that a "bottom-heavy" distribution could be reasonably defined as beginning once α > 0, but really α = 0 would lead to a median radius of 10.1 R ⊕ still (i.e.nearly a Jupiter radius).Truly what matters is whether the samples of distribution is dominated by potentially terrestrial planets, defined in this work as occurring at 2 R ⊕ .The transition from a sample dominated by terrestrials ("terrestrial-heavy") versus dominated by Neptunian and Jovian worlds ("terrestrial-light") in fact occurs at α = +0.9.Comparing this to our MCMC posteriors, we find that only 2.8% of samples correspond to α > 0.9 (2.2 σ ) and this is a more accurate assessment of the HBM significance.

Surprisingness of our Inferred Model
As was done in Section 3, we can take this model and ask how surprising it would be for WD 1856 b to be the first detection.Before attempting that calculation, we should expect the surprisingness to be small, since the model is of course conditioned upon the detection of WD 1856 b.
To accomplish this, we draw a random sample for α from our inferred posterior distribution to create a radius distribution realization (100 such realizations are shown in the bottom-right panel of Figure 2).We next draw a random radius from that distribution.As before, the period distribution is assumed to be log-uniform from 0.1 to 10 days.We then draw a random impact parameter for the planet between 0.01 and a/R ⋆ , where the latter is indeed the maximum physically possible impact parameter possible and the former is the smallest impact parameter over which our grid of completeness models were trained (see Section 2).Finally, we evaluate the detection probability using our model from Section 2. If a planet is transiting and detected, the resulting period, impact parameter and radius are saved are plausibly detectable objects.This is repeated until 10 5 such "detections" are accumulated.
The resulting density of realized planets is shown in the top-right panel of Figure 2. From this, and the marginalized version shown in the panel below, we can see that the vast majority of realizations are consistent with the first planet being as large as WD 1856 b, indeed 80% of realizations are such.However, the top-heavy nature of the inferred distributions means that just 5.2% of simulated planets are potentially terrestrial (< 2 R ⊙ ).
This analysis supports the conclusions implied by Section 3; the fact that WD 1856 b was the first detected transiting planet around a WD implies that terrestrial planets are rare around them.

Summary
To date, there is only one known transiting planet orbiting a white dwarf -WD 1856 b (Vanderburg et al. 2020).
The diminutive size of WDs has long been recognized as a pathway to detecting Earth-sized or smaller planets with existing facilities (Agol 2011).Further, there is an emerging view that Jupiter-sized planets represent the minority of the planet population (Wittenmyer et al. 2016;Hsu et al. 2018;Fulton et al. 2021).Thus, the fact that the first transiting planet detected around a WD was found to be a giant planet is somewhat surprising.
We have quantified this surpringness for two possible radius distribution models; one using the Kepler -derived broken power-law model of He, Ford, & Ragozzine (2019) and the other based off the radial velocity survey results of Fulton et al. ( 2021).The former is bottom-heavy (abundance of terrestrial planets) and is in apparent contradiction with the detection of WD 1856 b, with a p-value of 0.37%.However, WD 1856 b may have migrated inwards during the post main-sequence evolution and thus the Kepler demographics (probing < 1 AU) are not representative.Using the Fulton et al. (2021) mass distribution, converted into a radius distribution using forecaster (Chen & Kipping 2017), we find a top-heavy distribution more consistent with WD 1856 b.However, such a distribution implies a dearth of terrestrial planets, with just 0.016% of WDs possessing them in close-in orbits.
Rather than propose speculative radius distributions, we also directly fitted for a radius distribution using an HBM and again find that a top-heavy distribution is indicated.Using a single-power law model, we find that only 5% of WDs would have 0.1-2 Earth radii planets at periods of 0.1 to 10 days.
An implied conclusion from this work is that the great hope that transiting WD planets might prove a shortcut to biosignature detection (Kaltenegger et al. 2020) may be a false hope, at least for exoplanets.Specifically, the deduced hypothesis is that planets with radii 0.1-2 R ⊕ are much rarer (∼ an order of magnitude) than 2-20 R ⊕ around WDs for orbital periods of 0.1-10 days (encompassing most WD habitable zones; Agol 2011).However, we note that exomoons of WD giant planets, should they exist, would then potentially be a remaining source of optimism.

Other Evidence for Exoplanets around WDs
WD 1856 b may well be the only transiting exoplanet known, but other exoplanet candidates have been reported around single WDs, but we note that all of them are indeed Jupiter-mass or higher (Veras, McDonald, & Makarov 2020).MOA-2010-BLG-477L b is a Jupiter-mass exoplanet found at ∼3 AU and detected via microlensing (Blackman et al. 2021).At even an wider separation, WD 0806-661 b is a giant planet or sub-brown dwarf at 2500 AU detected via direct imaging (Luhman, Burgasser, & Bochanski 2011).Both of these are clearly at much wider separations than we have been focussed on in this paper, that is exoplanets close enough for plausible transit probabilities -as well being within the habitable zone (Kilic et al. 2013) to enable biosignature reconnassiance (Kaltenegger et al. 2020).
Perhaps the most useful exoplanet for comparison, although it remains only a candidate (Gaia Collaboration et al. 2023), is Gaia DR3 2098419251579450880 b5 .Discovered through astrometry in Gaia DR3, this is a ∼34 M J candidate companion, orbiting at 242 days.Thus, the candidate is highly unlikely to exhibit transits but perhaps represents the tail-end of the distribution we are seeing in the 0.1-10 days regime (where they are expected).In any case, it is almost certainly roughly Jupiter-sized and thus is fully consistent with the hypothesis proposed in this work.However, unlike the TESS data, the Gaia data is not sensitive to terrestrial planets for these stars (Perryman et al. 2014) and thus we cannot state that its existence adds greater weight to our argument, rather that it merely appear consistent with it.
Concordantly, none of the other known exoplanet candidates around WDs really speak to the claimed hypothesis.A more interesting observation is that many WDs show evidence for metal contaminated atmospheres indicative of accreting material (Koester, Gänsicke, & Farihi 2014).These observations have been interpreted as likely have been caused by tidally disrupted rocky bodies (Jura 2008;Debes et al. 2019) and thus would seem to imply the existence of close-in terrestrial planets afterall.However, this doesn't necessarily contradict this work's hypothesis, since the progenitor population could be one which, if scattered inwards, is consistently scattered inwards sufficient to lead to tidal disruption.
The case of WD 1145+017 (discovered with K2 photometry) perhaps serves as a possible example (Vanderburg et al. 2015), where an ensemble of co-orbiting transits were detected around the host star at approximately the Roche limit, indicative of fragments and dust from a planetary disruption.ZTF J0139+5245 is another such example discovered with ZTF (Vanderbosch et al. 2020), but a more extreme one, since the period is far beyond the circular Roche limit implying either a high eccentricity (e > 0.97) or rotational fission (Veras, McDonald, & Makarov 2020).Some additional possible candidate examples are reported in Guidry et al. (2021).
An alternative resolution is that the minor bodies indicated above are in fact dominated by progenitor moons rather than planets, as been suggested previously (Doyle, Desch, & Young 2021).

Other Possibilities
We can conceive of two other possible ways to resolve the evidence for WD minor bodies with the surprisingly large size of WD 1856 b.The first is that our assumed single power-law is wrong in enforcing a monotonically increasing probability density with respect to radius.Perhaps the distribution turns over at some radius, representing the most unlikely planetary radius, and then peaks back up.There's an infinite number of possible such distributions and thus we haven't attempted to propose specific speculations, but such a distribution would necessarily have to turn over at around an Earth radius in order to be consistent with the WD 1856 b first detection (else terrestrial planets would dominate the detected population again).
A second possibility is that WD 1856 b is simply a fluke.Perhaps there truly is a bottom-heavy distribution and it was indeed highly improbable that a WD 1856 b-sized exoplanet would be the first to be revealed in transit.Using the Kepler -like bottom-heavy distribution, we found the odds of this to be 0.37% earlier in Section 3. That's certainly interesting, but hardly overwhelming -in the history of astronomy, improbable events can and will occur given enough time.
For these reasons, we don't consider our hypothesis in any way established with conviction.It would certainly be premature to abort on-going and future efforts to look for terrestrial planets around WDs.However, it is a curious implication resulting from the discovery of WD 1856 b that deserves further attention.Fortunately, the hypothesis is em-inently testable through on-going searches and thus we can look forward to a firmer resolution in the future.

ACKNOWLEDGMENTS
This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST).We acknowledge the use of public TESS Alert data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center.

DATA AVAILABILITY
We make the two simulated detected populations from Section 3 publicly available, as well as the version for the HBM inferred model from Section 4 (all of which were presented in Figure 2) at this URL.Similarly, we make the HBM posteriors available in the same repository (described in Section 4 and presented in Figure 3) as well as our MCMC inference code.We also make the code underlying the model shown in Figure 1 available, as well as the resulting simulated radius samples.Finally, the simulated Sector 14 SNRs are made available, as well as the corresponding training grid, and the resulting detection contours calculated from them (as well as the code to do).

Figure 1 .
Figure1.Comparison of the original data presented in Figure1ofFulton et al. (2021) for the binned occurence rate of exoplanets in the range 1-5 AU (black), with our proposed asymmetric Gaussian model (red).

Figure 2 .
Figure2.Three models (each column) for the radius distribution of WD exoplanets and the resulting detection yields expected around a WD 1856-twin star with a sector of TESS data.The bottom-heavy Kepler -inspired model (left) produces a yield map (top row) that makes WD 1856 b an outlier.In contrast, the the top-heavy model (middle) naturally explains WD 1856 b but predicts very few terrestrial planets.The right column shows a simple power-law inference using an HBM, which favours to a top-heavy distribution.The middle row shows the marginalized radius dimension of the upper row, and the bottom row shows the intrinsic distributions which from they were drawn.

Figure 3 .
Figure 3. Corner plot of the posteriors from our HBM model.The parameters (b and R P ) and hyper-parameters (just α) are shown together with the priors highlighted with the dashed lines on the marginalized panels.