The First High-contrast Images of Near High-mass X-Ray Binaries with Keck/NIRC2

Although the study of X-ray binaries has led to major breakthroughs in high-energy astrophysics, their circumbinary environment at scales of ∼100–10,000 au has not been thoroughly investigated. In this paper, we undertake a novel and exploratory study by employing direct and high-contrast imaging techniques on a sample of X-ray binaries, using adaptive optics and the vortex coronagraph on Keck/NIRC2. High-contrast imaging opens up the possibility to search for exoplanets, brown dwarfs, circumbinary companion stars, and protoplanetary disks in these extreme systems. Here we present the first near-infrared high-contrast images of 13 high-mass X-ray binaries located within ∼2–3 kpc. The key results of this campaign involve the discovery of several candidate circumbinary companions ranging from substellar (brown dwarf) to stellar masses. By conducting an analysis based on Galactic population models, we discriminate sources that are likely background/foreground stars and isolate those that have a high probability (≳60%–99%) of being gravitationally bound to the X-ray binary. This paper seeks to establish a preliminary catalog for future analyses of proper motion and subsequent observations. With our preliminary results, we calculate the first estimate of the companion frequency and the multiplicity frequency for X-ray binaries: ≈0.6 and 1.8 ± 0.9, respectively, considering only the sources that are most likely bound to the X-ray binary. In addition to extending our comprehension of how brown dwarfs and stars can form and survive in such extreme systems, our study opens a new window to our understanding of the formation of X-ray binaries.

X-ray binaries are semidetached binary systems in which a compact object (white dwarf; WD, neutron star; NS, or stellar-mass black hole; BH) accretes material from a donor star.These systems undergo several extreme physical phenomena, such as processes acting predominantly in soft X-rays (e.g., Khargharia et al. 2010;Tetarenko et al. 2021), and detectable X-ray pulsations (e.g., Lutovinov et al. 2005).
The variations in physical processes among different X-ray binaries are directly linked to the mass of the donor star.Over 90% of these systems can be classified into two distinct categories: high-mass X-ray binaries (HMXB; M donor ≳ 8 M ⊙ ) and low-mass X-ray binaries (LMXB; M donor ≲ 1,5 M ⊙ ; e.g., Tauris & van den Heuvel 2006).LMXBs are relatively old systems (> 10 9 yr) harboring a K-M spectral type donor star, where the process of mass transfer is believed to be triggered by Roche-lobe overflow (RLO; e.g., Savonije 1978).RLO is triggered either by stellar evolution or by angular momentum loss (e.g., Paczyński 1967;Verbunt & Zwaan 1981;Strohmayer 2002;Justham et al. 2006;Chen & Podsiadlowski 2016;Seto 2018; Van et al. 2019).The transferred mass then agglomerates to form an accretion disk around the compact object, giving rise to transient accretion and X-ray emission (e.g., Charles & Coe 2006).
As for HMXBs, they are generally thought to be younger systems (≲ 10 7 yr) harboring a massive O-B spectral type donor star.The transferred and accreted matter is thought to predominantly come from the capture of a fraction of the stellar winds ejected from the donor star (e.g., Mohamed & Podsiadlowski 2007;Abate et al. 2013;El Mellah et al. 2019).There are two subcategories of HMXBs relevant to this work.Firstly, we emphasize Be/X-ray binaries (BeXRB), wherein the donor star is a fast-rotating Be star.In these systems, the X-ray emission is mainly triggered by the compact object passing through a diffuse and gaseous circumstellar disk surrounding the Be star (known as a decretion disk; e.g., Okazaki et al. 2002;Martin et al. 2011;Rímulo et al. 2018;Kravtsov et al. 2020).Secondly, we highlight supergiant fast X-ray transients (SFXT; Negueruela et al. 2006), characterized by the presence of a supergiant donor star, and by fast transient X-ray flaring activity within the system (likely induced by a NS; e.g., Sidoli 2013;Ducci et al. 2019).
X-ray binaries are important touchstone objects for high-energy phenomena in astrophysics.They have been widely used to study several high-energy astronomical phenomena, including accretion physics (e.g., Done et al. 2007;Kara et al. 2019) and outflow/jet processes (e.g., Markoff et al. 2001;Fender et al. 2004;Mooley et al. 2018).However, the immediate surroundings of X-ray binaries, at the scale of ∼100-10,000 astronomical units (hereafter au), have been poorly studied.This paper undertakes a pioneering exploration of the circumstellar environments of X-ray binaries through the application of adaptive optics (AO) and direct/high-contrast imaging techniques.The goal is to probe a variety of phenomena ranging from protoplanetary disks to debris disks and fallback disks, and particularly to search for wide-orbiting circumbinary companions (CBCs) -be they exoplanets, brown dwarfs, or stars.
Considering the discovery of planetary-mass CBCs orbiting both binary systems (e.g., Bakos et al. 2007;Desidera & Barbieri 2007;Eriksson et al. 2020) and compact objects (WD or pulsars; e.g., Wolszczan & Frail 1992;Sigurdsson et al. 2003;Spiewak et al. 2018; Vanderburg et al. 2020;Blackman et al. 2021), it is not unfounded to expect CBCs orbiting X-ray binaries.A recent study argued that X-ray binaries could host planetary systems in close orbits detectable via X-ray eclipses (Imara & Di Stefano 2018).In this paper, we explore wider orbits (∼100-10,000 au), as the increased number of interactions within the system could lead to the ejection of potential CBCs from the direct environment of the X-ray binary (e.g., Bonavita et al. 2016).
In Prasow-Émond et al. (2022), we presented the first set of observations from a pilot study aiming to survey all X-ray binaries amenable for direct imaging within ∼3 kpc.We first targeted a γ Cassiopeiae-like X-ray binary harboring a Be donor star, RX J1744.7−2713, for which we had observations from two different bands and two epochs.We unveiled the presence of three potential CBCs within this system, exhibiting a strong likelihood of being stellar-mass CBCs.Here, we present the first L ′ -band high-contrast images of 13 other systems and conduct a preliminary statistical analysis derived from the results of the first epochs of observations.
The paper is organized as follows: Section 2 presents the sample and how it was constructed.Section 3 presents the near-infrared observations, and the data reduction and processing.Section 5 presents the first high-contrast images of the observed X-ray binaries.In Section 6, we analyze the images and explore the nature of the detection.Finally, in Section 7, we discuss our results and their implications.au scales.We chose a distance limit of ∼3 kpc within our Galaxy, which enables the detection of structures and objects located within a couple of thousands of au from the X-ray binary.The outer limit (10,000 au) corresponds to the approximate limit of the Keck/NIRC2 field of view (foV).
2. Brightness and Adaptive optics.The donor star must be bright enough (I < 9 − 10 mag) for the AO loop to be closed.

3.
Age.We targeted young (≲ 100 Myr) X-ray binaries to favor the detection of sub-stellar CBCs, considering the steep decline in planet brightness with time (e.g., Burrows et al. 2001).This limited us to bright LMXBs, and HMXBs with a massive O/B donor star.
4. Visibility.We selected X-ray binaries visible from the W. M. Keck Observatory at the time of our observations (Keck Observatory Semesters 2017B and 2020A).
Applying these criteria narrowed down the initial list from ∼300 X-ray binaries to 19, out of which 14 were observed between 2017 and 2020 using Keck/NIRC2 (see Section 3.1 for more details).Our sample includes both HMXBs and LMXBs (e.g., MAXI J1820+070, V404 Cyg, 1A0620-00); however, we have only observed HMXBs to date due to observational constraints.Note that X-ray binary surveys are far from complete, as Xray binaries can be undetectable in quiescence (e.g., Bird et al. 2007;Belczynski & Ziolkowski 2009).New Xray binaries have also been discovered since we did our sampling (Avakyan et al. 2023;Neumann et al. 2023).Nonetheless, we assume that our sample is as complete as possible and that our results accurately represent the known population of X-ray binaries.
Figure 1 displays the position of the 14 observed targets onto the sky (Aitoff projection).It also indicates the distance from the observer (in kpc) and, if known, the nature of the system's compact object.In Section 4, we present a brief literature review for each of the observed X-ray binaries.Table 1 summarizes the known relevant physical properties of the systems, namely the X-ray binary type, X-ray emission class, donor star spectral type and compact object type.Additional relevant physical properties can be found in the Appendix (see Table 5).

Keck/NIRC2 Observations
On September 8, 2017, we observed four HMXBs from our survey sample using the Keck/NIRC2 vortex coronagraph (Mawet et al. 2005;Serabyn et al. 2016) in pupil-tracking mode in L ′ -band (λ = 3.776 µm, ∆λ = 0.700 µm; PI: Mawet), and with the narrow camera (plate scale of 9.971 ± 0.004 mas pixel −1 ; Service et al. 2016).On January 3, 2018, we observed three additional HMXBs using a similar setting.Due to the successful and promising preliminary results of this initial campaign, we were awarded three supplementary nights of observation on July 11, 12, and 13, 2020 (PI: Fogarty).On the first night, we observed one additional target and re-observed two targets (RX J1744.7−2713 and γ Cas) using a similar setting.However, due to saturation, we had to downscale the frame size of γ Cas from 1024×1024 pixel 2 to 512×512 pixel 2 .On July 12, 2020, we obtained data for three other HMXBs in L ′band, in addition to re-observing RX J1744.7−2713 in K s -band (λ = 2.146 µm, ∆λ = 0.311 µm).Finally, on July 13, 2020, we obtained data for three other HMXBs, totaling observations for 14 out of the 19 X-ray binaries in the sample.
During the observations, we used the Quadrant Analysis of Coronagraphic Images for Tip-tilt sensing (QAC-ITS; Huby et al. 2017) to make tip-tilt adjustments to maintain precise centering of the target on the vortex focal plane mask.The observations were AO-assisted using the Shack-Hartmann wavefront sensor (which performs wavefront sensing in R-band) in 2017, 2018 and the last night of 2020.For the first two nights of our 2020 observations, we opted for the Pyramid wavefront sensor (PyWFS) instead.It performs wavefront sensing in H (Wizinowich et al. 2000;Bond et al. 2018), which is more advantageous for the redder targets in our sample.
A summary of the observing log is presented in Table 2.

Data Reduction
Similarly to Prasow-Émond et al. (2022), we performed data reduction using the Vortex Image Processing (VIP) and NIRC2 Preprocessing packages (Gomez Gonzalez et al. 2017).To obtain a preprocessed data cube, we proceeded as follows: (1) flat-fielding of the frames, (2) bad pixel masking using the dark frames, (3) determination of the vortex center for each frame followed by cropping the science cube around the mean center, (4) removal of sky contribution via a Principal Component Analysis (PCA)-based technique, and (5) image registration to align the quasi-static speckle pattern across frames.
After acquiring the preprocessed data cube, we applied a PCA-based Angular Differential Imaging (ADI; Figure 1.Position of the 14 observed X-ray binaries on the sky with equatorial coordinates and Aitoff projection, color coded with distance from the observer in kpc.The nature of the compact object is illustrated with different markers: squares for white dwarfs, crosses for neutron stars, asterisks for black holes and circles for cases where the nature is either unavailable or uncertain.The red dotted line shows the approximate coordinates of the galactic plane.Marois et al. 2006) algorithm to obtain high-contrast images.Subsequently, we generated several images using two algorithms in VIP (annular PCA and full-frame PCA) for a broad range of principal components (from 1 to 50).This was done to ensure consistency in the detection of sources within the images (i.e., that the source was detected regardless of the number of principal components).We then listed sources with a signal-to-noise ratio (SNR) greater than 5 (or with a signal exceeding 4σ) and determined the optimal number of principal components (n comp ) that maximized the SNR for each source.
Due to the poor quality of the PSF in our 2020 observations, we used the synthetic 2D PSF for fitting processes related to these observations.However, we kept the original PSF for our 2017 and 2018 observations.We conducted tests using both synthetic and real PSFs on the 2017 and 2018 data, yielding consistent results.Consequently, the use of a synthetic PSF does not affect significantly the derived parameters.
The upcoming sections detail the acquisition of parameters used in the magnitude calculations.

Determining the Distance from the Observer
The distance from the observer is presented in the second column of Table 5 in the Appendix.Distances for 1H2202+501, 4U 2206+543, 4U 1700−377, IGR J17544−2619, Cyg X-1, X Per and Vela X-1 were ob- tained from Zhao et al. (2023), which uses parallax measurements to infer distances either using an inversion or Bayesian approach for a catalog of X-ray binaries.For the other X-ray binaries, the distance was estimated through a photogeometric calculation (Bailer- Jones et al. 2021).It was calculated using the parallax measurement and its uncertainty (geometric), as well as the G magnitude and the BP-RP color (photometric) from the Gaia Data Release 3 (DR3; De Angeli et al. 2023;Montegriffo et al. 2023).Table 3 presents the Gaia DR3 ID for each target.In cases where an object's distance was sourced from the literature, the respective reference is cited in the literature review of Section 2. As previously mentioned, this study targets X-ray binaries within ∼2-3 kpc accessible with Keck/NIRC2.

Determining the Apparent Magnitude of the Central X-Ray Binary
Observing X-ray binaries in the L ′ -band of Keck/NIRC2 is not standard practice, making the direct determination of the true apparent magnitudes of the central X-ray binaries (m XRB ) unfeasible.Nonetheless, we leveraged the W 1 filter from the Wide-field Infrared Survey Explorer (WISE), which has a central wavelength similar to Keck/NIRC2 L ′ -band (λ = 3.353 µm).Using the WISE Source Catalog (Wright et al. 2010), we approximated m XRB for all X-ray binaries.These values are presented in the fifth column of Table 5 in the Appendix.To account for the differences between the two filters, we considered an uncertainty of 0.5 mag on m XRB .

Determining the Age of the System
Most X-ray binaries in our sample lack age estimates in the literature, except for RX J1744.Lyubimkov et al. 1997).For the remaining X-ray binaries, we established an upper limit using basic evolutionary models and the spectral type of the donor star (see Table 1).We found the maximum age that the donor star can reach before exploding in supernovae.However, this approach provides only an approximate estimation of the age.These values are presented in the sixth column of Table 5 in the Appendix.

Determining the Errors
Errors in the fit parameters (θ, ρ, and f 1 ) were estimated using an injection/recovery approach (see Prasow-Émond et al. 2022).This approach relies on injecting synthetic sources with known parameters into the images.Subsequently, the same optimization method was applied (see Section 3.3), and the error was determined as the difference between estimated and known parameters.This method was employed across a range of parameter values to ensure consistency; it was observed that errors were more pronounced for smaller offset values (i.e., those closer to the central X-ray binary).Additional sources of uncertainty were taken into consideration in cases where the value of a parameter is expected to remain consistent between two sets of observations, specifically astrometric parameters in two different bands.The dominant source of uncertainty was defined as the total uncertainty.

γ Cassiopeiae
γ Cassiopeiae -also known as 2S 0053+604 (hereafter γ Cas) -harbors a well-studied optical component classified as a Be star (Moffat et al. 1973).Its X-ray luminosity (∼10 32 -10 33 erg s −1 ; Raguzova & Popov 2005) is higher than the typical luminosity for O/B stars (∼10 30 erg s −1 ), but too low to be a Be/NS system (∼10 34 erg s −1 in quiescence; Shrader et al. 2015).The nature of the system can be explained by two hypotheses: (1) the system is a HMXB, involving accretion onto a WD or a fast-spinning NS (Postnov et al. 2017); or (2) the excess of X-ray emission stems from physical processes in the high atmosphere of γ Cas (Kubo et al. 1998;Robinson & Smith 2000).Though γ Cas is not confirmed as being a HMXB, its resemblance to other sources, referred to as γ Cas analogs, warranted its inclusion in our sample.Located at a distance of 0.19 ± 0.02 kpc, it is the nearest system in our sample.Its proximity and proper motions (25.7 ± 0.5 mas/yr in RA, −3.9 ± 0.4 mas/yr in Dec; Perryman et al. 1997) allowed us to conduct a proper motion analysis within the interval between our two observation sets (see Section 6.2).

RX J1744.7-2713
RX J1744.7-2713 is classified as a BeXRB (Israel et al. 1997) and is composed of a B0.5 III-Ve star (Motch et al. 1997;Steele et al. 1999;Lopes de Oliveira et al. 2006) and a WD.However, the origin of X-ray emission is still uncertain and debated in the literature (Lopes de Oliveira et al. 2006).It is also known as a γ Cas analog, due to the similarities in their X-ray properties (e.g., Shrader et al. 2015).The first high-contrast images of this HMXB were presented in Prasow-Émond et al. (2022), in which more comprehensive information on the system can be found.

RX J2030.5+4751
RX J2030.5+4751, also known as BD+47 3129 and SAO 49725, was first discovered with ROSAT (Motch et al. 1997).It is identified as a BeXRB and γ Cas analog.It has a B0.5 III-Ve spectral type donor star and a maximum X-ray luminosity of ∼10 33 erg s −1 (Liu et al. 2006;Raguzova 2007).XMM-Newton/EPIC observations have revealed that RX J2030.5+4751 has a hard X-ray spectrum, suggesting the presence of a dense, large, and stable circumstellar environment surrounding the system (Lopes de Oliveira et al. 2006).On July 20, 2016, RX J2030.5+4751underwent type I (i.e., smaller and repetitive) bursts, reaching its maximum luminosity (progressive weakening since; Steele 2016a,b,c).Regarding the long-term variability, Servillat et al. (2012) reported a significant non-periodic variability of approximately one magnitude in the light curve of RX J2030.5+4751 over ∼100 years, likely caused by changes in the properties of the decretion disk.

γ Cas
As discussed in Section 4.1, observations of γ Cas in the L ′ -band were initially made on September 08, 2017, and subsequently revisited almost three years later on July 11, 2020.Figure 2 presents the L ′ -band highcontrast images of γ Cas for both epochs.A bright source, labeled B, was detected with a SNR ≫ 5.A much fainter source, labeled C, was also detected, with a SNR ∼ 3.In Section 6.2, we undertake a proper motion analysis to determine whether these sources are more likely to be bound CBCs or background stars.
Table 6 in the Appendix presents several physical properties of the detected sources, including the angular separation in mas, the position angle in degrees, the apparent magnitude in L ′ -band and the mass estimated from evolutionary models.
In the next sections, we analyze the nature of the detected sources and discuss the implications of the results.
6. ON THE NATURE OF THE DETECTED SOURCES

Background Contamination
We used TRILEGAL (Girardi et al. 2005) to discriminate background/foreground stars from gravitationallybound CBCs.TRILEGAL is a 3D model employed to simulate the photometric properties of star fields within the Galaxy (e.g., Chun et al. 2015;Dietrich & Ginski 2018;Jones et al. 2021;Williams et al. 2021).We compiled a list of predicted sources within a 1 × 1 arcmin 2 region surrounding the X-ray binary using its RA/Dec coordinates.Subsequently, we applied a geometric rescaling so that the TRILEGAL foV matches the foV of the high-contrast images (10.24 × 10.24 arcsec 2 ) for statistical consistency.We calculated the cumulative distribution of the expected number of sources in the foV below an apparent magnitude (m L ′ ), for the distribution of apparent magnitudes (L ′ ).It is denoted as n foV (L ′ ≤ m L ′ ).
Fig. 4 presents n foV (L ′ ≤ m L ′ ) as a function of m L ′ for both TRILEGAL and the detected sources in the high-contrast images.It includes all X-ray binaries with detected sources, except γ Cas (for which the confirmation of the nature of the CBCs relies primarily on proper motion analysis) and RX J1744.7−2713(presented in Prasow-Émond et al. 2022).In Table 6, we listed n foV (L ′ ≤ m L ′ source ) for each detected source.This represents the expected number of sources -based on TRILEGAL simulations -with apparent magnitudes (m L ′ ) below the magnitude of the corresponding source (m L ′ source ) while accounting for errors.If the calculated value, including the upper limit, was lower than the total number of sources detected below m L ′ source , we reject the hypothesis of background/foreground contamination for that particular source.
Many sources were not predicted by TRILEGAL.Specifically, all detected sources in 4U 1700−377, RX J2030.5+4751,Cygnus X-1, X Per and 1H2202+501 were not expected from the model given their magnitudes and would thus suggest that they are bound to the X-ray binary.The foVs of IGR J17544−2619 and SAX J1818.6−1703 were expected to be more populated than what we detected, suggesting that the detected sources are likely background or foreground contaminants.The disparity between the predicted and observed number of sources might result from the elimination of stationary sources during the ADI process or from an insufficient S/N for detection.In the case of IGR J18483−0311, sources B, C, and F might be CBC candidates, but their status remains uncertain in the current stage of our study.
We extended our analysis by calculating the probability of chance alignment for each detected source.This probability represents the likelihood that these sources are not associated with the X-ray binary system based on the angular separation and the density on the sky of unrelated objects.This method assumes that the distribution of unrelated sources across the area follows a Poisson distribution.Note that this method usually relies on sky surveys such as the 2MASS Point Source Catalog (e.g., Correia et al. 2006;Lafrenière et al. 2008Lafrenière et al. , 2014;;Prasow-Émond et al. 2022) to establish the distribution of unrelated objects.However, in this work, we used TRILEGAL as an alternative due to the absence of available K s -band observations.
To calculate the aforementioned probability, we first divided the cumulative distribution of the number of sources n TRILEGAL (L ′ ≤ m L ′ ) by the area from which the sources were retrieved (between 6 × 6 arcmin 2 and 15 × 15 arcmin 2 depending on the location of the Xray binary).This division enabled us to derive a surface density denoted as Σ.Using the angular separation Θ in arcsec, the probability of a source being drawn from the TRILEGAL distribution -thus indicating its lack of association with the central X-ray binary -is given by: In Table 6, we listed 1 − P unrelated (Σ, Θ) as percentages.Many sources have high probabilities (> 85%) of being associated with the central X-ray binary.However, some sources have lower probabilities (between 65% and 75%) which are not as statistically significant, but we nevertheless identified them as candidate CBCs given the early stage of the study.Sources with probabilities below 60% (not statistically significant) were rejected as candidate CBCs.

Proper Motion Analysis: γ Cas
Conducting a proper motion analysis is among the most robust methods for confirming the gravitational association of a source with the system (e.g., Bohn et al. 2021).Given the proximity of γ Cas and its proper motions (see Section 4.1), it was possible to conduct a statistically significant proper motion analysis between the two observed epochs (September 08, 2017, andJuly 11, 2020).Figure 5 presents the relative separations be-tween sources B, C, and γ Cas in RA and Dec, alongside the expected position of a stationary background star.The figure also displays the angular separation and position angle over time, along with the expected tracks for both comoving and stationary background objects.Source B's trajectory implies an underlying motion that necessitates additional epochs of observation for validation.In 2020, the angular separation data point aligns with the comoving track, but the position angle data point deviates by ∼ 3σ from the same track.This trajectory suggests that source B is more likely to be bound to γ Cas rather than an unrelated background or foreground object.Nonetheless, its motion suggests potential scenarios such as orbital motion, ejection from the system, or the presence of a non-stationary background or foreground object.Using the mass of γ Cas (∼ 13 M ⊙ ; Nemravová et al. 2012) and the radial separation of source B, we calculated the as v esc = 2GM/ρ ≈ 7680 m/s.Additionally, the projected velocity between 2017 and 2020 was determined, resulting in the vector v proj ≈ (1670 ρ − 17964 θ) m/s.Since the norm of the velocity vector is greater than the escape velocity, source B appears to be physically associated with γ Cas, but is not bound (as also suggested in Hutter et al. 2021).A comprehensive characterization of this motion using high-contrast imaging necessitates further epochs of observation.
Source C consistently follows the motion track of an unrelated object across all three plots.Thus, we excluded it from our list of candidate CBCs with a ≫ 3σ confidence level.

Frequency of Multiple Systems and Companion Frequency
Two key concepts are commonly defined in the stellar multiplicity literature (e.g., Duchêne & Kraus 2013): the frequency of multiple systems (MF) and the companion frequency (CF; i.e., the average number of companions per target).While a proper motion analysis is required to confirm most of the sources, we calculated a first estimation of MF and CF for the observed X-ray binaries in our sample.Among the total of 14 observed Xray binaries, we have identified candidate CBCs in eight systems: 4U 1700−377, RX J2030.5+4751,Cyg X-1, X Per, 1H2202+501, γ Cas, RX J1744.7−2713, and IGR J18483−0311.Thus, based on these numbers and at this stage of the study, MF for triple or higher-tier systems would be 8/14 ≈ 0.6 (∼ 60%).For CF, it is important to note that X-ray binaries inherently possess one companion, the donor star, which is included in our calculation of CF.Based on the status of the detected sources listed in Table 6, the calculated average number of com-  panions per compact object is 2.1 ± 1.1 (210 ± 110%).This means that, on average, the compact object has two companions (the donor star and one additional companion).However, this value reduces to 1.8±0.9(180±90%) when considering only the sources that are most likely to be gravitationally bound (1 − P unrelated (Σ, Θ) > 85%).These values are subject to change as new observations become available and further analyses are conducted.The discovery of candidate CBCs would imply that X-ray binaries can still be produced by multiple-star systems rather than exclusively binary systems.The total mass (i.e., compact object and donor star) of the HMXBs in our sample all exceeds ∼ 10 M ⊙ , with some reaching up to around 60 M ⊙ , placing them on the higher end of the mass spectrum.Stellar multiplicity is believed to be common in high-mass star systems (e.g., Chini et al. 2012;Duchêne & Kraus 2013), and high-order multiplicity is thought to increase with primary mass (e.g., Peter et al. 2012).However, surveys for high-mass stars remain incomplete.
For high-mass stars (≳ 16 M ⊙ ), MF and CF are estimated to be ≥ 80% and 130 ± 20%, respectively (Chini et al. 2012;Duchêne & Kraus 2013).The first estimation of MF for our sample (∼ 60%; calculated in Section 6.3) falls below this percentage.However, in this study, MF is constrained by the range of projected separations (up to ∼ 12, 000 au).This implies that increasing this limit could potentially lead to the discovery of more companions and hence increasing the estimation of MF.As for CF (210 ± 110%; calculated in Section 6.3), it currently exceeds the estimate obtained from the literature (130 ± 20%; Chini et al. 2012).Further observations will likely lead to the rejection of candidate CBCs we have detected, which would lower the sample's CF (along with the associated uncertainty range).In both cases, our preliminary estimations (MF and CF) seem to be broadly in line with current estimates in the literature.
For solar-type stars, the frequency N (n) of multiplicity n follows a geometric distribution N (n) ∼ β −n (up to n = 7 with β = 2.3 or 3.4; Eggleton & Tokovinin 2008).Such a relation has not been derived for massive stars.Thus, the results from this study could be used to derive one, allowing us to predict the frequency of multiplicity in X-ray binaries and high-mass systems.The current sample size may be insufficient to infer a statistically significant relationship, but these results can still contribute valuable data to future surveys of high-mass systems.
This pilot study, in addition to showing evidence for potential additional components in X-ray binaries, could contribute to stellar multiplicity surveys for massive stars.It could also enable us to probe stellar multiplicity at low mass ratios (below ∼ 0.1).

Stability in Wide Orbit
The projected separations within the scope of this study (from ∼ 350 to ∼ 12,000 au) would suggest that CBCs orbit at very large distances from the central Xray binary.Note that CBCs located closer to the X-ray binary within the 2-3 kpc distance range cannot currently be detected through direct imaging.To assess whether potential CBCs could be gravitationally bound to these systems, we calculated the binding energy E bind for every source likely to be a CBC.Assuming circular orbits, E bind is estimated using the following equation (e.g., Naud et al. 2014;Prasow-Émond et al. 2022): where G is the gravitational constant, M XRB is the total mass of the central X-ray binary, M comp is the mass of the CBC, r is the projected separation between the CBC and the X-ray binary and 1.27 is the average projection factor between r and the semi-major axis assuming a random viewing angle (e.g., Brandeker et al. 2006).
The binding energies for each candidate CBC range from ∼ −2 × 10 42 erg to ∼ −3 × 10 44 erg.A binding energy of ∼ −10 41 erg was obtained for the comoving exoplanet GU Psc b around a M3 spectral type star (∼ 0.46 M ⊙ ), with a mass of ∼ 9 − 13 M Jup and located at a distance of ∼ 2000 au (Naud et al. 2014).All candidate CBCs have binding energies largely exceeding this currently known lower limit.This suggests that these sources, if confirmed as CBCs, would fall within the gravitational binding range of the X-ray binary.
We must also consider dynamic stability for systems containing more than one candidate CBC.N -body simulations have shown that a configuration of two CBCs at the same projected separation from the central system can lead to dynamic instability (e.g., Kiseleva et al. 1994).Thus, in the case of 4U 1700−377, where B and C have similar projected separations, configurations B+D or C+D are more likely than B+C+D.As for RX J2030.5+4751 and IGR J18483−0311, the candidate CBCs have distant projected separations, suggesting that configurations B+C and B+C+D, respectively, are plausible.

Companion Formation and Capture Scenarios
If our findings are confirmed, we hypothesize that CBCs orbiting X-ray binaries could originate through two main mechanisms: (1) formation within the same environment as the central X-ray binary or (2) capture by the system.On the one hand, as detailed in Prasow-Émond et al. (2022), there are three scenarios in which CBCs could potentially form in the direct environment of X-ray binaries, and these scenarios may unfold at different times.Firstly, CBCs could have formed simultaneously with the initial stars that subsequently evolved to form the present X-ray binary, resulting from the direct fragmentation of a collapsing prestellar core (e.g., Bate 2012).Note that this scenario is unlikely because we would expect companions with a mass similar to that of the main X-ray binary.Secondly, their formation might have occurred within the circumbinary disk surrounding the initial binary system, prior to the supernova explosion of the now compact object (e.g., Kratter et al. 2010).Thirdly, CBCs could have formed within the supernova fallback disk arising from the explosion (e.g., Wolszczan & Frail 1992).Note that fallback disks do not contain enough mass to enable star formation, but they do possess enough mass for the formation of substellar objects.Additional observations are needed to detect and characterize such disks (see Section 7.4).
On the other hand, stellar and substellar CBCs could be gravitationally captured by the X-ray binary.This is analogous to the case of the PSR B1620-26 system, where a giant exoplanet is thought to have been captured by a binary system containing a neutron star and a white dwarf (e.g., Sigurdsson et al. 2003).Given that the HMXBs in our sample are massive (> 10 M ⊙ ) and located in the galactic plane (see Fig. 1), such events would not be unlikely.

Follow-Up and Complementary Observations
The findings of our study primarily consist of intermediate results, and additional observations are necessary to conduct further analyses and confirm that the CBCs are bound to the system.In this section, we provide a list of recommendations for follow-up and complementary observations.First, we recommend re-observing the systems in the same band (L ′ ) using the same instrument (NIRC2) at one or multiple additional epochs.This will allow us to conduct proper motion analysis (see Section 6.2) for all the candidate CBCs identified in this study.Multiepoch observations will also enable us to characterize the orbital motion of these CBCs.Since the HMXBs presented in this study are located within 2-3 kpc, the time interval between epochs ranges from a few months to a few years.Table 4 provides an estimate of the recom-mended year of re-observation for each system, ensuring a statistically significant (> 3σ) proper motion analysis.
Second, we recommend observing the remaining 5 sources, such as Scorpius X-1, 1A0620-00, and V404 Cyg, to complete the sample of all X-ray binaries within 2-3 kpc accessible with Keck/NIRC2.This would also enable us to incorporate LMXBs into the analysis and discussion.
Third, observations in other bands (e.g., K s ) would allow us to construct color-magnitude diagrams and employ evolutionary models to estimate the physical properties of the CBCs with greater constraints (e.g., Prasow-Émond et al. 2022).Similarly, obtaining the near-infrared spectrum of the candidate CBCs would enable us to characterize the nature of the source and extract additional physical properties.Finally, submillimeter observations of the continuum emission would allow us to detect and characterize potential AU-scale circumstellar or protoplanetary disks composed of dust and hot gas (e.g., Coleiro et al. 2013;Iyer & Paul 2017;Waisberg et al. 2019).
Contrast curves for all systems with CBCs can be found in the Appendix to assess the limit that has been reached during these observations.Table 4. Recommended year for follow-up observations for every X-ray binary with at least one candidate CBC -calculated using the distance, proper motions, and astrometric errors on the 2020 observations -enabling a proper motion analysis with at a > 3σ significance level.In this study, we presented the first L ′ -band highcontrast images of nearby, high-mass X-ray binaries using NIRC2 and the vortex coronagraph on the W. M. Keck Observatory.A total of 14 systems were observed from a sample of 19 X-ray binaries within ∼ 2-3 kpc and amenable for direct imaging.One or several sources with a SNR > 5 were found in 8 of the observed Xray binaries.To discern the nature of these sourceswhether unrelated objects or candidate CBCs -we employed Galactic population models for all systems, and proper motion analysis for γ Cas.We find that, if confirmed, these results would imply the presence of stellar and sub-stellar CBCs in the direct environment of Xray binaries (∼ 350 − 12, 000 au), which opens up the discussion on the binary nature of these systems.As a pilot study, this work presents a catalog of photometric and astrometric parameters for the first epochs of observations.Follow-up observations or additional characterization (e.g., infrared spectrum) will enable us to conduct proper motion analyses to discriminate more robustly background/foreground sources from comoving, gravitationally-bound CBCs.

APPENDIX
This Appendix presents additional data related to the main text.Figure 6 presents the 5σ contrast curves for all X-ray binaries with CBCs.Table 5 presents additional relevant physical properties for the 14 observed X-ray binaries.Table 6 presents the properties of the detected sources in the high-contrast image, including the status (background/foreground source or candidate CBC), optimization parameters, physical properties (astrometric and photometric parameters), and estimated mass.); see Section 6.1); 6.One minus the probability of being unrelated with the central X-ray binary (1 − P unrelated (Σ, Θ); see Section 6.1); 7. The radial separation from the X-ray binary in mas; 8.The position angle in the image in degrees; 9.The apparent magnitude in the L ′ -band (no extinction correction was applied); 10.The mass of the source, in M⊙ or MJup, estimated from evolutionary models (MIST or COND/DUSTY); 11.The projected separation from the X-ray binary in au.

Figure 2 .
Figure 2. Keck/NIRC2 L ′ -band high-contrast images of γ Cas acquired on 2017 Sept 08 (up) and 2020 July 11 (bottom), treated and reduced using a PCA annular ADI algorithm (using VIP; Gomez Gonzalez et al. 2017).The sources detected with SNR > 3 are labeled B and C. The white X symbol denotes the approximate position of γ Cas, masked by the coronagraph.The bottom-left insets show zoomedin high-contrast images with a focus on an annular region (obtained using the PCA annulus algorithm in VIP), highlighting the presence of the fainter source (labeled as C).

Figure 3 .
Figure 3. Keck/NIRC2 L ′ -band high-contrast images of all observed X-ray binaries -except RX J1744.7−2713(see Prasow-Émond et al. 2022) and γ Cas (see Fig 2) -acquired on September 8, 2017, January 3, 2018, and July 11-13, 2020.Images were treated and reduced using a PCA annular ADI algorithm (using VIP; Gomez Gonzalez et al. 2017).The sources detected with SNR > 5 when computing the signal-to-noise map are labeled.The white X symbol denotes the position of the X-ray binary masked by the coronagraph.North points upwards and East points to the left, as in Fig. 2.

Figure 5 .
Figure 5. Left: Relative separations between source B (top row), C (bottom row), and γ Cas in right ascension (α) and declination (δ).The first epoch astrometric point is plotted in blue (September 08, 2017) and the second epoch astrometric point is plotted in red (July 11, 2020).The expected position for a stationary background object is plotted in yellow, along with its proper motion track.Middle: Separation from γ Cas in mas as a function of time.A background object with zero proper motion would follow the orange track, while a bounded and comoving CBC would lie within the blue zone.Right: Same as middle, but for the position angle in degrees.

(
Physical properties of the detected sources in the high-contrast images.The columns are: 1.The target; 2. The label of the detected source (S/N > 5); 3. The current status of the source, either background/foreground object (bkg) or candidate CBC (cc); 4. The optimal number of principal components used for the fitting (ncomp; 5.The expected number of sources in the foV below the apparent magnitude (nfoV (L ′ ≤ m L ′

Table 3 .
Gaia DR3 ID for each target.

Table 5 .
Additional relevant physical properties for the 14 observed X-ray binaries.The columns are: 1. Name of the target; 2. Distance in kpc (see Section 3.3.1);3. Proper motion in right ascension in mas yr −1