© The Authors 2026

1 Introduction

Several thousand exoplanets are confirmed to exist around main-sequence stars¹, and on average each star in the Milky Way probably hosts more than one planet (Cassan et al. 2012). Virtually all known planet host stars will end their lives as white dwarfs (WDs) One might therefore expect large numbers of exoplanets to exist around WDs.

Although in specific cases the gravitational interactions between the planets can affect survival (Ronco et al. 2020), in general planets with orbits beyond a few astronomical units should mostly survive unscathed, while planets orbiting at separations closer than ~ 1 au should be engulfed during the giant phases of the host star (e.g., Villaver & Livio 2007; Veras 2016; Nordhaus & Spiegel 2013). Companions that survive the giantstar phases of their host stars are expected to migrate outward due to stellar mass loss.

Several surveys of WDs using various direct imaging techniques to find surviving massive giant planets have been performed without a clear detection (Burleigh et al. 2002; Debes et al. 2005; Hogan et al. 2009; Brandner et al. 2021). These results from high-contrast imaging campaigns are complemented by transit searches using Kepler data (van Sluijs & Van Eylen 2018) and surveys of WDs with infrared excess using Spitzer photometry (Farihi et al. 2008).

Despite these efforts, very few substellar companions to WDs have been detected so far. Two of the detected substellar objects, the hypothetical Neptune-like planet around WD J0914+1914 (Gänsicke et al. 2019) and the transiting Jupiter-like planet around WD 1856+534 (Vanderburg et al. 2020), orbit their central WDs with short orbital periods, which implies that they either must have been affected by other companions, or that they perhaps survived a common envelope phase (Muñoz & Petrovich 2020; Lagos et al. 2021; Maldonado et al. 2021; Chamandy et al. 2021; Stephan et al. 2021). In contrast, the 1.4 Jupiter-mass planet discovered by Blackman et al. (2021) orbits the central WD at a projected separation of 2.8 ± 0.5 au, which suggests that it might have been affected by the metamorphosis of its host star mostly through orbital expansion as a consequence of mass loss. The origin of the very distant (2500 au projected separation) substellar companion candidate to WD 0806-661 (Luhman et al. 2012) and the rocky planet candidates potentially orbiting close binaries at large separations (Thorsett et al. 1993; Sigurdsson et al. 2003; Luhman et al. 2012; Zhang et al. 2024) remains unclear. Brandner et al. (2021) argued that the low number of detected substellar objects around WDs is related to instrumental limitations. This situation might be changing as the James Webb Space Telescope (JWST) and new ground-based, high-contrast imaging instruments (e.g., ERIS at the VLT) have become available. Indeed, three candidate giant planets around WDs have recently been detected using the JWST (Limbach et al. 2024; Mullally et al. 2024), and complementary surveys with ERIS are underway. Furthermore, in the near future large telescopes such as the Extremely Large Telescope (ELT) will probably revolutionize the field.

While the prospect of providing a more robust observational census of planets orbiting WDs in the near future is promising, few theoretical predictions have been made. To the best of our knowledge, the only population-synthesis approach to predicting the occurrence rate of planets around WDs was published by Andryushin & Popov (2021). These authors assumed a rather generic planet population around WD progenitor stars and focused on predicting the fractions of planets that are engulfed or ejected or remain in orbit despite the metamorphosis of their host star. While this in general represents a useful exercise, it does not allow us to predict the relative number of WDs that is expected to host a planet and how this fraction depends on WD parameters such as age (effective temperature) and mass. In addition, by using a generic planet population, Andryushin & Popov (2021) did not take into account possible dependencies of the planet-occurrence rate on stellar mass and metallicity. Therefore, despite clearly being interesting, the results obtained by Andryushin & Popov (2021) do not permit direct comparison with the results of ongoing observational surveys dedicated to planets around WDs.

We followed a different population-synthesis approach to fill this gap. We generated a representative population of WD progenitors based on an age-metallicity relation and the initial mass function and assigned substellar companions with masses from 0.8 up to 30 M_J using occurrence rates derived from radialvelocity surveys of giant stars. We then simulated the evolution of the host stars and the planetary orbits and predict that the general fraction of white dwarfs hosting a substellar companion is below 3 ± 1.5% and that the vast majority of these companions are Jupiter-like gas-giant planets. This result is obtained assuming single planet systems. Taking into account additional companions that potentially expel part of the stellar envelope during the giant phases, this fraction could increase to ≲4.4 ± 2.4%. Furthermore, the predictions of our population synthesis significantly depend on the assumed age-metallicity relation. In what follows, we describe our model and present its main predictions, beginning with a description of the numerical codes used.

2 Survival of substellar objects

To determine the impact of stellar evolution on the orbits of substellar companions, we combined the Modules for Experiments in Stellar Astrophysics (MESA) code (e.g., Jermyn et al. 2023) with an algorithm solving the equations describing the time evolution of the orbit of companions affected by stellar mass loss and stellar tides. This algorithm is called Final exoplAnet orbiTal dEstination around evolved Stars (FATES) and was developed by Ronco et al. (2020). Both codes are described in detail in what follows.

2.1 Stellar evolution with MESA

The MESA code (we used version r-23.5.1) allows users to reproduce the details of stellar evolution, albeit at high computational cost. Consequently, we chose to employ MESA exclusively for stars that host a planet.

The MESA equation of state is a blend of the OPAL (Rogers & Nayfonov 2002), SCVH (Saumon et al. 1995), FreeEOS (Irwin 2004), HELM (Timmes & Swesty 2000), PC (Potekhin & Chabrier 2010), and Skye (Jermyn et al. 2021) EOSes. Radiative opacities are primarily from OPAL (Iglesias & Rogers 1993, 1996), with low-temperature data from Ferguson et al. (2005) and the high-temperature, Compton-scattering-dominated regime by Poutanen (2017). Electron conduction opacities are from Cassisi et al. (2007) and Blouin et al. (2020). Nuclear reaction rates are from JINA REACLIB (Cyburt et al. 2010), NACRE (Angulo et al. 1999) and additional tabulated weak reaction rates (Fuller et al. 1985; Oda et al. 1994; Langanke & Martínez-Pinedo 2000). Screening is included via the prescription of Chugunov et al. (2007). Thermal neutrino-loss rates are from Itoh et al. (1996).

Additional key parameters in MESA for the survival of substellar companions are the mass-loss parameters. During the red giant branch (RGB), we assumed the model proposed by Reimers (1975), and for stars on the asymptotic giant branch (AGB) we used the relation suggested by Bloecker (1995).

Regarding convective overshooting, in principle it should always be included when there is convection. However, while constraints on core overshooting have been derived from aster-oseismology, overshooting in the convective envelope remains poorly constrained. Therefore, we only included overshooting for convective cores following Schaller et al. (1992) and Herwig (2000).

2.2 Orbital evolution of the companions

To model the orbital evolution of the companions over the lifetime of their host star, we used the FATES code developed by Ronco et al. (2020), a Runge-Kutta-Fehlberg algorithm that integrates the equations describing the effects of mass loss and tidal forces on the orbit (see Zahn 1977 and Mustill & Villaver 2012 for details): $$\left(\frac{\dot{a}}{a}\right)=-\left(\frac{{\dot{M}}_{\star }}{M_{\star }+M_{\mathrm{p}}}\right)-{\left(\frac{\dot{a}}{a}\right)}_{\mathrm{t}},$$(1)

where M_p is the mass of the substellar object, M_* is the total mass of the star, and a is the semi-major axis. The first term describes the orbital expansion due to stellar mass loss, which is assumed to be isotropic. Mass loss from the companion is not included, as it is expected to be negligible compared to the stellar mass-loss rate (Villaver & Livio 2009; Schreiber et al. 2019). The term (ȧ/a)_t represents the change of the orbital distance caused by stellar tides and is given by $$\begin{array}{rl}{\left(\frac{\dot{a}}{a}\right)}_{\mathrm{t}}=& -\frac{1}{9{\tau }_{\text{conv}}}\frac{M_{\star }^{\mathrm{e}\mathrm{n}\mathrm{v}}}{M_{\star }}\frac{M_{\mathrm{p}}}{M_{\star }}(1+\frac{M_{\mathrm{p}}}{M_{\star }}){\left(\frac{R_{\star }}{a}\right)}^8\\ & \times [2p_2+e^2(\frac{7}{8}{p}_1-10p_2+\frac{441}{8}{p}_3)].\end{array}$$(2)

Similarly, the eccentricity of the orbit is changing as $$\begin{array}{rl}{\left(\frac{\dot{e}}{e}\right)}_{\mathrm{t}}=& -\frac{1}{36{\tau }_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}}}\frac{M_{\star }^{\mathrm{e}\mathrm{n}\mathrm{v}}}{M_{\star }}\frac{M_{\mathrm{p}}}{M_{\star }}(1+\frac{M_{\mathrm{p}}}{M_{\star }}){\left(\frac{R_{\star }}{a}\right)}^8\\ & \times [\frac{5}{4}{p}_1-2p_2+\frac{147}{4}{p}_3].\end{array}$$(3)

Here, M_*^env is the envelope mass, R_* is the radius of the star, e is the eccentricity of the orbit, τ_conv is the eddy turnover timescale, and p_i represents the frequency components of the tidal force according to Mustill & Villaver (2012), that is, $$p_i\approx \frac{9}{2}min[1,\left(\frac{4{\pi }^2{a}^3}{i^{2}G(M_{\star }+M_{\mathrm{p}}){\tau }_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}}^2}\right)],i=1,2,3.$$(4)

Given that the convective turnover time is uncertain, we considered two different expressions for τ_conv. The first one was presented in Rasio et al. (1996, hereafter R96) and is given by $${\tau }_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}}={\left(\frac{M_{\star }^{\mathrm{e}\mathrm{n}\mathrm{v}}(R_{\star }-R_{\star }^{\mathrm{e}\mathrm{n}\mathrm{v}})R_{\star }^{\mathrm{e}\mathrm{n}\mathrm{v}}}{3L_{\star }}\right)}^{1/3},$$(5)

while an alternative prescription was proposed more recently by Villaver & Livio (2009, hereafter VL09), that is, $${\tau }_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}}={\left(\frac{M_{\star }^{\mathrm{e}\mathrm{n}\mathrm{v}}(R_{\star }-R_{\star }^{\mathrm{e}\mathrm{n}\mathrm{v}}{)}^2}{3L_{\star }}\right)}^{1/3}.$$(6)

In both Eqs. (5) and (6), R_*^env is the radius at the base of the convective envelope and L_* is the luminosity of the star. The two prescriptions for the convective turnover time lead to different results, especially for extended convective envelopes. The turnover times from R96 are shorter, which implies stronger stellar tides for stars on the RGB and the AGB.

2.3 Impact of stellar evolution codes, stellar tides, mass-loss efficiency, and eccentricity

Before carrying out population synthesis for substellar companions to WDs, we tested the impact of stellar evolution and stellar tides on the survival of these companions. As we show here, the modeling of mass loss on the AGB phase and tides can be crucial.

As in most stellar evolution codes, we modeled mass loss on the first giant branch with the prescription by Reimers (1975) with an efficiency of η = 0.5. Mass loss on the AGB is more uncertain. The two preferred prescriptions are those proposed by Bloecker (1995) and Vassiliadis & Wood (1993). Originally, Bloecker (1995) proposed an efficiency on the order of one, but very different scaling factors have been suggested since then, including dependencies of the efficiency on mass, metallicity, or evolutionary state (Pignatari et al. 2016; Huscher et al. 2025). We reflect this general uncertainty by using the prescription from Bloecker (1995) with two different (but constant) values for the mass-loss efficiencies (η = 0.02 and η = 0.7). The lower value is probably more realistic (Ventura et al. 2012; Cinquegrana et al. 2022). We therefore employed the former value for our baseline simulations.

To illustrate the impact of different mass-loss prescriptions and stellar evolution models, we simulated the orbital evolution of a 1 M_J giant planet orbiting a 2 M_⊙ star with nearly solar metallicity (Z = 0.0187), considering a range of initial orbital separations. For simplicity, in this exercise all orbits were assumed to be circular. Stellar evolution was modeled with the simple single star evolution (SSE) code (Hurley et al. 2000) that does not include thermal pulses on the AGB and MESA assuming η = 0.02 and η = 0.7.

Figure 1 presents the evolution of the giant planet’s orbital separation calculated under these three scenarios, assuming the turnover timescale from VL09. Inspecting the two panels with MESA simulations, we note that a higher value for the mass-loss efficiency (middle panel) leads to fewer thermal pulses during the AGB, which in turn limits the maximum stellar radius. As a result, companions can survive at smaller initial separations. Conversely, lower mass-loss efficiencies lead to larger maximum radii, requiring significantly wider orbits to avoid engulfment. The results obtained using the SSE code predict a similar maximum radius to that obtained under high-efficiency mass loss in MESA. Most importantly, SSE does not describe the radius evolution of thermally pulsing AGB stars, which is key for determining the survival of companions.

A second uncertainty that needs to be considered concerns the prescriptions used for stellar tides. In Fig. 2, we highlight how different prescriptions for the convective turnover timescale, which impact the strength of stellar tides, affect the survival of substellar companions for different orbital eccentricities. Stellar tides are stronger assuming the approximation of the turnover time proposed by R96, which implies that substellar companions need to be located at larger initial separations to survive. In general, the larger the eccentricity of the orbit, the larger the semimajor axis needs to be for survival. This is because, for highly eccentric orbits, stellar tides can start to dominate during periastron and may ultimately cause engulfment.

The examples shown in Figs. 1 and 2 illustrate that the survival of substellar objects up to the WD phase depends on the assumptions made for stellar evolution on the AGB as well as on the assumed approximation for the convective turnover timescale. We took these findings into account when designing our population-synthesis approach.

3 Population synthesis

To theoretically predict the population of substellar objects orbiting WDs, we first generated a sample of potential WD progenitor stars taking into account the initial mass function and the age-metallicity relation. We further assumed a companion occurrence-rate distribution (depending on mass and metallic-ity of the host star) consistent with radial-velocity surveys of 0.8-4 M_⊙ giant stars (Wolthoff et al. 2022). We then simulated the future evolution of these planetary systems to estimate the present-day population of substellar companions to WDs. In the following sections, we detail our population-synthesis methodology.

Fig. 1 Orbital evolution of a 1 MJ gas-giant planet orbiting a 2 M⊙ star with Z = 0.0187, comparing calculations performed with SSE (left), the MESA default test suite (middle), for which the AGB mass-loss efficiency is set to η = 0.7, and MESA assuming a more realistic mass-loss efficiency (η = 0.02; right). The orbits are assumed to be circular, and initial separations range from 2 to 4.5 au, with a step size of 0.5 au. The red filled area corresponds to the stellar radius, while purple and green lines denote the orbital separation of the engulfed and surviving planets, respectively. The insets show zoomed-in views of the thermal pulse phase of the AGB and highlight its critical role in planetary engulfment. SSE does not account for the thermal pulses (left). For a large mass-loss efficiency (middle), most planets survive because the orbit expansion caused by stellar mass loss dominates. In the most realistic scenario assuming a small efficiency (right), the star evolves through more thermal pulses and reaches a larger radius, which causes most planets to be engulfed.

Fig. 2 Survival of a 10 MJ gas-giant planet orbiting a 2 M⊙ star for different initial semimajor axes and eccentricities. Red stars indicate engulfment, while blue circles stand for survival of the planet. We assumed two different approximations for the convective turnover time, τconv. Following R96 (right panel), many more planets are engulfed than when assuming the prescription suggested by VL09 (left panel). In both cases, the survival of substellar objects depends on the eccentricity of the orbit, with more eccentric orbits leading more frequently to engulfment. Due to stronger tides, this dependence is more pronounced assuming the approximation by R96. For initial semimajor axes ≥ 8 au, the effects of tidal forces are negligible except for very large eccentricities.

3.1 Initial distributions of potential WD progenitors

The population of potential WD progenitors was generated using the initial mass function from Kroupa (2001) for the mass range from 0.8 to 4 M_⊙. We excluded larger masses because the substellar companions’ occurrence rate is unconstrained. This also means that our predicted sample excludes WD masses exceeding ~0.90 M_⊙. We generated 243 617 stars.

We then assigned each star an age assuming a constant star formation rate using 10 Gyr as the age of the thin disk of the Milky Way. The metallicity of each star was incorporated taking into account an age-metallicity relation by interpolating the parameters defined in Carrillo et al. (2023, their Table 1). This age-metallicity relation was derived from evolutionary simulations of the Milky Way. Therefore our population synthesis does apply to global populations of WDs and not necessarily to local observed populations where the age-metallicity relation might have evolved differently (Fernández-Alvar et al. 2025).

3.2 Probability distributions for substellar objects around potential WD progenitors

The progenitor stars of the current WD population were intermediate-mass stars (i.e., -0.8-8 M_⊙). Detecting substellar companions around main-sequence stars with masses exceeding ~1.5 M_⊙ is challenging as the absence of a large number of absorption lines makes radial-velocity studies extremely difficult, and the larger contrast complicates the direct detection through imaging and transits.

Reliable constraints on the occurrence rate of giant planets and brown dwarfs at large orbital separations around stars that can evolve into a WD within the age of our Galaxy have been derived from radial-velocity studies of giant stars (Reffert et al. 2015; Wolthoff et al. 2022). We used their probability distributions which cover the companion mass range from 0.8 to 30 M_J to generate a representative population of substellar companions around the potential WD progenitors, that is, we assumed that the companion occurrence rate (f) depends on stellar mass and metallicity as $$f(M_{\ast },[\mathrm{F}\mathrm{e}/\mathrm{H}])=C\cdot \exp (-\frac{1}{2}{\left[\frac{M_{\ast }-\mu }{\sigma }\right]}^2)\cdot {10}^{\beta \cdot [\mathrm{F}\mathrm{e}/\mathrm{H}]},$$(7)

with

• C = (0.128 ± 0.027)

• μ = (1.68 ± 0.59) M_⊙

• σ = (0.70 ± 1.02) M_⊙

• β = (1.38 ± 0.36).

The orbital period distribution of the companions is given by (Wolthoff et al. 2022) $$f_{\mathrm{o}\mathrm{c}\mathrm{c}}\left(P\right)=\frac{A}{\sqrt{2\pi }\sigma \cdot P/\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}}\cdot \exp (-\frac{{(\ln (P/\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s})-\mu )}^2}{2{\sigma }^2}),$$(8)

with

• A = (9988 ± 4166)

• σ = (0.89 ± 0.19)

• μ = (7.46 ± 0.45).

We adopted the distribution recently obtained by Kane & Wittenmyer (2024) for the eccentricities of the orbits of the companions. Excluding planets with very short periods (P < 10 days), these authors showed that the eccentricity distribution of giant planets has a mean value of 0.29, a median of 0.23, and a 1σ RMS scatter of 0.23. The distribution is concentrated toward low eccentricities, but it extends up to e ≈ 0.9. These high eccentricities are thought to result from planet-planet scattering and/or the Kozai-Lidov effect (Alqasim et al. 2025).

Finally, we assumed a planet mass function derived from observations of companions to FGK main-sequence stars, which indicate that the mass distribution for companions in the range of 0.3 M_J ≤ M_p ≤ 10 M_J can be reasonably approximated by a power law ($\propto M_{\mathrm{p}}^{-1.3}$, Cumming et al. 2008; Mayor et al. 2011; Adams et al. 2021), and we extrapolated this distribution up to 30 M_J. While the number of companions detected by Wolthoff et al. (2022) was too small to derive a planet mass function, our approach seems to be roughly in agreement with their finding that the number of companions decreases with mass (see their Fig. 7). We furthermore note that by extrapolating the planet mass function up to the brown dwarf regime, we probably overestimated their occurrence. This is because the mass range between 13 and 30 M_J at separations of a few astronomical units encompasses the brown-dwarf desert, where such companions have been shown to be very rare (see, e.g., Grether & Lineweaver 2006; Zandt et al. 2025). However, we estimate that this discrepancy affects our results by ~10% at most, which is significantly smaller than the uncertainties associated with the other adopted distributions.

Using the distributions described above, we assigned a representative number of substellar companions to the population of potential WD progenitors (0.8-4 M_⊙) using a Monte Carlo sampling approach. The resulting fraction of these stars hosting a companion was 5.9 ± 3.4 %. To predict a present-day population of substellar objects around WDs, each star in the population we generated was evolved up to today, taking into account the impact of stellar evolution on the substellar companion orbit, if present.

3.3 Generating a representative WD sample

Equipped with a sample of stars that are potential progenitors of present-day WDs (according to their mass), knowing their metallicity and age, and having assigned a representative population of substellar companions, we simulated the evolution of these systems in order to predict a representative population of present-day WDs including substellar companions.

Given that the sample of stars that we generated was far too large to perform each simulation with MESA even in a supercluster, we decided to evolve all stars from their birth time to today using SSE. We found that of the initial sample of 243 617 stars (16 086 with substellar companions), 113 985 are predicted to be present-day WDs (5611 of them initially hosting a companion). While evolving the host stars with the fast stellar evolution code SSE should not introduce significant errors -the ages and masses of the resulting WDs are similar to those predicted with MESA- the survival of the substellar companions depends on the details of late AGB evolution, especially the thermal pulses and the mass loss. The details of both processes are not well described by SSE (see Fig. 1 for details), and we therefore used more detailed stellar evolution tracks to estimate the survival of the companions.

We randomly selected 1000 stars that have become WDs (according to our SSE simulations) and initially hosted a substellar companion, and we simulated the evolution of the host star with MESA. The selection of 1000 stars corresponds to a fraction of 17.82% of all stars that initially hosted a companion and are currently WDs. To establish our final representative sample of present-day WDs, we therefore randomly selected the same fraction of objects from the sample of stars that have become WDs but did not initially host a companion. This led to predicting 19 315 present day WDs. As the final step, we evaluated the survival of the 1000 substellar companions using the code FATES and the stellar evolution tracks calculated with MESA.

4 Results

Our final sample consists of 20 315 WDs whose parameters were calculated with SSE, and 1000 of these WDs initially hosted a companion whose final orbital parameters were calculated using the FATES code from Ronco et al. (2020) with stellar evolution tracks from MESA. We first examine the outcomes for the 1000 stars that initially hosted a companion.

4.1 Survival and initial companion characteristics

We simulated the orbital evolution of the 1000 substellar companions assuming the approximations of the convective turnover time from R96 and VL09. In the latter case, the resulting stellar tides are weaker, and more companions are therefore expected to survive.

The top panels of Fig. 3 show that the rate of surviving companions indeed depends on the assumed prescription for the convective turnover time. While more than a quarter of all substellar companions survived (26.2%) assuming the approximation by VL09, only 12.3% survived if the prescription for the turnover timescale from R96 was assumed. As more massive companions generate stronger tides, the median of the surviving companion mass distribution is shifted somewhat to smaller companion masses. Of the population of surviving companions, 93.5% would be classified as giant planets; i.e., they would have masses below ~13 M_J (Burrows et al. 1997) when weak tides are assumed, compared to 95.1% under strong tides.

Similarly, the initial semimajor axis has to be somewhat larger on average for companions to survive in the case of stronger tides (middle panels of Fig. 3). Assuming weak tides, the initial semimajor axis of some surviving companions is ~2.5 au (left), while the closest substellar companions that survive in the case of strong tides are initially located slightly farther away from the star (3 au).

Finally, companions with a larger initial eccentricity are more likely to be engulfed. This is because for a given orbital separation, the larger the eccentricity of the companion, the closer to the host star the companion is during periastron, which can lead to significantly increased tides. Therefore, surviving companions have smaller eccentricities on average. Again, this effect is more pronounced when assuming stronger tides (bottom panels of Fig. 3).

4.2 Characterizing WDs and their companions

Having inspected the relation between the initial distributions of companions that survive or are engulfed for both prescriptions of the convective turnover time, we now take a detailed look at the predicted final semimajor axis and eccentricity of the surviving companions and the properties of their host stars. Again, we present the results for weaker and stronger tides resulting from the different prescriptions for the convective turnover time.

In the top panels of Fig. 4, we show the distributions of the final semimajor axis and eccentricity of the surviving companions and compare them with the initial ones. Independently of the strength of the stellar tides, the semimajor axis of the surviving companions increases significantly, on average, due to the evolution of the host star into a WD. In addition, the distributions become much broader (top left panels of Fig. 4). This is easy to understand. For most surviving companions, mass loss dominates, and, consequently, the orbital separation generally increases. We note, however, that changes in semimajor axis depend not only on their initial value and the eccentricity, but also on the mass of the central star and the companion. In addition, for a few companions the semimajor axis increases only slightly, or even not at all, as the dominance of mass loss over stellar tides is much less pronounced. All these effects combined explain why the resulting final distributions are much broader than the initial one.

The strength of stellar tides also plays a role when comparing the initial and final eccentricities of the surviving companions. Eccentricities become significantly smaller due to stellar tides in the case of stronger tides. This effect mostly results from the fact that companions with initially intermediate eccentricities (-0.15-0.30) are reduced due to stellar tides. Companions with an initially very large eccentricity are either engulfed during the giant phases or have initial semimajor axes large enough that even during periastron passage stellar tides remain relatively inefficient (see Fig. 4, upper right panels). We therefore predict that a non-negligible fraction of the surviving companion population may retain eccentric orbits (especially under weak tides), which could contribute to the perturbation of smaller bodies, either delivering them into the WD or leading to their ejection (Frewen & Hansen 2014).

The masses of the predicted WD population with surviving companions are quite similar for stronger and weaker stellar tides. In both cases the distribution peaks at typical single-WD masses of ~0.6 M_⊙. This is largely because low-mass progenitor stars are more frequent and because the substellar companion occurrence rate peaks around 1.7 M_⊙. In both cases, very few massive WDs (~0.8 M_⊙) are predicted to host companions. These WDs must have formed from low-metallicity stars of -3-3.5 M_⊙ or even more massive (3.5-4M_⊙) high-metallicity stars (see Bauer et al. 2025, Fig. 3). We note that for this mass range the occurrence rate of companions is highly uncertain (see Wolthoff et al. 2022 for details).

In the bottom right panel of Fig. 4, we show the age distributions for predicted WDs with companions. In general, companions are more likely around young WDs because the age-metallicity relation implies a larger companion occurrence rate for younger stars. The total age of the companions hosts, however, peaks at around 5 Gyr, because younger stars in our potential progenitor sample had often not yet evolved into WDs.

Fig. 3 Survival of substellar companions as function of companion mass (top), initial semimajor axis (middle), and initial eccentricity (bottom). The left panels represent results when considering the tidal force model by VL09, while the right panels show results when using the tidal force model by R96. All distributions are normalized by the total number of stars (1000). The number of surviving companions is indicated in the figure for each model (blue histogram). The vertical lines correspond to the median value of each sample. The solid lines correspond to the cumulative distributions. In general, the initial distributions of surviving companions (orange) are shifted toward lower masses, larger semimajor axes, and smaller eccentricities compared to the initial values of all companions (blue). These effects are slightly more pronounced for stronger tides (right panels).

Fig. 4 Final and initial properties of systems that host a susbtellar companion. Each row shows a set of histograms comparing two WD samples: one modeled using VL09 (left) and one based on the R96 (right) approximation for stellar tides. All histograms are normalized to the total number of surviving companions for each simulation. The final semimajor axes are significantly larger than the initial ones for both tidal prescriptions (top left). For weak tides, the eccentricity distribution does not significantly change, while the distribution moves slightly toward smaller eccentricities for strong tides (top right). The bottom panels show the WD mass and age distributions for both approximations of tidal forces which appear to be rather similar.

4.3 Fractions of WDs with companions

The perhaps most important question we would like to answer through our population synthesis is how many WDs are expected to host substellar companions. In order to answer this question, we evolved a representative number of WD progenitors, independently of whether they host a companion, from their assigned birth time to the present day (using SSE). The survival and the final orbital parameters were determined using stellar evolution tracks from MESA. The resulting population of WDs allowed us to determine the relative number of WDs hosting substellar companions.

We find that in total 1.3% of the predicted WD population hosts a companion if weak tidal forces according to the turnover time approximation by VL09 are assumed. This number decreases to 0.6% if stronger tides resulting from the turnover times by R96 are assumed. Of the population of surviving companions, 93.5% would be classified as giant planets when weak tides are assumed, compared to 95.1% under strong tides. Therefore, the resulting global occurrence of giant planets around WDs would be approximately 1.2% for the weak-tide case and 0.6% for case of the strong tides. In other words, the vast majority of the predicted substellar companions are gas giant planets.

Given that the strength of stellar tides is uncertain and both prescriptions tested here could be poor descriptions of reality, as an additional test we simulated the evolution of the 1000 substellar companions without tidal forces to explore the maximum survival rate when the orbital evolution is governed solely by stellar mass loss. We find that the survival rate increases significantly (see third row of Table 1), yielding a global occurrence of 3.2%. This difference clearly shows how important stellar tides are in shaping some properties of the surviving companions population.

In addition to estimating the total fraction of WDs hosting substellar companions, we can evaluate how this fraction depends on the characteristics of the WD and its progenitor. In Fig. 5, we show the fraction of WDs with companions as a function of the WD mass and the total age of the system. For both assumed tidal forces, the predicted companion occurrence rate is largest for young (~1-6 Gyr) and relatively low-mass (~0.53-0.66 M_⊙) WDs. This result can be understood by the dependence of the companion occurrence rate on stellar mass (which peaks around ~1.7 M_⊙) and metallicity, which (according to our assumed metallicity-age relation) increases with age.

While the general tendencies seem to be independent of the efficiency of stellar tides, we note that the maximum occurrence rate slightly shifted toward older stars and smaller masses if strong tides were assumed (middle panel of Fig. 5). This small shift is related to the fact that for strong tides and large stellar masses only companions at very large semimajor axes survive (as also illustrated in Fig. 4, lower left panels).

Finally, in order to estimate the statistical uncertainty of the predicted fractions of WDs with companions, we generated 3000 populations of WDs (~ 100 000 WDs each) with potential companions, ignoring the possible effects of stellar evolution but varying the parameters in the occurrence-rate function (Eq. (7)) according to their uncertainties. We find that without taking into account the impact of stellar evolution (i.e., assuming that all companions survive), we can predict an occurrence rate of 4.4 ± 2.4% (1σ uncertainty). This value represents a strict upper limit. The relatively large uncertainty arises mainly from higher stellar masses (M_MS > 2.5 M_⊙), where the companion occurrence is highly uncertain (see Eq. (7)). This implies that the statistical error of the obtained fractions should generally be roughly 50%.

In summary, our simulations predict that the substellar companion occurrence rate around WDs, even if only weak tides are assumed, does not exceed -1.3 ± 0.7% and that the expected fraction is largest for low-mass WDs with total ages below ~6 Gyr. However, this result depends on several assumptions that we made in our population-synthesis simulations. In what follows, we investigate to what degree our results depend on these assumptions and what our results might imply for the observability of substellar companions around WDs.

Fig. 5 Fractions of WDs hosting substellar companions as function of total system age and WD Mass. The total fraction of WDs with substellar companions is 1.3% (assuming weak tides VL09, top left) and 0.6% for stronger tides according to R96. In both cases, companions are more likely around lower mass WDs and relatively young systems. In particular, for higher mass stars, stronger tides lead the radius expansion to dominate over orbital expansion, and therefore only planets around lower mass stars survive. The extension of the region with high probabilities is therefore very sensitive to the strengths of the tides. For example, further increasing the tides according to R96 by a factor of 8/3 would eliminate the remaining small high probability island in the right panel.

5 Discussion

In our base models presented in the previous section we found an upper limit of 1.3 ± 0.7% for WDs hosting substellar companions. The exact predicted fraction depends on the adopted prescriptions for tidal forces. The surviving companions are predicted to be located around intermediate-mass WDs with final semimajor axes covering a range of 3-24 au with a median of 11 au. These findings are based on the assumed initial companion distributions and the age-metallicity relation. In what follows, we address the reliability of our assumptions before briefly discussing the potential survival of engulfed companions and the observability of the predicted population of substellar objects around WDs.

5.1 Age-metallicity relation and star formation history

The substellar-companion occurrence rate depends significantly on the metallicity of the host star. We adopted an age-metallicity relation derived from Carrillo et al. (2023), who used simulations to predict a global age-metallicity relation. For most age ranges (from 10 up to 4 Gyr ago), they find a relatively smooth and well-known continuous trend: the older the star, the lower its metallicity. For stars younger than 4 Gyr the relation is essentially flat, clustering roughly around solar metallicity (with a relatively large scatter). This is consistent with classical approaches considering the chemical enrichment of the interstellar medium.

Alternative studies, however, have revealed that the Milky Way might have a more complex history. Ruiz-Lara et al. (2020) applied Gaia DR2 color-magnitude diagram (CMD) fitting for stars within 2 kpc of the Sun and found that the star formation rate has been slowly decreasing and perturbations to our galaxy from the Sagittarius dwarf galaxy triggered enhanced star formation rates, especially at three ages: 5.7, 1.9, and 1.0 Gyr. More recently, Gallart et al. (2024) and Fernández-Alvar et al. (2025) used improved CMD fitting to derive the star formation history and age-metallicity relation for different local samples. According to these works, the thin disk’s metallicity in the defined local regions was super-solar at the beginning (-8-10 Gyr ago), decreased until an age of ~3 Gyr, and then increased again to super-solar values (Fernández-Alvar et al. 2025, their Fig. 4).

In principle, these recent findings could directly affect the predicted fraction of substellar companions around WDs in the local volume. In order to estimate the potential impact of different star formation histories and age-metallicity relations, we performed two separate tests incorporating the age-metallicity relation (see right panel of Fig. 6) and star formation history from Fernández-Alvar et al. (2025). In both cases, we predicted the population of WDs with substellar companions ignoring the impact of stellar evolution (because of the tremendous computational costs of running a meaningful number of MESA simulations). The obtained results therefore have to be taken as strict upper limits and should be compared with the corresponding upper limits predicted by our default model (constant star formation rate and standard age-metallicity relation).

While changing the star formation history has very little impact on the predicted population of substellar companions around WDs (strict upper limit slightly increases from 4.4 ± 2.4 to 4.7 ± 2.6%), given the strong dependence of the companion occurrence rate on metallicity, a change in the age-metallicity relation affects the prediction significantly. In the top panels of Fig. 6, we compare the more local age-metallicity relation with that used in our population synthesis. Given that observational surveys will obviously target WDs close to the Sun, in this context the more local relation might provide more realistic predictions (right panel). Using the age-metallicity relation from Fernández-Alvar et al. (2025) increases the number of predicted companions by almost a factor of two (the strict upper limit - obtained by ignoring the fact that companions might be engulfed - increases from 4.4 ± 2.4 to 8.2 ± 4.4%). The properties of WDs hosting companions would also be affected, as considering a super-solar metallicity at very old ages results in the prediction of many more old WDs hosting substellar companions (see Fig. 6, bottom panels). In summary, due to variations in the age-metallicity relation in local volumes or the Milky Way, the fraction of WDs hosting planetary and brown dwarf companions may change by factors of -2. Also the characteristics of the WDs most likely hosting a companion significantly depends on the local age-metallicity relation.

Fig. 6 Age-iron abundance [Fe/H] relation for 243 617 stars generated by interpolating distributions obtained by Carrillo et al. (2023) (left panel) and for the estimate based on Fig. 4 from Fernández-Alvar et al. (2025) (right panel; assuming the same dispersion as our baseline age-metallicity relation). The age-metallicity relation plays an important role in our population synthesis as the occurrence of substellar companions depends strongly on stellar metallicity. The strict upper limit (assuming all companions survive stellar evolution) for the occurrence rate of companions to white dwarfs increases by a factor of -2 (from 4.4 ± 2.4 to 8.2 ± 4.4%), and old WDs are much more likely to host a companion when we assumed the age-metallicity relation from Fernández-Alvar et al. (2025).

5.2 Initial period distribution

The perhaps most fundamental assumption in our population synthesis is that the companion occurrence rates derived by Wolthoff et al. (2022) represent realistic approximations. Concerning the dependence of the occurrence rate on metallicity and stellar mass, this assumption appears to be well justified, as complementary studies broadly agree with the results obtained by Wolthoff et al. (2022). For example, Johnson et al. (2010) and Fulton et al. (2021) found a similar increase with metallic-ity and stellar mass up to ~2 M_⊙, and recent imaging surveys of FGK and A stars further confirm the trend that substellar objects are more frequent around more massive stars up to 1.5 M_⊙ (e.g., Vigan et al. 2017; Wagner et al. 2022).

With respect to the orbital period distribution of the companions, Wolthoff et al. (2022) obtained good fits to the observational data using a broken power-law or log-normal distribution. In general, the companion occurrence rate increases with period up to a peak (turnover) and then decreases for longer periods. The general shape of the companion period distribution is in agreement with previous findings for giant planets around lower mass main-sequence stars. However, the peak (or turnover) found by Wolthoff et al. (2022) is located at shorter periods (i.e., 718 ± 226 days for the broken power-law or 797 ± 455 days for the log-normal distribution) than derived from surveys of main-sequence stars that found the peak to be close to the snow line (Fernandes et al. 2019). This difference is somewhat surprising because, assuming the relationship between the peak of the occurrence rate and the snow line is a general trend, one would expect companions around more massive stars to be located farther away from the host star; instead, they seem to be found closer to it. Possible explanations for this discrepancy range from differences in the companion formation process, for example because the disk lifetimes are shorter around more massive stars (Ronco et al. 2024), to the impact of stellar evolution (where tidal evolution may have already played a role) or false positives (see also Wolthoff et al. 2022 for a discussion).

Indeed, it has previously been argued that pulsations of very massive and luminous giants can mimic the radial-velocity signal produced by companions (Saio et al. 2015). In particular, Döllinger & Hartmann (2021) identified several stellar properties associated with suspected radial velocity planet candidates around evolved stars; i.e., large radii, low metallicities, low effective temperatures, and low surface gravities. In particular, giants with R ≳ 21R_⊙ show an anomalously high incidence of radial velocity planet candidates, with orbital periods between 300 and 800 days. These last two criteria would compromise the validity of at least three of the 37 planets used by Wolthoff et al. (2022). In line with this, a recent study by Spaeth, Dane et al. (2025) demonstrates that one of these planet candidates (HIP 16335b) is in fact a false positive, because the radial velocity oscillations can be explained by intrinsic stellar processes.

If the period distribution by Wolthoff et al. (2022) is indeed affected by false positives or stellar evolution, our predictions would be directly affected, because we assumed their period distribution to generate a representative sample of potential WD progenitors with substellar companions. To explore how much a shift in the period distribution would affect our results, we moved the assumed period distribution to larger periods (simply by multiplying by a factor of 1500/800 = 1.875) to produce a peak at 1500 days, which is consistent with Fernandes et al. (2019). The results of the simulations are displayed in rows 4 and 5 in Table 1 and show an increase in the survival rate (f_surv increases by a factor of 2.2 for weak tides and 3.1 for strong tides). The final separation distribution follows the expected scaling, as the median increases by almost a factor of (1500/800)^2/32.

In summary, as expected, the predicted companion occurrence rate depends on the assumed period distribution, which remains uncertain. However, even by moving the period distribution to wider semimajor axes, we estimate that less than 3% of the current WD population hosts a substellar companion.

5.3 Eccentricity

In our population synthesis, we assumed the eccentricity distribution derived by Kane & Wittenmyer (2024). They used planets with masses exceeding that of Saturn from the NASA Exoplanet Archive (Akeson et al. 2013). Most of the host stars with known companions and determined eccentricity are solar or lower mass stars, i.e., smaller than those of WD progenitors. Moreover, a recent study by Alqasim et al. (2025) suggests that the planet eccentricity of transiting long-period giant planets depends on stellar metallicity, planet radius, and planet multiplicity. Such dependencies or a relation with any properties of the host star or the companions were not explored in detail by Kane & Wittenmyer (2024). A more complete and detailed eccentricity distribution for giant planets around WD progenitors has yet to be derived.

To very broadly explore the impact of the assumed eccentricity distribution, we simulated the survival of companions around WD progenitors assuming only circular orbits (sixth and seventh rows of Table 1). We found that the survival rate increases by more than a factor of 1.3 for both tidal prescriptions (close to a factor of 1.4 in the case of strong tides). Interestingly, under R96, the median of the final semimajor axis increases. Although more companions closer to the star are able to survive, at the same time companions that previously survived (assuming the eccentricity distribution from Kane & Wittenmyer 2024) are less affected by tides, and so their final semimajor axis is larger. The characteristics of the companion-hosting WDs are very similar to those predicted by the corresponding baseline simulations.

5.4 Impact of the mass-loss parameter

Another crucial factor affecting the results is the modeling of stellar mass loss. For the MESA simulations, we adopted massloss prescriptions from Reimers (1975), with an efficiency of η = 0.5 for the RGB, and the model from Bloecker (1995), with η = 0.02 for the AGB (Ventura et al. 2012). Although these mass-loss prescriptions are widely accepted and in general observationally rather well constrained, the efficiency of mass loss on the AGB is still somewhat uncertain (see Cinquegrana et al. 2022, their introduction). In general, a smaller value of η produces more thermal pulses and an extended AGB phase, while larger values truncate the late AGB phases much earlier. This clearly affects the survivability of companions.

To analyze how this mass-loss efficiency on the AGB impacts our results, we ran our simulations adopting a significantly larger efficiency parameter (η = 0.7) and found, as expected, that the survival rate increases, especially for stronger tides (see bottom two rows in Table 1). Noticeably, the median of the final semimajor axis becomes considerably larger under stronger tides (compared to the baseline model), because the star expands and pulsates less, and the surviving companions are therefore more affected by mass loss and much less by tides. Even with this change, we still predict that under 3% of the current WD population should host a substellar companion.

5.5 The potential impact of additional objects

One limitation in our predicted population is that we only considered single-companion systems; i.e., we ignored the potential impact of multiple companions. The related uncertainties are particularly important to discuss as multi-planetary systems could be more common than single-planet systems (Zhu et al. 2018; Zink et al. 2019). For instance, Mulders et al. (2018) reproduced Kepler data with simulated planet populations and found that roughly two thirds of all stars could host multi-planetary systems. However, it remains unclear to what extent multiplicity affects our predictions. First, multi-planetary systems seem to be less frequent around more massive stars. Second, how multiple planets affect survival has not been investigated in detail. Maldonado et al. (2021) and Maldonado et al. (2022) combined stellar evolution and N-Body simulations, but investigated the evolution during the WD phase rather than the survival during the giant branches. They found that the higher the number of planets is, the more likely instabilities in the planetary system are during the WD phase.

Furthermore, planets at low separations to WDs may result from a variety of dynamical processes involving the orbital perturbation induced by another companion, such as the eccentric Kozai-Lidov mechanism (Stephan et al. 2021; O’Connor et al. 2021) or planet-planet scattering (Maldonado et al. 2021). An additional companion located closer to the star might also truncate stellar evolution on the giant branches by generating common envelope evolution and thereby increase the chances of survival for the more distant companions (e.g., Lagos et al. 2021). A systematic population study including all these effects, however, is currently out of reach.

In general, it might also be that the survival of most companions would be little affected by multiplicity. Ronco et al. (2020) incorporated stellar tides, stellar evolution, and N-body simulations for two-planet systems. They find that the interactions between planets changed the final fate (engulfment or survival) of one of the companions in only roughly two per cent (2/99) of their simulations. This is because mean-motion resonances are required for this to happen. We therefore believe that our predictions are valid first-order estimates of the population of substellar companions around WDs.

5.6 Survival of engulfed companions

In this work, we considered that all substellar companions that enter the giant envelope are destroyed, although this may not always be the case. In fact, massive planets around less massive stars have been suggested to survive the common-envelope phase. If orbital energy is efficiently used to unbind the envelope, engulfed companions may spiral in toward the common center of mass, but can avoid being evaporated or disrupted at the Roche limit (Nelemans & Tauris 1998; Nordhaus & Spiegel 2013). However, this scenario most likely only applies to the most massive companions considered here. For instance, O’Connor et al. (2023) studied the potential survival of engulfed planets and found that for planets with masses of 1-10 M_J around stars with masses of 1.0-1.5M_⊙, survival was impossible.

This situation might change if several planets are engulfed. In fact, common envelope evolution can explain the existence of substellar companions very close to WDs, especially when invoking successive planet engulfment or contributions from additional energy sources such as recombination energy (Lagos et al. 2021; Chamandy et al. 2021; Merlov et al. 2021). We will investigate the number of companions that might survive engulfment in our population synthesis approach in a future work.

5.7 Observability implications

A key question that we did not answer in the present paper is how many of the predicted companions will be detectable with current and future observing facilities. A complete observational census of companions to WDs might allow us to constrain some of the uncertainties in our simulations, such as the strength of stellar tides, the mass-loss efficiency during the AGB phase, or maybe even the assumed initial companion distributions.

The potential to directly detect the predicted companions through imaging will depend on the age of the WD, the total age of the system, and the orbital separation, which converts to the required contrast and angular resolution of the chosen observing facility. Detecting the companions through an infrared excess in the spectral energy distribution will depend on the relative brightness of the WD and the companion. We will determine the detectability of the predicted companions and present the results in a future publication.

6 Conclusion

We generated a representative population of potential WD progenitor stars with substellar companions using initial distributions based on the results of radial-velocity surveys of 1-4 M_⊙ giant stars and a global age-metallicity relation. We then evolved these stars to the present day taking into account stellar and orbital evolution to predict the population of substellar companions around WDs. We find that under 3 ± 1.5% of WDs are expected to host a substellar companion, suggesting that substellar companions around WDs are intrinsically rare. Roughly ~95% of the predicted substellar companions are gas giant planets. The exact fraction of WDs with companions depends on initial distributions, stellar tides, and mass loss during the AGB.

To some extent, the properties of the surviving companion population are shaped by the adopted tidal prescription. The fraction of companions that survive the stellar evolution of their host star is roughly two times higher when adopting the weak-tides prescription of Villaver & Livio (2009) compared to the strong stellar tides suggested by Rasio et al. (1996). Stronger tides generally require companions to be initially located at somewhat larger semimajor axes and to orbit less massive stars in order to survive, which results in a distribution shifted toward slightly larger final semimajor axes and lower WD masses. The companion population is expected to be located at semimajor axes from around 3 to 24 au, with a median of 11 au, and is more likely to orbit lower mass white dwarfs (~0.53 M_⊙ to ~0.66M_⊙) with relatively low total system ages (~1-6 Gyr), where the occurrence reaches values of ≳3%.

Our general findings are largely independent of the assumed period and eccentricity distributions as well as the stellar massloss parameter. The local sample of white dwarfs close to the Sun might host more companions than we predict here, as the local age-metallicity relation implies a larger occurrence rate (≲8%). In a future article, we aim to derive concrete predictions concerning the observability of the predicted companions in order to inform ongoing surveys of substellar companions to white dwarfs.

Acknowledgements

AMS, MRS, DC, JP, CRJ, JV, MPR thank for support from FONDECYT (grant No. 1221059). MPR is partially supported by PICT-2021-I-INVI-00161 from ANPCyT, Argentina. FLV acknowledges support from ANID FONDECYT Postdoctoral Grant No. 3250464.

References

All Tables

All Figures

Fig. 1 Orbital evolution of a 1 MJ gas-giant planet orbiting a 2 M⊙ star with Z = 0.0187, comparing calculations performed with SSE (left), the MESA default test suite (middle), for which the AGB mass-loss efficiency is set to η = 0.7, and MESA assuming a more realistic mass-loss efficiency (η = 0.02; right). The orbits are assumed to be circular, and initial separations range from 2 to 4.5 au, with a step size of 0.5 au. The red filled area corresponds to the stellar radius, while purple and green lines denote the orbital separation of the engulfed and surviving planets, respectively. The insets show zoomed-in views of the thermal pulse phase of the AGB and highlight its critical role in planetary engulfment. SSE does not account for the thermal pulses (left). For a large mass-loss efficiency (middle), most planets survive because the orbit expansion caused by stellar mass loss dominates. In the most realistic scenario assuming a small efficiency (right), the star evolves through more thermal pulses and reaches a larger radius, which causes most planets to be engulfed.

In the text

Fig. 2 Survival of a 10 MJ gas-giant planet orbiting a 2 M⊙ star for different initial semimajor axes and eccentricities. Red stars indicate engulfment, while blue circles stand for survival of the planet. We assumed two different approximations for the convective turnover time, τconv. Following R96 (right panel), many more planets are engulfed than when assuming the prescription suggested by VL09 (left panel). In both cases, the survival of substellar objects depends on the eccentricity of the orbit, with more eccentric orbits leading more frequently to engulfment. Due to stronger tides, this dependence is more pronounced assuming the approximation by R96. For initial semimajor axes ≥ 8 au, the effects of tidal forces are negligible except for very large eccentricities.

In the text

Fig. 3 Survival of substellar companions as function of companion mass (top), initial semimajor axis (middle), and initial eccentricity (bottom). The left panels represent results when considering the tidal force model by VL09, while the right panels show results when using the tidal force model by R96. All distributions are normalized by the total number of stars (1000). The number of surviving companions is indicated in the figure for each model (blue histogram). The vertical lines correspond to the median value of each sample. The solid lines correspond to the cumulative distributions. In general, the initial distributions of surviving companions (orange) are shifted toward lower masses, larger semimajor axes, and smaller eccentricities compared to the initial values of all companions (blue). These effects are slightly more pronounced for stronger tides (right panels).

In the text

Fig. 4 Final and initial properties of systems that host a susbtellar companion. Each row shows a set of histograms comparing two WD samples: one modeled using VL09 (left) and one based on the R96 (right) approximation for stellar tides. All histograms are normalized to the total number of surviving companions for each simulation. The final semimajor axes are significantly larger than the initial ones for both tidal prescriptions (top left). For weak tides, the eccentricity distribution does not significantly change, while the distribution moves slightly toward smaller eccentricities for strong tides (top right). The bottom panels show the WD mass and age distributions for both approximations of tidal forces which appear to be rather similar.

In the text

Fig. 5 Fractions of WDs hosting substellar companions as function of total system age and WD Mass. The total fraction of WDs with substellar companions is 1.3% (assuming weak tides VL09, top left) and 0.6% for stronger tides according to R96. In both cases, companions are more likely around lower mass WDs and relatively young systems. In particular, for higher mass stars, stronger tides lead the radius expansion to dominate over orbital expansion, and therefore only planets around lower mass stars survive. The extension of the region with high probabilities is therefore very sensitive to the strengths of the tides. For example, further increasing the tides according to R96 by a factor of 8/3 would eliminate the remaining small high probability island in the right panel.

In the text

Fig. 6 Age-iron abundance [Fe/H] relation for 243 617 stars generated by interpolating distributions obtained by Carrillo et al. (2023) (left panel) and for the estimate based on Fig. 4 from Fernández-Alvar et al. (2025) (right panel; assuming the same dispersion as our baseline age-metallicity relation). The age-metallicity relation plays an important role in our population synthesis as the occurrence of substellar companions depends strongly on stellar metallicity. The strict upper limit (assuming all companions survive stellar evolution) for the occurrence rate of companions to white dwarfs increases by a factor of -2 (from 4.4 ± 2.4 to 8.2 ± 4.4%), and old WDs are much more likely to host a companion when we assumed the age-metallicity relation from Fernández-Alvar et al. (2025).

In the text