A Noninteracting Galactic Black Hole Candidate in a Binary System with a Main-sequence Star

We retrieved UV-IR photometry for this source from GALEX (Martin et al. 2005), SDSS (Fukugita et al. 1996), the APASS survey (Henden et al. 2016), Gaia DR3 (Riello et al. 2021), PanSTARRS (Tonry et al. 2012), SkyMapper (Wolf et al. 2018), 2MASS (Skrutskie et al. 2006), and WISE (Wright et al. 2010). Figure 2 displays the available photometry for this source, and Table 1 lists the measurements and associated errors in each band. We analyze the spectral energy distribution of this source later on and discuss our choice of photometry.

Zoom In Zoom Out Reset image size

2.2.1. Observations and Data Reduction

We obtained follow-up high-resolution spectroscopy of Gaia DR3 4373465352415301632 to determine its stellar parameters, search for signatures of a second luminous component in the binary, and measure the radial velocity as a function of time. We observed the star for five epochs with the Automated Planet Finder (APF) telescope (Vogt et al. 2014) at Lick Observatory, four epochs with the DEIMOS spectrograph (Faber et al. 2003) on the Keck II telescope, and obtained a high signal-to-noise ratio spectrum with the MIKE spectrograph (Bernstein et al. 2003) on the Magellan/Clay telescope (as well as a second, shorter MIKE spectrum) at Las Campanas Observatory.

The APF observations employed the Levy Spectrometer (Radovan et al. 2010) with a 2'' × 8'' slit, providing R = 80,000 spectra covering 3730–10200 Å. We observed Gaia DR3 4373465352415301632 on 2022 June 28, 2022 August 28, 2022 September 2, 2022 September 12, and 2022 October 4, obtaining three 1000 s exposures per epoch. The APF spectra were reduced with the California Planet Search pipeline (Butler et al. 1996; Howard et al. 2010; Fulton et al. 2015).

We observed Gaia DR3 4373465352415301632 with Magellan/MIKE for 2700 s on 2022 August 22, using a 07 slit to obtain an R ≈ 30,000 spectrum from 3300 to 9400 Å. We reduced the MIKE spectrum with the Carnegie Python data reduction package (Kelson 2003). To provide a comparison spectrum, we also obtained MIKE observations of the G1.5V radial velocity standard star HD 126053. Portions of the spectra of these two stars are displayed in Figure 3. We also observed Gaia DR3 4373465352415301632 for 745 s with MIKE during twilight on 2022 October 13 to obtain an additional velocity measurement. Figure 4 compares the two MIKE spectra of Gaia DR3 4373465352415301632, which represent the radial velocity extremes in our data set.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We observed Gaia DR3 4373465352415301632 with DEIMOS on 2022 September 19, 2022 September 24, 2022 September 27, and 2022 October 28 using the 1200G grating and a 07 long slit to obtain R ≈ 7000 spectra from 6500–9000 Å. A single 300 s exposure was obtained in twilight on each night. We reduced the data with a modified version of the DEEP2 reduction pipeline (Cooper et al. 2012; Newman et al. 2013; Kirby et al. 2015). The long-slit mask employed for these observations did not contain enough slits to adequately constrain the quadratic correction to the sky line wavelengths as a function of position on the mask described by Kirby et al. (2015). Instead, we examined the night sky emission line wavelengths in the extracted spectrum, and used the measured wavelength shifts near the A-band and Ca triplet lines (which ranged from 0 to 4 km s⁻¹) to correct the wavelength scale.

2.2.2. Velocity Measurements

All of the archival photometry of Gaia DR3 4373465352415301632 that we were able to locate, including measurements from GALEX (Morrissey et al. 2007), SDSS DR16 Ahumada et al. 2020, APASS DR10 (Henden 2019), Pan-STARRS DR2 (Flewelling et al. 2020), SkyMapper DR2 (Onken et al. 2019), 2MASS (Cutri et al. 2003), and WISE (Cutri et al. 2021) is displayed in Figure 2 and listed in Table 1. As is clear, the SkyMapper v, i, z bands are discrepant with other photometric measurements at comparable wavelengths, and we therefore exclude them from our SED fits. Using the publicly available SED fitting code ARIADNE (Vines & Jenkins 2022), we fit the observed photometry (excluding SkyMapper and the wide Gaia bands) to derive stellar parameters. Although Gaia offers extremely precise optical photometry of Gaia DR3 4373465352415301632, we note that there may be contamination in the G_BP and G_RP bands as a result of the star's location in the direction of the Galactic bulge. The phot_bp_rp_excess_factor for this source is 1.24 and the G_BP − G_RP color is 1.17, which indicates that the blending probability is ∼30% (Riello et al. 2021).

The resultant posterior distributions are shown in the corner plot for the stellar parameters, along with the best-fit extinction in Figure 5. ARIADNE fits broadband photometry to a variety of stellar atmosphere models. In this case, the best fit is obtained with the Phoenix models (Husser et al. 2013). Here, we have used a prior on the distance that is consistent with the measured Gaia parallax for this source. Including the Gaia G-band magnitude in the fit does not significantly change any of the derived parameters. We obtain similar values for the effective temperature and other derived parameters when only fitting to GALEX photometry, along with B and V bands from APASS, SDSS g, r, i, z, 2MASS J, H, K_s , and WISE photometry. Similar posterior distributions are also obtained if Pan-STARRS photometry is used instead of SDSS. The best-fit ARIADNE parameters are: ${T}_{\mathrm{eff}}={5796}_{-153}^{+93}$ K, $\mathrm{log}g={4.26}_{-0.27}^{+0.24}$, $R={1.02}_{-0.029}^{+0.031}$ R_⊙, and ${\rm{[Fe/H]}}=-{0.24}_{-0.08}^{+0.09}$. As the corner plot shows, the posterior distributions are tightly constrained. The SED fit, along with the residuals, is shown in the bottom panel of Figure 5.

Zoom In Zoom Out Reset image size

We have also used the BaSeL model spectral library (Lejeune et al. 1997, 1998) to fit the observed photometry with both a combined χ² and genetic algorithm optimization and a Markov Chain Monte Carlo (MCMC) approach. Here we assume a Fitzpatrick (1999) extinction law with R_V = 3.1 and reddening E(B − V) = 0.3 mag. The best-fit SEDs are shown in Figure 6(a) for a single source. The fit assuming two sources is shown in Figure 6(b). In the single-source case, we obtain a reduced χ² of 3.7, best-fit parameters T_eff = 5927 K, $\mathrm{log}g\,=\,4.0$, L = 1.08 L_⊙, and metallicity of Z = 0.009. For a two source fit, we obtain a reduced χ² of 5.5, with T_eff = 6241 K, $\mathrm{log}g=4.99$, L = 0.42 L_⊙, and Z = 0.0006 for the hotter component and T_eff = 5763 K, $\mathrm{log}g\,=\,4.01$, L = 0.67 L_⊙, and Z = 0.01 for the cooler component. Although the raw χ² value does improve upon adding a second component (17 for the composite spectrum versus 25 for the single SED fit), this is not enough to compensate for the increase in the number of model parameters in the reduced χ². Therefore, we conclude that there is no significant evidence in the SED for the presence of more than one luminous component.

Zoom In Zoom Out Reset image size

3.2.1. Full Spectral Fitting

3.2.2. Fundamental and Photospheric Stellar Parameters

Next, we determined the stellar parameters of Gaia DR3 4373465352415301632 more rigorously using a combination of the classical spectroscopic approach ¹⁷ and a comparison with theoretical isochrones to infer accurate, precise, and self-consistent parameters.

The inputs to this parameter inference procedure include the equivalent widths of Fe i and Fe ii atomic absorption lines, as well as multiwavelength photometry from Gaia DR3 G-band (Riello et al. 2021) and 2MASS (J, H, and Ks), the Gaia DR3 parallax (Lindegren et al. 2021b), and we adopt the V-band extinction from the three-dimensional Stilism reddening map (Lallement et al. 2014; Capitanio et al. 2017; Lallement et al. 2018). The absorption line data are based on the line lists presented in Meléndez et al. (2014) and Reggiani & Meléndez (2018). We measured the equivalent widths by fitting Gaussian profiles with the Spectroscopy Made Harder (smhr) software package to the continuum-normalized MIKE spectrum. The continuum normalization was performed using smhr. We assume solar abundances from Asplund et al. (2021) and follow the steps described in Reggiani et al. (2022b, 2022a) to obtain the fundamental and photospheric stellar parameters from a combination of spectral information and a fit to the MESA Isochrones and Stellar Tracks (Choi et al. 2016) stellar models. We do the fitting with the isochrones package ¹⁸ (Morton 2015), which uses MultiNest ¹⁹ (Feroz & Hobson 2008; Feroz et al. 2009, 2019) via PyMultinest (Buchner et al. 2014).

Our derived fundamental and photospheric parameters for Gaia DR3 4373465352415301632 are listed in Table 3. Its position relative to the MIST isochrones is illustrated in Figure 7. The effective temperature is in excellent agreement with the values determined from the single-source SED fit, as well as the full spectral fitting. According to the effective temperature calibration given by Gray & Corbally (2009), the spectral classification would be G0. The surface gravity that is required to enforce ionization balance on the Fe lines is somewhat higher than determined via other techniques, but the disagreement is only at the 1.6σ level. The spectroscopic metallicity of Gaia DR3 4373465352415301632 is also modestly lower than the photometric value. The overall good agreement on the stellar parameters that is achieved through multiple independent methods indicates that systematic uncertainties, e.g., related to the choice of stellar evolution tracks, are smaller than the adopted uncertainties (Table 3). We conclude that Gaia DR3 4373465352415301632 is an old, moderately metal-poor main-sequence star, which is slightly less massive than the Sun.

Zoom In Zoom Out Reset image size

Table 3. Derived Stellar Parameters

Property	Value	Unit
Effective temperature (Teff)	5972 ± 100	K
Surface gravity ($\mathrm{log}g$)	4.54 ± 0.15	$\mathrm{log}{{\rm{cm}}\,{\rm{s}}}^{-2}$
Metallicity ([Fe/H])	−0.30 ± 0.10
Microturbulent velocity (ξ)	1.10 ± 0.10	km s−1
Stellar mass (M*)	0.91 ± 0.10	M⊙
Stellar radius (R*)	1.003 ± 0.075	R⊙
Luminosity (L*)	1.186 ± 0.232	L⊙
Isochrone-based age (τiso)	7.1 ± 3.8	Gyr

Note. The parameters that are listed here are determined via the combined analysis of the high-resolution stellar spectrum along with theoretical isochrones, as described in Section 3.2.2.

Download table as:  ASCIITypeset image

3.2.3. Chemical Abundances

We used the photospheric stellar parameters from Table 3 to calculate the abundances of C i, O i, Na i, Mg i, Al i, Si i, S i, K i, Ca i, Sc ii, Ti i, Ti ii, V i, Cr i, Cr ii, Mn i, Fe i, Fe ii, Co i, Ni i, Cu i, Zn i, Sr ii, Y ii, Zr ii, Ba ii, and Ce ii, from the equivalent widths (EWs) of spectroscopic absorption lines. We measured equivalent widths from the continuum-normalized MIKE spectrum by fitting Gaussian profiles with smhr. We use the one-dimensional plane-parallel, α-enhanced, ATLAS9 model atmospheres (Castelli & Kurucz 2003) and the 2019 version of MOOG (Sneden 1973) to infer elemental abundances based on each equivalent width measurement, while initially assuming local thermodynamic equilibrium (LTE). Our calculations include hyperfine structure splitting for Sc ii (based on the Kurucz line list ²⁰ ), V i (Prochaska et al. 2000), Mn i (Prochaska & McWilliam 2000; Blackwell-Whitehead et al. 2005), Cu i (Prochaska et al. 2000, and Kurucz line lists), and Y ii (Kurucz line lists), and we account for isotopic splitting for Ba ii (McWilliam 1998). We report the adopted atomic data, equivalent width measurements, and individual line-based abundances in Table 4. Our final abundances are reported in Table 5.

Table 4. Atomic Data, Equivalent Widths and Line Abundances

Wavelength	Species	Excitation Potential	log(gf)	EW	${\mathrm{log}}_{\epsilon }(X)$
(Å)		(eV)		(mÅ)
6154.225	Na i	2.102	−1.547	20.80	5.993
6160.747	Na i	2.104	−1.246	31.52	5.927
4571.095	Mg i	0.000	−5.623	89.51	7.395
4730.040	Mg i	4.340	−2.389	43.79	7.381
5711.088	Mg i	4.345	−1.729	80.01	7.244
6318.717	Mg i	5.108	−1.945	26.75	7.284
5260.387	Ca i	2.521	−1.719	23.01	6.139
5512.980	Ca i	2.933	−0.464	63.20	5.966

Note. Full version online. This table will be published in its entirety as a machine-readable table. A portion is shown here for guidance regarding its form and content.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Table 5. Chemical Abundances

Species	N	log(X )	[X/Fe]	σ[X/Fe]
C i	4	8.01	−0.15	0.09
O i	3	8.65	0.26	0.03
Na i	2	5.96	0.04	0.03
Mg i	4	7.33	0.08	0.04
Al i	5	6.08	−0.05	0.02
Si i	11	7.27	0.06	0.03
S i	4	6.85	0.03	0.01
K i	1	5.18	0.42	0.01
Ca i	10	6.06	0.06	0.02
Sc ii	4	2.94	0.10	0.09
Ti i	9	4.60	−0.07	0.03
Ti ii	10	4.80	0.13	0.04
V i	3	3.51	−0.09	0.02
Cr i	8	5.26	−0.06	0.02
Cr ii	5	5.24	−0.07	0.03
Mn i	6	4.88	−0.23	0.05
Fe i	79	7.16	⋯	⋯
Fe ii	18	7.15	⋯	⋯
Ni i	14	5.94	0.04	0.03
Cu i	2	3.79	−0.09	0.02
Zn i	2	4.25	−0.01	0.07
Sr i	1	2.41	−0.12	0.00
Y ii	3	1.82	−0.08	0.05
Zr ii	1	2.29	−0.00	0.03
Ba ii	3	2.12	0.15	0.05
Ce ii	1	1.55	0.27	0.03
non-LTE Corrected Abundances
C i	4	8.07	−0.04	0.16
O i	3	8.57	0.24	0.03
Na i	2	5.85	−0.01	0.04
Al i	2	6.07	0.00	0.01
Si i	11	7.25	0.10	0.08
K i	1	4.74	0.03	⋯
Fe i	79	7.18	⋯	⋯
Fe ii	18	7.22	⋯	⋯
Ba ii	3	1.96	0.06	0.03

Download table as:  ASCIITypeset image

When possible, we update our elemental abundances derived under the assumptions of LTE to account for departures from LTE (i.e., non-LTE corrections) by interpolating published grids of non-LTE corrections for several elements. We make use of three-dimensional non-LTE corrections for carbon and oxygen (Amarsi et al. 2019), and one-dimensional non-LTE corrections for sodium (Lind et al. 2011), aluminum (Amarsi et al. 2020), silicon (Amarsi & Asplund 2017), potassium (Reggiani et al. 2019), calcium (Amarsi et al. 2020), iron (Amarsi et al. 2016), and barium (Amarsi et al. 2020). The NLTE-corrected abundances are listed in Table 5.

The abundance analysis shows that Gaia DR3 4373465352415301632 is a mildly α-enhanced star ([α/Fe]∼0.1), as expected for its metallicity. The iron-peak and neutron-capture elemental abundances are also fully compatible with the typical values for this metallicity range, based on galactic chemical evolution models (e.g., Kobayashi et al. 2020) and large spectroscopic surveys (e.g., Amarsi et al. 2020). We do not see any sign of the large enhancements in light elements and α elements that have been detected in some X-ray binary companions (e.g., Orosz et al. 2001; González Hernández et al. 2008; Suárez-Andrés et al. 2015; Casares et al. 2017).

The Gaia DR3 non_single_stars catalog provides orbital parameters for Gaia DR3 4373465352415301632 based on a fit to the individual astrometric measurements obtained by Gaia. However, the astrometric measurements themselves will not be available until the fourth Gaia data release ∼3 yr from the time of this writing. The parameters listed directly in the catalog include the period P, the time of periastron T₀, the eccentricity e, the secondary mass M₂, and the Thiele-Innes (TI) elements A, B, F, and G. Following the method of Binnendijk (1960), as presented by Halbwachs et al. (2022), we use the TI elements to compute the semimajor axis a, the position angle of the ascending node Ω, the longitude of periastron ω, and the inclination i. The values of these parameters are listed in Table 6. We note that these calculations are based on the results reported in the two body orbit catalog and the binary mass catalog from Gaia Collaboration et al. (2022a).

Table 6. Orbital Solutions for Gaia DR3 4373465352415301632for the Data Set Listed in Table 2

Parameter	Value
Astrometric fit
P	185.77 ± 0.31 day
T0	2457377.0 ± 6.3 day
e	0.49 ± 0.07
a	2.98 ± 0.22 mas
i	1213 ± 3.5
Ω	896 ± 3.8
ω	−106 ± 11.9
M2	${12.8}_{-2.3}^{+2.8}$ M⊙
RV fit
P	${184.28}_{-0.89}^{+0.75}$ day
T0	2459802.4 ± 3.0 day
e	${0.411}_{-0.021}^{+0.034}$
ω	$19\buildrel{\circ}\over{.} {1}_{-6.6}^{+6.0}$
K1	${60.6}_{-5.8}^{+10.0}$ km s−1
γ	${47.8}_{-1.4}^{+2.2}$ km s−1
f	${3.2}_{-0.9}^{+1.6}$ M⊙
Joint fit
P	${185.41}_{-0.087}^{+0.101}$ day
T0	${2459757.9}_{-1.5}^{+1.2}$ day
e	0.456 ± .023
a	${1.357}_{-0.061}^{+0.069}$ au
i	123.6 ± 1.6
Ω	$96\buildrel{\circ}\over{.} {2}_{-2.0}^{+2.3}$
ω	$10\buildrel{\circ}\over{.} {4}_{-2.0}^{+2.7}$
M2	${11.33}_{-1.32}^{+1.57}$ M⊙

Download table as:  ASCIITypeset image

The TI elements, together with the eccentricity, describe an orbit in the plane of the sky. We generate 50 random realizations of measurements drawn from the Gaia best-fit values and covariance matrix. We then plot the sky paths given by the eccentricities and TI elements. Figure 8 shows these sky paths about the barycenter; the best-fit orbit reported by Gaia DR3 is shown as a thicker black line. The orbit has a semimajor axis of 2.98 ± 0.22 mas, with the major axis oriented nearly east–west. ²¹

Zoom In Zoom Out Reset image size

We fit the radial velocity data using the adaptation of orvara(Brandt et al. 2021) described in Lipartito et al. (2021). This remains an MCMC-based fitting approach but differs from the main version of orvarain several ways:

• 1.  
  It fits only for the RV parameters of period, eccentricity, RV semi-amplitude, periastron time, argument of periastron, and RV zero-point;
• 2.  
  The argument of periastron refers to the star (as is the definition for Gaia astrometry); and
• 3.  
  Only the RV zero-point (barycenter RV of the system) is analytically marginalized off.

We use uniform priors on all of the RV parameters: period, semi-amplitude, eccentricity, periastron time, RV zero-point, and (stellar) argument of periastron. We find a highly multimodal posterior, with orbital solutions at periods of ∼142 days, ∼164 days, ∼184 days, ∼220 days, and additional possible periods out to at least 1000 days. The top panel of Figure 9 shows the wide range of RV orbits consistent with our Gaia DR3 4373465352415301632 velocities. The fits to the data are satisfactory, with a best-fit χ² of 3 for 13 data points with six free parameters, i.e., seven degrees of freedom.

Zoom In Zoom Out Reset image size

Because one of our main goals is to determine whether the Gaia astrometric orbital solution is correct, we restrict the subsequent analysis of our RV data to solutions with a period range within ±10 days of the Gaia orbital period of ∼185 days, discarding those with shorter or longer periods. A similar procedure was adopted by El-Badry et al. (2023b) in their analysis because their RV data (although more extensive than ours) also yields multiple period solutions. However, as we note below, the combination of our RV data and that of El-Badry et al. (2023b) gives a unimodal solution for the period, which is not possible using either data set alone. This result is illustrated in Figure 9b, where the combined RV data uniquely determine the orbital period without reference to the Gaia astrometric solution. The fit to the combined RV data set remains formally good, with a best-fit χ² of 24 for 52 total radial velocity measurements and six free parameters (46 degrees of freedom). The low χ² value suggests that, if anything, the RV errors and, by extension, our derived parameter uncertainties, may be overestimated.

The spectroscopic orbital parameters from the MCMC fit to our RV data set (the velocities listed in Table 2), restricted to periods between 175 and 195 days, are listed in the second section of Table 6. Figure 10 shows the corner plot of the posterior distribution. Using all of the available RVs gives the posterior distributions shown in Figure 11 and the best-fit parameters listed in the top section of Table 7, whereby the combination of both data sets yields very precise constraints on all parameters that can be determined with RV data. The radial velocity curve and the phase-folded radial velocity curve with residuals are displayed in Figure 12 for both our RV data alone and all available RV data. Although our RV data allow a wide range of possible orbits, especially at phases near 0 and 1, the inclusion of data from El-Badry et al. (2023b) restricts the range of possible orbits significantly. Just 0.4% of the points in the chain lie away from the mode at P ≈ 185 days, and are in a much sparser cluster at P ≈ 218 days. We obtain very similar results using different orbit-fitting packages such as TheJoker (Price-Whelan et al. 2017) or exoplanet (Foreman-Mackey et al. 2021a, 2021b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 7. Orbital Solutions for Gaia DR3 4373465352415301632 Using All Available Data

Parameter	Value
RV fit
P	${184.52}_{-0.80}^{+0.66}$ day
T0	2459801.16 ± 0.12 day
e	${0.4368}_{-0.0037}^{+0.0036}$
ω	1610 ± 0.63
K1	${65.58}_{-0.30}^{+0.31}$ km s−1
γ	${48.48}_{-0.36}^{+0.36}$ km s−1
f	${3.92}_{-0.042}^{+0.042}$ M⊙
Joint fit
P	${185.52}_{-0.08}^{+0.08}$ day
T0	2459759.3 ± 1.1 day
e	${0.439}_{-0.003}^{+0.003}$
a	${1.258}_{-0.010}^{+0.011}$ au
i	$126\buildrel{\circ}\over{.} {80}_{-0.63}^{+0.62}$
Ω	$98\buildrel{\circ}\over{.} {71}_{-2.01}^{+1.97}$
ω	1589 ± 0.61
M2	${9.326}_{-0.208}^{+0.216}$ M⊙

Download table as:  ASCIITypeset image

The mass function derived from our period-restricted spectroscopic orbital solution is $f={3.2}_{-0.9}^{+1.6}$ M_⊙. When combined with the primary mass determined from the spectral analysis, the corresponding minimum secondary mass is 4.6 M_⊙, which is consistent with the companion being a black hole. However, because our spectroscopic measurements missed the maximum radial velocity at the periastron of the orbit, our constraint on the velocity semi-amplitude K₁ is not very strong, as is apparent from the definition of the mass function:

$$\begin{eqnarray}&&f=\displaystyle \frac{{{PK}}_{1}^{3}}{2\pi G}{\left(1-{e}^{2}\right)}^{3/2}=\displaystyle \frac{{M}_{2}^{3}{\sin }^{3}i}{{\left({M}_{1}+{M}_{2}\right)}^{2},}\end{eqnarray} \tag{ 1 }$$

where P is the orbital period, and M₁ and M₂ are the masses of the two stars. The cubic dependence of the mass function on K₁ then results in a large uncertainty in the mass function and the mass of the secondary, such that with our radial velocities alone we cannot completely exclude masses below 4 M_⊙. The derived constraints on the spectroscopic mass function improve significantly when using all available RV data, while combining our measurements and those from El-Badry et al. (2023b). In this case, we obtain f = 3.9 ± 0.04 M_⊙, which implies a secondary mass of 9.2 M_⊙ (when the astrometric inclination is included).

Our joint fits rely on determining the set of orbits that simultaneously satisfy constraints from both the RV data and the Gaia astrometry. We begin with the MCMC chains derived from RVs alone and then condition these chains on the Gaia astrometry. We achieve this by first choosing random values for parameters unconstrained by RV: parallax, inclination, and position angle of the ascending node. With a random set of these parameters for each step of the chain, we can compute the TI constants and use the Gaia DR3 covariance matrix to compute a likelihood. We then re-weight each step of the chain by its corresponding likelihood. In practice, we improve our signal-to-noise ratio by choosing many possible inclinations, position angles, and parallaxes for each step of the RV chain, and computing the Gaia likelihood for each of them.

Our conditioning of the RV fits on the Gaia covariance matrix immediately reveals a problem. If we adopt the full RV data set—the union of our data and that of El-Badry et al. (2023b)—, then the goodness-of-fit is formally good for the RVs (χ² = 24 for 46 degrees of freedom). The step in this chain that agrees best with Gaia, even when freely choosing inclination, position angle, and parallax, has a χ² with respect to Gaia of 11 with five degrees of freedom. This step is, in turn, a much worse fit to the RVs, with ${\rm{\Delta }}{\chi }_{\mathrm{RV}}^{2}\approx 12$. The Gaia constraint provides an additional five degrees of freedom but increases the total χ² by about 20. This would be even worse if the RV errors are underestimated, as suggested by the low χ² for the RV-only fit.

The poor formal agreement of Gaia astrometry with the RVs renders any resulting confidence intervals dubious. The increases in χ² from conditioning on Gaia are about a factor of 4 larger what is expected; this suggests that the Gaia uncertainties are underestimated by a factor of ≈2. Figure 13 shows this tension in the five parameters constrained both by astrometry and radial velocities; the RV confidence intervals are from the entire RV data set. The 1σ confidence intervals do not overlap in any parameters apart from eccentricity. (We emphasize, though, that both the astrometry and RV data sets independently point to Gaia DR3 4373465352415301632 containing a ∼10 M_⊙ black hole with a ∼185 day period. The disagreement that is discussed here only relates to determining the exact orbital parameters of the system.)

Zoom In Zoom Out Reset image size

Based on the combined evidence from a satisfactory fit to the RV data and an unsatisfactory joint fit, we inflate the Gaia uncertainties by a factor of 2. We retain the form of the Gaia covariance matrix (i.e., the correlation coefficients between parameters) and simply multiply this matrix by a factor of 4. This increase in the uncertainties leads to acceptable joint fits, with χ² increasing by ≈8 when optimizing the best-fit orbit with the additional five degrees of freedom provided by astrometry.

We therefore adopt the following prescription to condition the RV fits on the Gaia astrometry:

• 1.  
  We select many possible sets of random values for parallax, position angle, and inclination for each step in an RV chain;
• 2.  
  We compute the TI constants and evaluate

  $$\begin{eqnarray}&&{\chi }^{2}={({\boldsymbol{p}}-\tilde{{\boldsymbol{p}}})}^{T}{{\boldsymbol{C}}}_{\mathrm{Gaia}}^{-1}({\boldsymbol{p}}-\tilde{{\boldsymbol{p}}}),\end{eqnarray} \tag{ 2 }$$

  where p is the set of parameters from the RV chain (including the additional parameters above), $\tilde{{\boldsymbol{p}}}$ are the best-fit Gaia values, and C _Gaia is the Gaia covariance matrix; and
• 3.  
  We weight each of these steps by ${e}^{(-{\chi }^{2}/8)}$, equivalent to multiplying the Gaia covariance matrix by 4 (or doubling the uncertainties). This produces a weighted chain.

Figure 14 presents another way of visualizing the resulting combined constraints and the tension of Gaia with the RVs. The top panel shows a random sampling of the astrometric orbits that are compatible with both the union of the RV data sets and the Gaia covariance matrix. It is more tightly constrained then the orbits shown in Figure 8 but less so than those shown in El-Badry et al. (2023b) due to our inflation of the Gaia uncertainties. We note that El-Badry et al. (2023b) did perform a check with inflated uncertainties but used the published Gaia covariance matrix for their baseline case. The lower two panels of Figure 14 show the tension between the astrometry and RVs. The black lines, denoting the steps in the RV chains that best match the Gaia astrometry, leave systematics in the RV residuals. This is true both when using only our RVs (middle panel) and using all RVs (lower panel). The light-blue lines are randomly drawn from the weighted chains computed as described earlier. These RV residuals can be compared to those shown in Figure 12 for the best fit to the RVs alone.

Zoom In Zoom Out Reset image size

Our joint constraints for Gaia DR3 4373465352415301632 using both our radial velocity data and the astrometry are listed in Table 6 (which also provides the astrometric fit and posteriors that result from fitting to the velocities only). Figure 15 displays the the resultant corner plot for the derived parameters of the combined fit using the RV data listed in Table 2. We find that the secondary mass of Gaia DR3 4373465352415301632 is ${11.64}_{-1.31}^{+1.52}$ M_⊙. Using all of the available radial velocity data gives the posterior distributions shown in Figure 16 and the parameters listed in the bottom section of Table 7, leading to a smaller secondary mass of ${9.326}_{-0.209}^{+0.216}$ M_⊙, which is in closer agreement with El-Badry et al. (2023b).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Prior to Gaia DR3, several studies calculated the expected population of black holes with luminous companions that could be detected with Gaia. The predicted Gaia yields range from tens to thousands of such binaries by the end of the mission (Breivik et al. 2017; Mashian & Loeb 2017; Wiktorowicz et al. 2020; Chawla et al. 2022; Janssens et al. 2022; Shikauchi et al. 2022; Wang et al. 2022). Although the quality cuts on orbital solutions described by Halbwachs et al. (2022) may have removed more binary systems that contain black holes, the identification of only two confident black hole-luminous companion binaries in DR3 so far (see El-Badry et al. 2023a, as well as the earlier identification of this source by Tanikawa et al. (2023) in their analysis of Gaia data) suggests that the upper end of those predictions is likely to be too optimistic. The properties of Gaia DR3 4373465352415301632 are quite similar to the most common binary parameters predicted by, e.g., Breivik et al. (2017), Yalinewich et al. (2018), and Shao & Li (2019) in black hole mass, companion mass, period, semimajor axis, and apparent magnitude. The only deviation from those expectations is the eccentricity, where the observed value differs from the most likely values of e < 0.1 or (in the case of significant kick velocities) e > 0.9. In this section, we speculate that the observed eccentricity of e ∼ 0.4 may have resulted from a possible architecture of a hierarchical triple system. Stellar triples are a common architecture and can affect the dynamical evolution of the system (Tokovinin et al. 2006; Tokovinin & Moe 2020; Chakrabarti et al. 2022; Lennon et al. 2022; Vigna-Gómez et al. 2022).

Although the characteristics of the Gaia DR3 4373465352415301632 binary system do not appear to be unusual with respect to population synthesis predictions for binaries detectable by Gaia, the particular combination of properties that we have determined may require fine-tuning to explain. To begin, the mass ratio of the system is quite large, M₂/M₁ = 10.3 ± 1.2. Although some low-mass companions to massive stars have been identified (e.g., GRAVITY Collaboration et al. 2018; Reggiani et al. 2022c), this does not seem to be the most common configuration for OB stars. The initial mass ratio of the system must have been even larger; for example, the initial-final mass relation of Raithel & Sukhbold (2018) suggests that a ∼10 M_⊙ black hole would have had a zero-age main-sequence mass of ∼15–50 M_⊙. Moreover, it is not clear that a binary with a separation of ∼1 au should have survived the evolution of the massive component. Prior to the formation of the black hole, the orbit must have been tighter because of the mass that the system lost when the black hole formed. (If the black hole received a natal kick, then that also would have heated the system, further reducing the inferred initial separation.) The maximum radius reached by a 15–50 M_⊙ star after crossing the Hertzsprung gap is ≳400 R_⊙(Martins & Palacios 2013), indicating that the system could not have avoided a common envelope phase unless the progenitor evolved quasi-homogeneously (e.g., Ramachandran et al. 2019; Gilkis et al. 2021; El-Badry et al. 2022b, suggesting that stars more massive than ∼30 M_⊙ may evolve quasi-homogeneously without respect to rotation.) Depending on the duration of that period and on the common envelope ejection efficiency, a merger might be expected to follow, which makes the survival of the G star until the present surprising.

One way to avoid this conundrum is if the system was initially a triple, consisting of an inner binary orbited by the G star. Massive stars are often found with close massive companions (Sana et al. 2012), so there would be nothing unusual about this architecture. Interactions between the two inner stars could then have proceeded either via mass loss or a merger in such a way as to prevent either star from expanding enough to engulf their wider companion. The question posed by this scenario is whether the inner binary ultimately merged, either before or after the formation of the black hole. If so, then that merger could have resulted in the production of the 9.3 M_⊙ black hole that we have identified. Meanwhile, if the stars did not merge, then an inner binary of two compact objects with a combined mass of 9.3 M_⊙ could still be present. Many such combinations are possible, and of course this suggestion is entirely speculative at the moment, but a binary with a combined mass of 9.3 M_⊙ would likely involve at least one object residing in the mass gap (e.g., Ö 1918), which would be of considerable interest. In the next section, we discuss whether there would be any observable consequences of the 9.3 M_⊙ dark mass being subdivided into a binary.

We have investigated with orvara the possibility that Gaia DR3 4373465352415301632 may be a triple system containing two unseen companions, but our present data set does not allow us to constrain the posterior distributions of a hierarchical triple. Extended radial velocity monitoring could make it possible to study the rich dynamics of a triple system, and we briefly speculate here on some of the possibilities that could be directly observed. Although the incidence of triples drops dramatically for periods greater than ∼10 days (Tokovinin et al. 2006; Tokovinin & Moe 2020), compact triple systems are not uncommon, and there is indeed a population of compact triples that manifest large eclipse transit variations in eclipsing binaries (Borkovits et al. 2015). Here, we focus on the more likely possibility of a compact inner binary with the G star as part of the outer system (Section 5.2). A similar scenario has been investigated in El-Badry et al. (2023b).

If this indeed is a hierarchical triple system, then the Kozai–Lidov mechanism would allow for an exchange between the eccentricity and the inclination (Chang 2009; Naoz et al. 2013; Suzuki et al. 2019). The Kozai–Lidov timescale (Antognini 2015), t_KL, is given by:

$$\begin{eqnarray}&&{t}_{\mathrm{KL}}\sim {P}_{\mathrm{in}}\displaystyle \frac{{m}_{1}+{m}_{2}}{{m}_{3}}{\left(\displaystyle \frac{{a}_{\mathrm{out}}}{{a}_{\mathrm{in}}}\right)}^{3}{\left(1-{e}_{\mathrm{out}}^{2}\right)}^{3/2},\end{eqnarray} \tag{ 3 }$$

where P_in is the orbital period of the inner binary, a_out and a_in are the semimajor axes of the outer and inner systems, e_out is the eccentricity of the outer system, and the masses of the three bodies are given by m₁, m₂, and m₃. For our specific case here, if we consider an inner binary of combined mass m₁ + m₂ ∼ 9.3 M_⊙ that has a third body with a semimajor axis that is 3(10) times that of the inner binary, then this gives a Kozai–Lidov timescale of ∼100(5000) yr. However, the shift in the velocity is ∼10⁻²(10⁻⁴) per orbital period. For this specific case, the velocity change in physical units translates to ∼1000(10) m s⁻¹ per orbit, which could in principle be observed, especially if the inner binary is not too close. This effect (the osculation of the orbit of the inner binary due to the Kozai–Lidov effect) is modified by a general relativistic correction term (Anderson et al. 2017). This is likely to be the dominant dynamical effect beyond that of Keplerian motion. Other dynamical effects for binary systems recently reviewed by Chakrabarti et al. (2022; including tidal effects and the general relativistic precession) are subdominant to this signal. Thus, the effects of the Kozai–Lidov effect may become manifest over several orbital periods for this system, even for a widely separated hierarchical system. In this context, it is advantageous that the companion star is a G star, with spectra that indicate that it has low stellar jitter (Wright 2005). For stars with low stellar jitter, high precision measurements should be able to capture this effect.

Earlier work on the report of the black hole binary system with a noninteracting giant companion by Thompson et al. (2019) indicated an unusually high [C/N] abundance ratio for its mass and evolutionary state. If the mass inferred for the red giant is correct, then this abundance pattern may have resulted from previous interaction with the black hole progenitor star (we note that van den Heuvel & Tauris (2020) argue that if the red giant's mass is ∼1M_⊙, then the companion could be a binary with two main-sequence stars). Low-mass X-ray binaries with black hole companions are also found to be metal-rich, with large abundances of α elements (Casares et al. 2017). For Gaia DR3 4373465352415301632, a comparison with synthetic spectra indicates an approximately solar [C/N] ratio, as expected for a main-sequence star of its metallicity. The [C/N] ratio of Gaia DR3 4373465352415301632 places it in the typical observed locus of main-sequence stars (Pinsonneault et al. 2018), which is consistent with its small Roche lobe to semimajor axis ratio. We find that the Roche limit to the semimajor axis for the G star is 0.2, which indicates that it is stable against tidal disintegration.

Given the current uncertainty in the formulation of the common envelope channel, we adopt a simpler approach and estimate the number of such systems that may exist in the galaxy. Here, we do not consider binary (and triple) evolution effects, which may be more significant than IMF dependencies. We assume a Salpeter IMF:

$$\begin{eqnarray}&&\displaystyle \frac{{dn}}{{dm}}={ \mathcal N }{m}^{-2.35},\end{eqnarray} \tag{ 4 }$$

where dn/dm is the number of stars per mass bin, and the normalization, ${ \mathcal N }$, can be set by integrating the number of stars per mass bin over the limits of the masses of stars in the Salpeter IMF, i.e., where ${m}_{\min }=0.1\,{M}_{\odot }$ and ${m}_{\max }=100\,{M}_{\odot }$ are the upper and lower limits, and equating to the observed number of stars in the galaxy, N_*. The binary fraction of high-mass stars is nearly unity (f_b ∼ 1) (Sana et al. 2009; Maíz Apellániz et al. 2016) and the companion masses, m_c , are taken from a flat distribution in mass ratio, q = m_c /m, e.g., dP/dq = 1 (Sana et al. 2012; Chulkov 2021).

$$\begin{eqnarray}\begin{array}{rcl}N & = & {f}_{b}{ \mathcal N }{\displaystyle \int }_{{m}_{\mathrm{lower}}}^{{m}_{\max }}\displaystyle \frac{{dn}}{{dm}}{dm}{\displaystyle \int }_{{q}_{1}}^{{q}_{2}}\displaystyle \frac{{dP}}{{dq}}{dq}\\ & = & 0.42{f}_{b}{ \mathcal N }\left({m}_{\mathrm{lower}}^{-2.35}-{m}_{\max }^{-2.35}\right)\left({m}_{2}-{m}_{1}\right)\end{array}\end{eqnarray} \tag{ 5 }$$

where m_lower is the minimum progenitor mass required to form a ∼10 M_⊙ black hole, which we take to be ∼20 M_⊙ (Sukhbold et al. 2016); and m₁ and m₂ are the lower and upper mass range of zero-age main-sequence masses of G-type stars, which we take to be 0.8 M_⊙ and 1.2 M_⊙, respectively. This gives a total number of binaries of ∼8 × 10⁵, which is comparable to the number of binaries with a black hole and a luminous companion in the zero-kick model by Breivik et al. (2017).

We integrate the orbit of Gaia DR3 4373465352415301632 backwards in time 500 Myr in the Galactic potential, as shown in Figure 17. Here, we use the potential derived from accelerations measured directly from pulsar timing (Chakrabarti et al. 2021), which provides the most direct probe of the Galactic mass distribution. Using potentials derived from kinematic assumptions (e.g., Bovy 2015) produces similar results to within a factor of ∼2. The star remains confined within a few hundred pc of the Galactic plane at all times, confirming that it formed as part of the thin disk population. Unless the natal kick of the black hole happened to be oriented within the disk of the galaxy, then the low scale height of the orbit indicates that the kick velocity must have been small. Given the low scale height of the orbit, ∼0.4 kpc, the vertical kick velocity is of order ∼10 km s⁻¹. This is similar to the negligible kick velocities for Cyg X-1 and VFTS 243 (Mirabel & Rodrigues 2003; Shenar et al. 2022). The orbit of Gaia DR3 4373465352415301632 around the Milky Way is somewhat eccentric, with a pericenter of ∼7 kpc and apocenter of ∼12 kpc, reaching substantially larger distances than the Sun does, which is consistent with its older age.

Zoom In Zoom Out Reset image size