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

Gliese 86 is a nearby K dwarf hosting a giant planet on a $\approx$16-day orbit and an outer white dwarf companion on a $\approx$century-long orbit. In this study we combine radial velocity data (including new measurements spanning more than a decade) with high angular resolution imaging and absolute astrometry from Hipparcos and Gaia to measure the current orbits and masses of both companions. We then simulate the evolution of the Gl 86 system to constrain its primordial orbit when both stars were on the main sequence; the closest approach between the two stars was then about $9\,$AU. Such a close separation limited the size of the protoplanetary disk of Gl 86 A and dynamically hindered the formation of the giant planet around it. Our measurements of Gl 86 B and Gl 86 Ab's orbits reveal Gl 86 as a system in which giant planet formation took place in a disk truncated at $\approx$2$\,$AU. Such a disk would be just big enough to harbor the dust mass and total mass needed to assemble Gl 86 Ab's core and envelope, assuming a high disk accretion rate and a low viscosity. Inefficient accretion of the disk onto Gl 86 Ab, however, would require a disk massive enough to approach the Toomre stability limit at its outer truncation radius. The orbital architecture of the Gl 86 system shows that giant planets can form even in severely truncated disks and provides an important benchmark for planet formation theory.


INTRODUCTION
Thousands of exoplanets are now known in a huge variety of systems, and in an enormous range of dynamical configurations (Luger et al. 2017;Shallue & Vanderburg 2018;Lam et al. 2020). These include hot Jupiters (Butler et al. 1997;Henry et al. 2000;Tinney et al. 2001), outer Jovian planets (Jones et al. 2010;Wittenmyer et al. 2014;Venner et al. 2021), smaller planets * NSF Graduate Research Fellow of all sizes and orbital distances (Barclay et al. 2013;Jenkins et al. 2015;Smith et al. 2018), planets around binaries (P-type systems; Welsh et al. 2012;Orosz et al. 2019;Kostov et al. 2021), and planets around individual stars within binaries (S-type systems; Teske et al. 2016;Hatzes et al. 2003;Ramm et al. 2016). This diversity suggests that planet formation is a robust, if not a universal, process accompanying star formation.
Planet formation in binaries is an especially important testbed for the planet formation process. The existence of the binary provides natural constraints on the properties of circumstellar and circumbinary disks, and therefore on the material available for planet formation. Both P-and S-type planets must form within a disk, but one that is dynamically interacting with the binary in an environment very different from the canonical Solar nebula. Su et al. (2021) conducted a statistical study of the S-type planetary systems detected from radial velocity (RV) surveys to generalize the characteristics of these systems. Table 1 of that paper summarizes the properties of 80 planet-hosting binaries; ten of them (HD 42936, HD 87646, HD 59686, HD 7449, γ Cep, HD 4113, HD 41004, 30 Ari, Gl 86, and HD 196885) have separations smaller than 30 AU. Jang-Condell (2015) argued that the frequent appearance of the planets in close binaries indicates the formation process is robust. However, that the binaries are close to each other limits the amount of materials in the circumstellar disk and significantly reduces the chance of forming planetary embryos. To better understand the planet formation under such conditions, we focus on Gl 86 and investigate its orbital parameters and possible planet formation scenarios in this paper.
Gl 86 is the second-closest planetary system containing a warm or hot Jupiter, after Gl 876, with a distance of 10.761 ± 0.005 pc (Lindegren et al. 2021) and an age of ≈10 Gyr (Fuhrmann et al. 2014). Queloz et al. (2000) discovered this system using the CORALIE echelle spectrograph and found an RV signal corresponding to an m sin i ≈ 4 M Jup planet with a 15.8-day orbital period as well as a distant and massive companion causing a long-term RV drift. Els et al. (2001) used the ADONIS adaptive optics system on the ESO 3.6-m Telescope at La Silla to observe Gl 86 and identified a wide companion that they inferred to be a brown dwarf causing the RV drift. Mugrauer & Neuhäuser (2005) performed additional high-contrast observations and found that this wide companion is, instead, a cool white dwarf, and they ruled out any additional stellar companions between 0. 1 and 2. 1, or 1-23 AU. Lagrange et al. (2006) used VLT/NACO to obtain photometric and astrometric measurements and confirmed that the companion is a white dwarf rather than a brown dwarf and inferred its mass to be 0.48 M ≤ m ≤ 0.62 M based on the amplitude of the RV trend observed by Queloz et al. (2000). Brandt et al. (2019) used all of the above measurements, as well as additional relative astrometry from Farihi et al. (2013) and the proper motion anomaly between Hipparcos and Gaia DR2 (Brandt 2018), to determine the mass and orbital parameters of the white dwarf. The Gl 86 system, with a white dwarf on a ≈20-AU orbit and a close-in, gas-giant exoplanet, challenges planet formation models. From a theoretical perspective, such close binary systems are expected to be hostile to the formation of giant planets due to disk truncation (e.g., Artymowicz & Lubow 1994) and destructive planetesimal collisions (e.g., Paardekooper & Leinhardt 2010;Rafikov & Silsbee 2015), and this is largely borne out by observations (e.g., Wang et al. 2014;Kraus et al. 2016). The Gl 86 system presents a further problem because when both stars were on the main sequence, the separation between them was even smaller and the stability and the feasibility of forming the inner planet becomes even more questionable. Both Lagrange et al. (2006) and Farihi et al. (2013) doubted the orbital stability of the inner planet since the semi-major axis of the primordial binary system was too small. In order to work out a theory regarding the formation of the warm Jupiter in the Gl 86 system, it is necessary to better constrain Gl 86's current and primordial orbital parameters, and the stellar masses.
In this paper, we perform a new fit to the masses and orbits of the Gl 86 system using absolute astrometry from the latest Gaia Data Release (EDR3, Gaia Collaboration et al. 2021;Lindegren et al. 2021), together with RV and relative astrometry data from the literature. We discuss the resulting orbital elements of the Gl 86 system in Section 2. In Section 3, we simulate the binary's evolution based on an N -body integrator program in order to constrain its primordial orbit. In Section 4, we discuss some implications that the primordial orbit has on the formation of the planet surrounding Gl 86 A and the challenges the planet faced by the time it was formed. (Because Gl 86 B starts with a higher mass and ends with a lower mass due to mass loss, we consistently call Gl 86 A the host and Gl 86 B the companion stars, instead of primary and secondary, to avoid confusion.) Finally, Section 5 summarizes our results.

THE CURRENT ORBIT OF GL 86
We use the open-source python package orvara  to fit for the current masses and orbits in the Gl 86 system. The program can fit any combination of RVs, relative, and absolute astrometry from Hipparcos and Gaia. We use all of these types of data to constrain Gl 86. In this section, we describe the input data and our resulting fit.

Data
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σ significance. We use relative astrometry from multiple literature sources. Table 1 lists our adopted relative astrometry for Gl 86 B, where ρ 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. 6220±0. 0040 with position angle 82. • 185±0. • 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.

Orbital Fit
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 A = 1.39 +0.24 −0.23 M by using the cross-calibrated Hipparcos and Gaia DR2 astrometry in a fit to Gl 86 B. Their prior was logflat, 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 A = 0.83 ± 0.05 M . 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, −2 ln L = χ 2 = χ 2 RV + χ 2 rel ast + χ 2 abs ast .
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/ √ 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 × 10 5 steps; we keep every 50th step and use the coldest chain for inference. We discard the first 1.25 × 10 5 steps as burn-in.

Results
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 ≈1 to the total χ 2 . Unfortunately, our best-fit χ 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 bestfit χ 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 ∼10 mas, which could be caused by a ∼20 M Jup companion on a ∼2 AU orbit. However, such companions are rare 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 χ 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 χ 2 values to an acceptable level (χ 2 relsep = 11.8 and χ 2 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 ≈25 AU to ≈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 100 µ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 sig-nal from the planet Gl 86 Ab, as it has a very short period and would obscure the overall trend.

THE PRIMORDIAL ORBIT OF GL 86
In the previous section we obtained the orbital parameters of the current Gl 86 system. Next, we work out the primordial orbit when both stars were on the main sequence. To account for the alteration of the or- ters, NGC 7789 and NGC 6819, plus data from the very old cluster NGC 6791, measuring the current masses of those white dwarfs, and calculating their corresponding progenitor masses. Their relation, is good down to M initial = 1.16 M . With this relation, the initial mass of Gl 86 B is predicted to be 1.36 ± 0.25 M , which is consistent with our MESA result. We adopt M B, initial = 1.39 M throughout this section.
In the remainder of this section, we first review analytic theoretical expectations for the evolution of semimajor axis and eccentricity and confirm that they are consistent with Mercury's results. We then evolve the orbit of a 1.39 M star around a 0.870 M star to infer the initial dynamical configuration of Gl 86 AB.

Mass loss
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, M ≈ −10 −8 M yr −1 ), the asymptotic giant branch (AGB,Ṁ ≈ −10 −8 to − 10 −4 M yr −1 ), and planetary nebula (Ṁ ≈ −10 −6 M 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 M in 0.602 million years, equivalent to a rate of 7.95×10 −7 M 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 ∼10 4 yr. This remains much longer than the orbital period of ≈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, settingv = 0 for isotropic mass loss, and take the orbit average: 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).

Binary evolution
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)). The star symbol at the origin stands for Gl 86 A. When both stars were on the AGB phase, where Gl 86 B mainly started to lose mass, it was at the innermost (yellow) orbit. As it shed mass, the semi-major axis increased, and it gradually spiraled out and ended up as a white dwarf at the outermost (purple) orbit.
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.

Mechanism of planet formation
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 ± 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 ≈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 ± 0.05 Gyr (Farihi et al. 2013), far too lit-tle 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, γ Cep A, HD 41004 A, HD 41004 B, HD 196885 A, and α Cen B. They tested a suite of eighteen different combinations of viscosity parameters α and accretion ratesṀ , given the binary parameters M host , M comp , a, and e, and counted how many (α,Ṁ ) 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 M ⊕ , 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 ≈10 M ⊕ 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 ⊕ and a minimum total mass of 5 M Jup in order to have any chance of forming Gl 86 Ab.

Total disk mass and dust mass
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 dustto-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.

Truncation radius
The protoplanetary disk around Gl 86 A could not have extended past the orbit of Gl 86 B, and in fact would have been truncated at a significantly smaller radius. The first criterion that needs to be satisfied to truncate a disk is that the resonant torque should be greater than the viscous stresses. A viscosity-dependent tidal distortion and resonant interactions exert torques on the disk with opposite directions. If the torque of the resonant interaction surpasses that of the tidal distortion, a gap would be opened and the disk would be truncated.
The magnitude of resonant torques varies at different resonant states (m, l), where m ≥ 0 is the azimuthal number, and l the time-harmonic number. At a given resonant state, there are three different types of resonances: the inner Lindblad resonance (ILR), the outer Lindblad resonance (OLR), and corotational resonance (CR). They occur at different positions according to Equations (9) and (10) of (Artymowicz & Lubow 1994, hereafter AL94): where µ equals M host /(M host + M comp ) for the circumstellar disk of M host , and M comp /(M host + M comp ) for the circumstellar disk of M comp , and 1 for the circumbinary disk. The minus and plus signs in Equation (10) correspond to the inner and outer LR, respectively. The ILR dominates in the circumstellar case, and OLR dominates in the circumbinary case. We are looking for the smallest possible radius at which the resonant torque is greater than the viscous stress or satisfies Equation (15) of AL94. We check the ILR first because it dominates in a circumstellar disk, and it has the smallest radius among the three types of resonant interactions given a (m, l) set. With the method elaborated in AL94, this first criterion breaks down as α 1/2 H r = a |φ ml | GM (πm) 1/2 (m ± 1) 1/6 |λ ∓ 2m| 2µ 2/3 l 2/3 , (11) where λ = m for the ILR of the circumstellar disk and λ = −(m+1) for the OLR of the circumbinary disk, and α is the viscosity parameter (Shakura & Sunyaev 1973). The LHS of Equation (11) indicates the magnitude of viscous stress. The larger the α, the stronger the viscous stress. We can also write the LHS as 1 √ Re in terms of Reynolds number Re, where With M 1 = 0.870 M and M 2 = 1.39 M , the left panel of Figure 9 shows the radial positions where resonant interaction occurs in units of semi-major axis versus the distribution of the eccentricity at which the torque is large enough to balance the viscous stress. As the eccentricity increases, the magnitude of the resonant torque gets larger, which means the disk can be truncated more easily. Squares connected by dotted lines share the same Reynolds number. The grey line in the figure indicates that the eccentricity of Gl 86 is 0.429. We can read the positions where the circumstellar disk is truncated, based on how the grey line intersects with the r vs. e curves. For Reynolds numbers equal to 10 3 , 10 4 , 10 5 , 10 6 , 10 8 , 10 11 , 10 14 , the truncation radii are 0.22a, 0.18a, 0.16a, 0.15a, 0.12a, 0.10a, 0.08a, or 3.3, 2.7, 2.3, 2.2, 1.8, 1.4, 1.2 AU, respectively. The right panel shows the resulting truncation radii with respect to different Reynolds numbers, ranging from 10 3 to 10 14 , at e = 0.429 and in terms of AU. The black dots are the truncation radii corresponding to a sequence of Reynolds numbers starting from 10 3 to 10 14 , incrementing with a factor of 10 for each step.
The second criterion is that gap opening time t open should last for a reasonably short time relative to the viscous timescale. The gap opening time, approximately the same as the viscous closing time, equals t open ≈ (∆r) 2 /ν, where ∆r is the radial extent of the gap, and ν the viscosity. According to the formula of viscosity, where c s the sound speed, and H the height of the disk. H = c s /Ω, where Ω is the Keplerian orbital angular frequency So, .
Since α ranges between 0.001 − 0.1, t open is approximately 1 P −100 P , which is short compared to the disk's lifetime on the main sequence.

Inner rim
The inner part of a protoplanetary disk can be divided into four components: a dust-free region, a dust halo, a condensation front, and an optically thick disk ( Figure  1 of Ueda et al. 2017, henceforth U17). The dust in this region functions as a feedback regulation. If the temperature of the disk is lower than T ev , then the dust condenses and heats up. Otherwise, the temperature goes down when the dust evaporates since the emissionto-absorption ratio of the dust is lower than that of the gas (U17). So, the temperature in this region is approximately high enough to evaporate dust, and the inner rim of the disk lies between the dust halo and condensation front. According to Flock et al. (2016) and U17's simulations, based on T * = 10000 K, M * = 2.5 M , and R * = 2.5 R models, the resulting radial profile of midplane temperature is step-like, with a temperature of T ev ≈ 1470 K in the dust halo region spanning from 0.35 − 0.45 AU from the star.
Such stellar parameters resemble those of Herbig Ae/Be stars, pre-main-sequence stars with a mass between 2 − 8 M (Natta et al. 2001). Kama et al. (2009) studied such stars and verified observational data by exploring the inner rim of their protoplanetary dust disks with a Monte Carlo radiative transfer simulation. They worked out an analytical expression of R rim in terms of dust, disk, and stellar properties. In terms of their result, high surface density, large silicate grains, and iron and corundum grains reduce the rim radius. With a typical surface density of 1 g cm −2 , ordinary 10 µm grains give a 0.4 AU rim, while grains as large as 100 µm give a 0.2 AU rim, and those in the minor extreme 0.1 µm give an inner rim of 2.2 AU.
Gl 86 A, however, is a Sun-like star. The disk structure of its inner rim is much different, as the dust sublimation temperature is reached closer to the star. With a specific accretion rate of 3.6 × 10 −9 M yr −1 , irradiated hydrostatic disk models show that the silicate sublimation front begins at around 0.08 AU, a curved dust rim exists between 0.08 AU and 0.15 AU, a small shadowed region between 0.2 − 0.3 AU and a flared disk beyond 0.3 AU (Flock et al. 2019). The simulation also shows a steep rise of gas density at 0.13 AU. Therefore, we set 0.13 AU as the start of the inner rim of Gl 86 A's protoplanetary disk. This is similar to Gl 86 Ab's current orbital radius of 0.11 AU.

Disk Mass
With the truncation radius and the inner rim worked out, the next step is to compute the total mass of the disk and the mass of solid material available to form the core of Gl 86 Ab. In this section we derive the expression of the disk surface density and then integrate it numerically.
The hydrostatic equilibrium equation readṡ where Σ G is the surface density of gas. We combine Equations (12), (14), (15) and (17) to get The integral of Equation (18) from the inner rim to the truncation radius is the total mass of a protoplanetary disk. This result is consistent with Table 1 of Jang-Condell & Sasselov (2003). In order to form Gl 86 Ab, the protoplanetary disk must have had enough solids to form the planetary core and enough mass to supply the gaseous envelope. The total dust mass is the disk mass multiplied by a dustto-gas ratio, which starts around 10 −10 near the star and increases rapidly near the inner edge of the disk. It reaches a plateau of 10 −2 when it reaches the optically thick disk (e.g., Figure 2 of U17). We adopt 10 −2 as the dust-to-gas ratio. A minimum 10 M ⊕ total dust mass implies a disk mass of 1000 M ⊕ , or ≈3 M Jup . The disk must also have supplied Gl 86 Ab's envelope, and must therefore have exceeded the planet's current mass of ≥4.3 M Jup (depending on the inclination). We adopt 5 M Jup as the minimum disk mass that could have possibly permitted the formation of Gl 86 Ab.
We ignore Reynolds numbers higher than 10 8 , as such Reynolds numbers rarely occur, and a typical Reynolds number is only ∼10 5 (AL94). For large Reynolds numbers, the truncation radii are so small that they are very close to the inner rim. However, as the disk mass is proportional to Re (see Equation (18)), it can still be unreasonably large when Re is large even if the disk mass is integrated from the inner rim to the truncation radius. Huge Reynolds numbers require extremely small viscous stresses α, and the resulting disk mass exceeds the current mass of Gl 86 A. Figure 10 shows the total disk mass in terms of Reynolds number and accretion rate. With accretion Figure 10. Total mass of the protoplanetary disk of Gl 86 A in terms of Reynolds number and accretion rate. The mass is negligible at the lower-left corner, with small Reynolds numbers and low accretion rates, and large at the upperright corner, with large Reynolds numbers and large accretion rates. The dashed white line indicates a disk mass of 1000 M⊕, or 3.0 × 10 −3 M , or a dust mass of 10 M⊕ assuming a dust-to-gas ratio of 10 −2 . The dash-dotted white line represents 5 MJup disk mass. The protoplanetary disk of Gl 86 A must have occupied the space to the right of these white lines in order to have possibly formed Gl 86 Ab. To the right of the pink line the disk is Toomre unstable (Q < 1) at its truncation radius. rates ranging from 10 −9 M yr −1 to 10 −7 M yr −1 and Reynolds numbers from 10 5 to 10 8 , the total disk mass goes from 10 −5 M to 1 M . The region to the right of the dashed (dash-dotted) white line indicates dust (disk) mass higher than 10 M ⊕ (5 M Jup ). To the right of both lines, the disk exceeds the minimum mass to supply the material to form Gl 86 Ab. Figure 10 shows that with a Reynolds number 10 7 , the protoplanetary disk can contain enough material to form Gl 86 Ab despite its truncation at ≈2 AU. This minimum mass corresponds to about twice the minimum mass Solar nebula (Hayashi 1981) integrated from 0.13 AU to 2 AU. Still, this requires a very high efficiency in converting disk mass to planet mass. At very high disk masses, which would allow for a less efficient conversion of disk mass to planet mass, the disk approaches the Toomre stability limit (Safronov 1960;Toomre 1964). The pink line on Figure 10 indicates that this stability limit has been breached at the disk's outer truncation radius.
Changes to our assumptions can present further difficulties in accounting for Gl 86 Ab. One example is our assumption of a Z = 0.010 metallicity for the system. Literature measurements of [Fe/H] for Gl 86 A vary from −0.30 (Ramírez et al. 2007) to −0.16 (Allende Prieto et al. 2004), while the majority fall between −0.25 and −0.20 (e.g. Santos et al. (2000Santos et al. ( , 2004; Ramírez et al. (2013)). We adopt −0.23 for [Fe/H] and 0.0134 for solar metallicity (Grevesse et al. 2012), which give Z = 0.0134 × 10 −0.23 ≈ 0.008. We approximate this value to be 0.01. Using Z = 0.01 and Figure 4 of Bitsch et al. (2015), which displays the relation between the initial formation distance, formation time, final distance, and final mass of a planet, if Gl 86 Ab is formed at 3 AU away from the host star (just beyond the truncation radius) and the final distance is ∼0.1 AU, the figure indicates that the final mass is only several hundred times the Earth's mass. This number is a factor of ≈5 smaller than the mass of Gl 86 Ab.
Another challenge is the reliability of the planet formation mechanisms, namely pebble accretion and planetesimal accretion. These mechanisms can be characterized by inefficient coagulation of the solid material onto the planetary core (Guillot et al. 2014). Inefficiencies in converting disk solids to a protoplanetary core, and subsequently accreting the disk onto the forming planet, would place the disk closer to the Toomre stability limit. Gl 86 Ab thus presents an important example of either a very massive, albeit truncated, disk, and/or efficient conversion of a protoplanetary disk's mass into a warm Jupiter.

CONCLUSIONS
In this study we fit for the orbital parameters of the Gl 86 system based on RV and astrometric data. We find the white dwarf secondary, Gl 86 B, to have a mass 0.5425 ± 0.0042 M and the inner planet Gl 86 Ab to have m sin i = 4.266 +0.11 −0.087 M Jup . We cannot constrain the inclination and orientation of the inner planet, but obtain good constraints on all other parameters. We obtain an eccentricity of 0.429±0.017 and a semimajor axis of 23.7 ± 0.3 AU for the white dwarf's current orbit. In order to obtain a satisfactory fit to HST astrometry, we require an inflation of ≈10 mas in the relative astrometry uncertainties. This could also point to an unseen massive companion orbiting Gl 86 B, which would make the system especially unique. The existence of such a companion might be confirmed or refuted with individual epoch astrometry from a future Gaia data release, since both Gl 86 A and Gl 86 B are detected in Gaia, or by further astrometric monitoring of the system.
Mass loss by the white dwarf Gl 86 B means that its current orbit differs from its primordial orbit. We combine the MESA initial-final mass relation with a current white dwarf mass of 0.5425±0.0042 M to infer an initial mass of 1.39 M . We then run a simulation backward in time to derive the primordial orbit of Gl 86 AB. We verify that the semimajor axis a satisfies M tot a = constant and eccentricity e satisfies e t = 0, in agreement with analytic theory. When both stars were on the main sequence, we infer a semimajor axis of 14.8 AU, and an eccentricity matching its current value of 0.429. This relation assumes isotropic, adiabatic mass loss, as the maximum size of the Gl 86 B was not large enough to initiate Roche lobe overflow.
Finally, we examine how the formation of Gl 86 Ab took place under such a dynamically challenging situation. We find a truncation radius of ≈2 AU for Gl 86 A's protoplanetary disk, with somewhat smaller values at higher Reynolds numbers, by balancing the viscous stress and inner Lindblad resonant (ILR) torque. We then derive an expression of the total disk mass and infer a total dust mass assuming a dust-to-gas ratio of 1%. Despite the short separation between the two stars when they were both on the main sequence, the disk mass is sufficient to provide the material to form Gl 86 Ab in the high accretion rate and large Reynolds number (low viscosity) range. This scenario requires an efficient conversion of dust to a planetary core. Inefficient conversion of material would require a more massive disk, which would then approach the Toomre stability limit at its outer truncation radius. The Gl 86 system, with a 4 M Jup planet formed within a disk truncated at ≈2 AU, demonstrates that giant planet formation is possible even in adverse circumstances. It shows that severely truncated disks around stars in binaries can birth super-Jovian exoplanets, and provides an important benchmark for planet formation theory.

ACKNOWLEDGMENTS
Some of this research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with program 14076. T.D.B. gratefully acknowledges support from the National Aeronautics and Space Administration (NASA) under grant #80NSSC18K0439. We acknowledge the traditional owners of the land on which the AAT stands, the Gamilaraay people, and pay our respects to elders past, present, and emerging.