High mass function ellipsoidal variables in the Gaia Focused Product Release: searching for black hole candidates in the binary zoo

The recent Gaia Focused Product Release contains radial velocity time-series for more than 9,000 Gaia long-period photometric variables. Here we search for binary systems with large radial velocity amplitudes to identify candidates with massive, unseen companions. Eight targets have binary mass function $f(M)>1\ M_\odot$, three of which are eclipsing binaries. The remaining five show evidence of ellipsoidal modulations. We fit spectroscopic orbit models to the Gaia radial velocities, and fit the spectral energy distributions of three targets. For the two systems most likely to host dark companions, J0946 and J1640, we use PHOEBE to fit the ASAS-SN light curves and Gaia radial velocities. The derived companion masses are $>3 M_\odot$, but the high Galactic dust extinctions towards these objects limit our ability to rule out main sequence companions or subgiants hotter than the photometric primaries. These systems are similar to other stellar-mass black hole impostors, notably the Unicorn (V723 Mon) and the Giraffe (2M04123153$+$6738486). While it is possible that J1640 and J0946 are similar examples of stripped giant star binaries, high-resolution spectra can be used to determine the nature of their companions.


INTRODUCTION
The mass distribution of stellar mass black holes (BHs) is directly related to the late-stage evolution of massive stars and core-collapse supernovae.By observing and characterizing the black hole population in the Milky Way, we stand to learn more about the initial-final mass relation (IFMR), which describes the connection between the pre-supernova mass of the star and the type of compact object produced.Various physical dependencies complicate the modeling of the IFMR, such as metallicity, mass-loss rate, and binary interactions (Sukhbold et al. 2016).An accurate census of compact objects is also needed to understand the "lower mass gap" between the most massive neutron stars (∼2.1  ⊙ , Antoniadis et al. 2016;Cromartie et al. 2020) and the least massive BHs in X-ray binaries (∼5  ⊙ , Özel et al. 2010;Farr et al. 2011;Kochanek 2015).
The Milky Way is expected to contain ∼10 8 stellar-mass BHs and ∼10 9 neutron stars (Timmes et al. 1996;Wiktorowicz et al. 2019).The compact objects detected so far are dominated by X-ray binaries (Neumann et al. 2023) and extra-Galactic gravitational wave mergers (e.g., Abbott et al. 2016).However, only a small number of Galactic BHs are expected to be in the tight binary orbits necessary to produce X-ray emission (Corral-Santana et al. 2016), and the fraction of systems that merge in gravitational wave events is also a strongly biased fraction of the overall population (e.g., Kruckow et al. 2018).The majority of Galactic BHs are instead most likely in non-interacting binaries that are not X-ray bright or are isolated free-floating systems.Discovering and characterizing these systems is crucial to understanding the population and improving models for the end states of stellar evolution.While isolated BHs can only be detected via microlensing surveys (Lam et al. 2020(Lam et al. , 2022;;Sahu et al. 2022), there are multiple methods that can be used to search for and characterize non-interacting BHs in binary systems.
Gaia astrometry has proven to be an effective tool to probe binary orbits.Population synthesis studies predict that Gaia astrometry will detect tens to thousands of BHs (Breivik et al. 2017;Yamaguchi et al. 2018), though these predictions require assumptions about complicated physical processes like common envelope evolution and neutron star/BH natal kicks.So far, only two strong BH candidates have been identified using Gaia astrometry (Gaia BH-1: El-Badry et al. 2023a; Chakrabarti et al. 2023;Gaia BH-2: Tanikawa et al. 2023;El-Badry et al. 2023b).The formation and evolutionary history of these systems are unclear, and there is some evidence that they formed from dynamical interactions in cluster environments (e.g., Rastello et al. 2023).
Spectroscopic surveys can be used to search for single-lined binaries (SB1s) with large amplitude periodic radial velocity (RV) variations.Systems identified through RVs can only yield a lower limit to the companion mass because, unlike astrometric binaries, the orbital inclination is not directly constrained.For example, Thompson et al. (2019) reported the detection of a possible stellar-mass black hole orbiting the spotted red giant 2MASS J05215658+4359220 using RVs from the Apache Point Observatory Galactic Evolution Experiment (APOGEE, Majewski et al. 2017).Some candidates have also been identified using spectroscopic data from the Large Sky Area Multi-object Fiber Spectroscopic Telescope (LAMOST, Cui et al. 2012;Gu et al. 2019), but additional observations are needed to confirm their orbits and to rule out luminous companions.More recently, Gaia Data Release 3 included spectroscopic orbits for > 180, 000 SB1s (Gaia Collaboration et al. 2023a,c), however no strong candidates have been identified (El-Badry & Rix 2022;Jayasinghe et al. 2023).Aside from large spectroscopic surveys, dedicated observations of targets in globular clusters have identified several candidates (Giesers et al. 2018(Giesers et al. , 2019)), but dynamical effects in clusters may drive the formation and evolution of these systems, making comparison to the field population challenging (e.g., Ryu et al. 2023).
For some spectroscopic binaries, photometric ellipsoidal variability can be used to place additional constraints on the mass ratio and inclination of the binary.Ellipsoidal modulations are caused by the tidal distortion of a star by its binary companion.Some X-ray binary systems are ellipsoidal variables (ELLs; e.g., Orosz et al. 2001), but the accretion and its variability can limit our ability to use ellipsoidal modulations to characterize these systems.There have been a number of searches to identify ELLs with high mass companions in photometric surveys (e.g., Rowan et al. 2021;Gomel et al. 2021Gomel et al. , 2023;;Green et al. 2023), but spectroscopic follow-up is required for these systems and no strong candidates have been identified (e.g., Nagarajan et al. 2023).
The presence of ellipsoidal modulations can also indicate a history of mass transfer that complicates the characterization of systems.For example, V723 Mon was initially reported as a red giant with a mass-gap black hole companion  ∼ 3  ⊙ (Jayasinghe et al. 2021b).The light curve shows clear ellipsoidal variability consistent with the RV orbital period and no strong evidence of a luminous companion in the spectra.Using spectral disentangling, El-Badry et al. (2022) showed that the system is instead better described by a very low mass stripped red giant  ≈ 0.4  ⊙ with a subgiant companion at a similar effective temperature.The subgiant is also rapidly rotating, blurring its spectral features.2M04123153+6738486 was similarly identified as a BH false-positive better explained by a stripped star scenario (Jayasinghe et al. 2022), demonstrating the complexities in ruling out luminous companions for ellipsoidal variables with large amplitude RV variations.
As astrometric, spectroscopic, and photometric surveys expand, it is becoming increasingly important to be able to promptly and accurately vet targets and identify the strongest non-interacting BH candidates.Gaia DR3 includes more than 800,000 non-single star solutions.The next Gaia data release, not expected before 2025, will include many more binary star orbits as well as epoch RVs and astrometric measurements.The recent Gaia Focused Product Release included time-series RV measurements for a sample of > 9, 000 stars identified as long-period variables (LPVs) in Gaia photometry (Gaia Collaboration et al. 2023d).Here, we use this catalog to identify ellipsoidal variables with high radial velocity amplitudes that could host compact object companions.
In Section §2, we describe how we select targets with high binary mass function and reject obvious eclipsing binaries.We fit spectroscopic orbit models in Section §3 and spectral energy distributions (SEDs) in Section §4.Using the SED fit results, we model the light curve and radial velocity curve using PHysics Of Eclipsing BinariEs (PHOEBE, Prša & Zwitter 2005;Conroy et al. 2020) in Section §5.Finally, we present possible interpretations of the systems and describe the additional observations that are necessary to better characterize these binaries in Section §6.

TARGET SELECTION
The details of the Gaia FPR catalog construction are described in Gaia Collaboration et al. (2023d).In brief, the process starts from the full sample of more than 1.7 million Gaia long-period variables (Lebzelter et al. 2023).Only ∼ 500, 000 of these have Gaia RV measurements.After making cuts on the magnitude ( RVS < 12 mag), the number of visibility periods (> 12), and on the ratio  RV / RV < 0.175 between the RV amplitude,  RV , and the RV uncertainty  RV , only ∼ 111, 000 targets remain.The photometric and RV time-series are then searched for periodic signals.Systems with long-term linear RV trends trends are also removed leaving only ∼ 44, 000 systems.Finally, targets are selected to have photometric periods equal to or twice the best-fit period of the RV time series.The final sample contains 9,614 long-period variables with epoch RVs.This is the second sample of epoch RVs from Gaia, following the samples of ∼ 2, 000 Cepheids (Ripepi et al. 2023) and RR Lyrae variables (Clementini et al. 2023) included in Gaia DR3.
Gaia DR3 also includes a large number of binary star orbits as part of the "non-single-stars" (NSS) tables.We crossmatch the Gaia FPR catalogs with the Gaia NSS tables and find a total of 862 matches, 854 of which are in the singlelined spectroscopic binary sample.Gaia Collaboration et al. (2023d) compares the orbital periods from the FPR analysis to the Gaia DR3 SB1 solution (their Figure 21) and find that the

✗
RV periods for the targets selected by Equation 1 agree within 10% for all but one system.Seven of the remaining matches are in the single-lined second degree trend sample.Only one target in the Gaia LPV FPR catalog, GDR3 2567779977831471232, has an astrometric orbit solution.Interestingly, the astrometric orbital period  = 927 ± 85 days for this system does not agree with the RV period in the Gaia FPR catalog of  = 695 days.Gaia Collaboration et al. (2023d) show how the 9,614 sources in the Gaia FPR can be separated into pulsators and binaries by comparing the amplitude of the photometric variability,   , to the amplitude of the radial velocity variability,  RV , using the criteria that systems with   < 0.3 mag and   < 0.25 (1) are binaries.Figure 1 shows the distribution of   and  RV and the 1,109 targets selected as binaries using Equation 1.
The binary mass function (2) is the absolute lower limit on the companion mass obtained in the limit of an edge-on inclination and a zero mass ( 1 = 0  ⊙ ) photometric primary.Figure 2 shows the distribution of mass functions for the binaries using the RV period and semi-amplitude from the Gaia FPR catalog and assuming the eccentricity is  = 0.There are 8 binaries with  () > 1.0  ⊙ which are listed in Table 1.

Rejecting Eclipsing Binaries
The targets in the Gaia FPR catalog were selected based on their Gaia photometric variability.The cadence and number of epochs in the Gaia light curves are not well suited for identifying short, narrow eclipsing features in long-period binaries.To that end, we use light curves from the All-Sky Automated Survey for Supernovae (ASAS-SN, Shappee et al. 2014;Kochanek et al. 2017;Hart et al. 2023) to identify eclipsing binaries in the sample of high mass function binaries.
Out of the eight, targets, four are ASAS-SN variables (J0946, J2017, J1640, and J1030) and all were classified as semiregular variables (Jayasinghe et al. 2021a;Christy et al. 2023).The ASAS-SN photometric period is also approximately half of the reported Gaia RV period, as expected for ellipsoidal variables.Figure 3 shows the ASAS-SN light curves of six of the eight high  () targets.We perform phase-dispersion minimization to refine the Gaia period by searching a narrow period window (0.05) around the Gaia period.
Two targets, J0812 and J1513 do not have ASAS-SN Skypatrol V2 light curves.J0812 has a nearby (< 5. ′′ 0) star, HD 68845, that is bright ( ∼ 9 mag) and blended with it in the ASAS-SN images.J1513 was omitted from the ASAS-SN Skypatrol V2 sample (Hart et al. 2023) because the ATLAS-REFCAT2 (Tonry et al. 2018) -band magnitude is  = 17.9 mag and the proximity statistic1 is < 20.′′ 0. We recomputed the ASAS-SN -band light curve using aperture photometry (Kochanek et al. 2017) and found that the target has a median  = 14.4 mag.However, we find no evidence of periodic photometric variability using a Lomb-Scargle periodogram (Lomb 1976;Scargle 1982).J2017 is  ∼ 10.2 mag while ASAS-SN begins to saturate for stars brighter than  ∼ 11 mag.For these three targets we instead show the Gaia light curves in Figure 3.
We remove J0131, J1030, and J2059 from consideration as non-interacting compact object binary candidates because of the presence of eclipses in their ASAS-SN light curves.The secondary eclipses in all three light curves are shallow or not visible, suggesting that the photometric secondary has a lower surface brightness.J0812, J0946, J1640, and J2017 have no evidence of eclipses and are likely ellipsoidal variables.J0812 also has uneven maxima in its Gaia light curve, which could be evidence for the presence of a spot on the giant star.The Gaia light curve of J1513 does have evidence for periodic photometric variability, but it is not a clear ellipsoidal variable.For all four targets, we find that the light curves phase with the Gaia RVs as expected for binary stars.
We also inspect the TESS light curves of our targets.As compared to the Gaia and ASAS-SN light curves, the TESS light curves are inherently less useful for characterizing these systems because of the 27-day TESS Sector length.Many of the targets are also near the Galactic plane, so crowding presents an additional challenge in interpreting light curves given the large 21" TESS pixels.Even with these limitations, TESS is still useful to search for low-amplitude features, such as shallow eclipses.We use three pipelines to retrieve TESS light curves: 1.The Quick-Look Pipeline (QLP, Huang et al. 2020a,b).These light curves are processed from the full-frame images (FFIs) and are available for targets brighter than  = 13.5 mag.Since the detrending process can remove real variability on timescales > 0.3 days, we use the raw, undetrended light curves.
2. The Science Processing Operations Center (SPOC, Caldwell et al. 2020).These light curves are also generated from the FFIs, but are only available for 10,000 targets per Sector for each CCD with a magnitude limit of  = 13.5 mag.

TESS-Gaia Light Curve (TGLC, Han & Brandt 2023).
This pipeline models the FFI with point-spread functions (PSFs) based on Gaia DR3 positions.We use the calibrated PSF flux to reduce contamination from nearby variable targets.
Since the long-period ELL or eclipsing binary (EB) signal produces a long-term trend in the TESS light curves for each Sector, we do a simple check for prominent variable signatures by-eye to identify any shallow eclipsing features.J0812, J0946, J1513, J1640, J2017, and J2059 show no evidence of variability aside from the primary ELL/EB signal.The primary eclipse of J1030 is detected in the Sector 64 light curve with a phasing consistent with the ASAS-SN light curve.No secondary eclipse or other variable signals are observed.Finally, for J0131, we find some evidence of short period variability and additional eclipsing features outside of the  = 145 day signal.We show this system's TESS light curve and discuss why it is likely a blended target in Appendix A.

RADIAL VELOCITY ORBITS
The Gaia FPR catalog reports the period and amplitude of the photometric variability but they do not fit a spectroscopic orbit model.Here, we fit the Gaia time-series RVs using a prior based on the photometric period.We fit a Keplerian orbit model of the form where  is the center-of-mass velocity,  is the radial velocity semiamplitude,  is the true anomaly, and  is the argument of periastron.The true anomaly,  , is related to the eccentric anomaly, , and the eccentricity,  by and the eccentric anomaly is where  is the period and  0 is the time of periastron.We start by using TheJoker (Price-Whelan et al. 2017) for rejection sampling, using a Gaussian prior on the period based on the ELL period from the light curve.We use the resulting solution to initialize walkers for an MCMC fit with emcee (Foreman-Mackey et al. 2013).Figure 4 shows the RV orbits for the five ellipsoidal variable targets and Table 2 reports the orbital parameters.Two of the ELLs, J0812 and J0946, were also included in the Gaia DR3 SB1 orbit catalog (Gaia Collaboration et al. 2023a;Katz et al. 2023).For J0812, the Gaia SB1 solution is  SB1 = 681 ± 2 days,  SB1 = 0.077 ± 0.008, and  SB1 = 29.9 ± 0.2 km/s, which  4. The ratio of the orbital motion to the parallax is calculated using Equation 6assuming an edge-on inclination.agrees with our orbital solution (Table 2).For J0946, the Gaia SB1 solution is  SB1 = 102.60 ± 0.07 days,  SB1 = 0.076 ± 0.019, and  SB1 = 56.4± 1.1 km/s, which also agrees with our orbital fit.Although it is not surprising that we find the same orbital solution since we are fitting the same RV measurements, some of the Gaia spectroscopic orbit solutions have been found to be spurious.Bashi et al. (2022) use radial velocities from large surveys to predict a "score", S, for each SB1 solution on a zero to one scale.J0812 and J0946 have S = 0.92 and S = 0.71, respectively, which are both well above the cutoff S > 0.587 that Bashi et al. ( 2022) use for their "clean" sample.
The RV orbit of J1513 is less well-constrained.We searched for orbital solutions at different periods using rejection sampling with TheJoker, but were unable to find any other solution.Gaia DR3 does not include an SB1 orbit for this target because the effective temperature of the RV template is 3500 K, below the cutoff of 3875 K (Gaia Collaboration et al. 2023c).Since the light curve of this target is also the least consistent with ellipsoidal modulations, additional RV and photometric observations are needed to better characterize this system.We do not list this target as an ELL in Table 1 and remove it from consideration as a non-interacting compact object candidate.

PROPERTIES OF THE GIANT STARS
To better understand the nature of the unseen companions we first need to characterize their photometric primaries.The Gaia LPV catalog (Lebzelter et al. 2023), the starting point for the Gaia FPR catalog, includes a color filtering criteria of  BP −  RP > 0.5 mag.The color-magnitude diagram (CMD) in Figure 5 shows that all of the targets in the Gaia FPR catalog are on the giant branch as expected.Absolute magnitudes are calculated using photogeometric distances from Bailer-Jones et al. ( 2021) and -band extinctions from the mwdust (Bovy et al. 2016) "Combined19" dust map (Drimmel et al. 2003;Marshall et al. 2006;Green et al. 2019).We convert the -band extinction to extinctions in the Gaia filters using the coefficients from Wang & Chen (2019).J1513 has a large fractional parallax error, /  = 2.3, so its CMD position is uncertain and it is not shown in Figure 5.
The Gaia parallax measurements could also be biased by the orbital motion of the binaries for these systems.Following Thompson et al. (2019), we compare the orbital motion of the binary to the measured parallax as Table 2 includes this ratio for each target computed using the MCMC posteriors for the RV orbit assuming an edge-on inclination.For J0812, J1513, and J2017, the expected astrometric motion from the binary orbit is larger than the measured parallax.This could bias the stellar luminosity and radii in either direction depending on the projection of the binary on the sky.None of these targets have an astrometric orbit solution in Gaia DR3.Only J2017 has a Gaia renormalized unit weight error RUWE > 2 (Table 1).
The three eclipsing targets (J0131, J1030, and J2059) are bluer than the majority of the Gaia FPR sample, which indicates that their companions are hotter subgiants or main sequence stars.The ELL J0946 also has extinction corrected color  BP −  RP ≈ 1.0 mag, which suggests the presence of a luminous companion.
Next, we fit the spectral energy distributions (SEDs) for three of the four remaining targets.We do not fit the SED of J0812 because it is blended with a nearby bright star.We use broad-band photometry from 2MASS (Cutri et al. 2003), WISE (Wright et al. 2010;Cutri & et al. 2012), and the lowresolution Gaia XP spectra (De Angeli et al. 2023).J1640 was observed by GALEX (Bianchi et al. 2017) as part of the All-Sky Imaging Survey (AIS) in the near ultraviolet (NUV) band with an exposure time of 64 seconds.There is no source detected at the Gaia position.We use gPhoton (Million et al. 2016) to derive a 5 upper limit NUV < 20.6 mag.
We fit the spectral energy distribution (SED) using the Castelli & Kurucz (2003) atmosphere models included in pystellibs2.We calculate synthetic photometry with pyphot3 and sample over stellar parameters with emcee (Foreman-Mackey et al. 2013) for 10000 iterations.We first fit the SED with the extinction from 3-dimensional dust maps "Combined19" dust map and use a Gaussian prior on the distance using the Bailer-Jones et al. (2021) photogeometric distance.Since J0946 and J1640 are in the southern sky, we use extinction maps from Lallement et al. (2022) rather than the Drimmel et al. (2003) map used by mwdust.Next, we re-fit the SED with the extinction as a free parameter.We also fit for  2.
a binary star SED model with fixed extinction and distance.
Finally, we consider a SED model with graphitic circumstellar dust using DUSTY (Ivezic & Elitzur 1997;Ivezic et al. 1999), following the procedure described in Neustadt et al. (2024), to model the IR excess for J0946 and J1640.Figure 6 shows the SED fits for the three targets and Tables 3-5 report the MCMC posteriors.

J0946: Gaia DR3 5405789050140488320
The single-star SED fits for J0946 predict a stellar radius  = 33±1  ⊙ and  eff = 4150 +50 −40 K.The single star fits (Figure 6a) show that a UV observation would rule out a massive main sequence companion of spectral type B4V ( = 5.1  ⊙ ), but lower mass main sequence stars with  ∼ 3  ⊙ can still hide below the luminosity of the giant star.Both single-star SED fits are poor fits in the IR and near-IR, especially for the WISE W3 and W4 filters, which we exclude when fitting the SED.
The two-star SED fit (Figure 6b) does not improve the fit in the IR, but does show that the SED could be fit by a slightly cooler giant ( 1 = 34.2+0.6 −0.7  ⊙ ,  eff,1 = 3980 +30 −50 K) and a hotter, smaller secondary ( 2 = 5 +3 −2  ⊙ ,  eff,1 = 6200 +800 −700 K).A fourth SED model including circumstellar dust is more consistent with the observed IR excess (Figure 6c).In this model, the giant star is slightly cooler ( eff = 4447 K) with a smaller radius ( = 28  ⊙ ), as compared to the singlestar fits.The interstellar extinction is lower for this model (  = 1.46 mag), and consequently, the stellar luminosity is lower,  = 282  ⊙ .The circumstellar dust around J0946 could be produced in a radiatively driven wind which would dominate the observed spectrum in the IR for dust optical TABLE 3 SED fit results for J0946.We fit the SED with four models.The "Free" model includes extinction as a free parameter while the "Fixed" model uses the value from 3-dimensional dust maps.We also fit a binary model with two stars.The extinction and distance are fixed for this model.Finally, the "Dust" model includes circumstellar dust with optical depth   and temperature  dust . dust is the distance from the giant to the dust shell.The core mass of the photometric primary is computed from Equation 7 using the luminosity from the SED model and assuming a Solar metallicity.SED fits are shown in Figure 6.We use the results of the model including circumstellar dust to set priors on the giant star in our PHOEBE models (Section §5).
First ascent late-type red giants can also have episodes of internal mixing that produce circumstellar dust associated with lithium enrichment (de La Reza et al. 1996;Jasniewicz et al. 1999).However, the timescale of these events is short (≲ 10 5 years), and only a few percent of giants are observed to be in this phase.Table 3 reports the best fitting dust temperature,  dust = 822 K, and the optical depth at  = 0.55 m is   = 0.18.Mallick et al. (2022) fit the SEDs of five Li-rich giants using a DUSTY model with silicate grains.All five of their targets have average dust temperatures < 200 K, much cooler than the best-fitting dust temperature for J0946.The optical depth in their systems is also smaller by a factor of 10.Hotter and more optically thick dust indicates that the dust is closer to the star and the mass loss episode occurred more recently or that the dust is forming in a wind.High-resolution spectra of J0946 could be used to measure the Li abundance and make more direct comparisons to the Li-rich giants with IR excesses.Some circumstellar dust could also be created during episodes of mass transfer (e.g., Debes et al. 2012;Lü et al. 2013).The large dust shell radius,  dust = 1400 ± 300  ⊙ , indicates that the dust is actually circumbinary, providing further evidence for a history of mass transfer in addition to the circular orbits and ellipsoidal variability.
It is also possible that the WISE photometry has been contaminated by a nearby bright star.Figure 7 shows the 2MASS   and WISE W1-W4 image cutouts.IRAS 09450−5056 is a bright (W1 = 5.4 mag) star ≈ 1. ′ 2 from J0946.The diffraction spike from IRAS 09450−5056 crosses J0946 in the W1 and W2 images.The W1 and W2 photometry measurements have an "A" quality flag, meaning that the source is detected with SNR>10, but J0946 is flagged with possible flux contamination from the diffraction spike in the W1-W3 bands.While the diffraction spike contamination is most relevant for variability analysis, and many targets are over-flagged (Cutri et al. 2013), systematic errors in the WISE photometry could explain the apparent IR excess.
If the mass transfer has stripped the photometric primary, the core mass represents an absolute lower limit on the primary mass.For giant stars with degenerate cores, the core mass   is related to the luminosity  by at Solar metallicity (Boothroyd & Sackmann 1988).Table 3 includes the core mass for each SED model luminosity.
If we assume that the photometric primary is filling its Roche-lobe,  1 =  Roche , we can use the Eggleton approximation for the volume-averaged Roche-lobe radius to find the minimum primary mass that would be consistent with the SED radius without Roche-lobe overflow.Starting from the Eggleton (1983) approximation for the Roche-lobe radius, we write the semimajor axis in terms of the Roche-lobe radius,  Roche , and the mass ratio,  =  2 / 1 , The orbital period is also known, so This can be solved to determine the range of mass ratios where  -Spectral energy distributions (SEDs) of J0946, J1640, and J2017.The blue points in each panel show synthetic photometry from the Gaia Synthetic Photometry Catalog (GSPC, Gaia Collaboration et al. 2023b).For J1640, we show the upper limit for the GALEX NUV band in panels (d) and (e).For J2017, we show the archival IUE spectrum in panels (g) and (h).The inset on panel (h) shows a zoomed-in view of the IUE spectrum.The top row shows the single-star SED fits.The red lines show the "Fixed" model where   is fixed at the values from the 3-dimensional dust maps and the black lines show the "Free" model where extinction is a free parameter.For comparison, the purple and pink lines show the spectrum of a B4V and a B9V star, respectively extinguished by the same amount of dust as the free-  model (solid lines) and fixed-  model (dashed lines).The middle row (panels b, e, and h) show binary star SED fits using the fixed extinction from the 3-dimensional dust maps.The red and blue lines show the contribution from each component and the combined model is in black.J0946 and J1640 show evidence of significant IR excesses in the WISE photometry.In panels (c) and (f) we show models including circumstellar dust computed using DUSTY (Ivezic et al. 1999).The combined model (black) is broken down into the intrinsic stellar flux (blue), the attenuated flux (dashed red), the scattered flux (dotted red), and dust emission (solid red).
Roche-lobe overflow would be expected for a given primary mass.Figure 8 shows this radius relative to the radius estimates  SED for the dusty SED model of J0946 and for the SED models of J1640 and J2017.For the J0946 system, we find that  SED <  Roche for all mass ratios, which means that the giant can be entirely stripped without overflowing its Rochelobe.

J1640: Gaia DR3 5968984160290449024
The SED of J1640 also shows a strong IR excess (Figure 6d).As with J0946, we exclude the W3 and W4 filters from the single-star and binary-star SED fits that do not include circumstellar dust (Panels d and e of Figure 6).As compared to J0946, J1640 is cooler and larger, though the estimated parameters differ significantly between SED models (Table 4).For the model where the extinction is fixed from the 3-dimensional dust maps, the giant is predicted to have radius  = 59 +4 −5  ⊙ and effective temperature  eff = 3630 +20 −20 K.The model with free extinction prefers substantially more interstellar extinction,   = 3.39 mag as compared to   = 2.44 mag, and has a larger stellar luminosity and higher temperature as a result.
Figure 8 shows the Roche-lobe radius relative to the SED radius for the fixed extinction model and the model with circumstellar dust.For the fixed extinction model, if the primary is entirely stripped, the stellar radius is larger than the maximum radius before Roche-lobe overflow for  < 1.The model including circumstellar dust predicts a smaller stellar radius, and the giant can be totally stripped without overflowing its Roche lobe for this model.
Even with a GALEX upper limit, the high extinction limits our ability to rule out massive main sequence companions.Figure 6d shows the SED of a B4V and B9V with the same amount of extinction found in the SED models.While we can rule out a B4V companion in the model with fixed extinction, a B4V companion could not be ruled out if the extinction is There is an IR bright star (IRAS 09450−5056, W1 = 5.4 mag) ≈ 1. ′ 2 away from J0946.The diffraction spike may affect the photmetry of J0946 in the W1 and W2 bands.For J1640, there is a nearby star with apparent W1 = 7.4 mag ≈ 10. ′′ 5 away marked with the blue crosshairs.The circles in each panel show the approximate PSF for each filter (Cutri et al. 2003;Wright et al. 2010).
= 3.39 mag.For either model, we find that a B9V star would be consistent with the upper limit.
The binary star solution for J1640 suggests a hotter, less luminous companion with a luminosity ratio < 0.1.A Swift UVM2 observation (effective wavelength  eff = 225 nm) could be used to place stronger UV limits and search for luminous main sequence/subgiant companions, though this is made challenging due to the high extinction towards this target (see Section §6).
As with J0946, only the model including circumstellar dust is able to describe the IR emission of J1640 (Figure 6f).The dust temperature is  dust = 630 +30 −20 K, and the optical depth is   = 0.68 +0.09 −0.08 , suggesting a moderate amount of dust absorption and scattering.The dust in this system is also circumbinary with  dust = 4000 ± 500  ⊙ .We use the stellar parameters from this model to set priors on our light curve fit in Section §5.
There is also evidence for systematic issues with the WISE photometry of J1640. Figure 7 shows the IR images of J1640.There is a second bright WISE target, ALLWISE J164023.34−411619.9 ≈ 10. ′′ 5 away from J1640.In the 2MASS -band image the two stars are clearly resolved but they are blended in the WISE W1-W4 photometry.The WISE catalog flags the W1 and W2 photometry as being potentially contaminated by this star.The star is also a long-period photometric variable in Gaia DR3 (GDR3 5968984164630392960, Lebzelter et al. 2023) with a photometric period of  = 753 days, but it was not included in the Gaia Focused Product release sample.While the two stars may be blended in the WISE photometry, the separation is large enough that no contamination is expected in the ASAS-SN light curves or the Gaia radial velocities.

J2017: Gaia DR3 2053893497434322048
J2017 (V2012 Cyg, HD 332077) is one of the most luminous targets in the Gaia FPR LPV catalog (Figure 5), and the singlestar SED models both suggest that the photometric primary has  3, but for J2017.Since this target does not have an IR excess we do not fit a model including circumstellar dust.Since this target has luminosity > 10 4  ⊙ , we compute the core mass using Equation 10 6g).The model where extinction is a free parameter predicts a higher extinction than the value from the mwdust extinction maps, and consequently predicts a higher stellar luminosity and larger radius.The binary star SED model for J2017 prefers a solution with a luminosity ratio near unity (Figure 6h).This model can be rejected because J2017 is observed as a single-lined spectroscopic binary.Unlike J0946 and J1640, J2017 has no IR excess, so we do not fit a model that includes circumstellar dust.J2017 has previously been characterized as an 'S-star' (Stephenson 1984), a transition object between M giants and carbon stars.Many S-stars show evidence of technetium (Tc) in their atmospheres, which is produced in the s-process and is short lived.However, no Tc is observed in the atmosphere -Roche-lobe radius relative to the radius from the SED,  SED , as a function of mass ratio, .Where the curve is above 1, the primary can have the given mass and SED radius without Roche-lobe overflow.Systems below the gray line have giants that overflow their Roche-lobes.of J2017 (Jorissen & Mayor 1992), and Tc-deficient S-stars have been shown to have a higher binary fraction (Brown et al. 1990).
Based on the SED fits, the luminosity of J2017 is > 10 4  ⊙ while equation 7 is only valid for lower mass stars with luminosities ≲ 2000  ⊙ .For red supergiants with carbon and oxygen cores, the core mass-luminosity relation is (Paczyński 1970).Hot bottom burning provides extra luminosity, breaking the linear relationship (e.g., Wagenhuber & Groenewegen 1998), but we use the Paczyński (1970) relation as a first order estimate.Table 5 reports the core mass of J2017 using this relation.
Jorissen & Mayor (1992) use the spectroscopic orbit and an upper limit from the International Ultraviolet Explorer (IUE, Boggess et al. 1978) to suggest that the companion to J2017 is an A-star.Follow-up low resolution spectra taken with the International Ultraviolet Explorer (IUE, Boggess et al. 1978) found evidence of an A-type star spectra from 173-198 nm, but Ake et al. (1992) note that a simultaneous fit to the long and short wavelength regions with a A-star plus giant star model is not successful.We include the IUE spectrum of the long wavelength region from May 1996 for comparison in Figure 6g and a zoomed-in view in the inset of Figure 6h.The flux from the IUE spectrum exceeds what is expected from our SED model, suggesting that there is a luminous companion.There is also evidence of MgII line emission at 280 nm, which is associated with chromospheric activity in FGK stars (Buccino & Mauas 2008).
A main sequence A6 companion with mass  2 = 1.8  ⊙ would not produce the observed radial velocity semi-amplitude if  1 > 12  ⊙ , as suggested by our SED fits.Jorissen & Mayor (1992) instead suggest that the giant is very low mass < 0.5  ⊙ with  1 < 245  ⊙ .This radius is inconsistent with our SED results (Table 5).
Our picture of a massive giant primary can be reconciled with the stripped star+A-star scenario proposed by Jorissen & Mayor (1992) if the Gaia parallax is biased by the binary motion.The Gaia parallax is  = 0.4930 ± 0.0535 mas after using the zero-point correction from Lindegren et al. (2021).The astrometric signal from the binary motion is then ≃ 0.801 mas/sin(i) (Equation 6).While a more detailed numerical simulation of the orbit incorporating the Gaia scanning law and observation times is needed to fully understand the parallax bias (e.g., Thompson et al. 2019), if the true Gaia parallax is larger than the reported value, the actual stellar radius will be smaller than the result from our SED fit.
New ultraviolet observations could be used to search for additional evidence of a luminous companion and validate the A-star signal identified by Ake et al. (1992).Our SED fits suggest that a main sequence A6V star would exceed the flux from the giant in the Swift UVM2 filter.However, the A-star companion would have a magnitude of ∼ 20.1 mag in the UVM2 band using the extinction and distance from the fixed extinction model (Table 5), which is approaching the sensitivity limit for the instrument (e.g., Molina et al. 2020) unless the binary motions have biased the parallax.While we are unable to provide a definitive characterization of this system here, it is unlikely that the companion is a compact object.We remove this target from further consideration.

MODELS OF ELLIPSOIDAL VARIABILITY
By modeling the ellipsoidal variable light curve and radial velocity curve simultaneously, we can constrain the orbital inclination and mass of the unseen companion.We used PHysics Of Eclipsing BinariEs Eclipsing Binary (PHOEBE, Prša & Zwitter 2005;Conroy et al. 2020) to fit the ASAS-SN light curves and Gaia radial velocities for J0946 and J1640.We exclude J2017 because of the evidence of an A-star companion in the UV spectrum (Ake et al. 1992).
For J0946 and J1640 we start by using PHOEBE to fit a "dark companion" model, where the secondary is fixed to be small  2 = 3 × 10 −6  ⊙ , and cold,  eff,2 = 300 K, and ignore effects of irradiation and reflection.We fix the mass of the photometric primary to be the core mass   computed using the SED model with circumstellar dust (Tables 3 and 4 and Equation 7) so that the resulting companion mass is the lower limit for a maximally stripped red giant.
To obtain an initial estimate for the fit, we use the PHOEBE differential evolution optimizer.For the MCMC fit, we set Gaussian priors on the effective temperature and photometric primary radius from the SED fit and a prior on the period from the Gaia FPR RV solution.We run the PHOEBE models for 5000 iterations and use 26 walkers.We visually inspect the walker distributions to select a burn-in period for each target.
We also run a second set of PHOEBE models where the companion is luminous and the effective temperature  eff,2 and radius  2 are allowed to vary.The photometric primary mass  1 is also a free parameter in this model and 32 walkers are used for the MCMC fit.

J0946: Gaia DR3 5405789050140488320
Figure 9 shows the corner plot for the dark companion PHOEBE model of J0946.Table 6 reports the MCMC posteriors.As compared the circumstellar dust SED model, the dark companion PHOEBE model predicts a smaller giant ( = 20.4+0.3 −0.3  ⊙ ) with an effective temperature that agrees with this SED within 1.The mass ratio is large  > 10, and the companion has  2 = 3.8 +0.3 −0.2  ⊙ .Based on the SEDs (Figure 6), ultraviolet observations could be used to rule out main sequence stars of this mass.
The PHOEBE model where the companion is treated as a luminous star finds a solution with two red giants.The orbital elements and masses are similar between the two models.The photometric primary in this model has  1 = 28 +1 −1  ⊙ and  eff,1 = 4430 +60 −150 K, which agrees with the circumstellar dust SED model.The companion is also large  2 = 22 +3 −2  ⊙ and hotter than the primary.The luminosity ratio this system is  2 / 1 = 0.85.While it seems likely that a companion with this luminosity ratio would be detected in the Gaia RVS spectra, the spectral lines of the smaller hotter companion could be broad and not detected if it is rapidly rotating.For comparison, V723 Mon, has a  2 / 1 = 0.74 and is only detected as a SB1 since the projected rotational velocity of the companion is high ( 2 sin  = 70 ± 10 km/s, Jayasinghe et al. 2021b;El-Badry et al. 2022).If the companion is also a luminous giant with this radius and effective temperature, its flux would exceed that of the photometric primary at wavelengths shorter than ∼ 6000 Å, so additional optical spectra could be used to search for evidence of rotationally broadened spectral lines.

J1640: Gaia DR3 5968984160290449024
Figure 10 shows the corner plot and light curve fits for the dark companion PHOEBE models of J1640.For both the dark The orbital inclination in both of the PHOEBE models for J1640 is predicted to be ∼ 50 • and the companion mass  2 = 3.7  ⊙ .As with J0946, if the unseen companion is actually a main sequence star with  2 = 3.7  ⊙ it could be searched for with ultraviolet observations, though this would be more challenging for J1640 due to the high extinction (  = 2.44 mag from the Lallement et al. (2022)

maps).
In the luminous companion model, the secondary is a cool  eff,2 = 3800 +400 −200 , smaller star  2 = 4 +4 −2  ⊙ , but the secondary radius is poorly constrained.The luminosity ratio is small,  2 / 1 = 0.006 but we note that the properties of such a star are not consistent with expectations from single star evolution.A star with  2 = 3.7  ⊙ would cool to 3700 K after it evolves off the main sequence, at which point its radius would be ≳ 10  ⊙ .
The PHOEBE models of J1640 predict near circular orbit  = 0.01 ( = 0.03 for the luminous companion model), which is significantly smaller than the ellipticity in our RV orbit model from Section 3 ( = 0.05).Figure 11 shows the RV orbit model from Section §3 compared to the PHOEBE result.While the ASAS-SN photometry is well-fit by the PHOEBE model (Figure 10), the PHOEBE RV model does not fit the Gaia RVs.This suggests that there is a small relative phase offset between the RV orbit and the ellipsoidal modulations.
Apsidal motion could produce a relative phase offset between the Gaia RVs and the ellispodial modulations.Since the PHOEBE model prefers a circular eccentricity, the argument of periastron, , is unconstrained in the PHOEBE model of J1640.
If we instead adjust  of the RV model from Section §3 to match the phasing of the PHOEBE model, we find a change in argument of periastron Δ ∼ 15 • is needed.The median time of the ASAS-SN -band light curve is ≈ 5.2 years later than the median time of the Gaia RVs, so Δ/Δ ≃ 3 deg year −1 .For a synchronized binary, Sterne (1939) and Rosu et al. (2020) write the rate of apsidal motion as, where  1 and  2 are the second order apsidal motion constants for the primary and the secondary, respectively, and (1 −  2 ) 5 , and There is also a contribution to  from general relativity, but for this system  GR ≪  and it can be neglected.The values of  1 and  2 depend on the interior structure of the stars.If the companion is a main sequence star or a compact object we effectively have  2 = 0. Using the period and eccentricity from the RV fit (Table 2), the values of  1 and  1 =   from the SED (Table 4), and setting the companion to be a 3.7  ⊙ BH, the value of  1 would need to be  1 ≃ 0.145 in order to have  ≃ 3 deg year −1 .This is roughly twice as large as the value of  1 estimated for the stripped star 2M04123153+6738486 (El-Badry et al. 2022), and is outside of the range of the theoretical grids from Claret (2019).While apsidal motion could produce some change in  in the time between the Gaia RVs and the ASAS-SN photometry, it is unlikely to be responsible for the observed phase offset.If J1640 is actually a hierarchical triple, light travel time effects (LTTE) could contribute to the effects of apsidal motion and increase the apparent phase offset between the RVs and the ellipsoidal variations.There are triple systems with apsidal motion (e.g., Bozkurt & Deǧirmenci 2007), but these are generally shortperiod ( < 10 day) inner binaries.
It could also be the case that one or more Gaia RVs are poorly constrained or have underestimated uncertainties.We test this by re-fitting the RV model from Section §3 with emcee dropping all possible combinations of one, two, or three RV measurements.The initial walker positions for all MCMC runs are set using the rejection sampling results from TheJoker using the full set of RVs.
For each trial, we phase the ASAS-SN light curve using the period and the time of superior conjunction.To determine the deviation from the expected phasing, we fit an analytic model of the form motivated by the analytic model of Morris & Naftilan (1993) and compute the reduced  2  of the ASAS-SN light curve.RVs.While the two models produce a similar fit for J0946 (top), there is a significant offset between the two models of J1640 (bottom).This could be due to apsidal motion or a subset of poor RV measurements (Section §5).
While the majority of these trials do not affect the phasing of the light curve significantly, we find that there are some combinations of RV measurements that, when removed, produce a better phasing agreement between the light curve and the radial velocity measurements.
Figure 12 shows an example where three RV points have been removed.The RV fit looks nearly identical to the fit using all of the Gaia RVs, and the MCMC posteriors show the results are consistent within 1.In this case, the orbital period differs by ∼ 0.9 days, which can produce an appreciable phase shift over the course of the five year gap between the Gaia RVs and the ASAS-SN light curve.
Additional PHOEBE models that exclude subsets of RV measurements could be used to better understand which RV measurements contribute most to the poor fit shown in Figure 11.Ground-based radial velocity follow-up could also help further constrain the RV orbit of this target.An extended baseline could also place stronger constraints on the possibility of apsidal motion.Ultimately, however, the large RV residuals for the PHOEBE model of J1640 do not alter the overall characterization of this system as a high mass function ellipsoidal variable on the giant branch.

DISCUSSION AND CONCLUSIONS
We use the Gaia DR3 light curves and RVs from the recent Gaia Focused Product Release (Gaia Collaboration et al. 2023d) to identify eight targets with large binary mass functions, all of which are on the giant branch.The presence of a massive, unseen companion suggests that these systems are 1.20  for BHs than for neutron stars (Repetto et al. 2012;Mandel & Müller 2020), we would expect these orbits to be moderately eccentric following the supernovae.Still, tidal effects would be expected to circularize the binaries within timescales short compared to the evolution of the giant.
The circularization timescale is where  1 is the mass of the giant primary,  env is the mass of the envelope of the giant, and  2 is the mass of the companion (Verbunt & Phinney 1995;Zahn 1977).The quantity  is a dimensionless factor of order unity that contains physics related to the convective and viscous properties of the star.For J0946, even if the primary mass is typical of a red giant ∼ 1  ⊙ with  env =  1 /2, the binary would be expected to circularize within ∼ 75, 000 years for companions ≲ 10  ⊙ .High-resolution spectra could be used to make more direct comparisons to V723 Mon (Jayasinghe et al. 2021b) and 2M04123153+6738486 (Jayasinghe et al. 2022).El-Badry et al. (2022) used spectral disentangling to identify compan-ions and rule out compact object companions in the systems.The stripped red giants in these systems dominate the combined SED redder than ∼ 5000 Å, and identifying the subgiant companions is made more challenging due to their rotationally broadened spectral lines from being tidally locked.For J0946 and J1640, the high extinction could make obtaining spectra with high enough signal-to-noise spectra for disentangling challenging.A subgiant companion would be expected to dominate the SED at  ≲ 5000 Å.
In the SDSS -band (effective wavelength  = 3608 Å), the SED models of J0946 and J1640 predict extinction-corrected magnitudes for the giants of  ∼ 17.9 mag and  ∼ 19.3 mag respectively.Even if a subgiant companion contributed four times as much flux at this wavelength, the magnitude of the binary would still be fainter than 16.1 mag and 17.6 mag for J0946 and J1640, respectively.Multi-band Swift UVOT or Hubble Space Telescope UV photometry may be a more promising method to detect any subgiant companions hotter than the giant stars.
Searching for and characterizing non-interacting BHs requires careful consideration of false positives originating from mass transfer.Although Gaia DR3 was predicted to detect dozens to hundreds of BHs (e.g., Breivik et al. 2017;Chawla et al. 2023), only two strong candidates have been detected.While this could indicate a need for revision of binary population synthesis models, both of the Gaia astrometric BHs are within 1.2 kpc, suggesting similar systems almost certainly exist all over the Galaxy.Spectroscopic and astrometric surveys are the most promising pathways to identify new candidates.Gaia DR3 included spectroscopic orbit fit results for ∼ 185, 000 binaries and astrometric orbit fits for ∼ 130, 000 binaries.The next Gaia data release is expected to include time-series RVs and astrometric measurements for a larger sample of systems.This Gaia Focused Product Release includes the first Gaia epoch RVs for binary systems and serves as a demonstration of the potential for future Gaia data releases to aid in the search for non-interacting black holes.

APPENDIX TESS LIGHT CURVE OF J0131
In Section §2.1 we used the QLP, SPOC, and TGLC pipelines to download TESS light curves for the eight high  () targets in Table 1.J0131 shows evidence for additional eclipses in the Sector 58 SPOC light curve shown in Figure 13.We use a Lomb-Scargle (LS) periodogram (Lomb 1976;Scargle 1982) to identify a signal at  LS = 3.30843 and a box-least-squares (BLS) periodogram (Kovács et al. 2002) to identify a second  BLS = 6.40193.Figure 13 shows the TESS light curve folded at both periods.
Neither of these signals are present in the QLP light curves.The Sector 24 and 25 TGLC light curves show some evidence for additional variability, but it is noisier than the SPOC Sector 58 light curve.No other Sectors are available from the SPOC pipeline.The large TESS pixels mean that this variable signal could be produced by a nearby variable star.Following Rowan et al. (2023), we cross-match with variable star catalogs to search for the contaminant with a 5' radius from the Gaia position.There are 5 variable stars in the Zwicky Transient Facility variable star catalog (Chen et al. 2020) within 5'.One target, J013146.49+621953.5 is 95.′′ 1 away from J0131 and Chen et al. (2020) reports a period  ZTF = 3.36842.Since  ZTF ≈  LS , and the -band light curve shape of J013146.49+621953.5 is similar to the middle panel of Figure 13, the variability seen in the TESS light curve of J0131 is likely due to contamination.
The ZTF light curve of J013146.49+621953.5 does not show evidence for the  BLS = 6.40193 eclipsing signal, and none of the other four nearby ZTF variables have a similar period.We find no matches to the ASAS-SN (Jayasinghe et al. 2021a;Christy et al. 2023) or ATLAS (Heinze et al. 2018) variable star catalogs.Although it is possible that this signal is intrinsic to the J0131, making it a triply-eclipsing binary, since the signal is only present in some TESS Sectors and with some reduction pipelines, it seems more likely to be contamination from a nearby blended variable star.
This paper was built using the Open Journal of Astrophysics L A T E X template.The OJA is a journal which provides fast and easy peer review for new papers in the astro-ph section of the arXiv, making the reviewing process simpler for authors and referees alike.Learn more at http://astro.theoj.org.
Fig. 1.-The amplitude of photometric and radial velocity variability for the full Gaia FPR LPV sample.The blue line shows the cut described by Gaia Collaboration et al. (2023d) that separates binary stars from pulsating variables.
Fig. 4.-Radial velocity orbits for the five ELL variables with high  (  ) identified in the Gaia FPR sample.The orbital parameters are reported in Table2.
Fig. 5.-Gaia color-magnitude diagram (CMD) for the Gaia FPR catalog.Systems classified as pulsators (binaries) are shown in gray (red) based on the Equation 1.The three high  (  ) eclipsing targets are marked as red crosses, and the four ELL targets are colored circles.J1513 is not shown because the Gaia parallax uncertainty is /  < 5.The black line shows a MIST isochrone (Dotter 2016; Choi et al. 2016) for an age 10 9.4 years.

Fig
Fig.6.-Spectralenergy distributions (SEDs) of J0946, J1640, and J2017.The blue points in each panel show synthetic photometry from the Gaia Synthetic Photometry Catalog (GSPC, GaiaCollaboration et al. 2023b).For J1640, we show the upper limit for the GALEX NUV band in panels (d) and (e).For J2017, we show the archival IUE spectrum in panels (g) and (h).The inset on panel (h) shows a zoomed-in view of the IUE spectrum.The top row shows the single-star SED fits.The red lines show the "Fixed" model where   is fixed at the values from the 3-dimensional dust maps and the black lines show the "Free" model where extinction is a free parameter.For comparison, the purple and pink lines show the spectrum of a B4V and a B9V star, respectively extinguished by the same amount of dust as the free-  model (solid lines) and fixed-  model (dashed lines).The middle row (panels b, e, and h) show binary star SED fits using the fixed extinction from the 3-dimensional dust maps.The red and blue lines show the contribution from each component and the combined model is in black.J0946 and J1640 show evidence of significant IR excesses in the WISE photometry.In panels (c) and (f) we show models including circumstellar dust computed using DUSTY(Ivezic et al. 1999).The combined model (black) is broken down into the intrinsic stellar flux (blue), the attenuated flux (dashed red), the scattered flux (dotted red), and dust emission (solid red).
Fig. 7.-Cutouts of the 2MASS -band and WISE infrared images for J0946 (top) and J1640 (bottom).There is an IR bright star (IRAS 09450−5056, W1 = 5.4 mag) ≈ 1. ′ 2 away from J0946.The diffraction spike may affect the photmetry of J0946 in the W1 and W2 bands.For J1640, there is a nearby star with apparent W1 = 7.4 mag ≈ 10. ′′ 5 away marked with the blue crosshairs.The circles in each panel show the approximate PSF for each filter(Cutri et al. 2003;Wright et al. 2010).

Fig. 10
Fig. 10.-MCMC posteriors and light curve for J1640.andluminous companions models, the radius of the giant primary agrees with the results from the circumstellar dust SED with within 1.The dark companion PHOEBE model predicts an effective temperature of the giant  eff,1 = 3830 +40 −30 K, which is within 2 of the SED result.The effective temperature of the giant in the luminous companion model agrees with the SED result at the 1 level.The orbital inclination in both of the PHOEBE models for J1640 is predicted to be ∼ 50 • and the companion mass  2 = 3.7  ⊙ .As with J0946, if the unseen companion is actually a main sequence star with  2 = 3.7  ⊙ it could be searched for with ultraviolet observations, though this would be more challenging for J1640 due to the high extinction (  = 2.44 mag from theLallement et al. (2022) maps).In the luminous companion model, the secondary is a cool  eff,2 = 3800 +400 −200 , smaller star  2 = 4 +4 −2  ⊙ , but the secondary radius is poorly constrained.The luminosity ratio is small, Fig. 11.-Comparison of the orbit fit from the Keplerian model (Section §3) and the joint PHOEBE model fit to the ASAS-SN photometry and GaiaRVs.While the two models produce a similar fit for J0946 (top), there is a significant offset between the two models of J1640 (bottom).This could be due to apsidal motion or a subset of poor RV measurements (Section §5).

Fig. 12 .
Fig. 12.-Comparison of the MCMC posteriors for the Keplerian orbit fit of J1640 using the full set of Gaia RVs (red) and one with three points removed (blue).The posteriors are nearly identical, but the slight change in the orbital period is enough to improve the phasing of the ASAS-SN light curve compared to the expectation from ellipsoidal variability.The RV orbits are shown in the upper right and the three dropped measurements are marked in green.

Fig
Fig. 13.-Sector 58 light curve of J0131 downloaded from the SPOC pipeline using the PDCSAP flux.The top panels shows the unfolded light curve.The middle and bottom panels show the light curve folded at the periods from the LS and BLS periodograms, respectively.The red points show the binned light curve.
SN -band (black)and Gaia -band (green) light curves of the eight targets with  (  ) > 1  ⊙ .Three targets (J0131, J1030, and J2059) show evidences of eclipses consistent with the spectroscopic period.The remaining targets are our candidate ellipsoidal variables.

TABLE 2
RV orbit fit results for the five ELL candidates.The RV orbits are shown in Figure

TABLE 4
Same as Table3, but for J1640.We use the stellar parameters of the SED model with circumstellar dust to set priors on our PHOEBE model in Section 5.