© The Authors 2025

1 Introduction

The number of directly imaged exoplanets has grown steadily over the past decade thanks to imaging surveys with ground-based facilities. About 40% of the imaged exoplanets are located at separations greater than 100 au from their host stars. For instance, AB Pic b is a 13.5 ± 0.5 M_J companion of a K0 star at a projected separation of 275 au (Chauvin et al. 2005); COCONUTS-2b is a 6.4 ± 2 M_J cool companion of a field M dwarf at a projected separation of 6471 au (Zhang et al. 2021); 2MASS J1155-7919 b is a 10 M_J co-moving companion of a M3 dwarf at a projected separation of 582 au (Dickson-Vandervelde et al. 2020); HD 203030b is an 8–15 M_J companion at a projected separation of 487 au orbiting around a G8 star (Miles-Páez et al. 2017); 2M 2236+4751 b is an 11–14 M_J companion at a projected separation of 230 au orbiting around a K7 star (Bowler et al. 2017); GU Psc b is a 9–13 M_J companion at a projected separation of 2000 au orbiting around a young M3 star (Naud et al. 2014). The discovery of these exoplanets challenges known formation theories of gas giant exoplanets.

It is unclear how these planets are formed based on our current understanding of how and where gas giant exoplanets are expected to form. There are two widely accepted scenarios for the formation of gas giant exoplanets: core accretion (e.g. Pollack et al. 1996) and gravitational instability (e.g. Boss 1997). Core accretion is based on the formation of a rocky embryo close to the ice line within a circumstellar disc. Once the embryo reaches a critical mass, runaway accretion of gas continues until the zone around the planet is cleared. However, core accretion is problematic for wide-separation planets. It takes more than 10⁷ yr for core formation by runaway planetesimal accretion beyond 5 au (Goldreich et al. 2004), which is longer than the gas dissipation time of the disc. Therefore, planet migration is required for wide-separation planets if they are formed by core accretion (Dodson-Robinson et al. 2009). Instead of core accretion by planetesimals, core growth by pebble accretion can form planets beyond 50 au much faster than via traditional core accretion (Lambrechts & Johansen 2012). However, this process still requires ejection or scattering processes to explain the observed extremely wide-separation planets. Gravitational instability is driven by the initial fragmentation of the protostellar gas cloud. This formation mechanism has a faster timescale than that of core accretion and forms planets typically at 20 to 200 au from the star, supporting in situ planet formation. However, for planets at even wider separations, it also requires planet migration. In addition, gravitational capture of free-floating planets is another possible formation mechanism for wide-separation planetary systems (Perets & Kouwenhoven 2012).

There are several direct imaging surveys to search for exoplanets, such as the Spectro-Polarimetric High-contrast Exoplanet REsearch (SPHERE, Beuzit et al. 2019) infrared survey for exoplanets (SHINE, Chauvin et al. 2017) and the Gemini Planet Imager (GPI) Exoplanet Survey (GPIES, Macintosh et al. 2015). Among these surveys, the Young Suns Exoplanet Survey has a relatively high planet detection yield (YSES, Bohn et al. 2019, 2020a,b, 2021). It is a direct-imaging snapshot survey of exoplanets targeting a homogeneous sample of 70 solar-mass pre-main-sequence stars in the Lower Centaurus Crux (LCC) subgroup of the Scorpius-Centaurus association (Sco-Cen, de Zeeuw et al. 1999; Pecaut & Mamajek 2016) with SPHERE. It has detected two wide-separated exoplanets orbiting around a K3 star so far: YSES 1bc, a 14 ± 3 M_J companion at a projected separation of 160 au and a 6 ± 1 M_J companion at a projected separation of 320 au (Bohn et al. 2020a,b).

In this work, we present the detection of a planetary-mass candidate companion of 2MASS J10065573-6352086 (hereafter: 2M1006) as part of YSES. The primary star is referred to as star A and the fainter star is referred to as star B. The observations of the system and data reduction are detailed in Sect. 2. In Sect. 3, we present the analysis of the central stars, a stellar system consisting of a young Sun analogue and an M dwarf. We find tension among all the data and speculate that there is at least one more stellar component in this system. The photometric and astrometric analysis of the planetary-mass candidate are presented in Sect. 4. We discuss the probability of the candidate being a free-floating object and co-moving planet in Sect. 5 and summarise our work in Sect. 6.

2 Observations and data reduction

We obtained six epochs of observations of 2M1006 between 2018 and 2024 with SPHERE and MagAO-X, including imaging and integral field unit (IFU) spectra. We summarise the imaging observation setup and conditions in Table 1.

2.1 SPHERE observations and data reduction

The SPHERE instrument is an extreme adaptive optics system and coronagraphic facility mounted on the Unit Telescope 3 of the ESO Very Large Telescope (VLT) (Beuzit et al. 2019). 2M1006 was observed with SPHERE the infrared dual-band imager and spectrograph (IRDIS, Dohlen et al. 2008) in five epochs on 2018-11-15, 2021-12-02, 2021-12-26, 2023-01-28, and 2023-03-03. Observations in the first, second, fourth, and fifth epoch were obtained in the classical imaging mode (CI) with the broad-band H and K_S filters, while observations in the third epoch were obtained in the dual-band imaging mode (DBI, Vigan et al. 2010) using the K12 and J23 narrow-band filter combinations. Only observations in the third epoch used the field-tracking mode, while observations in other epochs used the pupil-tracking mode. All the observations used a coronagraph (N_ALC_YJH_S) to mask the central star.

The central star was resolved to be composed of two stars in the first epoch. To measure the precise location of the primary star and fainter star, we performed a point spread function (PSF) fitting to the binary in the flux image where the central star is not masked. We built a PSF template from other stars in YSES in the broad H and K_S bands. We then fit two PSFs simultaneously to the binary using the NAUTILUS package (Lange 2023), which implements importance sampling and efficient space exploration using neural networks. The two stars were well resolved in the first epoch, marginally resolved in the second epoch, and moderately resolved in the fifth epoch. To better fit the location, we fixed the flux ratio obtained from the first epoch for the PSF fitting in the second and fifth epochs. For the third epoch in J2 and J3 we scaled the PSF template in the H band by the corresponding wavelength and performed the fitting. There are multiple non-coronagraphic images (flux images) of the central stars in each epoch. We fit the binary in each individual frame and took the mean and standard deviation as the fitting results and uncertainty of each epoch. The initial uncertainties of the position and the flux ratio are on the order of 10⁻² pixels and 10⁻³. We performed an astrometric calibration following the latest version of SPHERE User Manual¹ and Maire et al. (2021). To correct the systematic error between the PSF template and the observed images, we empirically included an additional 1 mas in RA and Dec in the final position uncertainty and the root mean square (rms) of the residuals in the final ratio uncertainty. This uncertainty agrees with the typical binary measurements in Bonavita et al. (2022), which do not use any PSF fitting. Table 2 lists the error budget for the binary fitting. The final uncertainty is calculated by the propagation of error.

In science images where the central star is blocked by the coronagraph, the binarity of the central source challenges image alignment between the flux images and science images. For coronagraphic observations of SPHERE/IRDIS, there are four waffle spots in the centred image, which are duplications of the PSF of the central star dispersed with wavelength (Langlois et al. 2014; Beuzit et al. 2019). They are used to determine the precise location of the central star behind the coronagraph by the intersection of the two pairs of spots. The positions of the four spots were fitted by a 2D-Gaussian model, which is effective when the central star is a single star. However, when the central star is a blend of two sources, it is problematic to use a single 2D-Gaussian fitting because the spots are dispersed images of the two sources. To measure the precise location of the brighter star in the coronagraphic images, we developed two methods to tackle this problem. The first one is to fit two 2D-Gaussian models to the spot simultaneously. The separation and flux ratio of the binary obtained from the PSF fitting of the flux image were taken as initial guesses for the two 2D-Gaussian fittings. After measuring the positions of the sources in the four waffle spots, we measured the positions of the brighter star by the intersection of the two pairs of the spots corresponding to the primary star. The second method is to perform cross-correlation between a simulated waffle image of the blended source and the observed image. We built an intensity profile of waffle spots of a single star and convolved it with the flux image. Then we cross-correlated the simulated waffle image with the observed waffle image. The peak position of the cross-correlation function is the offset between the two images.

The telescope was supposed to be pointed at the primary star. There are two sets of flux images taken before and in the middle of the science observations in the fifth epoch. We centred on the primary star and derotated the flux images to the orientation, where north is up and east is left. The fainter star is at the same position in the derotated flux images, verifying that the telescope was pointed at the primary star. After locating the primary star in the coronagraphic images, we derotated the coronagraphic images back to the north centred on the primary star. Then we stacked the derotated images by averaging. We applied angular differential imaging (ADI, Marois et al. 2006) and a principal component analysis (PCA, Soummer et al. 2012; Amara & Quanz 2012) to the de-rotated images. No objects were detected at separations smaller than $1_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}6$. Because the candidate is at a wide separation, we took the stacked image as the final reduced image.

2.2 MagAO-X observations and data reduction

The six-epoch observations were taken with the MagAO-X instrument (Males et al. 2022). MagAO-X is an extreme adaptive optics instrument for the Magellan Clay telescope at Las Campanas Observatory. MagAO-X has been specifically designed for direct imaging at visible wavelengths (Close et al. 2018). The system was observed on 2024-03-23 in the z′ filter with the EMCCD camera camscil on MagAO-X. The observations were made in direct imaging mode to perform follow-up astrometric measurements of the tight binary system. One minute of further observations were taken in pupil tracking mode. No significant field rotation occurred during this time and therefore the images were stacked before any processing was done. After stacking, basic detector calibration was performed by subtracting a dark frame. The data were then centred on the primary star. It was not possible to create a PSF model directly from the data due to the small separation of the two objects. YSES 2 was observed directly before the observations of this system. Those observations were used to create a reference PSF because of the targets’ similarity in elevation, seeing conditions and brightness during the MagAO-X observations. The YSES 2’s PSF was centred to create the reference PSF. Both components of 2M1006 were fitted with the reference PSF. A true north correction of $1_{\cdot }^{\circ }57$ ± $0_{\cdot }^{\circ }2$ (Long et al. 2025) was applied.

2.3 MagAO-X/VIS-X observations and data reduction

Spectroscopic data were taken with the MagAO-X system, which has an integral-field unit, the Visible Integral-field Spectrograph eXtreme (VIS-X; Haffert et al. 2022). VIS-X is a microlens-based IFU and it has been designed for high-spectral resolution (R = 13 500, Δv = 20 km s⁻¹) observations of the Hα emission line. The spectral range is limited to about 5 nm around the Hα line due to the microlens design. 2M1006 was observed with VIS-X to measure the radial velocity difference between the two stars. The VIS-X observations were taken on 2024-04-02. The total on-target time is 40 minutes, consisting of eight exposures of 5 minutes each. The seeing condition varied between $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}6$ and $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}75$.

The VIS-X data reduction follows a typical IFU data reduction pipeline (Haffert et al. in prep). First, the detector calibration steps are applied by subtracting a dark frame. Internal source calibrations were used to identify the spectral trace model for each spaxel. This model can then be used to extract the 1D spectra for all the spaxels. The shifts in the wavelength solution between spaxels were calibrated using standard star observations, which was the star Sirius in this case. The IFU calibration is completed after these steps. The second part of the data reduction follows a similar approach to the imaging data. The data cubes were derotated and a true north correction was applied. The spectra extraction is described in Sect. 3.4.

3 Stellar properties

2M1006 was identified to be a candidate member of the LCC subregion of Sco-Cen (Preibisch & Mamajek 2008; Pecaut & Mamajek 2016). Gaia DR3 reports a parallax of ϖ = 7.3241 ± 0.0689 mas (Gaia Collaboration 2021), which translates to an estimated distance of d = ${136.15}_{-1.58}^{+1.66}$ pc (Bailer-Jones et al. 2021). The first spectroscopic study of the star came from Torres et al. (2006) who reported it as a K0V(e) star with strong Li absorption (EW(Li I 6707) = 350 mÅ), EW(Hα) = 0.0 Å, fast projected rotation (v sin i = 77.0 ± 7.7 km s⁻¹), and a radial velocity of 14.0 km s⁻¹. Kiraga (2012) reported time-series photometry from the ASAS-South survey and reported variability with a period P = 0.7271 d. TESS observed it in Sectors 10, 11, 36, 37, 38, 63, and 64. It shows a distinct rotation modulation of 0.73 d in all these sectors. Fetherolf et al. (2023) reported a detection of a rotation period of P_rot = 0.725893 ± 0.00557 day in the TESS Sector 10 data in agreement with the variability period reported by the ASAS-South survey.

2M1006 was resolved as two stars in SPHERE observations in Fig. 1. Given their small separation, they most likely form a bound binary pair. After considering multiple lines of evidence (detailed in the following sections), there is also significant evidence for at least one more stellar companion in the system.

Fig. 1 Full image of this system of the first epoch in the H band. Left panel: central stars. The fainter component is in the south-east direction of the primary star on 2018-11-15. Middle panel: stacked coronagraphic image centred on the primary star. There are ten sources in the field of view in total. The planetary-mass candidate is the faint source between the two bright background stars in the south-west. Right panel: zoomed-in image of the candidate. All images are derotated to the orientation, where north is up and east is left.

3.1 Acceleration

A companion can affect the proper motion of a star through gravitational reflex motion. Therefore, a proper motion anomaly (PMA) is the difference between the proper motion measured in the long term and the proper motion measured in the short term is an indicative sign of a perturbing companion (Kervella et al. 2019). To search for any PMA of 2M1006, we checked the Hipparcos (Perryman et al. 1997), Tycho-2 (Høg et al. 2000), and UCAC5 catalogues (Zacharias et al. 2017). This star is not in the Hipparcos catalogue, but is in the Tycho-2 and UCAC5. The difference in proper motion between Tycho-2 and Gaia is within the large uncertainty of Tycho-2. The proper motion provided in the UCAC5 catalogue is a long-term proper motion, calculated as the difference in position between the mean UCAC epoch and the mean Gaia DR1 epoch, divided by the difference in time between those mean epochs. The long-term proper motion can be assumed to follow the barycentre of the system and to be relatively independent of perturbations from a companion (Kervella et al. 2019). The proper motion provided in the Gaia catalogue is a short-term proper motion, over the epochs covered in the Gaia DR chosen. The difference between the short-term proper motion and the long-term proper motion is the proper motion anomaly (PMA) and is an indicator of the degree to which an unseen companion causes the measured short-term proper motion to deviate from the barycentric proper motion. Thus, we constructed the Gaia-UCAC5 PMA as $$PMA={\mu }_{Gaia}-{\mu }_{UCAC5-Gaia},$$(1)

where μ_Gaia is the reported short-term Gaia proper motion and μ_UCAC5–Gaia is the reported long-term proper motion between the UCAC and Gaia DR1 epoch.

The star’s μ_α cos δ is −21.5 ± 1.2 mas yr⁻¹ and μ_δ is 8.1 ± 1.2 mas yr⁻¹ between UCAC5 and Gaia DR1. Its μ_α cos (δ) is – 17.606 ± 0.086 mas yr⁻¹ and μ_δ is 4.892 ± 0.079 mas yr⁻¹ in Gaia DR3. The proper motion difference is larger than 1σ, suggesting that a fainter companion is perturbing its on-sky motion. Moreover, although Gaia DR3 identifies 2M1006 as a single star, its renormalised unit weight error (RUWE) is 6.012, indicating a poor astrometric fit for a single-star model, which for single stars have values close to unity. Additionally, two flare events in Tess Sector 36 were observed in 2021 with an interval of ~14 d, which could be due to the fainter star, a likely M dwarf. There was another strong flare event at the end of the Sector 36 observations, but it was only captured at the very beginning of the flare, which increased the flux by 10%.

3.2 Orbital fitting

Table 3 lists the relative astrometry and flux ratio of the fainter star to the primary star measured in our observations. The flux ratio of the fainter star to the primary star varies between 0.18 and 0.34 from z′ to K_S, indicating a later spectral type for the fainter star. The images of the stars and residuals after fitted star subtraction can be found in Appendix A. The fainter star moves from the south-east to the northeast relative to the primary star from 2018 to 2024.

We used orbitize (Blunt et al. 2020) to find any plausible orbital solutions for the two stars by fitting the positions measured from direct imaging observations. Because the observations in Ks, K12, and J23 were taken very close to the observations in the H band, we only fit the positions in the three H-band epochs and the last z′-band epoch. We set a Gaussian prior for the parallax with 7.3 mas reported in the Gaia catalogue and inflated its uncertainty to 0.5 mas. We set a broad Gaussian prior for the total mass with 1.8 ± 0.8 M_⊙ with an initial guess from the flux ratio. We show the posterior distributions of orbital parameters in Fig. 2 and 100 randomly drawn orbits in Fig. 3. A high inclination of ${110}_{\cdot }^{\circ }4\pm 3_{\cdot }^{\circ }1$ is required to fit the orbit. The fitted total mass is 2.1 ± 0.5 M_⊙. The mass of the primary star is well-constrained to about 1.03 M_⊙ (see Sect. 3.6). This implies that the best-fitting mass of the fainter star is 1.06 M_⊙, resulting in a mass ratio between the fainter star and brighter star of 1. However, a mass around 1 M_⊙ for the fainter star contradicts the mass expected from the measured flux ratio between the two stars. From the measured flux ratio, we expect the fainter star to be considerably less massive than the brighter star.

Furthermore, using GAIAUNLIMITED (Castro-Ginard et al. 2024) which queries Gaia scanning laws and estimates Gaia DR3 detectability of unresolved binary systems, we calculated a RUWE of 1.0–2.4 for the orbits fitted by ORBITIZE depending on the initial phase of the orbits observed by Gaia, much smaller than 6.012. We also tried to fit the orbit of the fainter star with Stan (Wallace et al. 2025) combining RUWE and the direct imaging observations. Stan is developed to fit the mass and orbit of a hidden companion that can cause a high RUWE of a given star by fitting Gaia absolute astrometry including proper motion, parallax and RUWE and possibly combined with image observations using Bayesian inference. However, it is not able to find plausible orbit solutions for the two stars. While we can constrain the orbit very well by including RUWE and positions only between 2018 and 2023 in the orbital fitting, the observed position in 2024 does not agree with the predicted position of 2024, which should show an arc instead of still linear motion. It is therefore also problematic to explain the high RUWE value given these orbital solutions for the two stars. We calculated the PMA from the orbital fit posteriors following the method of Biller et al. (2022), adopting masses for the fainter star of 0.1, 0.3, 0.5, 0.7, and 0.9 M_⊙. Only secondary masses >0.7 M_⊙ can produce PMA on the order of the measured PMA for this system, however, such a high mass for the secondary is in tension with the measured flux ratio between the two stars.

Fig. 2 Posterior distributions of the orbital parameters by fitting the measured relative positions with ORBITIZE.

Fig. 3 100 orbits drawn from the posterior distributions of the binary fitted with ORBITIZE.

3.3 SED fitting

We performed a spectral energy distribution (SED) fitting with SPECIES (Stolker et al. 2020) to the two stars simultaneously using the multi-modal nested sampling algorithm (Feroz et al. 2009). We fit broad-band photometry (unresolved) of 2M1006 from GAIA, 2MASS, and WISE catalogues with the synthetic spectra of BT-Settl-CIFIST (Allard et al. 2013). We used the parallax from Gaia with an inflated error of 0.5 mas. We set uniform priors for the T_eff, log(g), and radii of both stars: 4000 K < T_{eff 1} < 7000 K, 3000 K < T_{eff 2} < 4000 K, 3 dex < log(g₁) < 6 dex, 3 dex < log(g₂) < 6 dex, 0.005 R_⊙< R₁ < 2 R_⊙ and 0.005 R_⊙ < R₂ < 2 R_⊙, and solar metallicity. To combine our direct image observations with broad-band photometry, we set Gaussian priors for the flux ratio of the fainter star to the primary star in SPHERE/IRDIS/B_H, SPHERE/IRDIS/B_ K_S, SPHERE/IRDIS/D_J23, and MagAO-X/camsci1/z′ bands. The optimised SED fitting is presented in Fig. 4. The posterior distributions of the fitted parameters are presented in Appendix B.

We derived T_eff = ${5196}_{-96}^{+162}$ K, R= 1.09 ± 0.05 R_⊙, log(g) = ${3.88}_{-0.50}^{+0.61}$ dex, log(L/L_bol) = −0.10 ± 0.05 dex for the primary star and T_eff = ${3458}_{-91}^{+110}$ K, R = 0.99 ± 0.06 R_⊙, log(g) = ${4.48}_{-0.43}^{+0.59}$ dex, and log(L/L_bol) = −0.90 ± 0.04 dex for the fainter star. With an estimation of T_eff for young stars from Pecaut & Mamajek (2013), the primary star should be a G8–K0 star. The radius of the primary star is constant with its T_eff for a 20 − 30 Myr 1 M_⊙ star from the BT-Settl evolutionary model (Baraffe et al. 2015). However, the radius of the fainter star is too large for an M dwarf with that T_eff. It has to be younger than 4 Myr to have a radius of 0.99 R_⊙ from the evolutionary models, but it should have a radius smaller than ~0.7 R_⊙ if it is in the same system as the primary star with the same age. There is likely some additional luminosity from a third unresolved object and the inflated radius of the fainter star in the SED fitting is accounting for that additional luminosity. The fitted A_V is smaller than 0.13 mag. No infrared excess was detected from the SED fitting.

Fig. 4 Simultaneous SED fitting of the primary star and fainter star to the unresolved broad-band photometric measurements of Gaia, 2MASS, and WISE with SPECIES. The black solid line is the combined best-fitted spectrum of the two stars. The black dashed line is the best-fitted spectrum of the primary star and the black dotted line is the best-fitted spectrum of the fainter star. Spectra drawn from 30 random samples of each star are shown in grey lines. The top panel shows the filter trans-mission profiles and the bottom panel shows the fitting residuals.

3.4 Radial velocity

The primary star is detected at high significance while the fainter star is detected with a lower signal-to-noise ratio (S/N) in the MagAO-X/VIS-X IFU data as shown in the left panel of Fig. 5. The spectrum of the primary star is extracted with an aperture of a radius of 3.5 spaxels centred on the star. For the fainter star, we extracted its spectra with the same aperture but shifted the aperture around the fainter star by 1 spaxel in either x or y direction each time. Starting from the position of the fainter star, x and y change between [−3, 3], respectively, with a total of 49 centring positions. This results in 49 spectra of the fainter star extracted from different places. This step is to include as many pixels containing the signal of the fainter star as possible while minimising the effect of possible bad pixels such as those affected by cosmic rays. Then we applied PCA to the 49 spectra to extract the spectra of the fainter star and also remove the contamination of the primary star. We compare the spectra of the primary star and the first three components of the 49 spectra in the middle panel of Fig. 5. We detected H_α absorption for the primary star, in agreement with the literature. We extracted the spectrum of the fainter star as the first principal component as the fainter star is the dominant source in these regions. We detected H_α emission in the first component, suggesting that the fainter star is likely a magnetically active M dwarf. The second principal component is featureless and likely originates from the background. The third principal component has the same shape as the primary star’s spectrum with the same H_α absorption feature. This verifies that this method successfully distinguishes the spectra of the two stars. We also injected artificial signals of the same brightness and blue shift as the fainter star at the opposite side at the same separation to the brighter star. We successfully retrieved the injected signal using the same method, verifying the reliability of this technique. The injection test is demonstrated in Appendix. C.

The right panel of Fig. 5 presents the smoothed spectra convolved with a Gaussian kernel of a standard deviation = 0.48 Å from one exposure. The H_α absorption and emission due to the two stars are easy to identify. We have three exposures and applied the same method to them. All of them present consistent H_α absorption in the spectrum of the primary star and H_α emission in the first principal component of the fainter star’s spectra. We fitted a Gaussian profile to the spectra and measured a blue shift of 2.11 ± 0.14 Å between the H_α emission line centre and the H_α absorption line centre. This corresponds to a radial velocity difference of 96.5 ± 6.6 km s⁻¹. The escape velocity of the two stars is calculated by: $$v_{esc}=\sqrt{\frac{2G(M_1+M_2)}{r}}.$$(2)

The total mass derived from orbital fitting is 2.1 ± 0.5 M_⊙ and the closest projected separation of the two stars is 32 ± 1 mas, which is 4.2 ± 0.1 au for the parallax of 7.32 ± 0.07 mas. Putting them into Eq. (2), the upper limit of the escape velocity is 29.6 ± 3.6 km s⁻¹. The radial velocity difference between the two stars measured in the H_α line is much larger than their escape velocity.

We are cautious about this comparison because H_α emission profiles of M dwarfs can vary due to chromospheric activity, although it is rare to produce such a large blue shift by chromospheric activity. An analysis of H_α emission lines of 72 M dwarfs finds that 94% of the stars have blue shifts smaller than 50 km s⁻¹ at H_α emission lines (Flasseur et al. 2018). We measured a full width at half maximum (FWHM) of 105 km s⁻¹ of the H_α emission line for the fainter star. However, all stars that have such large FWHM in Flasseur et al. (2018) have a shift of < −50 km s⁻¹ and more than half of them show redshifts. Therefore, while we could not ascertain that the two stars are not gravitationally bound from radial velocity measurements, it makes this two-body system more debatable.

Fig. 5 MagAO-X/IFU data of the binary. Left panel: IFU images of the binary stacked along the wavelength channel. The image is not rotated to place North up. Middle panel: spectrum of star A and the first three principal components of star B’s spectra. The grey area shows the Hα line band. Star A presents Hα absorption. The first component presents Hα emission from star B. The third component (K=3) is the contamination from star A. Right panel: smoothed version of the middle panel by a Gaussian kernel with a standard deviation of 0.48 Å. The grey area shows the Hα line band.

3.5 Possibility of the third companion

There are tensions in the data treating the two stars as a binary system:

1. The dynamical mass measured for the fainter star is too high to match the observed flux ratio between the two stars.

2. The fitted orbit from direct imaging does not agree with the acceleration measurement for this system.

3. The radius of the primary star from the SED fitting is too large if it has the same age as the primary star.

4. The radial velocity difference between the two stars is too large for them to remain gravitationally bound.

We used the equation from Brandner et al. (2000) to estimate the probability that the fainter star is a randomly aligned foreground or background star: $$P(\mathrm{\Theta },m)=1-e^{-\pi \rho (m){\mathrm{\Theta }}^2},$$(3)

where Θ is the angular separation of the object to the central star and ρ(m) is the cumulative surface density of background sources of limiting magnitude, m. We queried all sources within 10 arcmin of the primary star from the 2MASS catalogue (Skrutskie et al. 2006) and counted stars with magnitude down to 1 mag fainter than the faint star, which is 67 stars. Then we estimated the probability of the two stars being unrelated is 5.22 × 10⁻⁷.

Given that the probability of the two stars being unrelated is extremely low, we hypothesise that there is at least one additional unseen companion in the system that may explain the disagreements in our current analysis. Either the primary star or the fainter star can have an additional unresolved stellar component. If it is related to the brighter star, this could explain the PMA, RUWE, and SED and also bias the measured ΔRV. If it is related to the fainter star, this can increase the total mass of the fainter star. For example, if B is an equal mass binary, then all other SED parameters remain the same but the radius of each object decreases from 1.0 R_⊙ to 0.7 R_⊙ and the individual masses are ~0.4 M_⊙. Then the total mass of the fainter component is 0.8 M_⊙. The mass ratio is 0.8 and meets the required total mass of orbital fitting. While the PMA and ΔRV may also be explained by the multiplicity of the fainter star, the RUWE of A cannot be explained by the derived acceleration: with a mass ratio of 1, the derived highest RUWE is only 2.4.

Besides the astrometric perturbations caused by a companion, there are other means of producing a RUWE excess, such as variability (Belokurov et al. 2020). In TESS light curves of the system, the primary is variable at a level of ~2% with a period of 0.73 d. There are two flare events and one of them increases the total brightness by 4% and another very strong eruption event at the end of the light curve, which increases the total brightness by 10%. These flare events are irregular in period and intensity and are consistent with flaring events seen from low mass stars, which we attribute as coming from B. A few percentages of variability can lead to a photocentre shift of several mas, which can cause a large RUWE excess. Therefore, this additional component can either be related to the brighter star or the fainter star.

Fig. 6 Age estimation of 2M1006. Left panel: age estimation by Li depletion using BAFFLES (Stanford-Moore et al. 2020). Right panel: age estimation by combining the luminosity and Teff of the primary star from SED fitting with the isochrones and tracks of BT-Settl evolutionary model (Baraffe et al. 2015).

3.6 Mass and age estimation

The unresolved star is variable with a V-band amplitude of ~0.07 mag and rotation period, P_rot = 0.7271 d (Kiraga 2012). Adopting the colour conversions from Riello et al. (2021) and the Gaia DR3 G and BP − RP colour, we estimated V = 11.015 mag, which is in the middle of the previously published estimates. We used two methods to estimate the age of the star. The first one is to estimate the age by lithium depletion. We fit the equivalent width of lithium absorption (EW(Li)) using Bayesian Ages for Field Lower-Mass Stars (BAFFLES, Stanford-Moore et al. 2020). With EW(Li) = 350 mA (Torres et al. 2006) and B − V = 0.797 mag, BAFFLES estimates a median age of 54.5 Myr with a 95% confidence range of 2.95–396 Myr, shown in the left panel of Fig. 6. The second method is to place the primary star in the luminosity and T_eff diagram with isochrones and tracks from BT-Settl evolutionary model (Baraffe et al. 2015). With log(L/L_⊙) = −0.10 ± 0.05 dex and log T_eff = ${5196}_{-96}^{+162}$ K, we estimated a mass of 1.03 M_⊙ with a 68% confidence range of 0.99–1.08 M_⊙ and an age of 25 Myr with a 68% confidence range of 19–28 Myr, shown in the middle panel of Fig. 6. All the results suggest that the primary star is a young star. We adopted the age results from the second method. The unresolved photometric observations of the system from the literature and adopted parameters of the primary and fainter stars are presented in Table 4.

Fig. 7 Images of the candidate in five epochs of H, J23, K12, and KS bands. The candidate is the faint source between the two bright background stars highlighted by the white circles. It is only visible in the H band with marginal detection in J23 bands.

4 The planetary-mass candidate

Figure 1 shows all the sources and the central binary in the H band. The deepest images are in the H band, with ten sources in the field of view (FoV). Several sources are not detected in the other bands, including the candidate. We detected the candidate in the H band as shown in Fig. 7 with marginal detection in J23 bands.

4.1 Photometric analysis

We performed aperture photometry to measure the flux of the sources in the FoV and PSF fitting to measure the flux of the central two stars. We used the same aperture with a radius equal to the mean FWHM of all sources. For each source, we took the three σ clipped median of an annulus centred on the source, with radii ranging from 10 to 15 pixels as the background level per pixel, and subtracted it from the aperture photometry.

To measure S/N, we chose 14 random apertures of the same size as the aperture on the source adjacent to the source and took their standard deviation as the noise level. Table 5 summarises the contrast of the candidate to the primary star and S/N. We consistently detected the candidate in the H band with S/N > 5 but with only upper detection limits in K12 and K_S bands and marginal detection in J23 bands. We achieved the highest detection significance on 2023-03-03 with S/N = 8.8. Averaging the magnitude of the candidate in three measurements, its apparent magnitude in the H band is 22.04 ± 0.13 mag.

4.2 Astrometric analysis

We fitted a 2D Gaussian model to the candidate to measure its position. At first, we fit a 2D Gaussian model to the candidate with Astropy (Astropy Collaboration 2013, 2018) to obtain the initial guess of the position. Then we cut out a sub-image of 31 pixels × 31 pixels at the position of the candidate. We took three σ clipped median and standard deviation of an annulus with an inner radius of 4 pixels and an outer radius of 12 pixels from the central position of the candidate as the local sky background and uncertainty of each pixel in the sub-image, respectively. We subtracted the fitted model from the sub-image and minimise residuals by minimising the chi-square $\chi ={\sum }_{i=1}^N\frac{(Data-Model{)}^2}{Dataerro{r}^2}$ (N is the number of pixels) using the EMCEE package (Foreman-Mackey et al. 2013).

The uncertainty from the posterior distribution is tiny, and therefore we injected artificial stars using the PSF template and retrieved positions using the above algorithm. We injected 50 stars at the same separation of the object but at the opposite side of the image with one star per injection. We took three σ clipped mean of the deviations between the injected and retrieved positions of the 50 injections as the uncertainty of our position measurement algorithm. We also included uncertainties from the pixel scale (12.25 ± 0.021 mas per pixel), true north angle correction (−1.75 ± 0.08 deg), pupil offset (135.99 ± 0.11 deg), distortion (0.4 mas per 1 arcsec separation), and systematics such as the coronagraph centring uncertainty (we added an arbitrary 1 mas in RA and Dec) in the propagation of the error when calculating the separation between the primary star and the candidate. The error budget is demonstrated in Table 2.

Table 6 presents the separations and position angles of the candidate relative to the primary star from 2018 to 2023. The positions of the candidate relative to the primary star between 2018 and 2023 are shown in the left panel in Fig. 8. They match within 1σ, suggesting that this candidate’s proper motion is comparable to the proper motion of the primary star. Evolving from the position on 2018-11-15, it is more than 3σ away from the position of a static background source on 2023-03-03. The candidate is at a separation of 730 ± 10 au to the primary star assuming they are at the same distance. However, as noted in Section 3.5, the primary star is very likely a triple system; thus, the true barycentre of the system is unknown. Due to the unknown barycentre of the central stars, it is not possible to confirm this candidate as a co-moving source by common proper motion analysis. We discuss the plausible characterisation of this source in Sect. 5.

Fig. 8 Relative astrometry of the candidate in RA and Dec offsets. Left panel: relative astrometry of the candidate to the primary star (mass ratio = 0). Right panel: relative astrometry of the candidate to the barycentre if the fainter star and primary star are an equal mass binary (mass ratio = 1). The yellow and green regions are the candidate companion positions allowed by the escape velocity within 1σ on 2021-12-02 and 2023-03-03, respectively. The blue line is the trajectory of a static background star evolving from 2018-11-15 calculated with the proper motion and parallax reported in Gaia and the two points on the track are the positions on 2021-12-02 and 2023-03-03.

Fig. 9 Colour-magnitude diagram of brown dwarfs. The grey points are field brown dwarfs with high surface gravity. The blue points are young or low surface gravity objects. The brown points are low-mass stars. Several directly imaged exoplanets are labelled with black text. The red arrow shows the possible positions of the candidate assuming it is at the same distance as the central stars. If confirmed, it would be one of the coolest directly imaged exoplanets.

4.3 Colour-magnitude diagram

Assuming the contrast is the same between H_2MASS and H_{S PHERE}, we converted the candidate’s apparent magnitude to absolute magnitude with the primary star’s parallax. We placed the candidate in the colour-magnitude diagram (CMD) in Fig. 9. Assuming the same distance of the primary star (136 pc), the candidate has an absolute magnitude of 16.38 ± 0.13 mag. Due to the unknown magnitude in the J band, we could only derive a blue limit of its 2MASS J − H colour of 0.22, which is redder than most T-type objects of the same assumed H magnitude. We also present the colour of M dwarfs. Based on the current colour range, we cannot determine whether the candidate is a planetary-mass object, brown dwarf, or background M dwarf. If it is a planet, it would be one of the coolest directly imaged planets that reveals a similarity to 51 Eri b (Macintosh et al. 2015; Samland et al. 2017) and AF Lep b (De Rosa et al. 2023; Franson et al. 2023; Mesa et al. 2023). Its J − H colour is between these two confirmed planets. We derived a red limit to its H − K colour of 0.79, redder than most objects of the same H magnitude. Using the ATMO2020 evolutionary model (Phillips et al. 2020), we estimated a mass of 3–5 M_J under the assumption of chemical equilibrium.

5 Discussion

5.1 A co-moving planet or free-floating object

The discrepancies between the orbital fitting to the central two stars with their RUWE value, SED fitting, and radial velocity measurements suggest that there might be an additional low-mass stellar companion in the system. The uncertainty in the resulting barycentre of the system makes it challenging to confirm if the candidate is a co-moving companion by conventional common proper motion analysis. This candidate may be a comoving planet, a free-floating planetary-mass object, a brown dwarf, or a background star. We followed the method used in Nielsen et al. (2017) to quantify the probability ratio of the candidate being a background star to a bound planet of the barycentre of the central stars as a function of the mass ratio between the fainter star and the brighter star. Nielsen et al. (2017) defined three likelihood components of each scenario: the overall likelihood of the object being a background star P(bg) or planetary companion P(pl), the relative likelihood as a function of separation from the star P(ρ|bg) or P(ρ|pl), and the relative probability as a function of projected velocity P(v|ρ_bg) or P(v|ρ_pl): $$\frac{P_{\mathrm{b}\mathrm{g}}}{P_{\mathrm{p}\mathrm{l}}}=\frac{P(\mathrm{b}\mathrm{g})P(\rho \mid \mathrm{b}\mathrm{g})P\left(v\mid {\rho }_{\mathrm{b}\mathrm{g}}\right)}{P(\mathrm{p}\mathrm{l})P(\rho \mid \mathrm{p}\mathrm{l})P\left(v\mid {\rho }_{\mathrm{p}\mathrm{l}}\right)},$$(4)

where P is the probability, ρ is the projected separation, and v is the projected velocity. Below, we describe the calculation taking the mass ratio of 0 as an example, which is an extreme case assuming the barycentre is on the brighter star. We note that the positions of the candidate to the barycentre vary with the mass ratio and affect the corresponding terms.

Based on the hypothesis that the candidate is a background star, we used the Besançon galaxy model² (Czekaj et al. 2014) to simulate stellar populations in the direction of the stellar system. We simulated stars with distances from 0 to 50 kpc within a solid angle of one square degree. We selected stars with 2MAS S H magnitudes within 2σ of the candidate’s magnitude, 21.79 − 22.30 mag, which contains 4373 stars. The SPHERE image covers regions with separations from the bright central star from $0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}15$ to $5_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}5$; thus, the first term is $P_{\mathrm{b}}=4373\times \frac{\pi {\left(5_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}5\right)}^2-\pi {\left(5_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}15\right)}^2}{{\left({3600}^{\mathrm{\prime }\mathrm{\prime }}\right)}^2}=3.2\mathrm{\%}$. Then, P(ρ|b) is the ratio of the area of the 2σ separation to the whole image, which is $\frac{{\left(5_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}367\right)}^2-{\left(5_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}327\right)}^2}{{\left(5_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}5\right)}^2-{\left(0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}15\right)}^2}=1.4\mathrm{\%}$. For the third term P(v|ρ_bg), we assumed the candidate was free-floating and fitted its proper motion and parallax using EMCEE. We derived its values as μ_α cos δ = −13.86 ± 4.02 mas yr⁻¹ and μ_δ = 8.70 ± 3.94 mas yr⁻¹, as presented in Fig. 10. The parallax is not well constrained because two of the three astrometric measurements are almost exactly one year apart. We compare its proper motion with the primary star and simulated background stars in Fig. 11. The primary star is within the 2σ range of the candidate. Both of them deviate from the majority of the simulated background stars. Then P(v|ρ_bg) is the probability that a background star has a proper motion comparable to the measured proper motion of the candidate. This is calculated by the ratio of the number of background stars that have proper motions within the 2σ contour of the candidate to the total number, which is 7.6%.

In the hypothesis that the candidate is a co-moving planet bound to the central stellar system, we started with the companion occurrence rate in the SHINE survey. Vigan et al. (2021) estimated that 5.7% of FGK stars have companions of 1–75 M_J with a separation range of 5–300 au. Due to a lack of observed statistics for planets beyond 300 au, we used this rate as the lower limit of companion occurrence rate between 5 and 760 au. With an age range of 19–28 Myr and 21.79–22.30 mag, we estimated a mass range of 2.7–4.5 M_J and used the power-law distribution of planetary mass, dN/dM ∝ M^−1.31 (Cumming et al. 2008) to calculate P(p) = 5.7% × 0.15 = 0.8%. We used the power-law distribution of separation, dN/da ∝ a^−0.61 (Nielsen et al. 2017), and the 2σ range of the separation measured at the first epoch to calculate $P(\rho \mid \mathrm{p})=\frac{{\left(5_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}367\right)}^{0.39}-{\left(5_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}327\right)}^{0.39}}{(5_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}{5}^{0.39}-{\left(0_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}15\right)}^{0.39}}=0.4\mathrm{\%}$. The relative measured projected velocity between the candidate and star A (from 2018 to 2023) is 4.8 ± 3.8 mas yr⁻¹ (3.4 ± 3.9 mas yr⁻¹ in RA and 3.4 ± 3.6 mas yr⁻¹ in Dec). The escape velocity for a companion at a similar projected distance is 3.5 ± 0.5 mas yr⁻¹, where the uncertainty is due to the uncertainty in the total stellar mass. This velocity is the upper limit of the projected escape velocity due to the unknown orbital inclination and phase. We therefore calculated P(v|ρ_pl) by generating 10⁷ points, assuming Gaussian distributions for the measured projected velocity ~𝒩(4.8, 3.8²) and escape velocity ~𝒩(3.5, 0.5²). Here, the upper limit of the fraction of the points with measured projected velocity smaller than the escape velocity is 36.3%. Then, we have $$\frac{P_{\mathrm{b}}}{P_{\mathrm{p}}}\sim \frac{3.2\mathrm{\%}\times 1.4\mathrm{\%}\times 7.6\mathrm{\%}}{0.8\mathrm{\%}\times 0.4\mathrm{\%}\times 36.3\mathrm{\%}}=\mathrm{2.8.}$$(5)

The candidate is three times as likely to be an unbound object than a bound planet assuming that the mass ratio of B to A is zero. We then repeated this calculation increasing the mass ratio to unity and show the probability ratio in Fig. 12. The reason for this change is mainly due to two factors: the fitted proper motion of the candidate assuming it is free floating and the measured projected velocity of this candidate, assuming it is bound. It is 3–240 times more likely being a background star than a bound planet, with the central star mass ratio varying between 0 and 1.

Although the probability of the candidate being a bound planet is lower compared to that of being an unbound object, we note that there are several confirmed planets beyond 300 au, such as YSES 1c (Bohn et al. 2020b), b Cen b (Janson et al. 2021), and COCONUTS-2b (Zhang et al. 2021). Therefore we cannot definitely rule out the possibility of this candidate being a planet by statistical analysis. A spectroscopic confirmation is therefore crucial to determine if a source is a planetary-mass object, brown dwarf, or a background star.

In addition, the statistical analysis is based on the following set of assumptions:

1. Our knowledge of the wide-orbit planet occurrence rate is incomplete: there are no statistical results on the planet occurrence rate beyond 300 au. We used two approximations: a) the planet occurrence rate is based on separations within 300 au of 150 stars from multiple stellar regions (Vigan et al. 2021) while the planet frequency statistics may vary with regions; and b) the number of planets as a function of mass is taken from Cumming et al. (2008), which was calculated for the radial velocity (RV) detected planets with periods smaller than 2000 days (within 3 au for solar mass stars). These RV planets are likely formed by core accretion, while wide-orbit planets are likely formed by different mechanisms; therefore, the planet distribution may not follow this distribution.

2. The proper motion calculation of the candidate (assuming it is a free-floating object) relies on the proper motion of the barycentre, which is approximated by the proper motion of the primary star measured by Gaia. However, the proper motion of the primary star is likely affected by the fainter star and the potential hidden companion. The same is true for the proper motion of the barycentre.

3. Due to the unknown orbit of the planet, we adopted the escape velocity as the projected escape velocity, which is the upper limit of the projected escape velocity.

4. We used a galaxy population model to simulate the background star model in the direction of the central star, which is a good representative of the background star density on a large scale. However, there is some uncertainty regarding the local background density that has not been taken into account. This star is likely located in a crowded field given that there are a total of ten sources in the FoV besides the central stars.

Given these uncertainties, the probability ratio should not be taken as the deciding factor on the nature of this source. Regardless, the motion of the candidate is the most significant among the other nine sources in the images and is also significantly distinct from the other sources. We show the relative motion of all sources in two cases in Fig. D.1 in Appendix D. If it is a free-floating object, it might be a planetary-mass object, a brown dwarf, or a distant background star. We also show the positions of low-mass stars in Fig. 9 using the colour calculated in the evolutionary model of Baraffe et al. (2015), with an age of 1 and 5 Gyr. Their colour also falls into the colour range of the candidate. The colours of the candidate are compatible with a wide range of background M dwarfs. For instance, if the candidate were to be an M dwarf of 0.15 M_⊙, its absolute magnitude would be 8.7 mag in H. To match the apparent magnitude of the candidate, it would have to be 4.8 kpc away and have a very high proper motion within the Galaxy.

Fig. 10 Proper motion and parallax posterior distribution of the candidate, assuming it is a free-floating object when the mass ratio of B to A is 0. The proper motion is in the unit of mas yr−1 and the parallax is in the unit of mas.

Fig. 11 Proper motion of the candidate assuming it is free-floating when the barycentre is on the primary star. The blue populations are the stars from 0 to 50 kpc within 2σ magnitude of the candidate generated by the Besançon galaxy model. The 1 and 2σ contours are shown by the dashed blue circles. The red populations are the posterior distribution of the candidate’s proper motion with black circles showing the 1 and 2σ range. The purple point is the proper motion of the primary star reported in Gaia and the errorbar is smaller than the marker size. The green points are Sco-Cen candidate members within 20 angular degrees of the primary star on the sky.

Fig. 12 Probability ratio of the candidate being a free-floating object to a bound planet of the central stars as a function of the mass ratio of the fainter star to the brighter star. The red horizontal line is where the probability ratio = 1.

5.2 Group membership

The unresolved central star was classified as a member of LCC by Pecaut & Mamajek (2016). The age we estimated for the primary star is more than 19–28 Myr, older than the typical age of LCC members. LCC is a subregion in Sco-Cen and is revealed to be composed of several substructures of different ages, ranging from 6 to 15 Myr (Ratzenböck et al. 2023; Žerjal et al. 2023). This leads to a question about the membership of the star. Using the online association membership analysis tool BANYAN Σ (Gagné et al. 2018), the star is classified as a field star by 99.9% with Gaia proper motion and parallax. Furthermore, we selected Sco-Cen candidate members within 20 angular degrees of the star on the sky from Luhman (2022) and show them in Fig. 11. The primary star is located outside of the scattering region of these stars and we find that the closest Sco-Cen star to the primary star is at a separation of $8_{\cdot }^{\circ }9$. So the primary star is likely not a member of Sco-Cen. Detailed analysis is required to classify the star’s membership, such as including the resolved absolute radial velocity measurements, instead of the unresolved radial velocity from Gaia. However, this kind of analysis is beyond the scope of this work.

6 Summary

We report the detection of a planetary-mass candidate companion imaged around 2M1006. We consistently detected it in the H band with marginal detection in J23 bands and non-detection in K12 and K_S bands from 2018 to 2023. The central star is resolved to consist of a G8–K0 young Sun analogue with $T_{\text{eff\hspace{0.17em}}}\cong {5196}_{-96}^{+162}\mathrm{K}$ and an M dwarf for the first time. We estimated an age of 19–28 Myr for the primary star. We find tension among the acceleration data, orbital fitting, SED fitting, and radial velocity difference for the two stars if they are indeed part of a binary system. Therefore, we speculated that there is at least one additional low-mass companion in the central stellar system: either the primary star is not a single star or the fainter star is not a single star. Due to the unknown barycentre of the central stellar system, we are unable to confirm if the candidate is a co-moving planet of this system. Nevertheless, the candidate shows the most significant proper motion compared to other sources in the FoV and shares a common proper motion with the primary star. Multi-band colours and spectroscopic observations are necessary to identify the properties of this candidate. A narrow-band colour of H2–H3 may enable us to distinguish a low-mass planet with a clear atmosphere from a dusty planet, sub-stellar, and low-mass star as the latter three types of objects would have H2–H3 > 0 mag; meanwhile a clear-sky low-mass planet would have a very blue H2–H3 colour. A narrow-band colour J2–J3 can also help identify a clear-sky low-mass planet, as it would have J2–J3 > 1.5 mag when the mass is smaller than 10 M_J. From the marginal detection in 2021, its J2–J3 is 1.02 ± 0.51 mag, which cannot distinguish it between a low-mass planet, brown dwarf, and M dwarf due to the large uncertainty. Future high S/N detections in H23 and J23 may help verify whether it is indeed a clear-sky low-mass planet. Spectroscopic observations will fully confirm its nature. Due to its faintness, only JWST is capable of obtaining its spectra. Long-term imaging and radial velocity monitoring of the central stars will be needed to ascertain the kinematic nature of the two stars, which is crucial to confirm the companionship of this candidate via a common proper motion analysis. With an apparent magnitude of 22.04 ± 0.13 mag in the H band and the age of the primary star, if it is confirmed to be a planetary-mass companion, it would a planet with a mass of 3–5 M_J at a separation of 730 ± 10 au from the primary star. This would confirm it as a low-mass planet similar to 51 Eri b and AF Lep b and place it among known planets with the widest orbits imaged so far.

Acknowledgements

We thank Rico Landman, Bin Jia, Richelle van Capelleveen, Zephyr Penoyre and Anthony G.A. Brown for their valuable discussion. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes 0101.C-0153(A), 108.220Q. 001 and 110.23 NY.001. Part of this research was carried out in part at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). This research has used the SIMBAD database, operated at CDS, Strasbourg, France (Wenger et al. 2000). This research made use of ASTROPY, a community-developed core PYTHON package for Astronomy (Astropy Collaboration 2013, 2018), SCIPY (Virtanen et al. 2020), NUMPY (Harris et al. 2020) and MATPLOTLIB, a Python library for publication quality graphics (Hunter 2007). This research made use of LIGHTKURVE, a Python package for Kepler and TESS data analysis (Lightkurve Collaboration 2018). This work has benefitted from The UltracoolSheet, maintained by Will Best, Trent Dupuy, Michael Liu, Rob Siverd, and Zhoujian Zhang, and developed from compilations by Dupuy & Liu (2012); Dupuy & Kraus (2013); Liu et al. (2016); Best et al. (2018, 2020). This research made use of the Montreal Open Clusters and Associations (MOCA) database, operated at the Montréal Planétarium (J. Gagné et al., in prep) We thank the writers of these software packages for making their work available to the astronomical community.

Appendix A Central star images

The images of the central stars from 2018-11-15 to 2024-03-23 in H and z′ bands are presented in Fig. A.1, including residuals after model fitting and subtraction.

Fig. A.1 Images of the central stars from 2018-11-15 to 2024-03-23. All images are rotated to the direction where north is up and east is left. Upper panel: star images before subtraction. Bottom panel: residuals after model fitting and subtraction.

Appendix B Posterior distribution of SED fitting for the central two stars

The posterior distributions of SED parameter fitting of the central two stars are presented in Fig. B.1.

Fig. B.1 Posterior distributions of the SED fitting parameters for the central stars.

Appendix C Radial velocity injection test

We injected artificial signals in the VIS-X IFU data to verify the PCA method we developed in Sect. 3.4. We fit a 2D Gaussian model to the primary star and used the fitted model as the PSF template. Then we created a 1D Gaussian signal with an amplitude of 0.2 as the injected spectrum which is of similar brightness to the fainter star. We scaled the PSF template with the 1D Gaussian signal and injected it at the same separation and blue shift as the fainter star but in the opposite direction. We retrieved the signal successfully using the same PCA method we developed for the fainter star as demonstrated in Fig. C.1.

Fig. C.1 Injection test for radial velocity measurements in VIS-X IFU data. Upper panel: Original data. Bottom panel: Fake Gaussian signal of the same contrast and same blue shift as the fainter star injected at the same separation but in the opposite direction, shown by the cross marker in white. Before the injection, there is no signal in the first component (K=1) in the upper panel. After the injection, we detected the emission at the supposed blue shift and amplitude in the first component (K=1) in the lower panel.

Appendix D Relative astrometry of all sources

Figure. D.1 shows the relative astrometry of all sources in FoV when the mass ratio = 0 and 1. Regardless of the mass ratio of the central stars, the candidate has the most significant motion compared to other sources except for the central stars and it deviates from the proper motion cluster of other sources.

Fig. D.1 Relative astrometry all sources in RA and Dec offsets. Left panel: Relative astrometry of all sources to the primary star (mass ratio = 0). Right panel: Relative astrometry of all the sources to the barycentre when the mass ratio = 1. The positions of all sources of the first epoch are calibrated to the origin, so the positions of the later epochs are subtracted by the positions of the first epoch. The candidate is shown by the triangle and square markers with the epoch 2021-12-02 in orange and the epoch 2023-03-03 in green. The sources are labelled with the same number in Fig. 1, using the same colour convention for epochs. The blue line is the trajectory of a static background star evolving from 2018-11-15 and the two points on the track are the positions on 2021-12-02 and 2023-03-03 calculated by the proper motion and parallax reported in Gaia.

References

All Tables

All Figures

Fig. 1 Full image of this system of the first epoch in the H band. Left panel: central stars. The fainter component is in the south-east direction of the primary star on 2018-11-15. Middle panel: stacked coronagraphic image centred on the primary star. There are ten sources in the field of view in total. The planetary-mass candidate is the faint source between the two bright background stars in the south-west. Right panel: zoomed-in image of the candidate. All images are derotated to the orientation, where north is up and east is left.

In the text

	Fig. 2 Posterior distributions of the orbital parameters by fitting the measured relative positions with ORBITIZE.
In the text

	Fig. 3 100 orbits drawn from the posterior distributions of the binary fitted with ORBITIZE.
In the text

Fig. 4 Simultaneous SED fitting of the primary star and fainter star to the unresolved broad-band photometric measurements of Gaia, 2MASS, and WISE with SPECIES. The black solid line is the combined best-fitted spectrum of the two stars. The black dashed line is the best-fitted spectrum of the primary star and the black dotted line is the best-fitted spectrum of the fainter star. Spectra drawn from 30 random samples of each star are shown in grey lines. The top panel shows the filter trans-mission profiles and the bottom panel shows the fitting residuals.

In the text

Fig. 5 MagAO-X/IFU data of the binary. Left panel: IFU images of the binary stacked along the wavelength channel. The image is not rotated to place North up. Middle panel: spectrum of star A and the first three principal components of star B’s spectra. The grey area shows the Hα line band. Star A presents Hα absorption. The first component presents Hα emission from star B. The third component (K=3) is the contamination from star A. Right panel: smoothed version of the middle panel by a Gaussian kernel with a standard deviation of 0.48 Å. The grey area shows the Hα line band.

In the text

	Fig. 6 Age estimation of 2M1006. Left panel: age estimation by Li depletion using BAFFLES (Stanford-Moore et al. 2020). Right panel: age estimation by combining the luminosity and Teff of the primary star from SED fitting with the isochrones and tracks of BT-Settl evolutionary model (Baraffe et al. 2015).
In the text

	Fig. 7 Images of the candidate in five epochs of H, J23, K12, and KS bands. The candidate is the faint source between the two bright background stars highlighted by the white circles. It is only visible in the H band with marginal detection in J23 bands.
In the text

Fig. 8 Relative astrometry of the candidate in RA and Dec offsets. Left panel: relative astrometry of the candidate to the primary star (mass ratio = 0). Right panel: relative astrometry of the candidate to the barycentre if the fainter star and primary star are an equal mass binary (mass ratio = 1). The yellow and green regions are the candidate companion positions allowed by the escape velocity within 1σ on 2021-12-02 and 2023-03-03, respectively. The blue line is the trajectory of a static background star evolving from 2018-11-15 calculated with the proper motion and parallax reported in Gaia and the two points on the track are the positions on 2021-12-02 and 2023-03-03.

In the text

Fig. 9 Colour-magnitude diagram of brown dwarfs. The grey points are field brown dwarfs with high surface gravity. The blue points are young or low surface gravity objects. The brown points are low-mass stars. Several directly imaged exoplanets are labelled with black text. The red arrow shows the possible positions of the candidate assuming it is at the same distance as the central stars. If confirmed, it would be one of the coolest directly imaged exoplanets.

In the text

	Fig. 10 Proper motion and parallax posterior distribution of the candidate, assuming it is a free-floating object when the mass ratio of B to A is 0. The proper motion is in the unit of mas yr−1 and the parallax is in the unit of mas.
In the text

Fig. 11 Proper motion of the candidate assuming it is free-floating when the barycentre is on the primary star. The blue populations are the stars from 0 to 50 kpc within 2σ magnitude of the candidate generated by the Besançon galaxy model. The 1 and 2σ contours are shown by the dashed blue circles. The red populations are the posterior distribution of the candidate’s proper motion with black circles showing the 1 and 2σ range. The purple point is the proper motion of the primary star reported in Gaia and the errorbar is smaller than the marker size. The green points are Sco-Cen candidate members within 20 angular degrees of the primary star on the sky.

In the text

	Fig. 12 Probability ratio of the candidate being a free-floating object to a bound planet of the central stars as a function of the mass ratio of the fainter star to the brighter star. The red horizontal line is where the probability ratio = 1.
In the text

	Fig. A.1 Images of the central stars from 2018-11-15 to 2024-03-23. All images are rotated to the direction where north is up and east is left. Upper panel: star images before subtraction. Bottom panel: residuals after model fitting and subtraction.
In the text

	Fig. B.1 Posterior distributions of the SED fitting parameters for the central stars.
In the text

Fig. C.1 Injection test for radial velocity measurements in VIS-X IFU data. Upper panel: Original data. Bottom panel: Fake Gaussian signal of the same contrast and same blue shift as the fainter star injected at the same separation but in the opposite direction, shown by the cross marker in white. Before the injection, there is no signal in the first component (K=1) in the upper panel. After the injection, we detected the emission at the supposed blue shift and amplitude in the first component (K=1) in the lower panel.

In the text

Fig. D.1 Relative astrometry all sources in RA and Dec offsets. Left panel: Relative astrometry of all sources to the primary star (mass ratio = 0). Right panel: Relative astrometry of all the sources to the barycentre when the mass ratio = 1. The positions of all sources of the first epoch are calibrated to the origin, so the positions of the later epochs are subtracted by the positions of the first epoch. The candidate is shown by the triangle and square markers with the epoch 2021-12-02 in orange and the epoch 2023-03-03 in green. The sources are labelled with the same number in Fig. 1, using the same colour convention for epochs. The blue line is the trajectory of a static background star evolving from 2018-11-15 and the two points on the track are the positions on 2021-12-02 and 2023-03-03 calculated by the proper motion and parallax reported in Gaia.

In the text