Gaia-4b and 5b: Radial Velocity Confirmation of Gaia Astrometric Orbital Solutions Reveal a Massive Planet and a Brown Dwarf Orbiting Low-mass Stars

Since its launch in 2013 Gaia has been observing the entire sky as dictated by its ‘scanning law’ (Gaia Collaboration et al., 2016). As stars move across the focal plane of each of the two Gaia telescopes, their one-dimensional coordinates are measured with high precision, and the results are assembled into a global astrometric model. The DR3 results were based on 34 months of observations and about 40 scans per star, on average. As noted earlier, for a subset of the stars, the observed astrometric motion consisted not only of the usual parallax and proper motion effects but also two-body Keplerian orbital motion (Holl et al., 2023; Gaia Collaboration et al., 2023b). We retrieved the Gaia two-body solutions for all of the Gaia ASOI candidates from the Gaia archive²²2https://gea.esac.esa.int/archive/ and obtained the covariance matrix for the parameters using the nsstools³³3https://www.cosmos.esa.int/web/gaia/dr3-nss-tools code.

To detect potentially massive substellar companions around nearby M and K stars being unveiled by Gaia, we conducted spectroscopic observations of 28 northern-hemisphere Gaia ASOI candidates to characterize them and rule out false positives. These 28 systems constitute the northern M and K stars in the Gaia-ASOI list with a declination $\delta\gtrsim-20^{\circ}$ and $G_{\mathrm{Bp}}-G_{\mathrm{Rp}}>1.5$. Figure 1 gives an overview of the current Gaia ASOI list in a color-magnitude diagram. The 28 Gaia ASOIs observed as part of this work are highlighted with square markers, and Gaia-4b and Gaia-5b are shown with star markers. We highlight sources that are observed to be double-lined spectroscopic binaries (SB2) in red, systems that are still consistent with being single stars in blue, and systems which were previously known to host planets based on RV discoveries in green circles. We discuss the false positives and the designation in further detail in Section 6 and in Table 4 in the Appendix.

The Gaia two-body solution yields constraints on the following parameters:

where $A$, $B$, $F$, $G$ are the Thiele-Innes coefficients, $e$ is the eccentricity, $P$ is the orbital period, $t_{p}$ is the time of periastron, and $\varpi$ is the parallax. To convert between the Thiele-Innes coefficients and the Campbell orbital elements—the argument of periastron ($\omega$), longitude of ascending node ($\Omega$), angular semi-major axis ($a_{0}$) and the inclination ($i$)—we used the nsstools code. The relevant Campbell elements from the Gaia two-body solutions are listed in Table 2 and 3.

To create statistical samples of the posterior probability distribution based on the Gaia covariance matrix, we followed Winn (2022) and Fitzmaurice et al. (2023) by defining the likelihood function

where $\Theta$ is a vector of deviations between the values reported in the Gaia two-body solution and the sampled value, and $\mathcal{C}_{\text{scale}}$ is the covariance matrix from the Gaia two body solution including a multiplicative factor, $\sigma_{\text{scale}}$, used to optionally scale the original Gaia two-body solution covariance matrix $\mathcal{C}$,

Here, $\sigma_{\text{scale}}$ is a multiplicative factor applied to the tabulated uncertainties. It has the effect of increasing the overall level of uncertainty in the Gaia solution while preserving the relative uncertainties and correlations of the parameters.

For the RV only-fit, we sampled the following parameters:

where $K$ is the RV semi-amplitude, and $\gamma$ is an RV offset parameter (one for each instrument). We used the likelihood function

where $v_{i,\mathrm{obs}}$ is the $i$-th RV data point and $\sigma_{v,i}$ is the associated uncertainty, and $v_{i}$ is the model evaluated at the same time as the $i$-th data point.

To perform the fit, we first used the differential evolution algorithm implemented in PyDE (Parviainen, 2016) to find the maximum likelihood solution, placing broad uninformative priors on the MCMC parameters. We then use the emcee code (Foreman-Mackey et al., 2013), and initialized 100 walkers around the global maximum solution to perform Markov-Chain Monte Carlo (MCMC) sampling of the posteriors. We ensured that the resulting chains were fully converged and well mixed through visual inspection of the chains, removing initial chains as ‘burn-in’, verifying that the Gelman-Rubin statistic is within 1% of unity, and following the recommendation given in the emcee documentation to compute the autocorrelation length of the chains and verify that each chain has at least 50 effectively uncorrelated samples. For the final fits, we performed 50,000 MCMC steps, which was sufficient to ensure well-mixed chains.

To jointly sample the RVs and the Gaia two-body solution, we followed Winn (2022) and Fitzmaurice et al. (2023) and use the following parameters as jump parameters in the MCMC sampling:

where $M_{\star}$ is the mass of the host star, $m_{2}$ is the mass of the companion, $\cos i$ is the cosine of the orbital inclination, $\varpi$ is the parallax, $\epsilon$ is parameter denoting the flux ratio between the companion and the host star, and $\sigma_{\mathrm{scale}}$ is the uncertainty scaling factor for the Gaia covariance matrix (see Equation 3). We adopted the joint likelihood function

where $\log(\mathcal{L}_{g})$ and $\log(\mathcal{L}_{v})$ are given in Equation 2 and 5, respectively.

In the sampling of the joint likelihood, the $\sigma_{\text{scale}}$ parameter has the effect of scaling the importance of the two likelihoods relative to each other giving the RV and the astrometric datasets different weights. We considered different fits fixing $\sigma_{\mathrm{scale}}=1$, and fits allowing $\sigma_{\mathrm{scale}}$ to float, to test the agreement of the Gaia two body solution uncertainty estimates with the observed RVs. Due to the good agreement between the RV and the astrometric data for both Gaia-4b and Gaia-5b, we found that those runs resulted in fully consistent parameters. We elected to present the values from the $\sigma_{\mathrm{scale}}=1$ runs, given the simpler model.

To perform the joint sampling, we followed similar steps as the ‘RV-only’ fit discussed above. We used PyDE to find a global maximum likelihood solution with broad uninformative priors on the jump parameters, from which we used emcee to perform MCMC sampling around the global maximum solution. We assessed that the chains were well mixed with similar criteria as discussed in the ‘RV-only’ sampling.

Gaia-4b and Gaia-5b are substellar objects. We now ask the perennial question, are they ‘giant planets’ or ‘brown dwarfs’? The distinction of great interest to many researchers is in the formation of such objects (Burrows et al., 2001; Chabrier et al., 2014). One could use the term ‘brown dwarf’ for a substellar object that forms through instabilities or fragmentation of dense molecular gas either within a disk (e.g., Boss, 1997, 2006; Boss & Kanodia, 2023) or via collapse of molecular filaments (e.g., Bate et al., 2002; Bate, 2008), while reserving the term ‘giant planet’ for an object that forms via core-nucleated or pebble-assisted accretion within a disk (e.g., Pollack et al., 1996; Alibert et al., 2005; Mordasini et al., 2012). However, it is difficult to determine the formation history of a given object. Further, theoretical investigations into the different formation paradigms have shown that they produce overlapping ranges of companion masses and orbital separations. Gravitational disk instability can form objects as low in mass as Jupiter (Boss, 1997; Boss & Kanodia, 2023), while core-nucleated accretion can form objects as massive as 10 $M_{\mathrm{J}}$ under the right conditions (Mordasini et al., 2012), blurring the boundary.

An alternative way to differentiate between planets and brown dwarfs is based on whether deuterium burning occurs within the interior of the object (Chabrier & Baraffe, 2000; Burrows et al., 2001). However, deuterium burning cannot be observed directly. Instead, theoretical calculations are used to determine the lower limit on the mass of an object that undergoes deuterium burning. Current calculations lead to an expected lower mass limit for deuterium burning between 11-16$M_{\mathrm{J}}$ depending on the assumed composition (Mollière & Mordasini, 2012). A somewhat arbitrary and frequently used dividing line is $13M_{\mathrm{J}}$. From this perspective, Gaia-4b is a giant planet and Gaia-5b is a brown dwarf.

Given the blurred boundaries between massive planets and brown dwarfs, objects in these overlapping windows of massive planets and brown dwarfs present an intriguing opportunity for further study to gain further insights into their formation histories through observations of different observables including orbital distance, eccentricity, and the host star metallicity and mass.

From the perspective of the host star metallicity, core-accretion and gravitational instability predict different observable trends. Core-accretion theory predicts that stars with massive companions should preferentially be metal-rich, since a higher mass of solid material would facilitate the formation of massive companions (Fischer & Valenti, 2005; Santos et al., 2017; Osborn & Bayliss, 2020). In contrast, gravitational instability does not predict that giant-hosting stars should preferentially be metal-rich (Maldonado & Villaver, 2017; Maldonado et al., 2019). From the analysis of the HPF spectra, Gaia-4 has a metallicity of $\mathrm{[Fe/H]}=0.05\pm 0.13$ consistent with solar metallicity. This low metallicity for such a massive companion, could therefore favor formation via gravitational instability over core-accretion. For Gaia-5, we measure a sub-solar metallicity of $\mathrm{[Fe/H]}=-0.18\pm 0.13$, in line with tentative trends of solar and sub-solar metallicities measured for low-mass stars hosting brown dwarfs (Maldonado et al., 2019), also pointing towards formation via gravitational instability.

With respect to stellar host star mass, the occurrence of massive planets is known to decrease with decreasing stellar mass (Endl et al., 2006; Johnson et al., 2010; Maldonado et al., 2019; Sabotta et al., 2021; Schlecker et al., 2022; Gan et al., 2023; Bryant et al., 2023; Pass et al., 2023). This has been connected to the fact that less massive stars tend to have less massive protoplanetary disks (Pascucci et al., 2016; Ansdell et al., 2016; Manara et al., 2022), and a lower-mass protoplanetary disk has a smaller inventory of solid materials to nucleate runaway gas accretion. In contrast, brown dwarfs are observed to be rare across all spectral types with an occurrence rate $\lesssim$1% (Sahlmann et al., 2011; Carmichael et al., 2020; Barbato et al., 2023), which can be taken as evidence for a different formation pathway for brown dwarfs as compared to giant planets.

Figure 5 shows the masses of known planets and brown dwarfs orbiting stars with masses between 0.1 and 0.7 $M_{\odot}$. A continuum of planets and brown dwarfs has been found around the stars with masses between 0.3 and 0.7 $M_{\odot}$ (Figure 5). However, for stars less than $0.3M_{\odot}$, there is a lack of known Jupiter-mass planets ($1-10M_{\mathrm{J}}$) (Pass et al., 2023), highlighted in blue in Figure 5. It is unlikely that the gap is due to an observational bias against such systems, as recent RV surveys, such as Sabotta et al. (2021) and Pass et al. (2023) would have been sensitive to such massive planets around such low-mass stars. Could this ‘Jupiter desert’ reflect a dividing line between different formation channels for brown dwarfs and planets?

Core accretion struggles to form gas giant planets around the lowest mass stars due to the low inventory of dust available along with long formation timescales to form solid cores (e.g., Burn et al., 2021; Miguel et al., 2020; Stefánsson et al., 2023). To explain the formation of gas giant planets being discovered around low mass stars via photometry from TESS (e.g., Cañas et al., 2020; Kanodia et al., 2023; Bryant et al., 2024, Hotnisky et al., in prep.) and precise red-optical and near-infrared RV surveys (e.g., Morales et al., 2019; Quirrenbach et al., 2022) necessitates massive protoplanetary disks (Delamer et al., 2024). In such massive disks, formation via gravitational disk instability is an attractive alternative, requiring protostellar disks that undergo rapid gravitational collapse if they are massive and cool enough to satisfy the Toomre criterion (Toomre, 1964; Boss & Kanodia, 2023). As shown in Boss & Kanodia (2023), the disk-to-star mass ratio required for the protoplanetary disk to undergo collapse increases for lower stellar masses: $\sim$25% at 0.1 $M_{\odot}$ compared to $\sim$10% $\geq$0.3$M_{\odot}$, suggesting that more massive and rarer disks are required for gravitational collapse to occur around the lowest mass stars. For increasing mass ratios, the likelihood of the cloud fragmenting into separate filaments from which the formation of a star and a massive companion proceeds increases (e.g., Bate et al., 2002; Bate, 2008), representing a favorable mechanism to explain the highest mass ratio brown dwarfs around the lowest mass stars. Surveys such as the ongoing Searching for Giant Exoplanets around M Stars survey (Kanodia et al., 2024), and MANGOs (M dwarfs Accompanied by close-iN Giant Orbiters with SPECULOOS) (e.g., Triaud et al., 2023), along with new astrometric detections from Gaia such as those discussed here will shed further light into the occurrence and rarity of objects within this ‘Jupiter desert’.

To visualize the eccentricity distribution of massive planets and brown dwarfs comparable to Gaia-4b and Gaia-5b, Figure 6 shows the measured eccentricities of massive planets ($1M_{\mathrm{J}}<M<13M_{\mathrm{J}}$) and brown dwarfs ($>13M_{\mathrm{J}}$) with periods shorter than $10^{4}\,\mathrm{days}$ around low-mass stars. Gaia-5b has the most eccentric orbit of the brown dwarfs plotted in this figure. Given the orbital separations of both Gaia-4b and Gaia-5b, neither is likely to experience significant tidal circularization. As such, neither system requires a special explanation (such as an additional body) to maintain a high eccentricity. However, how did these objects obtain their eccentricities? Possible scenarios include companion-companion gravitational scattering (e.g., Rasio & Ford, 1996; Ford & Rasio, 2008; Petrovich & Tremaine, 2016), secular Kozai-Lidov perturbations with a massive outer companion (e.g., Naoz, 2016), or companion-disk interactions during formation (Goldreich & Sari, 2003). From the currently available data, we see no obvious evidence of additional massive companions in the systems. Even so, it still could be a possibility that other companions are present in the systems at distant orbits. To test for potential dynamical influence of a binary star, we checked the Gaia binary star catalog in El-Badry et al. (2021) for any co-moving stars, and we found no evidence for such stars.

From Figure 6, both the planet and the brown dwarf sample show similarly distributed tails of eccentric systems, while the planet sample shows a larger number of circular systems. This could hint at a two-component underlying distribution in the planet sample of a close-to an isotropic distribution in high eccentricities combined with a well-aligned component. However, as this is not a homogeneously selected sample, we did not attempt to further characterize the underlying distributions. We expect that in its all sky volume-limited survey, and due to the large number of expected detections of planets and brown dwarfs around low-mass stars, Gaia will enable the creation of such homogeneous samples to shed further light into both the eccentricity distribution between massive planets and brown dwarfs at intermediate orbits.

Including Gaia-4 and Gaia-5, we obtained observations of 28 Gaia-ASOI stars with HPF and NEID orbiting nearby M and K stars. Figure 7 shows broadening functions from HPF and CCFs from NEID of those systems. We highlight systems that are consistent with being single stars (i.e., showing only a single set of lines in any of the epochs), as well as systems that show clear evidence of double lines in the spectra. In all cases, when double lines are seen, they have approximately the same height in the BFs or CCFs, confirming that the primary false positive scenario is composed of double-lined nearly equal mass binaries that create small photocenter motions masquerading as planets, similar to the false positive scenario discussed in Marcussen & Albrecht (2023). We note that an observation of a single line is not sufficient to demonstrate that a system is definitely a single star, as it could be that the velocity separation and flux ratios would be such that at the epoch of the observation that the lines would overlap. However, as additional observations are obtained as a function of phase, the possible parameter space for a binary being simultaneously consistent with the astrometric data and the spectroscopic data rapidly shrink. In this work, we present the systems that still show single lines, and we leave the further analysis of them to future work. Additionally, in Figure 7 we highlight Gaia-ASOI-016 as a single-lined system that is likely an SB2 system, as although it is formally single lined, the broadening function profiles are very broad and show clear line-width variations suggestive of contaminating light from another star, although additional work would be needed to definitely separate out the two lines.

In Table 4 in the Appendix, along with Figure 1 in Section 2 we provide an overview of the targets observed and all of the designations (consistent with single, confirmed sub-stellar companion, double lined SB2s) from this work and the literature. Of the 28 systems we observed, 21 are double lined spectroscopic binary systems (SB2s), and six do not show evidence of double lines yet and are thus still compatible with being single stars and are promising candidates for hosting substellar companions consistent with the Gaia astrometric solutions. The last system is Gaia-ASOI-016, which is formally single lined but shows broad time-varying line profiles likely suggestive of an SB2, but we still formally count it in the group of ’single lined’ systems. For the promising single-lined systems, additional observations are needed to resolve the full RV curve to independently confirm the RV curve which can then be used to assess the validity of the expected astrometric orbit under the assumption of a single dark companion as we have done for Gaia-4b and Gaia-5b. Such observations are underway. Of the 72 original ASOI candidates, with the confirmation of Gaia-4b and Gaia-5b, in total, 12 are now confirmed to be substellar companions, of which 11 are exoplanets, where nine of the systems were already known RV detections before Gaia. One of the systems is the recently discovered Gaia-3b (Holl et al., 2023; Marcussen & Albrecht, 2023; Winn, 2022; Sozzetti et al., 2023) where the ASOI solution is withdrawn, but RV followup suggests a giant planet on an eccentric orbit, and the other two systems are Gaia-4b and Gaia-5b, discussed in this work.

From these observations we see that the false positive rate is between 23/72=32% and 60/72=83%, where the former only counts the confirmed SB2s across the full sample, and the latter assumes all non-confirmed exoplanets would be false positives. Restricting our sample of 28 to the 26 systems that did not have any spectroscopic observations before to the best of our knowledge¹⁰¹⁰10The two systems that did have previous observations are G-ASOI-001 (Gaia-3; single lined with likely planet) and G-ASOI-054 (confirmed SB2)., then the false positive rate is at minimum 20/26 = 77%, with the remainder being compatible with being single stars, but only 2/26=8% confirmed substellar companions.

Given the complex selection function steps that went into construction of the ASOI candidate sample, we do not attempt to interpret these numbers as reflective of astrophysical occurrence rates, which would need to wait for robust quantification of biases and sensitivity analysis only possible when the time-series astrometry comes available as part of Gaia DR4. In the meantime, these observations highlight the importance of spectroscopic observations in effectively and efficiently discerning between false positive scenarios. Additional RV observations are ongoing to quantify the agreement with the RV-only solution with the Gaia solution for the other systems that are still compatible with being single.

We note that three systems in the ASOI table, have had their astrometric solution retracted due to a software bug. This includes the polluted white dwarf system WD 0141-675 (Gaia DR3 4698424845771339520) which has been discussed in Rogers et al. (2024), that had a now withdrawn putative substellar candidate with a period of 33.65d. Gaia DR3 1712614124767394816 (Gaia-ASOI-1) is a system discussed by Holl et al. (2023); Unger et al. (2023); Winn (2022); Gaia Collaboration et al. (2023b); Marcussen & Albrecht (2023). This system was discussed in Winn (2022) and Marcussen & Albrecht (2023) to show a discrepancy in the RV semiamplitude. Most recently, Sozzetti et al. (2023) provided a detailed study in the system, showing that the system has a very high eccentricity of $e=0.95$. They provided simulations suggesting that the Gaia data indicated that the companion would be compatible with being an exoplanet. However, with the formal astrometric solution retracted, the impact on that conclusion has been made unclear. Finally, to the best of our knowledge, Gaia DR3 522135261462534528 (Gaia-ASOI-3) has not been observed spectroscopically.