Revealing the companion of NGC 1850-BH1

Introduction

A significant fraction of high mass O- and B-type stars are born in binaries and higher-order multiples (e.g. Sana and Evans [1]; Moe and Di Stefano [2]; Offner et al. 2023), and interactions between the components are thought to affect a majority of them throughout their lifetime (e.g. Sana et al. 2012). This can shape the physical and chemical evolution of one or both components, through the exchange of mass and angular momentum (e.g. de Mink et al. 2013).

The most massive among these stars are likely to end their lives as black holes (BHs) or neutron stars (NSs) following core-collapse supernovae. If the binary survives the explosion, the compact object may remain in a bound system with a stellar companion. If the compact object is a BH and X-ray quiet, it is commonly referred to as a “dormant Black Hole” (dBH). Further evolution of the system may lead to an X-ray binary, a second detonation, and/or a binary compact merger (e.g. Abbott et al. 2016).

This intermediate phase when the compact object is bound to a luminous stellar companion is thought to be comparatively long-lasting, and thus models suggest the presence of $O(10^{5})$ of such systems in the Milky Way Galaxy (e.g. Breivik et al. 2017; Langer et al. 2020; Nagarajan et al. 2025; Wagg et al. 2026). Their observational signatures are, however, challenging to detect, since the compact object must be inferred from the motion it induces on the luminous companion, and systems containing two non-compact stars can produce similar observables (e.g. Bodensteiner et al. 2020; El-Badry et al. 2022; El-Badry and Quataert 2020; Shenar et al. 2020). These “imposter” systems contain two stellar non-compact objects with one of them evolving though a short-lived bloated phase following substantial mass loss. In this phase, the star appears much brighter than expected given its mass and can outshine its companion, leading to the binary being incorrectly identified as a BH-star or NS-star system.

Despite ongoing efforts, only a few dBHs have been discovered in the Galactic field (e.g. El-Badry et al. 2023a,b; Gaia Collaboration 2024) and globular clusters (Giesers et al. 2018, 2019; Whitaker et al. 2026). Among massive-star binaries, the only known OB+dBH systems VFTS 243 in the LMC (Shenar et al. 2022) and HD 130298 in the Milky Way (Mahy et al. 2022) demonstrate that such systems can be identified through spectroscopic and orbital signatures.

For the Magellanic Clouds, the smaller overall stellar mass implies a lower total number of dBHs, although their lower metallicity leads to weaker stellar winds and hence favours the formation of more massive BH remnants (e.g. Mapelli et al. 2013). Additionally, the Magellanic Clouds host a substantial population of young massive clusters, which may provide favourable formation sites for dBHs, e.g. via dynamical interactions. Interestingly, the most massive dBH known in the Milky Way, Gaia-BH3 (33 $M_\odot$, Gaia Collaboration 2024), likely formed via a dynamical interaction in the low-metallicity (and now dissolved) cluster ED-2 (Marín Pina et al. 2024), and thus poses an example of the massive BH population that can result from dynamical low-metallicity environments. Thus, the presence of young massive clusters, and the lower metallicity in the Clouds make them promising targets for dBH searches.

NGC 1850 – BH1 was initially proposed by Saracino et al. (2022) as a system containing a main sequence turn-off primary orbiting an unseen secondary. Situated in the LMC cluster NGC 1850, the object’s location in the CMD and the cluster’s well-constrained age placed initial estimates of the primary mass at around $5\,M_\odot$, with simultaneous modelling of the radial velocities (RVs) and lightcurve suggesting a secondary mass around $11\,M_\odot$. A non-compact object of this mass would be expected to vastly outshine the proposed $5\,M_\odot$ primary, and due to the cluster’s age, no stars of $11\,M_\odot$ could be present in NGC 1850. Thus, the companion was labelled a BH.

Soon thereafter El-Badry and Burdge (2022) argued that the system’s tight orbit of around 5 days makes the assumption of a normal single-star evolutionary track questionable. As the primary mass inferred by Saracino et al. (2022) was based on the system’s location in the CMD, a divergence from single-star evolution could also affect the proposed mass of the primary, and thus the inferred companion mass. El-Badry and Burdge (2022) suggested constraining the mass via the star’s Roche Lobe (RL) geometry and orbit, deriving its mean density, and thus a mass of about $1\,M_\odot$. This implies a much lower mass companion, which could feasibly be hidden in the light of the puffed up primary. They further note that such a star could viably be the result of stripping via binary mass transfer, observed in the aforementioned bloated evolutionary stage, after which the primary is expected to contract and heat up again. Stevance et al. (2022) used BPASS binary population models to predict that a cluster of the mass and age of NGC 1850 should contain a few of such post-interaction systems, though none of the predicted systems matched the turn-off location of NGC1850-BH1.

In a subsequent work, Saracino et al. (2023) updated the optically derived RVs, thus increasing the implied secondary mass. Further, they placed an upper limit of 10% on the flux contribution of the unseen companion in the MUSE wavelength range, based on spectral disentangling of the MUSE data around the Paschen series. The analysis suggested that even in the case of a substantially stripped primary star, as proposed by El-Badry and Burdge (2022), the implied secondary mass of $\sim5\,M_\odot$ necessitates either a luminous object hidden in an optically thick disk, or a compact companion.

Here, we apply the spectral disentangling technique discussed in Seebergu et al. (2024), based on Simon and Sturm (1994), to analyse new multi-epoch FLAMES observations of the system in question. These data have a higher spectral resolution than the existing MUSE data and reach further into the blue optical wavelength regime, where the contribution of a luminous companion is expected to be stronger. Using these data, we aim to constrain the nature of NGC1850-BH1 and propose a potential evolutionary history.

Observations and data reduction

FLAMES-ARGUS spectroscopy

NGC 1850-BH1 was observed with the FLAMES-GIRAFFE spectrograph (Pasquini et al. 2002) installed at UT2 of the ESO Very Large Telescope (VLT) in programme gram 113.26S5.001 (PI: Saracino). The ARGUS integral-field unit (IFU), providing continuous spatial coverage across a field of view of $12 \times7$ arcsec with a sampling of $0.52$ arcsec/spaxel, was centred on NGC 1850-BH1 to enable us to deblend its spectrum from contributions by nearby stars during the data analysis (cf., § 3 below). Observations were carried out using the LR02 grating, which covers the wavelength range from 3964 Å to 4567 Å at a spectral resolution of $R=8\,500$. In Table 1, we provide an overview of the observations, carried out across 9 individual nights between October and December 2024. Within each night, a single observing block (OB) of two exposures of 1 350 s exposure time were recorded.

Date	OB start	OB end	Grade
4-Oct-2024	08:13:14	09:13:50	A
2-Nov-2024	06:42:40	07:43:44	A
24-Nov-2024	05:45:19	06:44:34	A
29-Nov-2024	05:36:47	06:32:52	A
1-Dec-2024	05:23:52	06:23:59	A
3-Dec-2024	03:58:17	04:55:11	A
4-Dec-2024	06:23:51	07:21:14	B
6-Dec-2024	05:38:30	06:37:39	A
25-Dec-2024	03:59:43	05:00:17	A

We reduced the data using P3D (Sandin et al. 2010), a multi-instrument software package for reducing the data of fibre-fed integral-field spectrographs written in idl. P3D performs the different reduction steps (bias subtraction, fibre tracing & extraction, wavelength calibration, and flat-fielding) on the individual exposures to create a set of raw-stacked spectra (RSS) files. Each RSS file contains a 1D spectrum per spaxel, with all spectra having been resampled to the same wavelength solution. To clean our science exposures from cosmic rays, we utilized the PyCosmic routine (Husemann et al. 2012), which is available within P3D.

In Figure 1, we show a reconstructed 2D image from one reduced exposure. To create the image, we integrated over the central part of the wavelength range and assigned each spaxel to the $(x, y)$-location in the focal plane indicated in the ARGUS calibration files available from ESO. The four corners of the ARGUS field of view do not contain fibres. Furthermore, there are three spaxels for which the signal is not recorded by the CCD. Those are visible to the top right of Figure 1.

HST Photometry

The Hubble space telescope (HST) photometry used in this work was published by Kamann et al. (2023). It is based on archival WFC3-UVIS observations collected during programs 14069 (PI: Bastian) and 14174 (PI: Goudfrooij) and includes exposures in the filters F275W, F336W, F343N, F438W, F656N, and F814W. The data were analysed with Dolphot (Dolphin 2016) following the approach described in Gossage et al. (2019) to create a single photometric catalogue. We matched the final catalogue to Gaia data release 3 (Gaia Collaboration 2023), resulting in the recovery of few bright sources that were missing in the HST photometry because they were saturated in all available exposures.

Following the publication of the catalogue in Kamann et al. (2023), Saracino et al. (2023) complemented it with WFC3-IR imaging from program 11193 (PI: MacKenty), resulting in near-infrared photometry in filters F105W, F110W, F125W, F140W, and F160W. The image analysis was analogous to the WFC-UVIS data.

Data analysis

Spectrum extraction

To extract stellar spectra from the individual ARGUS exposures, we used PampelMUSE (Kamann et al. (2013)), a software designed to deblend individual stellar spectra from integral-field spectroscopic observations of crowded stellar fields. PampelMUSE requires a reference catalogue of sources within the field of view, which it uses to identify the sources bright and isolated enough to extract their spectra. In our case, we used the HST catalogue described in § 2.2 above. To estimate the brightnesses of the stars in the ARGUS data, we adopted the F438W magnitudes available in the catalogue, given the filter’s good overlap with the wavelength range of our spectroscopy. In Figure 1, we indicate the locations of the sources in the HST catalogue relative to the ARGUS field of view as blue circles. Sources for which spectra were extracted are highlighted in orange, while the location of NGC 1850-BH1 is indicated by a red square.

Following the source selection stage, PampelMUSE determines the wavelength-dependent point-spread function (PSF) of the ARGUS data and the wavelength-dependent coordinate transformation from the catalogue to the IFU. The former is modelled using a Moffat profile, while the latter uses a six-parameter transformation accounting for any translation, scaling, rotation, and global distortion along either of the spatial axes. With this information in hand, the last step of the analysis is the extraction of the individual spectra.

As discussed in Saracino et al. (2022), NGC 1850-BH1 has a neighbour that is only $\sim0.5$ mag fainter, visible as an orange circle to the bottom right in Figure 1. While we were able to deblend the spectra of the two stars in the ARGUS data, there is a second neighbour, which is about 1.7 mag fainter than NGC 1850-BH1 (visible as a blue circle to the bottom left). This star is too faint and too close to NGC 1850-BH1 to be resolved in the ARGUS data.¹ Hence, it is likely that some of the flux from this faint neighbour will contaminate the extracted spectra of NGC 1850-BH1. We will discuss any potential impacts of such contamination on our results in § 5 below.

Following the extraction, we combined the spectra of NGC 1850-BH1 from the two exposures performed within each OB. This resulted in a final sample of nine spectra that were processed further.

Spectral disentangling

To establish the nature of a potential luminous secondary companion (or place upper limits on its maximum flux), spectral disentangling (Simon and Sturm 1994; Shenar et al. 2020) is an ideal tool. It uses multi-epoch composite spectra and decomposes them into the individual components, red- and blue-shifted according to their orbital motions. This is particularly useful when the spectra in question are in some way unusual and/or poorly matched by available templates, such as in post-interaction systems.

For a detailed description of the pipeline employed here, see Seeburger et al. (2024). In short, spectral disentangling assumes that the spectra of the components are time-invariant in their shape (i.e., lines are always the same shape and depth, and are at the same rest wavelength), and that the observed changes in the composite spectrum arise solely from RV shifts of the components. Then, the problem can be framed as asking “which two spectra can be shifted relative to each other such that they reproduce the observations at all epochs?”. Hence, for a set of trial primary velocities and mass ratios, the secondary velocities can be obtained and the two spectra computed.

Importantly, we can assess the quality of a disentangling result by comparing the observed data to a composite reconstruction from the disentangled components at the relevant shifts for each epoch and computing their $\chi^2$. A small $\chi^2$ suggests good agreement (within the errors) between observation and reconstruction, and thus a successful disentangling.

To successfully disentangle the spectra, a number of input parameters are needed. These are, at their most basic, the velocities of the primary and secondary component at each epoch. If they are in orbit, then their velocities are related via the mass ratio $q = \frac{M_{\mathrm{B}}}{M_{\mathrm{A}}}$. Thus, the set of input parameters reduces to the primary velocities and the mass ratio.

An initial guess of the primary velocities can be acquired by cross-correlation with a suitable template, or from an existing orbital solution. Here, it does not matter significantly if the template-data match is imperfect, as velocities will be refined in a later step. We use the orbital solution from Saracino et al. (2023) to compute the primary RVs at the dates our spectra were taken.

Holding these fixed, we can now determine an initial guess of the mass ratio. To do so, we trial different values of $q$ and perform spectral disentangling, recording the $\chi^2$ achieved by the reconstructed spectra for each $q$. Our starting guess is then the $q$ with the lowest associated $\chi^2$.

We show this in Figure 2, for a total of 300 orbital realisations. These realisations were generated by sampling from the orbital fit by Saracino et al. (2023), assuming Gaussian-distributed uncertainties in period, eccentricity, RV semi-amplitude and systemic velocity. Thus, for each realisation we can compute a set of predicted primary RVs. We hold these fixed, and attempt disentangling at a range of values of $q$, recording each solution’s $\chi^2$. The top panels show the behaviour of $\chi^2$ as a function of $\log q$, with each line representing one realisation. Note that we normalise the $\chi^2$ curves to a common value of 1 at $\log q = 1.1$ to account for the fact that not all orbit realisations fit the data equally well. We split the range into $0 < \log q < 0.7$ and $0.7 < \log q < 1.1$, as initial explorations revealed two minima in the $\chi^2$ curves, one within each of the stated ranges. The implications of this $\chi^2$ behaviour are twofold. First, presence of defined minima suggests that the analysis recovers two luminous companions in the data. Second, the presence of several minima indicates that the mass ratio of the two components is not well constrained. To account for the latter, we select the $\log q$ with the smallest associated $\chi^2$ within

either sub-range included in Figure 2. The bottom panels show, for each sub-range, a 2D histogram of these best-fit $\log q$ and their associated $\chi^2$. The titles at the top give the mean and standard deviation of the samples, weighted according to their associated $\chi^2$ value. Subsequently, we demonstrate the analysis using $q = 3.3$, but perform the whole pipeline using both values, with a comparative discussion of both values in § 5.

After determining initial guesses for RVs and $q$, we next perform $\chi^2$-minimisation, allowing all primary velocities and the mass ratio to vary. This should return the set of parameters that allow the disentangler to most accurately decompose the spectra into their components. This step leads only to small changes in the velocities and mass ratio, of the order of a few percent.

The initial results of this process are shown in Figure 3. The top panel shows the composite for one epoch (observation in black, reconstruction in grey), the middle panels the disentangled components, and the bottom panel the residual. We note the most obvious result, namely the clear presence of two luminous components in the two middle panels. We see the Balmer series present in both stars, identifying them as likely of BAF spectral type. Additionally, we observe MgII at 4481 Å in both the magenta and cyan spectra, as well as Helium-I lines at 4009, 4026, 4121, 4144, 4388 and 4471 Å in the cyan spectrum. Collectively, these place the magenta spectrum as a likely late B- or early A-type star, and the cyan as a potential mid B-type star. Further, the absence of obvious lines or other strange features in the residual indicates that disentangling has likely been successful, and the residual consists primarily of noise.

We can also inspect the component spectra and associated reconstruction at different epochs, shown in Figure 4. Here, each panel shows a different epoch, with the observed composite spectrum in black and the disentangled components in magenta (primary) and cyan (secondary). Co-adding them yields the reconstructed spectrum in grey (dotted line). Generally good agreement between the black observation and grey reconstruction increases our confidence in the disentangling result, especially considering the substantial velocity shifts of the components, visible by eye, between the two epochs.

An additional result of the disentangling is a dynamical estimate of the mass ratio $q$, as the secondary velocities at each epoch are calculated from the primary velocities and the value of $q$ which minimises $\chi^2$ determined during the disentangling process. Here, we find $q = 3.3$, i.e. a secondary that is 3.3 times heavier than the primary. Usually, in binary stars, we use “primary” to refer to the more massive, more luminous component of the system - here, the luminosity ratio between the two components is wavelength dependent. We thus elect to refer to the initially more massive star as the primary.

Template matching

Once the components of the binary have been identified, it is much simpler to match them to available template spectra to derive stellar parameters such as effective temperatures, surface gravities and projected rotation rates. At the same time, we can also determine the flux ratio $f$ between the templates in the relevant wavelength range. This works because the disentangled component spectra have not (yet) been scaled according to a known flux ratio. We thus perform the template fitting by selecting a section of the disentangled component spectra, and computing the total $\chi^2$ between sets of two templates, scaled by a trial flux ratio. For this step, it is important to note that we expect disentangled spectra to be in the rest frame, allowing consistent comparison with templates.

Here, we use synthetic spectra from Allende Prieto et al. (2018). The spectra are split into 5 libraries; relevant for this work in particular are libraries 3 and 4. For library 3, $T_{\rm eff}$ ranges from 7000 K to 12000 K in steps of 500 K, while $\log g$ ranges from 2 to 5 in steps of 0.5. For library 4, $T_{\rm eff}$ ranges from 10000 K to 20000 K in steps of 1000 K, and $\log g$ ranges from 3 to 5 in steps of 0.5. Due to the known location of the system in NGC 1850, and to reduce the parameter space, we fix [Fe/H] at -0.5, [$\alpha$/Fe] at 0 and $\log\xi$, the logarithm of the microturbulence, at 0.3. We produce our own set of rotation-broadened spectra using RotBroadInt (Carvalho and Johns-Krull 2023), with $v \sin(i)$ ranging from 0 km/s to 300 km/s in steps of 10 km/s for each spectrum in the aforementioned libraries.

The synthetic spectra are provided at a resolution of 10000, which we artificially degrade to FLAMES’ resolution of $\sim8000$ by convolving with a Gaussian. Additionally, we resample the spectra onto the FLAMES wavelength grid, and perform normalisation the same way as on the observed spectra to ensure coherence.

To ascertain uncertainties and increase the robustness of our results, we perform the fitting using Monte Carlo. We do so by creating multiple noise realisations of the templates and repeating the fitting procedure for each realisation. From the set of best fit templates for each realisation, we get the stellar parameters of each component, as well as the flux ratio between them. We select the wavelength ranges and spectral lines that are most sensitive to the parameters we are trying to fit and refine the fit in turn as explained in the following paragraphs. For clarity, we will show the procedure for one noise realisation as an example.

As the ratios between metal and H/He lines are strong diagnostics of stellar effective temperature, we use almost the entire disentangled spectra to ascertain this. Edge effects can cause spurious solutions here, so we trim off 50 Å on the left and right ends before fitting. Figure 5 shows the disentangled spectra of both components, as well as their closest template match (magenta and cyan, respectively). Visual inspection suggests reasonably good fits in both cases; most lines present in the disentangled spectra are accounted for by the templates, providing constraints on the temperatures². The primary (magenta) shows a star with $T_{\rm eff} \simeq11$ kK, and the secondary (cyan) with

$T_{\rm eff} \simeq19$ kK, in agreement with our previous by-eye assessment of the components as a late- and mid-B type stars respectively.

To reliably determine projected rotational velocities, we want to select lines that are minimally affected by shape changes caused by varying surface gravities, disqualifying the Balmer series. Instead, we use Mg II at 4481 Å, present in both stars, to determine $v \sin(i)$. For the secondary, there is a nearby HeI line at 4471 Å; we consider this in our fit as well.

Next, we hold the effective temperatures of both components fixed at the previously determined values, allowing only $v \sin(i)$ to vary and selecting the template with the smallest associated $\chi^2$. Our results for this realisation are shown in Figure 6. Here, we find $v \sin(i) \simeq60$ km/s for the primary, and $v \sin(i) \simeq280$ km/s for the secondary.

Finally, we wish to determine $\log g$ for both components. For this, the Balmer series is suitable; we perform simultaneous $\chi^2$-minimisation for all three Balmer lines present in our disentangled spectra (H$\gamma$, H$\delta$ and H$\epsilon$). The surface gravity primarily affects the shape of the wings of the lines; we determined $v \sin(i)$ in a previous step and now hold it constant to avoid degeneracies.

In Figure 7 we show the result of our $\log g$ fit. For the primary, we determine $\log g_{\rm A} \simeq3$, and for the secondary, $\log g_{\rm B} \simeq4$.

As a by-product of the template matching process, we also get a likely flux ratio between the components from the parameter the templates were scaled by. For the lines considered in this work, this was consistently within the errors of a 1:1 flux ratio, i.e. both stars contribute roughly the same amount of flux in the FLAMES wavelength range.

By repeating this process multiple times with different noise realisations, we can take the mean and standard deviation of each assessed parameter to determine likely stellar parameters as well as an error range. These are summarised in the upper part of Table 2.

Method	Parameter	Primary	Secondary
Disentangling	$T_{\rm eff}$ [K]	$10650 \pm250$	$18660 \pm650$
	$\log g$	$3.0 \pm0.1$	$3.8 \pm0.2$
	$v \sin(i)$ [km/s]	$68 \pm14$	$275 \pm16$
	$q$		$3.3 \pm0.7$
SED fit	Radius [$R_{\odot}$]	$6.7 \pm0.16$	$4.1 \pm0.15$
	$T_{\rm eff}$ [K]	$10830 \pm210$	$18470 \pm490$
Combined	Mass [$M_{\odot}$]	$1.7 \pm0.4$	$5.2 \pm2.6$
	$q$		$3.3 \pm1.8$
	$R_L$ [$R_{\odot}$]		$6.8 \pm0.5$

Spectral Energy Distribution fitting

In addition to spectroscopic data, broad-band photometry can help us further constrain our system parameters and confirm our results. We used the Hubble photometry presented in § 2.2 to fit the components’ spectral energy distributions (SEDs). Each star’s SED is determined chiefly by its temperature and radius, and attenuated by dust extinction. Assuming we know the distance and E(B-V), we can create a model composite SED consisting of two individual SEDs where we fit effective temperatures and radii of each component until the observations are well matched by our model. We use pystellibs (Fouesneau 2016) to generate stellar SEDs, the extinction package to apply dust attenuation, and pyphot (Fouesneau 2024) to generate mock filter photometry, which we compare to the HST observations.

To perform the fit, we use the Markov chain Monte Carlo sampling package emcee (Foreman-Mackey et al. 2013). We set Gaussian priors on the effective temperatures, centred on the values obtained during the template matching of the disentangled spectra described in Sec. 3.3. We hold the parallax fixed at its literature value, and place a flat prior on the radii of both components.

By considering our spectroscopic solutions again, we can place further constraints on the SED fit. Specifically, we can compare the spectroscopic flux ratio from disentangling with the flux ratio between the two model SEDs in the wavelength range used for disentangling. This helps break degeneracies in the fitting process and ensures a spectroscopically and photometrically consistent solution. We do this by adding a term to our likelihood function which emcee seeks to maximise:

Here, $m_{\rm model}$ is the mock photometry, $m_{\rm HST}$ the HST observations, $m_{\rm error}$ their associated errors. $f_{\rm don,model}$ is the ratio of $f_{\rm donor}$ over $f_{\rm total}$ in the simulated model fluxes over the FLAMES wavelength range, and $f_{\rm don,data}$ the flux ratio determined from disentangled spectra. The weight $w$ is chosen to ensure the flux-ratio term does not dominate the likelihood function.

Figure 8 shows the cornerplot of the parameters (effective temperatures and radii) fit using the SED method. Expectedly, we see a slight anti-correlation between a component’s temperature and radius, as both a larger and hotter star lead to a stronger flux, albeit with a different distribution. Overall, the posterior looks well-sampled, strengthening our confidence in our results.

We show the best-fit SED in Figure 9, noting good agreement between model (black squares) and observations (green dots). Comparing the model flux ratio to the spectroscopically obtained result also suggests agreement within the errors. The SED fit suggests a primary radius of $R_A = 6.7$ $R_{\odot}$, and a secondary radius of $R_B = 4.1$ $R_{\odot}$. SED temperatures are consistent within the errors with spectroscopic temperatures - this is not particularly surprising as the spectroscopic temperatures acted as priors on the fit. However, a large discrepancy here would potentially be indicative of an error in the analysis or some other issue, so the congruence is reassuring.

These results can be combined with the previously determined surface gravities of both components to infer stellar masses, providing an additional independent measure of the mass ratio $q$. Here, a primary mass of $M_A = 1.7$ $M_{\odot}$ and secondary mass of $M_B = 5.2$ $M_{\odot}$ are implied, suggesting a mass ratio of $q = 3.3$. This is solely computed using effective temperatures and surface gravities from analysing the disentangled spectra, and stellar radii from the SED fit, omitting any dynamical analysis.

With the two component masses and the period as measured by Saracino et al. (2023), we can compute the primary’s Roche lobe using the approximation from Eggleton (1983):

Here, $A$ is the orbital separation, $q$ the mass ratio, and $R_{\mathrm{L}}$ the radius of a sphere with the same volume as the Roche lobe. Following this, we get $R_{\mathrm{L}} \sim6.8 \pm0.5 R_{\odot}$, implying a primary which is close to filling or fills its Roche lobe. Saracino et al. (2022) and El-Badry and Burdge (2022) both note that the system is a double-periodic variable (DPV; e.g. Mennickent et al. (2003)). OGLE photometry shows both a $P = 5$ d and $P = 156$ d periodicity, with the shorter variability likely due to orbital motion, and the longer related to mass transfer, though its precise origin is not clear. DPVs frequently do not show strong emission lines, despite ongoing mass transfer (Mennickent et al. (2005)).

The determination of these stellar parameters allows us to shed a light onto the evolutionary history of the system, which we discuss in § 5 below.

Comparison with previous results

This work marks an additional analysis of the target NGC 1850-BH1, and the first such attempt including high-resolution, blue, multi-epoch spectroscopy. We compare and contrast our findings with those of previous studies in Table 3. The properties of the binary components for which this work provides tighter constraints (e.g. masses and flux contributions) are discussed in detail below.

	Saracino et al. (2022)	El-Badry and Burdge (2022)	Saracino et al. (2023)	This work
$M_{\mathrm{A}}$	$\sim5\,M_{\odot}$	$\lesssim1.5\,M_{\odot}$		$1.7 \pm0.4\,M_{\odot}$
$M_{\mathrm{A}}$ method	CMD/MSTO	RL geometry, mean density		$\log g_{\mathrm{A}}$ & $R_{\mathrm{A}}$
A State	Normal MSTO star	Stripped, RL-filling	Stripped	Stripped, near-RL-filling
$M_{\mathrm{B}}$	$\sim11\,M_{\odot}$	$2.5$–$5\,M_{\odot}$	$\sim5\,M_{\odot}$	$\sim5\,M_{\odot}$
$M_{\mathrm{B}}$ method	$K_1$, $P$, $i$			$\log g_{\mathrm{B}}$ & $R_{\mathrm{B}}$
B nature	BH	normal MS star	unclear	accretor
$f_{\mathrm{B}}$	dark	hidden by primary	$\lesssim10\%$ of MUSE flux	visible, $\sim50\%$ of FLAMES flux
$i$	$\sim38^\circ$	Higher, diluted ellipsoidal signal	$\lesssim67^\circ$ (no eclipses)	N/A

Masses

Saracino et al. (2022) derived a primary mass of $\sim5 M_{\odot}$ from the targets location in the CMD. However, this mass estimate did not take into account the possibility of binary interaction, which this work has revealed most likely took place. Mass transfer fundamentally alters the properties of stars, and makes it impossible to use single-star based isochrones to accurately determine stellar parameters. An alternative scenario considering a (partially) stripped donor was originally introduced by El-Badry and Burdge (2022), lowering the primary mass estimate to $\lesssim1.5 M_{\odot}$ (El-Badry and Burdge (2022)). In this work, we find $M_{\mathrm{A}} \sim1.7 M_{\odot}$ consistent with the upper limit of El-Badry and Burdge (2022), with the main uncertainty related to the donor $\log g$.

The higher primary mass from Saracino et al. (2022), combined with the system’s large RV variability and short period lead to a large implied secondary mass of $\sim11 M_{\odot}$. As the El-Badry and Burdge (2022) interpretation significantly decreases the assumed primary mass, this, in turn, lowers the secondary mass to $2.5 - 5 M_{\odot}$. This picture is consistent with this work, where we find a secondary mass of $\sim5.2 \pm2.6 M_{\odot}$ with the uncertainty being dominated by the secondary’s $\log g$.

Fluxes

In Saracino et al. (2022), the authors conclude that the secondary is dark, and thus contributes no flux to the system, as no MS star could reasonably fulfill the mass constraints on the secondary while remaining hidden. With the updated masses from El-Badry and Burdge (2022), a luminous secondary “hiding” in the light of the puffed-up primary becomes viable. The updated velocities and thus binary mass function of Saracino et al. (2023) imply a secondary that, if it is luminous, contributes at most 10% of the light in the MUSE wavelength range. They infer that such a secondary must be hidden in an optically thick disc to evade detection in the MUSE disentangling. This work identifies the signature of the luminous secondary, and places its flux contribution at about 50 % in the FLAMES wavelength range. SED fits paired with spectroscopic temperatures imply a secondary flux contribution between 35-50 % in the MUSE data, which is in conflict with the 10 % limit from Saracino et al. (2023). Possible explanations for this discrepancy include:

• An optically thick disc, as potentially implied by IR observations and suggested by Saracino et al. (2023), which hides (some of) the flux of the secondary in the redder regime. This may also account for the 156 d secondary period. However, the SED fit computed in this work does not necessitate a further component to explain the observed photometry. Additionally, circumstellar material is not expected to selectively absorb at longer optical wavelengths.

• As mentioned, NGC 1850-BH1 has a 1.7 mag fainter neighbour, that is too faint and too close to NGC 1850-BH1 to be resolved in the FLAMES data, and thus might contaminate the extracted spectra. While this is possible, and may contribute to the preference for the disentangler for a large mass ratio (assuming this star is single, and thus its RV does change, at least not on the timescale of the NGC 1850-BH1 orbit), we argue that our recovered secondary spectrum is due to NGC 1850-BH1’s companion, rather than this source. First, the neighbour is about 5 times dimmer, while we find our recovered secondary makes up about half the received flux in the FLAMES wavelength range. Second, the neighbour is unlikely to show He-lines, as it is cooler than our recovered secondary. It may, potentially, more strongly affect the MUSE spectra, precisely due to it being cooler and thus more

visible in the red. However, as the star could be deblended and a spectrum extracted in the MUSE data, this should not be the case.

• Similarly, an unresolved tertiary to the system might contaminate the spectra, as well as induce the additional periodicity.

• The accretor as recovered in this work shows only very few, and very broadened lines. In the MUSE wavelength range, we would mainly expect to see H and He lines, where the H lines are also present in the donor star, and thus hard to disentangle. Applying the code used in this work to the MUSE spectra, we find evidence of a secondary with He-absorption at 5016 Å, but do not manage to recover any other He-lines. Disentangling in the Paschen-line dominated wavelength regime did not reveal a clear signature of a companion, in agreement with the analysis performed by Saracino et al. (2023) in the same regime. One may speculate that the properties of the MUSE spectra (low spectral resolution, non-Gaussian noise properties, only broad Paschen lines) pose a more substantial challenge to spectral disentangling codes than initial tests suggested.

Inclination

Saracino et al. (2022)’s initial work constrained the inclination of the to around $38^\circ$, based on lightcurve data. Subsequent work by El-Badry and Burdge (2022) called this into question, specifically the assumption of the visible star being a normal $5\,M_{\odot}$ MSTO star which fills its Roche lobe, and the system’s photometric variability being undiluted by a secondary. Allowing for larger inclination removes the need for a very high mass secondary. The updated analysis in Saracino et al. (2023) contributes two important points to this discussion. First, the lack of visible eclipses constrains the inclination to $i \lesssim67^\circ$. And second, the newly measured velocities imply a larger velocity semi-amplitude of the primary, and thus a larger binary mass function of $f_{\mathrm{bin}} = 2.83$. Here, we explore the consequences of these updates.

Figure 10 follows the idea presented in El-Badry and Burdge (2022), showing the implied secondary mass for a range of primary masses as a function of inclination, using the binary mass function. Here, we use the larger binary mass function from Saracino et al. (2023), and show the mean $\pm1\sigma$ for the primary mass found in this work (coloured lines). The grey shaded region is the $1\sigma$ range of secondary masses as found here from spectral and SED fits. Considering the constraint from Saracino et al. (2023), that $i \lesssim67^\circ$, we note that our combination of masses falls essentially just inside the “allowed” region. Tighter constraints on the $\log g$ for both components using e.g. higher resolution spectroscopy would be helpful in further exploring the here proposed solution and its agreement with constraints found in other work.

Discussion

Evolutionary history— As spectral disentangling was able to determine two bright components in NGC-1850 BH1, we can confidently rule out a dark secondary in the system. Nevertheless, the current state and evolutionary history of the binary can help shed light on some of the most poorly understood processes in binary evolution, namely binary mass transfer, angular momentum transport, and stellar rotation.

Qualitatively, the primary likely started off as the initially more massive star, expanding as it reached the MSTO. As it did so, it filled its Roche Lobe, and material overflowed onto the secondary. The secondary thus gained both mass and angular momentum, increasing its rotational velocity. After reversal of the mass ratio (the secondary becoming more massive than the primary), further stable mass transfer led to a widening of the orbit, and eventually the primary no longer filled its Roche Lobe, and mass transfer stopped. We infer the cessation of mass transfer from the lack of emission visible in Hα as seen in the MUSE spectra in Saracino et al. (2022). Ongoing or recently ended mass transfer would be expected to be accompanied by an accretion disk around the secondary, with emission wings around H-α, such as in some objects in Seeburger et al. (2026).

Previously, it was assumed that the donor star in such a mass-transfer scenario would contract quickly, evolving into a subdwarf. Recently, however, multiple findings of “bloated” or “puffed-up” stripped stars (e.g. Shenar et al. 2020; Bodensteiner et al. 2020; El-Badry et al. 2022) have caused this to be reconsidered. This intermediary stage, before the stripped star settles into a hot, compact subdwarf, is now expected to last much longer, about 10% of the He-MS lifetime (Dutta and Klencki 2024).

NGC-1850 BH1 thus belongs to a class of binary systems which masquerade as containing dormant black holes (e.g. Rivinius et al. 2020; Jayasinghe et al. 2022, 2021; Liu et al. 2019), but were all subsequently found to be post-interaction binaries with two luminous components (e.g. Shenar et al. 2020; Bodensteiner et al. 2020; El-Badry et al. 2022). In many of these cases, spectral disentangling has provided the definitive proof for a bright companion, as in this work.

The mass ratio— Interestingly, initial disentangling work showed two plausible minima in the residual for different values of $q$: $\sim3.3$ and $\sim8.6$. We rule out the larger value by appealing to stellar evolution: The system’s position at the MSTO constrains the initial mass of the primary to around 5 $M_\odot$, with an initially less massive companion. For the system to reach the mass ratio of 8.6, assuming initial masses of 5 and 4.5 $M_\odot$, the primary must have lost over 4 $M_\odot$ of material which the secondary accreted. This is questionable, both as we expect the core of a 5 $M_\odot$ star to be approximately 1 $M_\odot$, and as it would require mass transfer to be nearly fully conservative. Additionally, even if one were to assume $q = 8.6$ and proceed through the analysis with this value, the resulting values of surface abundances and radii inferred from template and SED fits, and the thus inferred masses, are inconsistent with the assumed value of $q$. This is shown in Appendix A.

Finally, appealing to simple equations of angular momentum transfer from Soberman et al. (1997), we can compute the evolution of the period as a function of mass ratio as we go back in time, as in Bodensteiner et al. (2020). This approach, which is based on the idea that the orbit evolution of a binary depends on the properties of the mass transfer and hence the orbit evolution of the binary can be reconstructed for a given mass-transfer prescription, is summarized in Figure 11. Black and grey stars show the present-day period and mass ratio for the two minima found by the disentangler. Corresponding solid and dashed lines illustrate the previous evolution of the period under the assumption of different isotropic re-emission parameters $\beta$. The minimum period is shown with horizontal shading, depending on the initial mass of the donor. Notably, for the larger mass ratio case, all lines, including the one indicating fully conservative mass transfer, dip below this minimum. While it is important to highlight that the Soberman et al. (1997) equations make a number of simplifying assumptions that mean this does not necessarily rule out the $q = 8.6$ case, it is another supporting piece for the $q = 3.3$ scenario. In any case, the fact that even in the $q = 3.3$ scenario, the lines dip below the minimum allowed period, suggests that future studies of NGC1850-BH1 can be used to stress-test assumptions of binary mass transfer.

Remaining issues and uncertainties — At this point, it is important to highlight some of the remaining problems and uncertainties with the analysis.

• The uncertainties associated with the log $g$ values of both components (cf. Table 2) are important limitations of our analysis, as the component masses are very sensitive to them. Higher resolution/SNR observations and a more tightly constrained orbital solution might help reduce the errors and better constrain the other parameters in the cascade (masses and mass ratios).

• Our templates assume solar-scaled abundances, with a “normal” surface hydrogen fraction. However, due to the previous mass transfer, we expect the chemical composition of both components to be altered. Initial exploration (Dutta et al., in prep, using CMFGEN (Hillier and Miller (1998))) suggests that a changing surface H-fraction has only small effects on the resultant spectra down to very strong H-stripping, but this may regardless have a knock-on effect on measured $T_{\mathrm{eff}}$, log $g$, $v\sin(i)$ and flux ratios.

Conclusion

Spectral disentangling has demonstrated that the system known as NGC 1850 BH1 does in fact not contain a compact companion, but falls into the class of “imposter” systems with two bright components. The system presently consists of a luminous, comparatively cool ($T_{\mathrm{eff,A}} \simeq10650$ K), puffy (log $g_{\mathrm{A}} \simeq3.0$) primary, and a more massive, hotter ($T_{\mathrm{eff,B}} \simeq18660$ K) and denser (log $g_{\mathrm{B}} \simeq3.8$) secondary. This atypical configuration is most readily explained by binary interaction.

This work provides additional analysis, building on that of Saracino et al. (2022), El-Badry and Burdge (2022) and Saracino et al. (2023). While initial work classified the primary as a 5 M$_\odot$ MSTO star and inferred a mass of 11 M$_\odot$ for the (then unseen) companion, subsequent analysis (including this work) has lowered the mass estimate of the primary, revealed the secondary, and placed constraints on its nature, as well as the system’s evolutionary history. There are, however, some open questions.

In light of the poorly constrained mass ratio, a full characterisation of the system is still not possible. However, the location of NGC 1850-BH1 inside a cluster of well-known age, metallicity, and distance suggests the system could reveal important insights into our understanding of binary interactions. Therefore, pinning down the mass ratio through observations at higher spectral resolution and/or SNR appears like an important next step. Given the location of the system in a crowded stellar field and the relative faintness of the two components, such observations can currently only be performed from space (e.g., HST-STIS) or with blue-efficient long-slit (e.g., X-Shooter) or integral-field spectrographs at 8m-class ground-based telescopes.

These data might also aid in explaining the discrepancy between the 10% flux limit on the secondary in the MUSE wavelength range found by Saracino et al. (2023), and the here implied 30–50% flux contribution in the same range. Whether this inconsistency is due to the limited MUSE spectral resolution and similarity of the two spectra in this wavelength regime, dilution from an unrelated star or a tertiary companion, or contamination from a circumbinary disc can only be determined with high-resolution multi-epoch data in the red and/or NIR. Such data may also be useful in exploring the origin of the 156 d secondary period, and provide additional constraints on the parameters of the system and its components.

Acknowledgments

Based on observations collected at the European Southern Observatory under ESO programme 113.26S5.001. RS and SK gratefully acknowledge funding from UKRI through a Future Leaders Fellowship (grant no. MR/Y034147/1). SS acknowledges funding from the European Union under the grant ERC-2022-AdG, ‘StarDance: the non-canonical evolution of stars in clusters’, Grant Agreement 101093572, PI: E. Pancino. TS acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement 101164755/METAL) and the Israel Science Foundation (ISF) under grant number 0603225041. JB is supported by an NWO Veni fellowship (VI.Veni.242.199).

We thank Ylva Göttberg, Hugues Sana, Matti Dorsch and Ulrich Heber for helpful discussion of the target, as well as support with the original FLAMES proposal.

The q = 8.6 case

In Figure 12, we show the different results from disentangling when using the same primary and centre of mass velocities, but different $q$. Qualitatively, the disentangled spectra using $q = 8.6$ are similar to those achieved using $q = 3.3$ (top two panels, darker, dotted and lighter, solid lines respectively). Differences become more apparent when looking at the residuals, shown in the bottom panels. Particularly around the Balmer series, the line shape changes slightly, leading to deeper lines in the primary for the $q = 3.3$ case over the $q = 8.6$ case. The secondary lines are affected in the opposite way, with the Balmer series deeper for larger mass ratios. This likely happens as a larger mass ratio leads to less computed secondary motion, and thus more of the line depth gets ascribed to the secondary during disentangling. These differences are subdominant to noise when fitting templates, leading to similar (within the errors) stellar parameters, shown in Table 4. We note that the parameters generally fall within the errors of those acquired in the $q = 3.3$ analysis. The biggest change occurs in the masses and thus the combined-method (spectral $\log g$ and SED radii) mass ratio. Even here, however, due to the large errors owing to the comparatively large uncertainty on $\log g$, the results between the two analyses are consistent with each other. However, we note that the value of $q = 3.5 \pm2.4$ derived from a combination of SED-based radii and spectroscopically determined $\log g$ is inconsistent with the value of 8.6 obtained from disentangling.

Method	Parameter	Primary	Secondary
Disentangling	$T_{\mathrm{eff}}$ [K]	10500 $\pm$ 280	18860 $\pm$ 400
	$\log g$	3.0 $\pm$ 0.2	3.9 $\pm$ 0.2
	$v\sin(i)$ [km/s]	56 $\pm$ 12	278 $\pm$ 16
	$q$		8.6 $\pm$ 1.0
SED fit	Radius [$\mathrm{R}_{\odot}$]	6.8 $\pm$ 0.18	4.1 $\pm$ 0.11
	$T_{\mathrm{eff}}$ [K]	10720 $\pm$ 200	18550 $\pm$ 340
Combined	Mass [$\mathrm{M}_{\odot}$]	1.9 $\pm$ 1.0	5.4 $\pm$ 2.6
	$q$		3.5 $\pm$ 2.4
	$R_{\mathrm{L}}$ [$\mathrm{R}_{\odot}$]		6.8 $\pm$ 0.5