TTV constraints on additional planets in the WD 1856+534 system

ABSTRACT

WD 1856+534 b (or WD 1856 b for short) is the first known transiting planet candidate around a white dwarf star. WD 1856 b is about the size of Jupiter, has a mass less than about 12 Jupiter masses, and orbits at a distance of about 2  per cent of an astronomical unit. The formation and migration history of this object is still a mystery. Here, we present constraints on the presence of long-period companions (where we explored eccentricity, inclination, mass, and period for the possible companion) in the WD 1856+534 planetary system from transit timing variations. We show that existing transit observations can rule out planets with orbital periods less than about 500 d. With additional transit observations over the next decade, it will be possible to test whether WD 1856 also hosts additional long-period planets that could have perturbed WD 1856 b into its current close-in orbit.

1 INTRODUCTION

The past three decades have ushered in an exoplanet revolution with the discovery of over 5000 exoplanets (Akeson et al. 2013). Astronomers have learned that exoplanets are common in our galaxy (e.g. Kunimoto & Matthews 2020) and show a greater diversity in sizes (Barclay et al. 2013; Weiss et al. 2013; Zhou et al. 2017), compositions (Masuda 2014; Santerne et al. 2018), and architectures (Lissauer et al. 2011; Becker et al. 2015; Bourrier et al. 2021) than the planets in our own Solar system. However, almost all of the exoplanets known today orbit main-sequence stars that are still actively burning hydrogen into helium. Eventually, each of these main-sequence planet-hosting stars will undergo a dramatic transformation at the end of their lifetimes that will have major impacts on their orbiting planetary systems.

A white dwarf is the remnant core of a low-mass (⁠|$M_\star \lesssim 8 \, \mathrm{M}_{\odot }$|⁠) star at the end of its life, after it runs out of hydrogen fuel. When a star transitions from burning hydrogen on the main sequence into an inert white dwarf, it first puffs up into a red giant, more than a hundred times its initial size. As the red giant expands, it is expected to engulf and evaporate planetary material (except for the most massive objects; Passy, Mac Low & De Marco 2012) and carve out a zone without surviving planets in the inner region of the solar system (Nordhaus & Spiegel 2013). It will also likely destroy many small asteroids within a few au of the host star (Jura 2008).

As a result of the red giant phase of evolution, exoplanets orbiting close to white dwarf stars are expected to be rare, which is why the discovery of a close planetary-mass companion to the white dwarf WD 1856+534 was surprising (Vanderburg et al. 2020). WD 1856 is a white dwarf star that was found to host a transiting companion in data from the Transiting Exoplanet Survey Satellite (TESS) mission (Ricker et al. 2015). Vanderburg et al. (2020) detected periodic dips in the TESS light curve of WD 1856, indicating a close-in companion orbiting the star every 1.4 d called WD 1856 b. Follow-up observations from ground-based telescopes and the Spitzer space telescope showed that the companion has a size similar to Jupiter and a mass less than about 12 Jupiter masses. This planetary system is intriguing because the closely orbiting Jupiter-sized planet could not have existed in its current location while the host star was on the main sequence, or it would have been engulfed. Instead, the planet must have orbited farther away from the star, and migrated into its current location either during or after the star’s evolution into a white dwarf.

The evolution of WD 1856 b’s orbit is still uncertain and the subject of significant debate. There are two main mechanisms that could plausibly lead to the planet’s current orbit: (1) WD 1856 b either could have survived engulfment by the red giant and emerge at its current orbit intact (Chamandy et al. 2021; Lagos et al. 2021; Merlov, Bear & Soker 2021) or (2) it could have migrated via tidal circularization after being driven to a high-eccentricity orbit by other objects in the system [either other planets (Maldonado et al. 2021; O’Connor, Teyssandier & Lai 2022) or a pair of distant M-dwarf companions (Muñoz & Petrovich 2020; O’Connor, Liu & Lai 2021; Stephan, Naoz & Gaudi 2021)].

The first option is called ‘common envelope evolution’. In this scenario, the expanding star grows large enough to engulf a lower mass binary companion, and friction causes it to rapidly spiral inwards while depositing orbital energy into the giant star’s envelope. If the companion is able to deposit enough gravitational potential energy, then the envelope is ejected and the companion’s inward migration is halted. If there is insufficient gravitational potential energy to unbind the envelope, then the companion and the stars merge. WD 1856 b has a low mass and long orbital period compared to other post-common-envelope binaries (Vanderburg et al. 2020), which makes it challenging to explain its formation by this mechanism, but it could be possible with additional energy sources such as recombination energy in the envelope or energy deposited by other planets in the system (Chamandy et al. 2021; Lagos et al. 2021; Merlov et al. 2021).

The second scenario that could produce WD 1856 b’s current orbit is high-eccentricity migration, one of the main ways by which hot Jupiters around main-sequence stars are believed to have arrived in their current orbits (Dawson & Johnson 2018). In this scenario, gravitational interactions between WD 1856 b and other bodies in the system increase WD 1856 b’s orbital eccentricity. At the planet’s closest approach to the star, it is distorted by tidal forces, which dissipates a small amount of orbital energy and decreases the orbital semimajor axis. Over the course of ∼10⁸–10⁹ yr of tidal dissipation, the planet could slowly migrate into its present-day close-in orbit. One way the planet could reach such high eccentricities is via the von Zeipel–Kozai–Lidov (ZKL) mechanism (von Zeipel 1910; Kozai 1962; Lidov 1962), a secular dynamical process that causes oscillations in orbital inclination and eccentricity in the presence of distant, orbitally misaligned companions. WD 1856 is orbited by two M-dwarfs at a distance of ≈1500 au (Vanderburg et al. 2020) that could trigger ZKL oscillations (Muñoz & Petrovich 2020; O’Connor et al. 2021; Stephan et al. 2021). Another way WD 1856 b could reach high eccentricity is through interactions with other planets in the system, if they exist. In the presence of additional Jupiter-mass planets, WD 1856 b could plausibly reach high enough eccentricities to tidally circularize, but it is unknown whether WD 1856 b has had enough time to do this (Vanderburg et al. 2020; Maldonado et al. 2021; O’Connor et al. 2022).

Currently, astronomers are planning and performing follow-up observations of the WD 1856 system to better understand the system’s formation. The main goal of most of the follow-up observations planned so far is to detect spectroscopic features in the planet’s atmosphere, and use these features to constrain the planet’s mass. A precise measurement of the planet’s mass would help differentiate among the different formation scenarios, since common envelope evolution likely requires more massive companions than high-eccentricity migration (Xu et al. 2021). Recently, Alonso et al. (2021) and Xu et al. (2021) reported transmission spectroscopy observations to constrain the planet’s mass. These observations, carried out on 10-m-class ground-based telescopes, detected no statistically significant features in the planet’s transmission spectrum, and yielded a lower limit on the mass of about 0.84 Jupiter mass (Xu et al. 2021). Meanwhile, two programmes with the JWST have been approved for Cycle 1 observations (MacDonald et al. 2021; Vanderburg et al. 2021) to help constrain the planet’s mass and study its atmosphere.

In this work, we take a different approach to try to shed light on WD 1856 b’s present orbit: searching for additional planets in the system. If the system reached its current orbit via high-eccentricity migration excited by gravitational interactions with other planets, there should be other gas giant planets still orbiting in the system (unless they were ejected in the ensuing dynamical chaos). We search for planets by measuring the times of many transits of WD 1856 b to see whether the planet sometimes transits slightly earlier or later than expected. Such changes are called transit timing variations¹ (TTVs) and have been used to detect previously unknown planets (e.g. Ballard et al. 2011) and measure the masses of known planets (Carter et al. 2012). TTV measurements have also been used to detect planet candidates around eclipsing white dwarf/M-dwarf binary stars (e.g. Marsh et al. 2014), but it is challenging to interpret these variations because magnetic fields in the M-dwarf companions (Applegate 1992) often seem to drive changes in the binary orbit that can mimic orbital variations caused by additional planets (Bours et al. 2016; Pulley et al. 2022). Since WD 1856 b is much lower in mass than the M-dwarf companions that often show orbital variations, it may be possible to circumvent this problem, although any detected signals must still be subjected to the highest scrutiny for these reasons.

Our paper is organized as follows: In Section 2, we describe how we obtained the transit timing data for WD 1856 b. In Section 3, we discuss our process for modelling the TTVs and placing limits on the presence of planets. In Section 4, we summarize the results of our analysis, and in Section 5 we discuss its implications and conclude.

2 OBSERVATIONS

For our analysis in this work, we use a combination of previously published transit times and newly measured transit times. The transit times we used in this paper are listed in Table 1 and their deviations from a perfectly periodic model are shown in Fig. 2. In this section, we document the sources of the transit times we use in our analysis.

2.1 Literature observations

Since its discovery, WD 1856 b has been observed by several different groups. We collected observations from the literature from three different papers in particular.

First, we used transit observations collected by Vanderburg et al. (2020). They observed WD 1856 b with a number of different telescopes, including TESS, the Gran Telescopio Canarias (GTC) on the island of La Palma, Spain, the Telescopio Carlos Sánchez (TCS) on the island of Tenerife, Spain, the Spitzer Space Telescope, and three privately owned telescopes in Arizona, USA: one at Raemor Vista Observatory (RVO), one at Junk Bond Observatory (JBO), and one at Hereford Arizona Observatory (HAO). We used all of these observations except for the TESS data, which have much lower signal to noise per transit than the other observations, and the TCS observation, which was conducted simultaneously with the GTC observation, but with significantly lower signal to noise. Vanderburg et al. (2020) did not report individual transit times for their observations, so we used the transit times derived by Alonso et al. (2021) for the GTC and Spitzer data, and we report transit times for the Raemor Vista, Junk Bond, and Hereford observations following the procedure described in Section 2.2.

We also used transit times measured and reported by Alonso et al. (2021). They observed transits of WD 1856 b on several occasions with the GTC and derived transit times. Their observations were conducted using two different instruments and in different photometric bandpasses to support their primary science objective of measuring WD 1856 b’s transmission spectrum. In particular, Alonso et al. (2021) collected one photometric observation with the GTC’s OSIRIS instrument on 2020 March 2, one infrared observation with the GTC’s EMIR instrument on 2020 June 13, and three observations with OSIRIS in a spectroscopic mode on 2020 May 3, 2020 May 10, and 2020 May 20. Alonso et al. (2021) also derived transit times from archival transit observations collected by Vanderburg et al. (2020) from the GTC and Spitzer. In total, they presented seven precisely measured transit times, all of which have uncertainties smaller than 6 s, and five of which have uncertainties smaller than 1 s.

We used two observations of WD 1856 b collected by Mallonn (2020) using the 1.2-m STELLar Activity (STELLA) telescope on the island of Tenerife, Spain. These observations were conducted using the Wide Field STELLA Imaging Photometer imager on the 1.2-m STELLA in Sloan i band. We use the transit times measured by Mallonn (2020) in our analysis, which have uncertainties of 3–5 s.

Finally, we used transit times measured by Xu et al. (2021) from the Gemini North telescope on Mauna Kea, Hawaii, USA. Xu et al. (2021) observed 10 transits of WD 1856 b using the GMOS multi-object spectrograph, also with the primary science goal of measuring the planet’s transmission spectrum. These observations were taken on June 15, July 2, July 30, August 23, August 30, September 2, September 9, September 23, October 3, and October 17, all in 2020. Xu et al. (2021) derived transit parameters, including times, for each of the individual observations. For their transmission spectrum, Xu et al. (2021) excluded data from two transits that were taken in worse weather conditions, but we found that the transit times derived from these observations were still useful, albeit with somewhat larger uncertainties than the rest of the observations. All 10 Gemini observations have transit time uncertainties less than 1 s, and all but the two excluded by Xu et al. (2021) have uncertainties smaller than 0.25 s.

2.2 New observations from Hereford Arizona Observatory

We obtained our own observations of WD 1856 b using the 16-arcsec telescope at HAO, which is operated by Bruce Gary. We observed WD 1856 b on 21 nights between 2020 May 18 and 2020 December 21 (in addition to the two HAO observations from October 2019 published by Vanderburg et al. 2020). The data collection and photometric analysis to produce light curves were performed as described by Rappaport et al. (2016). From the light curves, we derived transit times by fitting a simple model light curve to our observations. The transit dip is defined by five straight line segments: two steeply sloping segments at ingress/egress, two more gradual slopes nearing the transit centre, and a flat bottom.² Six parameters define the transit model. Three of these parameters define the shape of the transit: the relative duration and depth of the steep segments and the more shallow segments, and the duration of the flat bottom. Three other parameters vary the overall depth, duration, and mid-transit time. After creating the piecewise linear transit model, we smooth that model with a running-average (boxcar) smoothing over five time bins, each 49 s long. While the resulting transit shape model does not perfectly describe the shape of high signal-to-noise observations like those from the GTC, it is certainly sufficient to describe the much noisier observations from the HAO 16-arcsec telescope.

The high signal-to-noise GTC light curve from Vanderburg et al. (2020) was used for establishing the shape parameters of the transits (the relative depth and duration of the ingress/egress and more gradual slopes, and the duration of the flat bottom). Once we optimized the shape parameters based on the GTC light curve, we measured transit times by varying the overall depth, overall transit duration, and mid-point time of the individual transits. We performed the optimization using a spreadsheet that calculated the sum of χ² between the observed and model light curves for many different parameter values and identifying the value that minimized χ².

We found no evidence for variations in any of the three parameters we optimized for each HAO observation. Unsurprisingly, the uncertainties on the transit times measured from the 16-arcsec HAO telescope are significantly (1–2 orders of magnitude) larger than the best uncertainties we measured with 10-m-class telescopes. However, the HAO observations could still be important because they have better time sampling, and can rule out large variations in the transit times when the oversubscribed 10-m telescopes are unable to observe.

3 ANALYSIS

3.1 TTV modelling

We analysed our transit times to assess the evidence for the presence of planets using a forward modelling approach. In brief, we used the ttvfast software (Deck et al. 2014) to calculate the expected transit times for WD 1856 b in the presence of many different hypothetical planetary companions. We then compared the expected transit times for each planetary companion to the transit times we actually observed to determine the likelihood of each set of calculated transit times, and therefore identify which of the hypothetical orbits for planetary companions are preferred, plausible, or unlikely.

The heart of our analysis is a likelihood function that, given a list of measured transit times and input planetary and orbital parameters for a hypothetical second planet in the WD 1856 system, returns the relative Bayesian likelihood of that particular system configuration. Our likelihood function, which we implemented in the python programming language, performs the following steps:

• First, the function identifies the ‘transit number’, or the number of planetary orbits that had taken place since a reference transit, of each measured transit time. For our reference orbit, we took the transit observed by Vanderburg et al. (2020) with the GTC on 2019 October 22. We calculated the transit number of each time by calculating the time difference between each measured transit and the GTC observation, dividing by the orbital period reported by Vanderburg et al. (2020), and rounding to the nearest integer. The transit number for each observation we used is listed in Table 1.

• We then calculated the expected transit times for WD 1856 b for a system given arbitrary input parameters using the ttvfast code (Deck et al. 2014). We initialized the model using the orbital parameters for WD 1856 b with an orbital period and time of transit as reported by Vanderburg et al. (2020), an inclination and longitude of the ascending node of 0° (defining the reference plane of the system to be the plane of WD 1856 b’s orbit), and eccentricity of 0. We chose to only explore circular orbits for WD 1856 b itself because the tidal circularization time-scale for roughly Jupiter-mass companions (2 Myr) is much less than the age of the system (6 Gyr; Vanderburg et al. 2020), and because adding two additional variables to our parameter study would be prohibitively expensive. We define the reference frame for our calculations such that WD 1856 b’s inclination is 90° exactly, and we measure the inclination of the companion with respect to this value. We note that this value is slightly different from the actual inclination of WD 1856 b’s orbit from the plane of the sky. For these experiments, we assumed that WD 1856 b has a mass of 10 M_J, but because the planet’s mass is much less than that of the star, this choice has little effect on our results. We allowed the planetary parameters of the hypothetical second planet in the system to vary. We experimented with different time-steps for the numerical integration and found that using a step of 0.02 d (1.5 per cent the orbital period of the innermost planet WD 1856 b) ensured that errors in the numerical integration were at least an order of magnitude smaller than the smallest uncertainties in our transit times. In each call to the likelihood function, we calculate the transit times of WD 1856 b over the course of 600 d, longer than the baseline (431 d) of our observations.

  Because ttvfast was originally designed to quickly model TTVs of low-mass planets near mean motion resonances, it makes two approximations that are not well justified for our experiments with the WD 1856 system. First, when ttvfast converts from orbital elements to positions and velocities, it assumes that the planets are massless. Since our assumed mass for WD 1856 b of 10 M_J is much smaller than the star’s mass (about 1.5  per cent) but not completely negligible, this assumption causes the orbital period of long-period companions in the ttvfast simulations to differ slightly from their expected values. We worked around this issue by initializing the ttvfast simulations with a planetary mass for WD 1856 b of 0, and a stellar mass equal to the sum of its measured value (0.576 M_⊙; Xu et al. 2021) and our 10 M_J assumed value for WD 1856 b. The second approximation made by ttvfast is that it ignores the component of TTVs caused by the Rømer delay, which is also known as the light traveltime effect (Deck, private communication). In most known multiplanet systems, timing variations due to the Rømer delay are negligible compared to TTVs caused by orbital variations, but in the case of long-period companions orbiting a close inner pair like in the WD 1856 system, this effect can dominate (see, for example, fig. 7 of Rappaport et al. 2013). We therefore performed our own analytical calculation of the Rømer delay (following equations 6 and 7 in Rappaport et al. 2013) and added the resulting curve to models returned by ttvfast.

• With these transit times from ttvfast (and corrected for the Rømer delay) in hand, we performed some post-processing to adjust for slight differences between the orbital period we used to calculate the model transit times and the true orbital period of WD 1856 b. In systems with TTVs, the true mean orbital period of each planet can differ significantly from the instantaneous orbital periods that might be measured by observations over a short time baseline. To account for this, it is necessary to allow the true orbital period and transit time of the transiting planets to vary, often as free parameters in a Markov Chain Monte Carlo analysis. However, this would require adding two additional free parameters to our analysis, so we take a different approach to account for this effect. We instead calculated a weighted linear least-squares fit between the transit times and transit number for both our observed transits and the model transit times from ttvfast where the weights were set equal to the inverse square of the uncertainties on the measured transit times. We then subtracted the best-fitting line from both the observed and model transit times, so we could compare the residuals – i.e. the TTVs themselves – when calculating the likelihood of each model. Optimizing the mean orbital period and time of transit within each model calculation using this polynomial fit allows our model to efficiently tweak the model’s assumed orbital period to match the observations without explicitly adding it as a free parameter.

• Finally, we calculated the likelihood of each model. We used a simple χ² likelihood function where the likelihood |$\mathcal {L}$| is given by

  $$\begin{eqnarray} \mathcal {L} = \mathrm{ e}^{-\chi ^2/2} , \end{eqnarray}$$

  (1)

  and

  $$\begin{eqnarray} \chi ^2 = \sum _i{\frac{(y_i-m_i)^2}{\sigma _i^2}} , \end{eqnarray}$$

  (2)

  where e is the base of the natural logarithm, y_i are the individual TTVs (after subtracting the best-fitting linear model), σ_i are the uncertainties on each transit time, and m_i are the model-predicted TTVs (also after subtracting the best-fitting linear model).

3.2 Parameter study

With our likelihood function, we have the tool we need to assess the evidence for additional planets in the WD 1856 system in hand; we simply need to calculate the likelihood of many different orbital configurations for a hypothetical second planet in the system. Often, astronomers use sampling techniques like Markov Chain Monte Carlo to explore parameter spaces like this one, but because it is computationally expensive to calculate TTVs using numerical integration, and the space of possible companion orbits is very broad, we opt instead to explore slices in parameter space using grid searches.

For our grid searches, we focused primarily on the planetary parameters that would likely change WD 1856 b’s TTVs the most. In particular, we focused on the hypothetical second planet’s mass, period, eccentricity, inclination, mean anomaly at the start of the integration, and argument of periastron. We treat two of these parameters, the mean anomaly and argument of periastron of the second planet’s orbit, as nuisance parameters, and present our results in terms of the mass, period, eccentricity, and inclination.

We created a three-dimensional grid of different values of mass, period, and mean anomaly, and calculated the likelihood of a hypothetical companion to WD 1856 b with the parameters at each grid point.The grid had 300 mass points, 100 period points, and 100 mean anomaly points. We collapsed the grid along the mean anomaly axis by selecting the values for these nuisance parameters that maximized the likelihood, leaving grids showing only the most important physical parameters: mass/period and eccentricity/inclination. As a point of comparison, we calculated the likelihood of a single-planet model (i.e. no TTVs). Fig. 4 shows our results for a grid of outer companion orbits with different companion masses and orbital periods, with eccentricity fixed at 0, and the inclination of the outer planet fixed at 90 (perfectly coplanar³ with the WD 1856 b). The colour represents the ratio between the maximum calculated likelihood within each grid point and the likelihood of a scenario with no additional planets (that is, no TTVs). In this case, we are able to rule out companions with large masses and short orbital periods in the purple region of the plot, while companions with low mass and long orbital periods are permitted, and in some cases, even weakly favoured (although not significantly; see Section 4.1).

To explore the relationships between the mass, orbital period, and two additional parameters (eccentricity and inclination of the companion’s orbit), we generated more grids slicing through a six-dimensional likelihood space. These will be available upon publication. We repeated our calculations to produce grids of the relative likelihood of different values of mass/period and eccentricity/inclination while fixing the other parameters to different values. In cases where we explored circular orbits only, we used three-dimensional grids as described in the previous paragraph, but when we explored non-zero values for eccentricity, we calculated the likelihood function along a four-dimensional grid that also included the argument of periastron for the outer companion’s orbit.⁴ In cases where eccentricity was non-zero, we collapsed our results into a two-dimensional grid by maximizing the likelihood over both the mean anomaly and argument of periastron. These additional grid calculations give us a better understanding of the impact the other variables had on the final likelihood and the effect different combinations of these parameters have on which planetary orbits are plausible and which can be ruled out.

We discuss the results of these calculations in Section 4.2. We performed these computationally intensive calculations using the MIT SuperCloud computer cluster (Reuther et al. 2018). Additionally, we attempted to reproduce the results from our numerical experiments using ttvfast with an analytical approach. We used the expressions for the physical and Rømer delays from Rappaport et al. (2013) to estimate TTV amplitudes over the baseline of our observations for similar grids of parameters. While the analytical expressions from Rappaport et al. (2013) are approximate and ignore higher order terms, they are very computationally efficient to calculate. In general, we found good agreement between the numerical and analytical approaches for relatively low mutual inclinations and eccentricities, but at higher inclinations and eccentricities the analytical approximations were less accurate. The analytical results gave good qualitative agreement with the numerical calculations, increasing our confidence in these results.

4 RESULTS

4.1 No evidence for additional planets orbiting WD 1856

The first result of our analysis is that we found no evidence for TTVs, and therefore no evidence for additional planets orbiting WD 1856 b, at least to the sensitivity limits of our data set. We identified the model with the highest likelihood from our entire parameter study and compared it to the likelihood of a model without a second planet. We consider parameter space ruled out when the likelihood is a factor of 10⁴ times worse than a 1 planet model. Although the likelihood of the best-fitting model from our entire study was higher than the likelihood of a single-planet model, this could be because a two-planet model has significantly more degrees of freedom than a single-planet model. To account for the increased model flexibility, we calculated the Bayesian information criterion (BIC) to determine whether the higher likelihood of our best two-planet model than a single-planet model is strong evidence for the existence of another planet. The BIC is an approximate method for comparing two models; the model with the lowest BIC is favoured, and one model is only strongly preferred over another when the difference in BIC (or ΔBIC) is more than 10. We found for the single-planet scenario that the BIC was 64.2, while for the best two-planet model the BIC was 78.1. Therefore, the single-planet model is strongly favoured over even our best two-planet model, so we find no convincing evidence of another planet in the system. We note that the parameters for the best two-planet model were 2 Jupiter masses, 50 d, 0.284 eccentricity, 70.5° inclination, 301.2° argument of periastron, and 88.2° mean anomaly.

In the absence of TTVs from additional planets, we can measure the best-fitting orbital period and time of transit for WD 1856 b. We performed a linear fit to the transit times (allowing for some excess noise in the transit time measurements) using a Differential Evolution Markov Chain Monte Carlo algorithm (Ter Braak 2006). We measured an orbital period of 1.407 939 217 ± 0.000 000 016 d and a time of mid-transit of 2459038.4358981 ± 0.00000114, in BJD, chosen at the epoch that minimizes the covariance between the period and transit time. We summarize this information in Table 2 as well. This ephemeris has high enough precision to predict transit times with an uncertainty better than 40 s for the next century.

4.2 Limits on the presence of additional planets orbiting WD 1856

Although we find no affirmative evidence for additional planets in the WD 1856 b system with the current data, it is still possible that such planets could exist but that the current data are not sensitive to them. To determine the parameter space in which planetary companions have and have not been solidly excluded based on the current data, we constructed portraits of the computed model likelihood as it depends on the orbital parameters of the additional planet. We examined the likelihood for various sets of orbital parameters (the companion’s orbital period P_c, eccentricity e_c, and inclination i_c) and physical parameters (companion planet mass m_c) by constructing parameter grids and computing the model likelihood for a discrete set of companion parameters.

In the portraits that follow, in general, regions of low probability (which tend to have shorter orbital period and larger companion planet mass) are ruled out because they would produce TTVs significantly larger than the observed values. On the other hand, the current data do not allow us to rule out companions in all of the studied parameter space; companions whose likelihoods are commensurate with (or even slightly higher than) the raw likelihood of the single-planet model cannot be excluded with the present data, even though the BIC indicates that the data quality does not support the added complexity of the two-planet solution generally.

4.2.1 The effect of orbital eccentricity

In Fig. 5, we show the result of one such analysis, where the likelihood is computed as a function of m_c and P_c for four slices of companion planet eccentricity e_c. The colour represents the ratio between the likelihood of a model at each grid point in mass/period and the likelihood of a model with no additional planets. The general trend persists that larger companion planet masses and shorter companion orbital periods (which together would tend to exert larger TTVs on WD 1856 b) are disfavoured. Within the most favourable area of parameter space (lower m_c and higher P_c), an additional band of apparent increased likelihood emerges for the planetary companions that induce small but non-zero TTVs to match the fluctuations in the data. Higher companion eccentricity results in a smaller range of possible orbital periods, as a closer periastron results in more significant planet–planet interactions (and larger TTVs).

In addition to those general trends with eccentricity, a secondary structure in likelihood emerges, which takes the form of horizontal spikes of allowed regions at particular orbital periods, and in some cases extending to high masses. We think that these spikes are caused by a combination of poor sampling of transit times at those periods, and occasional outliers in our transit timing data. In regard to sampling, these spikes appear at difficult periods for us to measure. We often see them at periods close to 1 or 1/2 yr since we are unable to observe WD 1856 b when its location in the sky is too close to our Sun, and the star is only overhead during daytime. This causes gaps in our observations each year, which prevent us from ruling out large timing variations at those times. We find that the spikes in the likelihood graphs are more prominent at high eccentricity because higher eccentricity introduces faster changes into the TTV curve, so it is possible that one or two outlying points can strongly affect the likelihood and favour a fast change to the TTV curve.

4.2.2 The effect of orbital inclination

In Fig. 6, we show the result of another analysis, where the likelihood is computed as a function of m_c and P_c for four slices of companion planet inclination i_c. The general trend persists that larger companion planet masses and shorter companion orbital periods are disfavoured in these inclination slices. In Fig. 6, we label the inclination of the companion with respect to the plane of the sky – that is, face-on orbits have an inclination of 0. Since in our reference frame the inner planet, WD 1856 b, is transiting with an inclination of 90°, and we set the longitudes of ascending nodes for the two planets to be equal, the mutual inclination between the two orbits is 90 − i_c.

In general, we find that we are able to rule out lower mass and shorter period companion planets when the mutual inclination between the planets is low. The data allow highly inclined companions with relatively large masses and short orbital periods. We are able to rule out more parameter space when the companion’s inclination is near 90° (low mutual inclination) because of the Rømer delay, which is proportional to sin(i_c). For inclinations near 0, we see additional effects at shorter periods. These effects are caused by the physical delay, or actual changes in the orbit of the potential planet. Unlike the Rømer delay, the amplitude of the physical delay is actually larger at high mutual inclinations (see equations 8–10 of Rappaport et al. 2013).

4.2.3 The effect of orbital period on an eccentricity/inclination slice

In Fig. 7, we show the result of an analysis wherein the likelihood is computed as a function of the e_c and i_c for four slices of companion planet period in days P_c. In this analysis, the mass of the planet in our parameter search was set to 2 M_J. In general, as we increase the period of our parameter slices we are able to rule out less parameter space.

At 50-d periods, we are able to completely rule out companions in a significant fraction of parameter space. This causes a particular pattern to emerge in the likelihood slice, where only very particular combinations of eccentricity and inclination can produce TTVs as small as we measure. At short periods, the TTV curves are dominated by the physical delay, which produces large TTVs at high mutual inclination and high eccentricity. We also see the contribution of the Rømer delay at low mutual inclination, but this effect is subdominant (see fig. 7 of Rappaport et al. 2013).

At 200-d periods, we see a crossover to where the Rømer delay dominates over the physical delay. This fundamentally changes the pattern of which eccentricities and inclinations are allowed or ruled out. Even though the Rømer delay is strongest at low mutual inclination, we find that we cannot rule out this space when the eccentricity of the planet is low. We can, however, rule out parameter space with high eccentricity and low mutual inclination, because as the eccentricity increases, the amplitude of both the Rømer delay and the physical delay increase. At even longer periods, ranging from 400 to 800 d, we are able to rule out relatively little parameter space for planets of this mass (2 M_J).

4.2.4 The effect of mass on an eccentricity/inclination slice

Finally, in Fig. 8, we show the result of an analysis where we calculate likelihood as a function of eccentricity and inclination of the companion’s orbit, while holding the orbital period fixed at 130 d, and using various values of the planet’s mass. Because the mass of the planet linearly scales both the TTV amplitude from the Rømer delay and the physical delay, the pattern we see (which parts of parameter space are ruled out and which are allowed) in the likelihood slices does not change significantly at the different masses. Instead, changing the mass mostly affects the ‘stretch’ of the images and enhances or decreases the contrast between the regions. In general, we find that for low planet masses, m_c, the companion can exist in almost any orbit, and larger m_c rule out more parameter space with these specific configurations. Near 90°, companions are largely ruled out due to the Rømer delay, since that was where it has the strongest effects. On the other hand, near 0° we rule out parameter space because of the physical delay, which is due to the prominence of this effect at lower inclinations (higher mutual inclinations). At middling inclinations around 45°, both of these effects are attenuated, so we can rule out relatively few companion orbits in this region of parameter space.

5 DISCUSSION AND CONCLUSIONS

5.1 Implications for WD 1856 b’s formation

The main goal of our study was to determine whether additional planets, beyond WD 1856 b, reside in the WD 1856 system, and subsequently to determine whether they may have had a role in WD 1856 b’s migration. If additional planets orbit WD 1856 at farther distances from the star, they could plausibly have acted as dynamical perturbers to send WD 1856 b into its current short-period orbit. However, we found no evidence for such planets in the transit time measurements that we analysed. Although there were some combinations of parameters that produced TTV curves that yielded slightly higher likelihoods than a model with no additional planets, these models were not favoured in a comparison of the BIC for the two scenarios. We therefore find no evidence supporting the hypothesis that WD 1856 b migrated to its current orbit as the result of dynamical interactions with additional planets in the system.

However, even though we find no evidence for additional planets orbiting WD 1856, this does not rule out dynamical interactions with other planets as the cause of WD 1856 b’s inward migration. Our data were collected over the course of just two observing seasons, and even in the most favourable circumstances, we can only rule out planetary-mass companions at periods less than about 1200 d (or 4 yr). Any planets massive enough to dynamically perturb WD 1856 b into a high-eccentricity orbit must have originally orbited at distances greater than 2–4 au in order to survive post-main-sequence evolution, and therefore have orbital periods of at least a decade (Nordhaus & Spiegel 2013). Therefore, even if WD 1856 b did migrate to its current location as the result of dynamical interactions with additional planets in the system, we would not expect to detect them given the current limits of our observations. If we find no evidence for additional planets with transit timing measurements, there may be other ways to detect other planets in the WD 1856 b system, including astrometry (Sanderson, Bonsor & Mustill 2022), direct imaging (Xu et al. 2015), and thermal emission (Limbach et al. 2022).

While we currently find no evidence of additional planets in the WD 1856 system, more data could tell a different story. Additional transit observations with large telescopes would allow updates to the likelihood maps presented in this work and increase our sensitivity to long-period planets. Occasional observations over the course of the next 5–10 yr would likely extend our sensitivity to planetary companions orbiting where we expect dynamical perturbers of WD 1856 b might reside. So, while our findings can eliminate some parameter spaces, more transit data are needed to fully understand the geometry and population of the WD 1856 system and identify any possible neighbours that it has that could explain the movement of WD 1856 b.

5.2 The value of transit times from small telescopes

One novel component of our work is that we combine transit times measured from small (≲1-m diameter) telescopes with those measured from large (≳8-m diameter), professional observatories. Although the observations from large, professional telescopes yield much more precise transit time measurements than smaller observatories, it is more difficult to secure observing time on these facilities, and relatively few observations are possible. Given these trade-offs, we investigated whether data from the small telescopes significantly improved our constraints compared to only using observations from large, professional telescopes.

To test how much value is added by the inclusion of transit times from small, privately owned telescopes, we performed a comparison of how much parameter space we were able to rule out when using our full transit timing data set as opposed to only the sparse, but precise, timing observations from 8–10 m-class telescopes. We repeated our analysis for two likelihood slices in mass/period space (with coplanar orbits, and eccentricity fixed at 0 and 0.5) using only data from 8–10 m-class telescopes, and compared them with our nominal analyses using our full data set. Fig. 9 summarizes the differences between the two analyses. The figure highlights regions of parameter space where the addition of data from small telescopes either marginally increases or marginally decreases the likelihood ratio between a multiplanet model and a single-planet model across a threshold of 10⁻⁴. In general, the constraints obtained in the two analyses are very similar, with only slight differences on the margins (where the small telescopes may help provide slightly stronger constraints for planets with eccentric orbits). This suggests that the Gemini and GTC data dominate the combined data set, and transit observations from small telescopes do not contribute significantly.

We also compared the constraints we were able to obtain when combining the data set with only the less precise transit times from HAO. We repeated our analysis for one likelihood slice in mass/period space, with coplanar circular orbits for the companion. Fig. 10 shows the results of this analysis. Qualitatively, the parameter space ruled out by the HAO observations is similar to that ruled out by the full data set, with strong constraints on the presence of short-period massive companions and weaker constraints for low-mass companions at longer orbital periods. However, the lower precision on the transit times measured by HAO shifts the constraints towards higher masses, and in the best cases HAO observations can only rule out masses larger than about 30 M_J. Evidently, HAO data alone are not competitive with the sparse but precise timing measurements from 8–10 m-class telescopes.

5.3 Prospects for future observations

In this work, we demonstrated that it is possible to measure transit times of WD 1856 b with sufficient precision to detect (or rule out) distant super-Jovian planets in the system. Currently, our observational baseline is too short to detect planets with the expected orbital periods for planets that would have survived WD 1856’s red giant phase, but future observations will increase our sensitivity to long-period planets. It would be highly beneficial to collect more transits from the WD 1856 system over the coming years to extend our sensitivity to detect or rule out longer period planets.

Some highly precise transit observations are already planned; the JWST will observe multiple transits of WD 1856 b during its first observing cycle (MacDonald et al. 2021; Vanderburg et al. 2021). Although these observations are aimed at studying the planet’s atmosphere, they will likely yield highly precise transit time measurements, with uncertainties similar to the GTC and Gemini observations taken from the ground. The JWST data will extend the transit timing baseline with large, professional telescopes to more than 3 yr. Our current observational baseline of about 300 d gives sensitivity to planets at orbital periods up to 1200 d, so increasing our baseline to over 1300 d could give us sensitivity to planets with periods of up to a decade. Over that time, it likely will not be necessary to maintain the same high cadence of transit observations we present here, since we already strongly constrain short-period planets. Instead, we may be able to increase the spacing of observations exponentially over time as we probe longer orbital periods.

6 CONCLUSIONS

In this work, we have presented the first constraints on additional planets orbiting WD 1856+534 from transit timing measurements. Our conclusions are summarized as follows:

• We find no compelling evidence for additional planets orbiting WD 1856. Our observations are sensitive to two Jupiter-mass planets in orbits out to 500-d periods and 10 Jupiter-mass planets in periods as long as 1000 d.

• We are most sensitive to planets with short orbital periods and high masses, as expected from analytical arguments (e.g. Rappaport et al. 2013). Our constraints depend in more complex ways on the orbital eccentricity and inclination of hypothetical companion planets.

• We find that our constraints are dominated by highly precise transit timing measurements from 8–10 m-class telescopes, and even relatively few and sparsely sampled observations provide much stronger constraints than larger numbers of more frequent but less precise transit measurements from small telescopes.

• Our constraints on the existence of additional planets do not yet extend to long enough orbital periods to inform theories of WD 1856 b’s orbit. Future observations with 8–10 m-class telescopes over the next decade should achieve the sensitivity needed to detect long-period planets that might have perturbed WD 1856 b into its current orbit.

While we cannot yet firmly constrain the process by which WD 1856 b came to its current orbit, our work shows that there is a path forward towards addressing this question. Although our best data so far indicate that WD 1856 b may be lonely, we will have to wait for more transits to determine whether it has any friends that could explain its migration.

SUPPORTING INFORMATION

Supplementary data are available at https://zenodo.org/record/7682979#.ZBnfhezMK3K.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

Acknowledgement

We thank the anonymous referee for a very careful and helpful report. We thank Tom Kaye for making his observations available for our analysis. This research has made use of NASA’s Astrophysics Data System, the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Programme, and the SIMBAD data base, operated at CDS, Strasbourg, France. The authors acknowledge the MIT SuperCloud and Lincoln Laboratory Supercomputing Center for providing high-performance computing resources that have contributed to the research results reported within this paper.

SX was supported by the International Gemini Observatory, a programme of NSF’s NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation, on behalf of the Gemini partnership of Argentina, Brazil, Canada, Chile, the Republic of Korea, and the United States of America.

SK would like to thank her first computer, Freddie 1, for its service and ultimate death in trying to finish this research. Its memory lives on in her new computer, Freddie 2. She would also like to thank her dog, Chevi, for the joy she brought when coding became frustrating.

DATA AVAILABILITY

Light curves from HAO are available upon request to B. Gary.⁵ Fits files containing the likelihood maps are available as supplementary materials on Zenodo⁶ (Kubiak 2023). We have created a jupyter notebook demonstrating how to use these files and how to recreate the plots that we used in this paper.⁷

Footnotes

References

Author notes