The Gliese 86 Binary System: A Warm Jupiter Formed in a Disk Truncated at $\approx$2 AU

We use the absolute astrometry from Hipparcos (ESA, 1997; van Leeuwen, 2007) and Gaia’s latest data release (Gaia Collaboration et al., 2021) as cross-calibrated by Brandt (2021). The Hipparcos-Gaia Catalog of Accelerations (HGCA, Brandt, 2018, 2021) provides three proper motions on the Gaia EDR3 reference frame; differences between them indicate astrometric acceleration. For Gl 86, the two most precise proper motions are the one computed from the position difference between Hipparcos and Gaia and the Gaia EDR3 proper motion. These two measurements are inconsistent with constant proper motion at nearly $300\sigma$ significance.

We use relative astrometry from multiple literature sources. Table 1 lists our adopted relative astrometry for Gl 86 B, where $\rho$ is the separation between the two stars and PA is the position angle (east of north). The last data point in Table 1 was imaged on 2016 November 10 using the Space Telescope Imaging Spectrograph (STIS) as part of program GO-14076 (PI Gänsicke). The imaging was performed using the narrow-band filter F28X50OII (central wavelength 3738 Å, full width at half maximum 57 Å) with a series of four two-second exposures in a standard dither pattern. This filter has no red leak, and thus the bright host star remains unsaturated and in the linear response regime. The white dwarf companion is also detected in all four frames, at approximate signal-to-noise ratios between 13 and 19. This set of images was used to robustly measure the separation of the binary, where the companion star was found at offset $2.\!\!^{\prime\prime}6220\pm 0.\!\!^{\prime\prime}0040$ with position angle $82.\!\!^{\circ}185\pm 0.\!\!^{\circ}098$ under the J2000 frame.

We note that the relative position measurement of the two stars in Gaia EDR3 is less straightforward than the other relative astrometry measurements. It is the difference between the 2016.0 positions of five-parameter astrometric solutions to each star in the binary. The formal uncertainties are tiny, but are subject to possible systematics from the proximity of the two stars and from their straddling of the $G=13$ mag boundary where the window function changes (Gaia Collaboration et al., 2021; Cantat-Gaudin & Brandt, 2021). We treat the measurement as instantaneous and, somewhat arbitrarily, adopt uncertainties similar to the HST uncertainties. This avoids having Gaia solely drive the result and mitigates the possible impact of systematics.

The RV data come from the UCLES échelle spectrograph (Diego et al., 1990) on the Anglo-Australian Telescope. Those results spanning 1998 to 2005 were published in Butler et al. (2006). We also include a further 34 previously unpublished UCLES RVs spanning 2006 to 2015 (Table 2), for a total time baseline of 17.8 years. The RV data are processed through the same pipeline, but they have some discrepancy with Butler et al. (2006) due to minor pipeline tweaks and the fact that they have the mean stellar RV subtracted. Thanks to Gl 86 A’s very large acceleration, the mean RV has changed appreciably with an additional nine years of data.

We use orvara to fit a superposition of Keplerian orbits to the Gl 86 A astrometry and RVs and relative astrometry. We use log-uniform priors for semimajor axis and companion mass, a geometric prior on inclination, and uniform priors on the remaining orbital parameters. We adopt the Gaia EDR3 parallax as our parallax prior; orvara analytically marginalizes parallax out of the likelihood. We use a log-uniform prior on RV jitter.

We adopt an informative prior on the mass of Gl 86 A. Brandt et al. (2019) obtained $M_{\rm A}=1.39_{-0.23}^{+0.24}\,\mbox{$M_{\sun}$}$ by using the cross-calibrated Hipparcos and Gaia DR2 astrometry in a fit to Gl 86 B. Their prior was log-flat, but stellar evolution allows a much narrower prior. Fuhrmann et al. (2014) modelled Gl 86 A’s atmosphere using high-resolution spectroscopy and concluded $M_{\rm A}=0.83\pm 0.05\,\mbox{$M_{\sun}$}$. We adopt this as our prior on Gl 86 A’s mass.

The log likelihood function consists of three parts: the chi-squared values of RV (including a penalty term for RV jitter), relative astrometry, and absolute astrometry. To maximize the likelihood is equivalent to minimizing the negative log likelihood,

In addition, the RV zero point, parallax, and the proper motion of the system’s barycenter are marginalized out as nuisance parameters. We refer readers to the orvara paper for more details on the formulas and techniques used.

We use Markov Chain Monte Carlo (MCMC) to explore the posterior probability distribution with the emcee (Foreman-Mackey et al., 2013) and ptemcee (Vousden et al., 2016) packages. Parallel tempering MCMC uses walkers at many temperatures, each multiplying an extra factor of $1/\sqrt{T}$ to the exponent of the likelihood. A larger temperature means a more flattened out posterior probability distribution and enables hotter temperature walkers to explore more parameter space. Temperature swaps can happen periodically while preserving detailed balance. Parallel-tempering helps to explore multimodal posteriors and avoid getting stuck at some local minimum. We use $100$ walkers, $30$ temperatures, and $2\times 10^{5}$ steps; we keep every 50th step and use the coldest chain for inference. We discard the first $1.25\times 10^{5}$ steps as burn-in.

Our first step with our resulting chains is to test whether they include formally well-fitting orbits. A satisfactory fit will have each data point contribute $\approx$1 to the total $\chi^{2}$. Unfortunately, our best-fit $\chi^{2}$ of relative separation is $59.4$, and that of position angle is $50.5$; both are much too large for 8 data points. Our high best-fit $\chi^{2}$ values show that either we have underestimated uncertainties, or there is an additional component in the system. Any additional component massive enough to affect the astrometry would have to orbit Gl 86 B to avoid detection in the precision RVs and direct imaging of Gl 86 A. Bringing the relative astrometry into agreement requires a perturbation of $\sim$10 mas, which could be caused by a $\sim$20 $M_{\rm Jup}$ companion on a $\sim$2 AU orbit. However, such companions are rare (Marcy & Butler, 2000; Halbwachs et al., 2000), and we have just two relative astrometry measurements at mas precision. This is insufficient to constrain the mass and orbit of a hypothetical substellar companion to Gl 86 B. We provisionally attribute the high $\chi^{2}$ values to a combination of underestimated uncertainties and systematics in the data.

For our final orbit analysis, we inflate the uncertainties in the absolute astrometry by a factor of $2$, add $10$ mas to our relative separation uncertainties and $0.05$ degrees to our position angle uncertainties, both in quadrature, in addition to the error inflation used by Brandt et al. (2019). This brings the $\chi^{2}$ values to an acceptable level ($\chi^{2}_{\rm relsep}=11.8$ and $\chi^{2}_{\rm PA}=11.3$), and has only a minor impact on our derived parameters. The mass of Gl 86 B changes by just 0.5%, while the best-fit semimajor axis decreases from $\approx$25 AU to $\approx$24 AU and the best-fit eccentricity increases from 0.38 to 0.43. Table 3 lists the full set of orbital parameters.

The MCMC results are summarized in Figures 1 and 2. The parameters for the white dwarf companion Gl 86 B are well constrained, but some parameters are strongly correlated. For example, its semi-major axis is anti-correlated with the eccentricity, and its inclination is also anti-correlated with its longitude of ascending node. The inner planet, Gl 86 Ab, has an orbital period much shorter than either the Hipparcos or Gaia mission baseline. This, combined with the planet’s low $m\sin{i}$, means that we have almost no constraint on its orbital inclination or orientation. Even epoch astrometry from Gaia DR4 might not detect the $\lesssim$100 $\mu$as orbit of Gl 86 A about its barycenter with Gl 86 Ab.

Figure 3 shows the relative astrometric orbit of Gl 86 AB, and Figure 4 shows the radial velocity, separation, position angle, and proper motions as a function of time. For display purposes, we have removed the signal from the planet Gl 86 Ab, as it has a very short period and would obscure the overall trend.

After a star evolves off the main sequence and passes the sub-giant branch, its radius expands and its envelope becomes loosely bound. The radiation pressure due to photon flux expels the envelope. The star will experience mass loss during the red giant branch (RGB, $\dot{M}\approx-10^{-8}\,\mbox{$\mbox{$M_{\sun}$}$\,yr${}^{-1}$}$), the asymptotic giant branch (AGB, $\dot{M}\approx-10^{-8}\leavevmode\nobreak\ \text{to}-10^{-4}\,\mbox{$\mbox{$M_{\sun}$}$\,yr${}^{-1}$}$), and planetary nebula ($\dot{M}\approx-10^{-6}\,\mbox{$\mbox{$M_{\sun}$}$\,yr${}^{-1}$}$) phases, during which mass is cast away in the form of an isotropic stellar wind.

MESA computes the (isotropic) mass loss of Gl 86 B as a function of time. Figure 5 shows the mass of Gl 86 B versus the star age. Although Gl 86 B lost mass most rapidly by the end of the AGB, it lost $0.478\,\mbox{$M_{\sun}$}$ in $0.602$ million years, equivalent to a rate of $7.95\times 10^{-7}\,\mbox{$\mbox{$M_{\sun}$}$\,yr${}^{-1}$}$, if we approximate the average mass loss as linear. Even during the final thermal pulses that expel the envelope, mass loss proceeds on a timescale of $\sim$10⁴ yr. This remains much longer than the orbital period of $\approx$100 yr, making the mass loss very nearly adiabatic.

We can apply this in the relation between speed $v$ and separation $r$ for a Keplerian orbit with semimajor axis $a$,

We differentiate with respect to time, setting $\dot{v}=0$ for isotropic mass loss, and take the orbit average:

where

and $P$ is the orbital period. Equations (5) and (4) yield

which recovers Equation (16) of Jeans (1924).

For the secular evolution of the eccentricity, Hadjidemetriou (1963) presents the evolution of eccentricity in a binary system with slow mass loss,

where $f$ is the true anomaly. Taking the time average and using Equation (13) of Dosopoulou & Kalogera (2016), the secular result is then

Equations (6) and (8) show that isotropic, adiabatic mass loss will expand the orbit, but at fixed eccentricity. This still holds with anisotropic mass loss, as long as we keep the adiabatic assumption and further assume that the anisotropic wind and jets are symmetric with respect to the equator (Veras et al., 2013; Dosopoulou & Kalogera, 2016).

We next consider mass transfer due to the ejection of Gl 86 B’s envelope and its possible accretion onto Gl 86 A through the Roche lobe. Once a star is large enough to fill out its Roche lobe, it transfers mass to its companion via the first Lagrangian point $L_{1}$. We can tell whether a star passes its Roche lobe or not using a byproduct of the MESA simulation, the radius of Gl 86 B as a function of time, and a calculator for Roche lobe properties (Leahy & Leahy, 2015), to show the 3D illustration of the Roche lobe of Gl 86 B (Figure 6). The result shows that the Roche lobe is very spacious, and the surface of Gl 86 B is far from touching that of the Roche lobe when Gl 86 B reaches its maximum size during the AGB stage. This suggests that there was no mass transfer/ejection throughout the evolution of the binary, when the primordial semi-major axis was $14.8$ AU (see Section 3.2).

In Section 3.1, we have argued that Gl 86 B’s mass loss was mostly isotropic. We further assume that the mass loss was adiabatic and linear for the whole duration when Gl 86 B lost mass, because the fractional mass loss rate is small on an orbital time scale. Our simplifying assumption of a constant mass loss rate does not change the evolution’s overall physical results (see Equation (6)).

We integrate the orbit of Gl 86 B around Gl 86 A using Mercury with the general Bulirsch-Stoer algorithm. It is slow but accurate in most situations and can dynamically adjust the step size. Our Mercury integration reproduces the orbit of the companion with respect to the host (Figure 7) and enables us to visualize its slow and steady expansion throughout time.

Mercury’s simulation agrees with the secular evolution of semi-major axis and eccentricity in Equation (6) and Equation (8); these are shown in the two panels of Figure 8.

There are two main phases of the formation of a Jovian planet by core accretion. The first phase is the formation of the planetary core, for which two main theoretical mechanisms have been advanced: planetesimal accretion and pebble accretion. Planetesimal accretion, once the dominant hypothesis (Pollack et al., 1996), has difficulty forming a core and accreting a gaseous envelope before the protoplanetary disk dissipates (Lambrechts & Johansen, 2012). The theory of pebble accretion has become increasingly popular as a solution to this problem (Bitsch et al., 2015; Johansen & Lambrechts, 2017; Bitsch et al., 2019). Although planetesimal accretion still dominates during the early stage of the growth of the planetary embryo, pebble accretion becomes significant when the core is several hundred kilometers in size. Pebbles passing by the core experience the drag force of the gas in the protoplanetary disk, undergo Bondi-Hill accretion, and quickly spiral in (Johansen & Lambrechts, 2017). This drastically reduces the time scale and enables the formation of the planet within the lifetime of the protoplanetary disk.

The second phase starts after the core reaches the isolation mass, during which the core carves a shallow gap and stops the drift of pebbles from accumulating to the core. Then, the core starts to accrete from the gaseous envelope. This process is slow until runaway gas contraction initiates when the mass of the envelope exceeds that of the core (Bitsch et al., 2015; Bitsch et al., 2018).

In rare cases, a planet could start as a circumbinary planet and be tidally captured by one of the stars, ending up as a circumstellar planet (Gong & Ji, 2018). It can also form around one star and be captured by the other (Kratter & Perets, 2012). However, both mechanisms typically produce highly eccentric planets, in conflict with Gl 86 Ab’s low observed eccentricity of $0.0478\pm 0.0024$. It is possible that Gl 86 Ab was captured as an eccentric planet and tidally circularized. However, its orbital period is 15.8 days, while the tidal circularization period cutoff for the 4 Gyr-old open cluster M67 has been measured to be $\approx$12 days (Meibom & Mathieu, 2005; Geller et al., 2021). Planets on longer periods remain eccentric; they need more than 4 Gyr to circularize. The cooling age of the white dwarf Gl 86 B is just $1.25\pm 0.05$ Gyr (Farihi et al., 2013), far too little time for tidal circularization of Gl 86 Ab. In the rest of the paper, we consider only the scenario in which Gl 86 Ab formed around its current host star.

The mechanism of planetary formation in specific close binaries has been investigated by Jang-Condell (2015). They studied the feasibility of planet formation in HD 188753 A, $\gamma$ Cep A, HD 41004 A, HD 41004 B, HD 196885 A, and $\alpha$ Cen B. They tested a suite of eighteen different combinations of viscosity parameters $\alpha$ and accretion rates $\dot{M}$, given the binary parameters $M_{\rm host}$, $M_{\rm comp}$, $a$, and $e$, and counted how many ($\alpha$, $\dot{M}$) models satisfied core accretion or disk instability mechanism, respectively. The conventional criterion for core accretion to occur is whether the total solid mass in the disk exceeds $10\,\mbox{$M_{\oplus}$}$, the least amount of mass to create a rocky core and initiate gas accretion. Jang-Condell (2015) found that except for HD 188753 A, where none of the models fit the observed system properties, core accretion was overwhelmingly more likely than disk instability for the formation of the planet. Jang-Condell et al. (2008) also pointed out that disk instability can only occur in the most massive disks with extremely high accretion rates.

In order to form Gl 86 Ab, we must then satisfy two requirements. First, we must have enough solid material in the disk to assemble a $\approx$10 $M_{\oplus}$ core. Second, the disk itself must exceed the current mass of Gl 86 Ab in order to supply the gaseous envelope. In the following section we will work out the total mass of Gl 86’s disk under different models. We will assume a minimum dust mass of 10 $M_{\oplus}$ and a minimum total mass of 5 $M_{\rm Jup}$ in order to have any chance of forming Gl 86 Ab.

In this section we compute both the total mass of the disk and its mass in dust suitable for forming the core of Gl 86 Ab. We take the total dust mass to be a dust-to-gas ratio multiplied by the total disk mass, which is the integral of the disk’s surface density from the inner rim to the outer truncation radius. We will compute the truncation radius due to the stellar companion first and the inner rim next.