© V. Squicciarini et al. 2022

1 Introduction

Prior to the discovery of exoplanets, the Solar System appeared to be the inevitable outcome of a smooth process of planetary formation: a process going from the primordial protoplanetary disk to the present ordered sequence of planets with nearly circular orbits and a tendency for a decrease in density with distance, such that the rocky planets are followed by giant planets as a consequence of the initial temperature gradient in the disk itself (Hayashi et al. 1985). However, the detection of the first exoplanet around a solar-type star (Mayor & Queloz 1995) thoroughly challenged this view. With a mass of 0.5M_J, 51 Pegasi b lies so close to its star that it completes one orbit in 4.2 days. This planet became the archetype of a whole new family of planets, the Hot Jupiters. In the following 27 yr, the discovery of exoplanets with properties markedly different from the ones in our system – that is to say extremely eccentric orbits, orbits perpendicular to the equatorial plane of their star, or evaporating atmospheres – has profoundly changed our view of planet formation and evolution (Winn & Fabrycky 2015). Indeed, the current view is that of a complex process involving mechanisms such as mutual resonances, orbital migration, and scattering (see, e.g., Bonomo et al. 2017), which emphasizes in turn the need for larger surveys to investigate the diversity we see systematically.

Despite advancing at a steady pace, laying the foundations of (exoplanet) demographics, a comprehensive and unbiased census of the planetary population is still missing, since every detection technique is more sensitive to certain regions of the parameter space¹ (Gaudi et al. 2021). Most previous and ongoing surveys have focused on stars as massive as (or less massive than) the Sun², and about 90% of the 5000 known exoplanets lie closer to their host stars than the Earth is to the Sun. For a few years, this strong observational bias has started being alleviated by direct imaging which, in contrast to the prevailing radial velocity and transit methods, is more sensitive to thermal emission from giant planets orbiting far from their stars. Since at young ages (≲ 1 Gyr) giant planets are mainly heated by gravitational contraction rather than by stellar radiation, they fade in brightness over time: the best targets are therefore young and nearby stars.

The few previous direct imaging and radial velocity surveys targeting more massive stars than the Sun (e.g., Johnson et al. 2010; Reffert et al. 2015; Wagner et al. 2022) have shown that gas giant planets are more frequent around more massive stars, likely due to the increased dust and gas reservoir in protoplanetary disks (Mordasini et al. 2012). According to radial velocity studies, however, the occurrence frequency has a turnover at about 2 M_⊙ and goes down to zero at M > 3 M_⊙ (Reffert et al. 2015), at least for those orbiting within a few au of their star. This is in line with theoretical expectations from the most popular scenario for planet formation, the core accretion model (CA; Pollack et al. 1996). This mechanism, which is thought to be responsible for the formation of the giant planets in the Solar System, requires the growth of solid planetary cores from dust, a slow process requiring a few million years. Therefore, due to a more rapid dispersal of the protoplanetary disk around heavier stars, giant planet formation should be increasingly hampered and eventually halted around more massive stars (Mordasini et al. 2012).

Radial velocity surveys, though, are strongly biased toward close-in planets. Unlike CA, gravitational instability (GI; Boss 2003) may well form massive planets and substellar objects at a wide separation on a timescale shorter (~ 10⁴ yr) than that the lifetime of the disk. In addition, the ratio between the mass of the disk and that of the star is known to increase with stellar mass (Pascucci et al. 2016), favoring the onset of GI (Cadman et al. 2020). If GI is a viable formation pathway for giant planets, the dearth of planetary-like companions around massive stars would be due to selection effects inherent to radial velocity searches.

On the other hand, as already mentioned, massive planets at large separations can be detected through direct imaging provided that they are young enough. Past high-contrast imaging surveys (e.g., Rameau et al. 2013; Bowler 2016; Chauvin et al. 2017; Nielsen et al. 2019; Vigan et al. 2021) have focused on stars of spectral type A or later (M ≲ 2.4 M_⊙); to fill in the gap, in 2018 we started the B-star Exoplanet Abundance STudy (BEAST) survey (Janson et al. 2021b), exploiting the capabilities of the Spectro-Polarimetric High-contrast Exoplanet Research (SPHERE) instrument (Beuzit et al. 2019) at the Very Large Telescope (VLT): its scientific goal is to look for exoplanets around a sample of 85 B stars in the young (5-20 million years) Scorpius-Centaurus OB association (Sco-Cen), the nearest region of ongoing stellar formation. While the survey is still in progress, we recently reported the discovery of a 10.9 ± 1.6 M_J planet around the stellar binary b Cen, whose total mass reaches 6−10 M_⊙ (Janson et al. 2021a).

Here we present a second planetary system surrounding the massive star μ² Scorpii (μ² Sco), composed of a robustly-detected 14.4 ± 0.8 M_J companion and a probable closer 18.5 ± 1.5 M_J companion, and discuss its implication for the current scenarios of planetary formation. We present an extensive discussion of stellar properties in Sect. 2. Section 3 details the observations and Sect. 4 discusses the data reduction methods. In Sect. 5 we identify the companion candidates, treating with special care the innermost one. In Sect. 6 we discuss the possibility that the two apparently comoving candidates are not physically bound to the star. Aftering ruling out this possibility, we derive their properties and constrain their orbits. In Sect. 7 we discuss the properties of these low-mass ratio companions in a more general context and examine the implications of their detection on the origin of this class of companions around B stars.

2 The Star

μ² Sco (also known as HR 6252, HD 151985, HIP 82545, Pipirima³) is a naked-eye star belonging to the young Scorpius-Centaurus association (de Zeeuw et al. 1999). The main astro-metric, kinematic, and photometric properties of the star are reported in Table 1.

Given that several stellar parameters (notably distance, age and mass) are of utmost importance for the characterization of directly-imaged exoplanets, our primary goal was to reduce the large uncertainty on distance that had historically limited a self-consistent physical analysis of μ² Sco. After exploiting for the first time kinematic information to indirectly constrain its distance, we combined this new information with data from the literature to perform a Monte Carlo analysis that determines posterior distributions for stellar mass, age, radius and effective temperature. A discussion on the adopted priors is provided in Sects. 2.1–2.6, and the final derived parameters are reported in Sect. 2.7. Complete details on the derivation are provided in Appendix A.

2.1 Distance and Membership to the Eastern Lower Scorpius Group

The Scorpius-Centaurus association, the nearest extended region with ongoing stellar formation to the Sun (de Zeeuw et al. 1999), is classically divided into three subregions: Upper Scorpius (US), Upper Centaurus-Lupus (UCL) and Lower Centaurus-Crux (LCC; de Zeeuw et al. 1999). Each region in turn has a complex morphology with a high degree of substructuring. μ² Sco, in particular, appears to lie within a small clump of stars (a region centered on galactic coordinates (l, b) = (343.1°, 4.7°) with a radius of about 2° that, while classically assigned to UCL (de Zeeuw et al. 1999), has recently started being recognized in its own right as an independent entity (Röser et al. 2018). We refer to this group as Lower Scorpius (LS).

Starting from the Gaia DR2-based catalog of bona fide Sco-Cen sources by Damiani et al. (2019), we create a census of LS stars (575 members) through the criteria shown in Table 2. We also crop from the same catalog an UCL sample, to be used for comparison purposes, using the classical boundaries by de Zeeuw et al. (1999) and excluding the sources already assigned to LS.

Inspection of the Gaia EDR3 catalog (Gaia Collaboration 2021) reveals an inconsistency of the star’s proper motions with those of LS and a parallax that, although consistent with LS, suffers from a large uncertainty. The values of Gaia EDR3, however, are also inconsistent with those from Gaia DR2 and Hipparcos (Perryman et al. 1997) (Table B.1). The inaccuracy of the astrometric solution in Gaia for such a bright star (which saturates the detector) is expected without invoking the presence of an unresolved stellar companion (Appendix B). Adopting the robust long-term proper motion by Kervella et al. (2019), computed as the difference between the astrometry of Hipparcos (J1991.25) and Gaia DR2 (J2015.5), the star is fully consistent with membership to LS (Fig. 1).

We notice that our target (together with its sibling μ¹ Sco) actually lies in a peripheral area of LS, ~2° eastward of its core: using the more sensitive Gaia EDR3 to expand the census of LS to fainter stars, we constructed a catalog of members of this small clump that includes 73 stars, that we call Eastern Lower Scorpius (ELS; see Table 2, Fig. 2). ELS appears somewhat closer to the Sun than the whole LS: we adopt its parallax distribution (modeled as a normal, π = 5.9 ± 0.2 mas) as our parallax prior (Fig. 3).

2.2 Mass and Radius

The mass estimate in the first BEAST paper (Janson et al. 2021b) was based on an underestimated distance from Gaia DR2: we therefore expect a higher mass than previously assumed. We conservatively start from a uniform prior distribution, M ∈ [5,12] M_⊙. Likewise, we do not favor a particular initial value for the radius: R ∈ [3.5,10.0] R_⊙. For a B2V star on the main sequence we expect R ~ 4 R_⊙ (Pecaut & Mamajek 2013), but the classification of the star as a B2IV (Hiltner et al. 1969) argues for a larger radius.

Fig. 1 Proper motion components of LS (blue) and ELS (orange) members. Both →1 Sco and →2 Sco are consistent with membership to ELS.

2.3 Age

As already mentioned, the derivation of a reliable stellar age is crucial in the context of the characterization of directly-imaged planets, since it is needed to turn the observed fluxes into mass estimates. Direct isochronal age determinations for B stars are especially difficult: their luminosity has only a weak dependence on age because they are already on the main sequence after ~ 10⁵ yr, and deviations from theoretical expectations might arise due to unresolved companions, strong magnetic fields and/or rapid rotation (Barrado 2016).

An indirect dating technique, based on isochronal analysis of comoving stars rather than the target itself, can be envisioned, in a similar way as in Squicciarini et al. (2021): using the MADYS tool (Squicciarini et al., in prep.), we collected photometry from Gala DR2/EDR3 and 2MASS (Cutri et al. 2003) for both LS and ELS and compared it to pre-MS isochrones to simultaneously derive individual mass and age estimates for their members (Squicciarini et al. 2021). To avoid known systematic uncertainties on absolute ages with this method (see discussion in Squicciarini et al. 2021), we decided to employ the whole UCL (~4000 stars) as a control sample, comparing the derived age distributions with each other.

The age distribution of ELS stars, derived by comparison with the BHAC15 isochrones (Baraffe et al. 2015) of solar metallicity, is consistent with that of the whole UCL according to a Kolmogorov-Smirnov test (α = 0.05)⁴.

In the absence of tighter constraints, we adopt the age of UCL (t = 16 ± 7 Myr: Pecaut & Mamajek 2016) as our age prior.

We highlight that an independent age estimate could come in the next years from asteroseismology: based on its light curve from TESS (Ricker et al. 2015), the star shows evidence for β Cephei type variability and in particular for that of a slowly pulsating B-star (SPB) β Cep hybrid. A quantitative analysis of the observed wealth of pulsational modes – especially if complemented by new observations in different photometric bands – might constitute the subject of a future work aimed at estimating the core hydrogen fraction, closely related to stellar age.

Fig. 2 Sky coordinates of bona fide Sco-Cen members from Damiani et al. (2019), shown as black dots. Classical boundaries for UCL are indicated by the red box. A zoom of LS (blue) and ELS (orange) region is drawn in the upper-right corner; →2 Sco is indicated by the black star.

Fig. 3 Parallax distribution of ELS stars. We adopt → = 5.9 ± 0.2 mas as parallax prior for →2 Sco.

2.4 Reddening

To estimate the reddening E(B − V) in the direction of μ² Sco in a robust way, we considered six different determinations:

• Starting from UBV photometry from Mermilliod (2006) and uυbyβ photometry from Hauck & Mermilliod (1997), we employ the Q-method of Johnson & Morgan (1953) to deredden OB stars with the modern calibration of Pecaut & Mamajek (2013): E(B − V) = 0.022 mag;

• dereddening the uυbyβ photometry through the calibration of Castelli (1991): E(B − V) = 0.019 mag;

• dereddening the uυbyβ photometry with the older, independent calibration of Shobbrook (1983): E(B − V) = 0.020 mag;

• from Lyman α observations, Savage & Panek (1974) estimate the interstellar column density of HI toward μ² Sco to be N(H) = 2.5 x 10²⁰ cm⁻². Adopting the recent relation between interstellar reddening and hydrogen column densities by Lenz et al. (2017), valid in the low column-density regime as in our case: E(B − V) = 0.031 mag;

• starting from the absorption EW(D1) of interstellar neutral sodium (NaI) D1 line toward μ² Sco by Hobbs (1978), and using the reddening vs. EW(D1) correlation by Poznanski et al. (2012): E(B − V) = 0.019 mag;

• integrating the STILISM 3D reddening map (Lallement et al. 2018) along the line of sight up to a distance d = 165 pc: E(B − V) = 0.025 ± 0.024 mag.

By averaging these estimates, we obtain E(B − V) = 0.022 ± 0.006 mag. Adopting a total-to-selective extinction ratio R_V = A_V/E(B − V) = 3.05, appropriate for early B-type stars (McCall 2004), we estimate A(V) = 0.068 ± 0.015 mag and A(K) = 0.0062 ± 0.0014. Being negligible with respect to photometric errors, extinction values will be from this moment on treated as constants.

2.5 Effective Temperature

Concerning the stellar effective temperature, μ² Sco – persistently classified as a B2IV star (Hiltner et al. 1969) – has a nearly identical combination of colors and reddening (B − V = - 0.214, U − B = −0.844, E(B − V) = 0.022) to the Morgan-Keenan B2IV standard star δ Cet (B − V = −0.219, U − B = −0.850, E(B − V) = 0.018) which has a median T_eff in the recent literature of T_eff ≃ 21 600 K (Cardiel et al. 2021). Hence we do not expect the effective temperature of μ² Sco to differ too much from that of δ Cet.

In Table 3 we list several new T_eff estimates for μ² Sco based on photometric data from the literature. The estimates are derived in the following way:

1. photometry from Hauck & Mermilliod (1997), dereddened through Castelli (1991) using Balona (1994) T_eff calibration;

2. photometry from Hauck & Mermilliod (1997), T_eff/H_β relation by Balona (1984);

3. employing the (U - B)_O vs. T_eff trend based on B2IV standard stars Pecaut & Mamajek (2013);

4. photometry from Hauck & Mermilliod (1997), dereddened through Castelli (1991) using Napiwotzki et al. (1993) T_eff calibration;

5. photometry from Hauck & Mermilliod (1997), using [c1] index adopting Nieva (2013) T_eff scale;

6. comparing IUE spectrophotometry⁵ in the wavelength range [110-195] nm with the grid of model atmospheres by Castelli & Kurucz (2003);

7. photometry from Hauck & Mermilliod (1997), using [u-b] index adopting de Geus et al. (1989) T_eff scale.

We adopted the averaged value (T_eff = 21 900 ± 1000 K) as our T_eff prior.

2.6 Surface Gravity

The stellar surface gravity was not used as a free parameter, but rather as one of the observational constraints for the optimization tests. Based on literature estimates (Table 4), we adopted a value of log g ∊ [3.6,4.0].

2.7 Final Stellar Parameters

While we refer to Appendix A for extensive details on the optimization analysis, we list in Table 5 the final stellar parameters obtained for μ² Sco.

With a spectral type of B2IV, corresponding to a best-fit effective temperature of 21 700 ± 900 K, μ² Sco is one of the brightest stars in the BEAST sample (Janson et al. 2021b). Indeed, its fitted mass of M = 9.1 ± 0.3 M_⊙ qualifies the star, from a stellar evolution standpoint, as a “massive star” (M > 8 M_⊙): a star that could explode as a supernova in the next 10−20 million yr (Nomoto & Leung 2017), in particular as an electron-capture supernova (Nomoto 1984). As a 20 ± 4 Myr star, μ² Sco is currently evolving off the main sequence, consistently with its spectral classification.

We note that in the specific case of μ² Sco, the use of model isochrones -properly constrained to reduce degeneracies- should be considered reliable, given that the star does not show any evidence for unresolved companions (see Appendix B), strong magnetic field or rapid rotation (see Appendix C).

3 Observations

High-contrast near-infrared images and spectra of μ² Sco were acquired by means of SPHERE’s dual-band imager (IRDIS; Dohlen et al. 2008) and integral field spectrograph (IFS; Claudi et al. 2008) at the ESO Very Large Telescope (Beuzit et al. 2019). The first epoch was obtained on April 24th, 2018, the second one on June 4th, 2021 (Table 6). Both observations were carried out in the so-called IRDIFS-EXT pupil-tracking mode and are diffraction limited. In this mode, IRDIS was used in the dual-band mode ( see Vigan et al. 2010) over a field of view (FoV) of 11'' × 11'', while IFS collected images – which can be combined to build low-resolution spectra (R = 30) – in the λ ∊ [0.96 − 1.64] μm wavelength range for each position over a field of about 1.7'' × 1.7'', Nyquist sampled at the shortest wavelength. Once extracted with the DRH data reduction pipeline (Pavlov et al. 2008), each spectrum is 39 pixel long. The APLC2 coronagraph (Carbillet et al. 2011) was used, masking the star out to a radius of 92.5 mas. As for all observations in the BEAST survey, calibration observations were taken together with science exposures: they included a flux calibration, allowing normalization to the peak of the star image, obtained by offsetting the star position off the coronagraphic mask with a suitable neutral density filter to avoid saturation of the image; a center calibration, obtained by imprinting a sinusoidal pattern to the deformable mirror, providing satellite images of the star; and empty sky exposures centered a few arcsec from the star position. Details of the observations are provided in Table 6.

4 Data Reduction

The raw data were reduced at the SPHERE Data Center by means of the SpeCal pipeline (Delorme et al. 2017; Galicher et al. 2018). For IRDIS data, the standard Template Locally Optimized Combination of Images (TLOCI) reduction technique (Lafrenière et al. 2007) was employed, as in previous BEAST publications (Janson et al. 2021b,a) and an additional PACO simultaneous Angular and Spectral Differential Imaging (ASDI) analysis (Flasseur et al. 2020) was carried on at a later stage only with regard to one specific companion candidate (Sect. 5.2). For IFS we used three data analysis methods: ASDI and monochromatic Principal Component Analysis (PCA; Mesa et al. 2015), PACO ASDI (Flasseur et al. 2020) and TRAP (Samland et al. 2021). The calibration of true north and pixel scale employs observations of far compact stellar cluster as in (Maire et al. 2016), while the recovery of astrometry and photometry of the companion candidates relied on the injection of negative planets on their position and minimization of residuals.

The final IRDIS and IFS images are shown in Figs. 4 and 5, respectively.

5 Companion candidates

5.1 Identification of companion candidates

Two companion candidates (CCs) were identified inside the FoV of SPHERE’s Integral Field Spectrograph (IFS), while 34 sources were identified inside the wider field of the Infra-Red Dual-beam Imaging and Spectroscopy (IRDIS). The quality of the final IRDIS image is somewhat better for the second epoch, which is on average ~0.25 mag deeper than the first one: this explains why some of the faintest CCs in the IRDIS FoV are only detected in the second epoch. The derived properties of all the CCs within the IFS and IRDIS FoV, labeled through a unique numeric ID, can be found in Table E.1.

5.2 The closest companion candidate: CC0

Unlike the other CCs, the proximity of CC0 to the edge of the coronagraphic mask and its faintness invokes special caution to rule out the possibility of a systematic artifact. In order to assess the reliability of the detection of CC0, the analysis of IFS data was performed – as already mentioned above – by means of three different methods:

• ASDI-PCA consists in the subtraction of the first n principal components from the 4D data cube, obtained by combining the 39 individual monochromatic images – corresponding to the different DITs on the detector – after radially scaling them with wavelength; to avoid the data cube to be too large, we averaged each consecutive DIT; after dropping a few poor images, the final data cube includes 90 χ 39 = 3510 images. We considered n= 25 modes, a value suited for small separations, and averaged the resulting data cube according to the expected contrast spectrum for a late M star. We then constructed a signal-to-noise ratio (S/N) map where the noise is evaluated along rings at constant separation and corrected it for the low number statistics penalty factor given by Mawet et al. (2014). In this way, we detect CC0 in the first epoch as a source with peak S/N ~ 8.0, separation = 128 ± 2 mas and PA = 99.5 ± 1.1°. Although not significant (S /N = 2.5), a source at a similar location (separation = 123 ± 7 mas, PA = 114 ± 4°) is tentatively spotted in the second epoch too.

• PACO-ASDI is based on local learning of patch covariances, in order to capture the spectral and temporal fluctuations of background structures; the statistical modeling is exploited to provide a detection algorithm and a spectrum estimation method. The modeling of spectral correlations is useful in reducing detection artifacts. Using this algorithm, we detected the same source found with ASDI-PCA in both epochs but with higher S/N: S/N =10.6, separation ~126 ± 2, PA = 99.9 ± 0.8° for the first epoch, S/N = 5.0, separation = 125 ± 3, PA = 111.0 ± 1.6° for the second epoch. These S/N estimation include the contribution of spectral correlations (Flasseur et al. 2020) and are performed locally (at a scale of a few pixels). This method tends to be more efficient in detecting very faint signals but it also tends to produce detection maps with slightly higher residual correlations. We also checked the detection maps obtained with the less local (at a scale of a few dozen pixels) algorithm, given S/N slightly above 5 for 2018 data with and without priors. The distribution of both detection criterion is (approximately) Gaussian in the absence of sources. While the first epoch detection employed a flat spectrum (i.e., no prior) and a late-M spectral template, the second epoch detection was obtained only when using the latter prior;

• TRAP is a data-driven approach to modeling the temporal behavior of stellar light contamination, rather than the spatial distribution of the speckle halo. Due to the field-of-view rotation of the image sequence, each pixel affected by planet signal sees a distinctive light curve which can be modeled and thereby distinguished from other systematic temporal trends of the data. This is achieved by a causal regression model, trained on other pixel light curves at a similar separation from the host star and fit simultaneously with the planet model for each pixel. This approach works well at small separations: no training data is lost due to insufficient field rotation, as the model is not trained locally. Default parameters, as described in the original paper, were used. The resulting contrast maps for each channel were then fitted with a L-type, T-type and flat contrast spectral template, assuming a BT-NextGen stellar model (Allard et al. 2011) of the companion. In addition we used the extracted PACO-ASDI spectrum of CC0 as a template to derive contrast limits for the object. Given the measured contrast by PACO-ASDI, a S/N between 2.5 and 3.5 would be expected at the position, but no signal above S/N = 0.5 was detected.

Assuming Gaussian noise distributions and considering the number of independent points sampled in an IFS image (~2 × 10⁴), the S/N obtained with PCA-ASDI and PACO-ASDI corresponds to values of the False Alarm Probability (FAP) equal to ~10⁻¹¹, that are extremely low even accounting for the fact that we have observed at least once 75 stars so far in the BEAST survey. However, these values underestimate the real FAP because the distributions are actually not exactly Gaussian⁶ (likely because of the edge effects introduced by the corona-graphic mask) and the FAP is therefore higher and not well determined.

The physical nature of CC0 is reinforced by the observation that, at both epochs, CC0 is the peak of the IFS S/N map, and that the second-epoch PA is found within ~ 11° from the first-epoch PA, the separation being the same (d_CC0 « 127 mas) within errors. The probability p_RPA of this fact happening by chance was estimated in an empirical way. We expect the astrometric shift between the two epochs ∆r not to be larger than the fraction Ar_max of the face-on circular orbit with r = d_CC0 covered over the timespan At ≈ 3.1 yr separating the two epochs: (1)

where yr is the period of the circular orbit with s = 1 au d_CC0/parallax ≈ 21 au around μ² Sco (M_*= 9.1M_⊙). As a comparison, the observed shift of CC0 is ∆r_obs = 26 ± 4 mas.

Due to the presence of the coronagraph at ~100 mas, we conveniently define as “interesting area” the annular sector, centered on CC0-s first-epoch PA, with separation 100 mas < d < 200 mas and semiangular width .

We then examined all BEAST images and verified that – after deleting those with clear companion candidates – in 19 out of 99 IFS images the S/N peak lies within 200 mas from the central star. The probability p_RPA of a nonphysical peak of the S/N map (e.g., a speckle) to be found within the interesting area defined above is given by pRPA ~ 19/99 · 35°/360° ~ 4%.

Finally, the separation of CC0 is constant with wavelength at both epochs (Fig. 6), further supporting the detection, since separation of speckles from the center of the image is expected to be proportional to wavelength. Intriguingly, a bright spot exactly at the location of CC0 is also detected with TLOCI in IRDIS K1 and K₂ data at first epoch - although at a very low S/N ~ 2 –However, given the very low S/N of these IRDIS detections, we consider them as upper limits in the following discussion.

On the other hand, CC0 was not detected using TRAP. The limiting magnitude obtained with this reduction is brighter than the PCA-ASDI and PACO-ASDI detections (because the wavelength dependence of speckle position is not exploited by TRAP), but should result in a marginal detection; it is at present unclear if the nondetection is due to lacking sensitivity as no Spectral Differential Imaging is used, or because the ASDI techniques underestimate errors and false alarm probabilities. For the time being we suggest that caution must be taken about the actual detection of CC0, that we consider as a probable detection but not a definite one.

Fig. 4 IRDIS images of →2 Sco obtained with TLOCI (left panel: first epoch; right panel: second epoch). The star, artificially obscured by the coronagraphic mask, is at the center of the image. Several background sources can be easily seen as bright point sources. →2 Sco b is the source inside the white circle.

Fig. 5 IFS images of →2 Sco obtained with ASDI-PCA with 25 modes (left panel: first epoch; right panel: second epoch). The star is at the center of the image but not visible due to the presence of the coronagraphic mask – which artificially masks the innermost ~90 mas of the image – and the aggressive data analysis method used. The probable companion CC0 is the source inside the white circle. A background source is visible on the lower left of the first-epoch image.

Fig. 6 Position of CC0 at different wavelengths. (A) J-band image, first epoch; (B) H-band image, first epoch; (C) J-band image, second epoch; (B) H-band image, second epoch. The Y-band image is not shown due to the extreme faintness of the source. White circles indicate the best-fit positions from Table E.1; circle radii have been enlarged to ~10 times the uncertainties on best-fit positions to enhance visibility. Unlike a speckle (red circle), the separation of CC0 does not increase with wavelength.

6 Analysis of Companion Candidates

6.1 Companion Confirmation

The confirmation of a directly-imaged companion candidate is usually performed by checking if the proper motion is similar to that of the target star and significantly different from that of the background sources, that is stars that are located far behind the target but appear projected close to it. Taking the first epoch as a reference, sources with a null proper motion are shifted in the second epoch as a reflex motion due to the target’s proper motion over the elapsed time.

As expected from the low galactic latitude (b = 3.86°) of the target, the final images are abundant in detected sources: two companion candidates are seen inside the IFS FoV Fig. 5, while in the wider IRDIS FoV there are 46 objects (Fig. 4) detected either in the first or in the second epoch, with the majority – 34 – being observed in both epochs. Excluding CC1 – seen in both epochs by IFS – we are left with 11 dim CCs only seen by IRDIS in the second epoch.

Fig. 7 (K1 → K2, K) color magnitude diagram of known substellar companions and field objects. IRDIS CCs only seen in one epoch are shown as red dots. The position of CC0 ( = →2 Sco b) is consistent with that expected for a young object.

6.1.1 Photometric Analysis

Before analyzing the astrometric shifts of the CCs seen at both epochs, we tried to assess whether 9 dim CCs that are only seen in the second epoch could be confidently labeled as background contaminants based on their colors (Fig. 7). However, owing to the large color uncertainty, their position in the (K₁ − K₂, K₁) color-magnitude diagram (CMD) does not provide us with a definitive answer. The same applies for CC29 and CC30, only seen in K₁ in one epoch. Given that the geometric probability of a random alignment between a background source and a target star with separation d is ∝ d², we are only able to argue that the presence of several secure background sources that have a smaller separation to μ² Sco than each of these 11 sources⁷, is strongly suggestive of a background nature for the 11 sources as well. We excluded all these CCs from the following analysis, in absence of any evidence for a physical association to our target.

6.1.2 Astrometric Analysis

With regard to the IRDIS CCs seen in both epochs, we are able to confidently label 33 out of 34 objects as background interlopers through astrometric analysis (see Fig. 8).

On the other hand, CC0 and CC2 show a completely different motion compared to the remaining sources. The mean astrometric shifts of this cloud of sources with respect to the position expected for a static source are −5.5 ± 1.2 mas (rms: 6.7 mas) along right ascension, −15.2 ± 1.8 mas (rms: 10.2 mas) along declination. The astrometric shift of CC0 is −44.4 ± 3.5 mas along RA, −96.1 ± 4.0 mas along declination; that of CC2 is −40.7 ± 2.1 mas along RA, −77.9 ± 2.1 mas along declination. Hence, the observed shifts of CC0 and CC2 are ~8σ and ~9σ away from the cloud of background sources (we call it “CC cloud”). The conclusion is robust to systematic astrometric offsets between the two epochs, equally affecting all the CCs. CC0 and CC2 are then plausible candidates as physical companion. We examine their cases in more detail in Sect. 6.2.

Fig. 8 Astrometric motion of CCs present in both epochs. As the plot shows the difference in separation between the two epochs, a comoving source should be close to the origin (labeled by “2018”). A background source with null proper motion will move according to the dashed curve as a reflection of the target’s motion, ending in the blue circle labeled as ‘2021’ at the second epoch. The motion of CC0 (labeled by '0') and CC2 (labeled by ‘2’) is distinct from the cloud of background contaminants.

6.1.3 Motion of Background Sources

As shown by Fig. 8, the CC cloud is shifted on average by approximately–16 mas with respect to the position expected on the basis of the motion of μ² Sco. After carefully checking the centering of our images to exclude a large systematic effect, we are left with three possible causes (or a combination of them):

1. the adopted proper motion is not correct, the error being approximately 3 mas yr⁻¹ in both right ascension and declination;

2. the photocenter of μ² Sco is offset with respect to the barycenter of the system by a different quantity at the two epochs. This might or might not be related to point 1;

3. the field stars have on average a proper motion that is not null. This might be a reflection of the galactic rotation curve (it is important to note that the line of sight toward μ² Sco passes at about 15 degrees from the galactic center and the background interlopers are expected to be stars of the inner part of the galactic disk at a few kpc from the Sun).

Scenario 1 seems unlikely: the proper motion considered for μ² Sco is very similar to that of its neighbor μ¹ Sco (μα = −11.867 ± 0.043 masyr−1 in RA and μδ = –22.511 ± 0.035 mas yr−1 in Dec) and to the average value for ELS (μα = −11.57 ± 0.48 mas yr−1 and μδ = −22.12 ± 0.70 mas yr−1); the differences are too small to explain the observed residuals.

If the correct scenario were 2, it would apply to both field stars and physical companions; indeed, the motion of CC0 and CC2 appears similar to the offset of the field stars with respect to the prediction for null motion. This would require a substantial motion of μ2 Sco caused by an (unseen) companion with an orbital period less than 24 yr; a massive short-period companion would likely cause large variations in the radial velocities, which instead are rather constant within a few hundreds of m/s over ~10 yr (see Appendix B).

To investigate scenario 3, we looked for Gaia EDR3 data of field stars near μ2 Sco. Given that none of the IRDIS background stars is bright enough to be present in Gaia, we searched for field stars within 5 arcmin from our target (we label it “Gaia bg sample”). The main properties of the Gaia bg sample, composed of 6286 stars, are:

The mean proper motion of the CC cloud, that is simply the ratio between the mean astrometric shift and the time baseline, is instead −1.77 ± 0.37 mas yr−1 (rms: 2.17 mas yr−1) along RA and −4.87 ± 0.56 mas yr−1 (rms: 3.29 mas yr−1) along Dec. We anticipate that a small systematic offset, of about −0.6 ± 0.4 mas along RA and + 1.7 ± 0.4 mas along Dec, is likely to affect the proper motions of our CCs (see Sect. 6.2.1). With this caveat in mind, we see that the proper motion of the Gaia bg sample is indeed similar to that of the CC cloud; the latter is likely the M-type part of the same population probed, at brighter magnitudes, by the former, as indicated by the following qualitative argument. From the median parallax of the Gaia bg sample we infer a distance of ~3 kpc, and absolute magnitudes MG ~ 7.5, that correspond to late K stars; since we are neglecting reddening, these objects may be even intrinsically brighter. This G magnitude translates to MK ~ 4.8 according to the tables by Pecaut & Mamajek (2013). The corresponding apparent magnitude K ~ 17.4 gives a contrast with respect to μ2 Sco of dK = 13.1 mag; these stars are roughly two magnitudes brighter than the CC cloud (median contrast: dK ~ 15 mag). If we assume that IRDIS CCs are on average at the same distance of the Gaia bg sample, their absolute K magnitudes are MK ~ 6.7; according to the tables by Pecaut & Mamajek (2013) this corresponds to M3V stars that have an absolute MG ~ 10.0, and to an apparent magnitude G ~ 22.6. IRDIS CCs are too faint to be detected by Gaia (that has a limiting magnitude of G ~ 21), and conversely we expect about a couple of Gaia stars (with a density of 1 star every 45 square arcsec) within the IRDIS FoV; the relative frequency agrees with expectation for a reasonable mass distribution (e.g., Salpeter-like) if they belong to the same parent population.

6.2 Confirmation of Physical Companion(s) to μ² Sco

In this subsection, we examine the possibility that CC2 is not a bound companion of μ2 Sco considering two possible alternatives:

• that CC2 is a high proper motion background star;

• that CC2 is a brown dwarf (BD) member of Scorpius-Centaurus that appears projected close to μ2 Sco.

We note that while the following considerations refer to CC2, they equally apply to CC0 which has a higher proper motion and a smaller separation than CC2, hence intrinsically lower false alarm probabilities (FAP).

6.2.1 CC2 is Not a Background Star

Having shown that the cloud of sources seen in the IRDIS images is made of background interlopers, it is necessary to estimate the probability that an object drawn from the same population could have been considered an “interesting” companion candidate to μ2 Sco because of a “remarkable” astrometric shift relative to μ2 Sco. Two factors must be taken into account: the probability of finding an object within certain boundaries of and μδ; the fact that we have fully8 observed 25 stars in the survey.

Starting from the Gaia bg sample, we identify as interesting the stars with mas yr−1 and μδ < –21.1 mas yr−1 9. After excluding a few objects that, despite passing the test, have parallax >5 mas (which would correspond to 4 ≲ dK ≲ 10 mag) and are thus so bright that they would have been disguised as stellar – and not substellar – companions, we obtain that 5 out of 6286 stars satisfy this criterion. This implies a fraction of interesting background objects of fG = 1.67 × 10−5. The FAP of finding one star with these features in our observations of μ2 Sco (that is, within the IRDIS field of view), given that we observe 36 CCs in both epochs, is then 6 × 10−4.

A similar argument can be made by creating a synthetic sample of background stars by means of the Besançon Galaxy Model interface10 (Czekaj et al. 2014). We selected a sample of stars with distance d ∈ [0,50] kpc, apparent K magnitude similar to that observed for CC0 (K ∈ [15,17] mag), radius ρ = 5 arcmin around the position of μ2 Sco. Out of 6595 stars, just 4 (corresponding to fB = 1.34 × 10−5) pass the proper motion test described above; this yields a FAP -given 36 CCs- of 4.8 × 410−4.

Thus, the Besançon-based test and the Gaia-based test give very similar results. We point out that, by defining an unbound region of the proper motion space of ‘fast-moving stars’, we are actually overestimating the FAP: a source with an unusually large absolute value of the proper motion would not have been consistent with being bound to the star, and would have been discarded.

The probability of having seen at least one background source with these features throughout the entire survey can be estimated in this way. The median proper motion of BEAST targets is -in absolute value- larger than that of μ2 Sco, but we assume for simplicity a strict equality; if we further assume that the fractions fB and fG computed above do not depend on the sky coordinates, we can extend the above reasoning to the whole survey. Given that we have fully observed 25 targets, and that on average we see ~11 sources per observation (Janson et al. 2021b), we obtain a FAP of having seen at least one background object disguised as a possible companion within the entire survey of 4.6 × 10−3 (Gaia-based test) and 3.7 × 10−3 (Besançon-based test). We conclude that CC2 (and even more likely CC0) is not a background object at a high level of confidence.

6.2.2 CC2 is Not a Free-floating Brown Dwarf

Whilst confirming that the sources are comoving, the astrometric analysis does not strictly allow exclusion of an alternative scenario: namely, that the sources are free-floating UCL substellar objects that happen by chance to be close to the line of sight of μ2 Sco. In order to quantify this false alarm probability, three points need to be considered: the probability of finding an “interesting” (M = 1–75 MJ) object; the probability of finding it within the IRDIS field of view; the fact that we have fully observed 25 stars in the survey.

Assessing the number of free-floating UCL members necessarily requires some assumptions on the initial mass function (IMF) of the association11.

Recently, Miret-Roig et al. (2022) have uncovered a rich population of free-floating planets and brown dwarfs in Upper Scorpius, extending the IMF of the association down to 0.005 M⊙. As US is a subregion of Sco-Cen, we expect that this IMF can be safely adopted for UCL too. By normalizing their IMF we obtain a probability density function (PDF); we are thus able to compute the fraction of objects in the mass range [5 MJ, 75 MJ]. We do not consider objects below 5 MJ not only because they are not covered by their data, but also, more importantly, because this is approximately the lower mass to which we are sensitive in BEAST. Integration of the PDF yields 0.205, meaning that one out of five objects in Sco-Cen is expected to belong to this mass range. This PDF should be multiplied by the projected density of objects in the region around μ2 Sco. To estimate this quantity, we take again the list of bona fide members compiled by Damiani et al. (2019). The faintest members of UCL have apparent G = 18.5–19.5, and the corresponding masses reach down to 0.013 M⊙ (see “Stellar system analysis”). Since Gaia is complete within this magnitude range (its limiting magnitude is G ~ 21), we assume a sharp transition between 100% completeness above 15 MJ to 0% completeness below 15 MJ. We rescale the number of sources we see in UCL (3842), ELS (575) and LS (73) by dividing by the integral of the PDF above 0.015 M⊙, obtaining a complete census of approximately ~4021 sources for UCL, ~602 for LS, ~76 for ELS.

Multiplying these numbers by the integral of the PDF from 5 MJ to 75 MJ, we obtain the expected number of interesting objects Ñ:

In order to turn the expected number of objects into a projected density (i.e. the number of expected BD interlopers per arcsec2), the areas of LS, ELS and UCL must be evaluated. For UCL, coordinate boundaries (l, b) = [313°, 343°] × [2°, 28°] as in Pecaut & Mamajek (2016) were used. For LS and ELS we computed the mean 〈l〈, 〈b〈 and the related standard deviations σl and and defined the areas as [〈l〉 – 2σ_l, 〈l〉 + 2σ_l] × [〈b〉 –2σ_b 〈b〈 + 2σ_b,]. We get:

so that the projected densities are:

The mean density of sources in the environment of BEAST stars is usually not as high as in LS. The value that is more representative of the median behavior of the sample is probably that of UCL; we retain the value for LS as a (very high) upper limit to this source of contamination. It is possible now to estimate the number of BD interlopers expected in the IRDIS FoV (a square of 11″ × 11″):

and finally the false alarm probability of having seen at least one of these objects across the 25 targets that we have observed at least twice using a binomial distribution is:

The probability of a chance alignment by free-floating objects is negligible for CC2 (and even more for CC0).

Finally, the position of both CC0 and CC2 is constant with wavelength, as expected for physical sources. We are therefore able to confirm in a robust way the source CC2 as being a substellar companion to μ2 Sco (μ2 Sco b). As regards CC0, based on the reliability arguments presented in the Sect. 5.2, we consider it to be a probable detection that would make it a second physical companion to μ2 Sco.

6.3 Companion Characterization

The star μ2 Sco appears to be surrounded by one, and possibly two, physical companions. We recall that, while the existence of CC2 is firmly established (S /N =80–120), the detection of CC0, owing to its extreme proximity to the target, is more sensitive to the reduction method employed, stretching from a robust detection (S/N = 10.6 in the first epoch, S/N = 5.7 in the second epoch) in one case to a nondetection where a marginal detection (S/N = 2.5-3.5) was expected in another case.

Given the not unambiguous outcome of the different reductions and the subtleties in the derivation of confidence levels and false alarm probabilities, we highlight the need for follow-up observations to definitely confirm or rule out its existence. With this caveat in mind, throughout the text, we continue to refer to both the robust μ2 Sco b and the probable CC0, deriving the properties of the latter from the results of our successful detections.

6.3.1 Spectra and Photometry

The PACO-ASDI algorithm (Flasseur et al. 2018) provides a spectrum for the probable companion CC0 for both epochs. As expected, the spectrum relative to the second epoch is characterized by a significantly lower S/N than the first-epoch one12. From this moment on we always refer to the first-epoch spectrum, the only one which allows a tentative characterization of CC0.

After combining the spectrum from the first epoch with upper limits in the K1 and K2 band provided by the same algorithm applied to the IRDIS first-epoch dataset (see Fig. 9), we fit the spectrum with those provided by the AMES-Dusty models (Allard et al. 2001) – suitable for substellar objects and low-mass stars – assuming our age and distance estimates for μ2 Sco. The best-fit model has a mass of Mc = 19.5 ± 0.9 MJ, a temperature of Teff = 2262 ± 28 K and log L/L⊙ = –3.08 ± 0.03. The comparison is generally fairly good (reduced;χ2 = 2.55); most of the contribution to ;χ2 is due to the region around the J-band, where the observed spectrum is much lower than expected. As no strong molecular band is expected in the region around 1.25 μm according to any realistic atmospheric model (compare, e.g., Marley et al. 2021), we attribute this discrepancy to residual speckle noise.

A similar fit with BHAC15 models yields Mc = 17.4 ± 0.9 MJ, Teff = 2274 ± 28 K and log L/L⊙ = –3.04 ± 0.03. We account for theoretical uncertainties by averaging the two estimates to obtain a final mass of Mc = 18.5 ± 1.5 MJ.

We highlight that, despite the high level of irradiation at its expected orbital distance (21 ± 1 au), radiation from the probable companion would be almost completely due to its self-luminosity rather than reflected light from the star. In fact, for an albedo similar to the Earth, the equilibrium temperature of CC0 is about 900 K; comparing it with our best-fit effective temperature of about 2270 K, we estimate that stellar irradiation contributes only about 2% to the total luminosity of the object, a value well within observational uncertainties.

As regards μ2 Sco b, for which we could not extract a spectrum as it lies outside the IFS FoV, the available photometric information is limited to the measured contrasts in the K1 and K2 band; combining them with the K magnitude of the primary, we derive absolute (K1, K2) magnitudes. No significant photometric variation is observed between the two epochs. The position of the companion in the (K1 – K2, K) color-magnitude diagram confirms its compatibility with a substellar object lying at the very beginning of the sequence of L dwarfs (L0-L2 type; see Fig. 7).

A comparison of the photometry with theoretical magnitudes from the same two models used for CC0 provides us with two mass estimates, which again we average to get a final mass Mb = 14.4 ± 0.8 MJ.

As regards CC0, the reported magnitudes are computed by collapsing IFS spectral channels 12–21 (1.159–1.333 μm) and 30–38 (1.504–1.636 μm) from the first epoch to build the two bands JIFS = 1.246 μm (band width = 0.174 μm) and HIFS = 1.570 μm (band width = 0.132 μm); the contrasts in the two bands are translated into absolute magnitudes by adding 2MASS J and H magnitudes of the primary.

Fig. 9 Contrast spectrum for CC0 obtained with PACO-ASDI from 2018 IFS data (black dots) and from PACO-ASDI for IRDIS data acquired at the same epoch (5→ upper limits; red squares). The second-epoch ASDI-PCA spectrum (gray squares) is shown for reference. The solid line is the contrast spectrum expected for a 20 Myr old, 19.5 MJ brown dwarf using AMES-Dusty models. The dot-dashed line is a black body curve with the same temperature and radius as this model.

6.3.2 Astrometry and Orbits

The position of μ2 Sco b and CC0 relative to the star is not exactly the same in the two epochs, showing a significant displacement of ~10 mas and ~20 mas, respectively: a possible hint of orbital motion around the star. Before analyzing this aspect quantitatively, we quantified the accuracy of the relative astrometry provided by SPHERE observations exploiting the large amount of background CCs present in the images.

Starting from the Gaia bg sample defined in Sect. 6.1.3 (6286 sources), we select the 2741 sources with parallax 0.1 mas < π < 5 mas (i.e., 0.2 kpc < d < 10 kpc). The mean proper motion of this “restricted Gaia bg sample” is –2.51 ± 0.07 mas yr−1 (rms = 3.53 mas yr−1) along right ascension, and –3.82 ± 0.07 mas yr−1 (rms = 3.76 mas yr−1) along declination.

To test whether the two samples are drawn from the same parent distribution, we performed two independent Kolmogorov-Smirnov (KS) tests (one for , one for μδ) at level α = 0.05. While the null hypothesis could not be rejected for (p = 0.13 > α), a significant difference exists with respect to μδ (p = 0.008 < α). In particular, the CC cloud appears to have a somewhat higher and lower μδ than the restricted Gaia bg sample. To quantify this, we identified the range of and μδ shifts to be solidly applied to the whole all CC cloud so that the KS test is passed at level α = 0.05.

The range of shifts needed for mas yr−1, while for μδ it is mas yr−1. The mean values within these ranges are –0.61 mas yr−1 and 1.68 mas yr−1 , respectively. Multiplying by the temporal baseline, and equally splitting the correction between the two epochs, we get Δα = –1.34 mas, Δδ = 3.73 mas, which can be attributed to a not perfect centering of the star in one or both epochs (note that this is much smaller than what cited at Sect. 6.2.1 and well within typical uncertainties of star centering in SPHERE astrometry; see Maire et al. 2016).

Instead of using these values to fix the astrometry of μ2 Sco b and CC0, we conservatively opted for treating Δα and Δδ as an additional source of random uncertainty on their relative astrometry, and propagate it to derive final uncertainties σd and σPA on separation and PA, respectively, that are somewhat broadened with respect to those in Table E.1: for CC0,(σd,σPA) = (2.5 mas, 2.0°) and (3.3 mas, 2.2°) for the first and second epoch, respectively; for CC2 = μ2 Sco b, (σd, σPA) = (3.1 mas, 0.12°) and (2.8 mas, 0.11°) for the first and second epoch, respectively.

Starting from the measured separations and PA and their broadened uncertainties, orbital parameters were estimated separately for the probable companion CC0 and for the robust one μ2 Sco b using the orbitize! Python package13, run with the recommended parameters for reliable convergence, for a total of 2 × 106 orbits (Blunt et al. 2020). The priors used for the MCMC were stellar mass and parallax, taken from Table 5, as well as the starplanet separations and position angles, taken from Table E.1. A subsample of suitable orbits, as well as the posterior distributions of semimajor axis, eccentricity, and inclination, are shown in Figs.10–13.

To understand whether the best-fit orbital configuration of the system can be dynamically stable, we refer to the Hill criterion. Let M* be the mass of the primary star, (m1, a1, e1) and (m2, a2, e2) the mass, semimajor axis and eccentricity of the inner and outer companion, respectively. Under the hypothesis that m1 ≪ M* and m2 ≪ M*, a system is Hill stable, meaning that the two companions will avoid close approaches at all times, if:

(Gladman 1993). Here .

Using nominal values for orbital parameters and masses, the orbits of b and CC0 are Hill stable; more accurately, taking the uncertainty on the orbital parameters into account, the system is Hill stable about 70% of the time. Likewise, the high nominal eccentricity values make it difficult to have additional stable orbits over a wide range of semimajor axes (from ~5 au up to ~800 au), though this cannot be excluded at present due to the large uncertainties still existing on the orbital parameters.

The final parameters of μ2 Sco b and the candidate CC0 derived throughout this section are reported in Table 7.

Fig. 10 Subset of possible orbits for the robust companion candidate →2 Sco b. (A) Sky-projected orbital fit (orange star icon) based on the 2018 and 2021 epochs which overlap at this scale, showing 100 randomly drawn orbits from the orbitize! MCMC chains. The position of →2 Sco is indicated by the yellow star. (B) Evolution over time of the planet-star separation (→); measured points and the corresponding error bars are shown in purple. (C) Evolution over time of the position angle (PA); again, measured points and the corresponding error bars are shown in purple.

7 Discussion

The two objects, having Mb = 14.4 ± 0.8 MJ and Mc = 18.5 ± 1.5 MJ, are just above the deuterium-burning limit (M ~ 13 Mj) that is classically used to mark the transition between planets and brown dwarfs. However, given the host mass, the mass ratios of these companions to the star are two of the lowest for objects discovered in imaging (qb = 0.0015, qc = 0.0019) and comparable to that of Jupiter to the Sun (q = 0.00095). This qualifies them as planet-like from this point of view. Very interestingly, the mean irradiation that the two planet-like companions of μ2 Sco receive from their parent star is similar to those of two Solar System planets: while the outer one has a mean irradiation similar to that of Jupiter, the inner one should have an irradiation comparable to that of Mercury, that would make it, if confirmed, the most irradiated substellar companion discovered by direct imaging so far (Fig. 14). The μ2 Sco system appears to be in many regards a scaled-up version of the Solar System.

Leaving aside the case to be confirmed of CC0, the robust discovery of a physical companion to μ2 Sco with a mass ratio similar to that of Jupiter to the Sun is not only noteworthy as an individual case, but acquires its greater significance when coupled to the recent discovery, within the same survey, of a similar object around the binary b Cen system. In that case the primary star has a mass of ~6 M⊙ and the overall mass is similar to that of µ2 Sco (Janson et al. 2021a). The possibility of capture of an object originally formed elsewhere, owing to the very low density of the Sco-Cen association to which both targets belong, was already unlikely for a single system14; the discovery of similar companions in two different systems firmly argues for in-situ formation, pointing toward a whole range of questions about how these objects came into being.

µ² Sco b and b Cen b are added to the growing population of directly-imaged substellar companions (Nielsen et al. 2019; Vigan et al. 2021); as shown by Fig. 14, this population is able to probe a unique niche of the parameter space, building upon the sensitivity of direct imaging to young massive companions at large separations from their stars.

Companions with 1 MJ ≲ M ≲ 40 MJ are more frequently found around stars more massive than the Sun in radial velocity studies (Reffert et al. 2015; Wolthoff et al. 2022), and the trend has been independently confirmed by direct imaging surveys15 (Bowler 2016; Nielsen et al. 2019; Vigan et al. 2021). As mentioned in Sect. 1, radial velocity studies indicate a turnover in this trend at about 2 M⊙; but the increasing separation of the snow line – beyond which most giant planets are thought to form – with stellar mass (Kennedy & Kenyon 2008), coupled with the scarce sensitivity of RV to separations larger than a few au, do not allow to conclude whether a wide-orbit population around more massive stars exists. Direct imaging is therefore the only technique able to investigate the extreme limits of planet formation. Early results from the BEAST survey are showing that giant exoplanets or small brown dwarfs can form even around B stars, but the completion of the survey is necessary to insert these qualitative findings into a robust statistical framework; a full comparison of companion mass and separation distributions with those from other large direct imaging surveys like SHINE (Chauvin et al. 2017) or GPIES (Nielsen et al. 2019) is expected to quantify the relative contribution of different formation pathways, clarifying the role of stellar mass in shaping planetary systems.

We must ask ourselves at this point what we mean by “planet” and “planetary system”, and if the μ2 Sco system should be seen as a planetary system or rather as a multiple stellar system composed of a massive star and one or two brown dwarfs. Just as the transition from substellar to stellar objects at 75–80 Mj is determined by the possibility to ignite hydrogen burning in their core, the transition from giant planets to brown dwarfs is often set to ~ 13 MJ, the so-called deuterium-burning limit (DBL). According to this definition, the μ2 Sco system should be considered a multiple stellar system, while b Cen b would be a (circumbinary) giant planet. The clear similarities between the two systems, though, highlight how a distinction that is based uniquely on a process happening in the companion core might not be adequate in every circumstance, and not necessarily related to different formation pathways.

On the other hand, a second distinction can be operated between objects being formed “like stars”, that is through turbulent core fragmentation, and objects being formed “like planets”, that is within a protoplanetary disk. This definition generally agrees with the previous one for solar-like stars: a correlation exists between the occurrence of giant planets with M < 4 Mj and stellar metallicity, hint of formation within a protoplanetary disk; more massive (M ~ 10 MJ) objects do not show this correlation (Schlaufman 2018), pointing toward a star-like formation process. However, when analyzing lower or higher stellar masses, the situation becomes more and more blurred: it is here that the mass ratio comes into play. On the one hand, a few known giant planet companions to very low-mass primary stars or brown dwarfs, likely outcome of turbulent fragmentation within the natal molecular cloud (e.g., Fontanive et al. 2020), should be considered “star-like”. On the other hand, the two brown dwarfs companions (M sin i = 22 MJ and 25 MJ) to the evolved 2.7 M⊙ star ν Oph show a 6:1 mean motion resonance of their orbits, a clue of a formation within a protoplanetary disk, and should be labeled as “planet-like” (Quirrenbach et al. 2019). This issue has been recently discussed within the IAU Commission F2, and led to a revised version of the definition of planets: although no explicit distinction based on the formation pathway (which is not easily related to physical properties) was set, an upper limit was established to the planet-to-star mass ratio of q < 0.04 (Lecavelier des Etangs & Lissauer 2022) in addition to the DBL mass limit at M < 13 Mj. While this revision is appropriate for systems around low-mass stars, a general observation is that, due to a scaling of disk mass with star mass (see, e.g., Pascucci et al. 2016), one could expect that more massive stars – originally surrounded by more massive protoplanetary disks – can form companions within a disk that are more massive than the deuterium-burning limit. We notice that the mass ratios q of µ2 Sco b, the candidate CC0 and b Cen b are consistent with one another and – within a factor two – with that of Jupiter to the Sun (Fig. 14); the mass ratios are smaller than those in the ν Oph system and the recently-adopted revised IAU definition of planet. As protoplanetary disks around massive stars can extend to thousands of au (Beltrán & de Wit 2016), μ2 Sco b, the candidate CC0 and b Cen lie well within the primordial extent of the protoplanetary disk of their host stars. In this sense, μ2 Sco b and b Cen b should both be considered as planets.

A third way to assess whether μ2 Sco b and the candidate CC0 are to be considered “planets” is based on a purely empirical basis. Giant planets around less massive stars than μ2 Sco show a bottom-up mass distribution, that is, a larger occurrence of less massive planets; this is not true for stellar companions, which instead show a top-down distribution that favors larger values of q. This stellar companion population, by construction, encompasses every possible pathway for multiple star formation, ranging from turbulent core fragmentation, to fragmentation of pseudo-disks, to GI within a protoplanetary disk16. One might wonder if objects like μ2 Sco b and b Cen b constitute the highmass tail of the bottom-up planet-like population (PP), or rather the low-mass tail of the top-down star-like population (SP).

A tentative comparison between a SP and a PP – based in turn on previous multiplicity and direct imaging studies, respectively – can be made in a similar way as in Janson et al. (2021a): for the SP, we adopt a log-normal separation distribution as in (De Rosa et al. 2014) – suited for A stars – and a mass ratio distribution as in Reggiani & Meyer (2013): (3)

normalizing it to the median frequency of 1 – 75 MJ companions at [5–300] AU for BA hosts taken from Vigan et al. (2021). For the PP, we adopt the parametric model from Vigan et al. (2021) and Meyer et al. (2021) with the set of parameters referring to BA stars: (4)

The expected number of companions around one star with a e [100, 1000] au, q ∈ [0.0005,0.0030] is 1.4 × 10−4 for the SP scenario, 9.0 × 10−3 for the PP scenario. Taking into account that we have observed 25 targets in the survey, the probability of finding at least one companion within these ranges is 3.4 × 10−3 for the SP scenario and 0.20 for the PP scenario. If we further assume that – as suggested for stars with M < 2.5 M⊙ (Vigan et al. 2021) – the peak of the orbital distribution shifts to larger separations with stellar mass, and we employ irradiation levels as our scaling factor, the probability under the PP scenario rises to 0.79. In this last case, the probability of finding at least two companions – as μ2 Sco b and b Cen b – is fairly high too (0.44). We stress that the comparison is based on a naive extrapolation of the known frequencies of planetary and stellar companions from A stars to ~9 M⊙ stars, and that a companion mass ratio distribution getting steeper at larger separations (De Rosa et al. 2014; Moe & Di Stefano 2017) might mitigate – at least partially – such a strong prevalence for the PP scenario.

From a formation standpoint, the bottom-up PP distribution is naturally associated with CA. CA assembles giant planets by building a solid core from dust present in the protoplanetary disk; as the core attains a critical mass (~10 M⊙), it rapidly starts accreting hydrogen and helium to become a gas giant (Pollack et al. 1996). The shift of the orbital peak of the PP distribution would be a natural consequence of the wider distance of the snow line, and the increased reservoir of gas and dust in the protoplanetary disk would make it possible to form an object as massive as µ2 Sco b; in other words, the PP scenario would explain the similarity of q ratio and irradiation with CA planets like Jupiter.

On the other hand, CA requires a few million years to operate, a timescale comparable to the disk lifetime around low-mass stars (Gorti et al. 2009; Gorti & Hollenbach 2009). The survival of a protoplanetary disk is mainly dictated by the strength of ultraviolet (FUV: 6 eV < hv < 13.6 eV; EUV: 13.6 eV < hv < 0.1 keV) and X-ray (0.1 keV < hν < 2 keV) irradiation from the central star. While disk lifetimes do not vary much for stellar masses in the range 0.3–3 M⊙, more massive stars (>7 M⊙) are expected to lose their disks rapidly (in ~105 yr) due to extremely high EUV and FUV fields (Gorti & Hollenbach 2009). Indeed, only very feeble disk remnants have been observed around O stars, while disks around B stars survive for a few hundred thousand years, leaving not enough time for CA to operate (Gorti et al. 2009).

In this regard, it is interesting to qualitatively compare the disk survival timescales of the binary b Cen (MA ~ 5.5 M⊙, MB ~ 3. 5 M⊙) and the single 2 Sco ( M ~ 9.1 M⊙). While a naive comparison between b Cen A – expected to emit nearly all the X and UV flux in the b Cen system and treated thus as a single star – and µ2 Sco would yield a three times shorter disk lifetime around µ2 Sco (see Eq. (7) by Gorti & Hollenbach 2009), the actual ratio should be much larger due to an initial disk mass for b Cen related to the total system mass, hence comparable to that of μ2 Sco. Combining Eq. (11) from Cesaroni et al. (2007) for the photoevaporation outflow rate with the expected ionizing photon flux Φ for the three stars, the ratio17 between the two disk survival timescales should be around ~20. While the impact of this on planet formation is difficult to be properly assessed, the presence of (at least) a companion with q ~ 0.002 around µ² Sco looks much more challenging than that around b Cen in the framework of a CA scenario. Recent updates of the classical CA model, such as pebble accretion, have been indicated as a possible solution to the conundrum (Lambrechts & Johansen 2012)¹⁸.

On the other hand, giant planets may possibly form inside massive protoplanetary disks on a much shorter timescale of ~10⁴ yr by means of GI. While this mechanism is not likely to be at the origin of most planets observed around solar-type stars, it might be considered for stars much more massive than the Sun; there is observational evidence of disk fragmentation around 10 M_⊙ stars (Cadman et al. 2020; Suri et al. 2021). GI preferentially produces massive planets in wide orbits, although rapid migration can force some of them to move much closer to the star (Malik et al. 2015). In this regard, it is interesting to notice that our orbital analysis shows a preference for large eccentricities for both the confirmed µ² Sco b and the probable CC0. A recent analysis of 27 directly-imaged giant planets and brown dwarfs in the separation range 5–100 au around a wide range of stellar hosts (0.2–2.8 M_⊙) has hinted at an interesting eccentricity dichotomy between the two populations: while the former usually has low eccentricity values , the latter is characterized by a flat distribution over the range 0 < e < 1 (Bowler et al. 2020). If this dichotomy continues at larger separations and stellar masses, it might favor an in-situ GI scenario for µ² Sco b and CC0 (see Appendix D); alternatively, the large eccentricities might be simply a result of a strong dynamical evolution after formation, possibly causing the migration of µ² Sco b and CC0 into their current orbits (bringing a CA scenario back into the game; see, e.g., Marleau et al. 2019). We stress that any definitive conclusion on this point cannot be reached without a follow-up of the system in the next few years, aimed at both constraining the orbit of b and establishing if CC0 exists.

According to the parametric model described above, GI should be considered as one of the top-down processes underlying the SP scenario, preferentially producing massive BD and low-mass stars rather than giant planets (Kratter et al. 2010); indeed, extensive population synthesis simulations, that take into account migration and N-body interactions, find that 90–95% of the surviving objects have masses above the DBL already around stars with 0.8–1.2 M_⊙ (Forgan & Rice 2013; Forgan et al. 2018). The corresponding typical mass ratio q = 0.01–0.1 is one or two orders of magnitude larger than that of b Cen b and µ² Sco b and CC0, implying an unusually low conversion efficiency of disk mass into companion mass for BEAST companions. Finally, although the completion of the survey is needed before drawing robust conclusions, the presence of already two or three companions near the deuterium-burning limit is not expected according to the mass distribution of objects generated by current GI models.

In this regard, recent 3D radiative magneto-hydrodynamic simulations have shown that, in a GI disk, the magnetic field acts to reduce the fragmentation scale by 1–2 orders of magnitude, possibly turning GI into a viable formation path for Jupiter-analogs around B stars (Boley et al. 2010; Deng et al. 2021).

Irrespective of the formation scenario invoked, theoretical efforts toward the modeling of the exoplanet population have often been restricted to solar-like stellar hosts. First steps are being taken to extend toward lower masses (as a result of the attention drawn by M stars as potential hosts of habitable planets), but very little effort is going toward larger masses. This is mostly due to the lack of detections and the challenges of developing models without observational constraints. Upcoming BEAST observations will provide these crucial constraints, both by delivering a new population of companions and through nondetections, paving the way to a new suite of planet formation models around intermediate and massive stars.

The µ² Sco system presented in this work appears remarkable in many regards: taken in isolation, it is the most massive star hosting a planetary system, pushing for the first time the frontier of exoplanetary studies well into the massive star regime; coupled with the b Cen system, it hints at an emerging companion population resembling a scaled-up version of the known giant planet population found around Sun-like stars; a population that, intriguingly, is not yet understood by the known mechanism of planet formation, urging for the development of new models. Follow-up of this unique system will be required to confirm the existence of CC0, and to better constrain the orbital configuration of the system, which in turn might give crucial clues on its formation history.

Fig. 11 Posterior distributions of orbital parameters for the robust companion µ2 Sco b. Corner plot showing the posterior distributions from orbitize! for semi-major axis (a), eccentricity (e) and inclination (i). The 16th (left dashed line), 50th (middle dashed line) and 84th (right dashed line) percentiles of each parameter distribution are indicated, and the 50th percentile (median) value is listed with 1→ uncertainties derived from the lower and upper percentiles. One percent of all chains, representing long-tail values, have been excluded from the corner plots for clarity, but are still considered for the percentile calculations. No priors or constraints have been given for any of the parameters.

Fig. 12 Subset of possible orbits for the probable CC0. (A) Sky-projected orbital fit (orange star icons) based on the 2018 and 2021 epochs, showing 100 randomly drawn orbits from the orbitize! MCMC chains. The position of →2 Sco is indicated by the yellow star. (B) Evolution over time of the planet-star separation (→); measured points and the corresponding error bars are shown in purple. (C) Evolution over time of the position angle (PA); again, measured points and the corresponding error bars are shown in purple.

Fig. 13 Posterior distributions of orbital parameters for the probable CC0. Corner plot showing the posterior distributions from orbitize! for semi-major axis (a), eccentricity (e) and inclination (i) derived in the same way as in Fig. 11.

Fig. 14 Mass ratio vs. irradiation for known exoplanets. Only exoplanets whose stellar host mass is known to a precision of at least 30% are shown. Each planet is labeled according to its detection method: transits in green, radial velocity in red, microlensing in orange and direct imaging in blue. BEAST discoveries are overplotted with larger blue symbols, and circular orbits with radius equal to the observed projected distance are assumed; Solar System planets (images from NASA) are also shown for reference. µ2 Sco b can be considered a Jupiter analog both in terms of irradiation and mass ratio (similar to the directly imaged 51 Eridani b; Macintosh et al. 2015, obscured here by the icon for Jupiter) while the irradiation received by probable CC0 is similar to Mercury’s. Sources: NASA Exoplanet Archive (https://exoplanetarchive.ipac.caltech.edu/), the Extrasolar Planets Encyclopaedia (http://exoplanet.eu/).

8 Conclusions

The BEAST program is devoted to the search of substellar – possibly planetary – companions at separations of tens to hundreds au from B stars (M ≳ 2.4 M_⊙) in the young (5–20 Myr old) nearby association Scorpius-Centaurus. The final aim of the program is to understand if the actual dearth of planets known around massive stars is due to an early dispersal of the protoplanetary disk or simply to a selection effect – the planets being too far from the star to be discovered through RV searches. After the discovery of a planet companion to b Cen (Janson et al. 2021a), in this paper we have presented the intriguing case of µ² Sco. After a careful re-evaluation of the kinematic properties of this star, we found it to belong to a small group of stars within Upper Centaurus-Lupus (UCL) that we have labeled Eastern Lower Scorpius (ELS). The star – a slowly pulsating B-star β Cep hybrid – appears to be a slow rotator, and shows no evidence of a strong magnetic field; the combination of the two factors with a new distance estimate and suitable priors derived from the literature allowed us to determine the properties of the star. With an age of 20 ± 4 Myr and a mass of 9.1 ± 0.3 M_⊙, µ² Sco appears to be a massive main-sequence star that will likely end its life as an electron-capture supernova.

Thanks to high-contrast imaging observations with SPHERE@VLT, we have detected a comoving source in IRDIS images with measured contrasts ΔK₁ ≈ 11.7 mag and ΔK₂ ≈ 11.3 mag. We examined the possibility that this object is a background interloper or small mass member of the Sco-Cen association unrelated to µ² Sco and found both hypothesis to be highly improbable. We have therefore concluded that this object – µ² Sco b – is a real physical substellar companion to the star. With a best-fit M = 14.4 ± 0.8 M_J, µ² Sco b is the first “planet-like”, in the sense defined in Sect. 7, companion to a supernova progenitor.

A second 18.5 ± 1.5 M_J substellar companion, CC0, is tentatively detected in IFS images, much closer to the star, but we highlight the need for follow-up observations to definitely confirm its existence. The object would be the most irradiated substellar companion ever discovered through direct imaging.

CC0 and µ² Sco b have angular separations of 0.128 ± 0.002 arcsec and 1.709 ± 0.005 arcsec from the star that translate, at a distance d = 169 ± 6 pc, into projected separations of 21 ± 1 au and 290 ± 10 au, respectively. The small astrometric displacements over the three years separating the two epochs (Fig. 8) are a consequence of their orbital motion around µ² Sco: even though our observations are not enough to completely determine the orbits, we could constrain the set of orbital parameters that are consistent with observations. A preference for large inclinations and/or large eccentricities emerges; owing to the large b-c distance, the two-companion system should nonetheless be dynamically stable (Figs. l0–13).

µ² Sco b (and possibly CC0) are naturally added to b Cen b, a 10.9 ± 1.6M_J companion recently discovered around the binary b Cen system in the same BEAST survey (Janson et al. 2021a). While the survey is still ongoing and a sound statistical analysis is not yet possible, the discovery of already two (possibly three) companions with very small mass ratios – similar to Jupiter’s – around the first 25 stars for which observations have been completed, dwelling at separations within the expected size of the protoplanetary disks, appears to be more compatible with a bottom-up “planet-like” formation scenario rather than a top-down “binary-like” scenario. The possible scaling of planetary system sizes with stellar mass would be in tension with the current definition of planets adopted by the IAU, which would classify the companion(s) of µ² Sco as brown dwarfs and b Cen b as a giant planet.

The completion of the BEAST survey will allow a more definite answer to this issue, providing a census of companions around B stars over a wide range of masses. In addition, new observations of µ² Sco with high-contrast imaging and near-IR interferometry are needed both to confirm the existence of CC0 and for a better characterization of the properties of µ² Sco b; coupled with asteroseismologic follow-up of the primary star, they will enable us to shed light on the fascinating history of this unique system.

Appendix A Optimization Tests for Stellar Parameters

Our optimization tests for the parameters of µ² Sco start from the prior distributions of (M, R, T_eff, π) described in Section 2. We create a synthetic sample of 10⁷ stars, each one described by a quadruplet (M, R, T_eff, π), where every parameter is randomly drawn from its prior distribution. The empirical table of intrinsic colors and temperatures of 5-30 Myr old stars by Pecaut & Mamajek (2013) allows a straightforward conversion of bolometric luminosities into V magnitudes by means of calibrated bolometric corrections. Starting from the stellar luminosity expected for main sequence stars: (A.1)

we apply the interpolated bolometric correction as a function of T_eff, the fixed extinction and the derived distance modulus to get synthetic apparent V magnitudes. A filter was then applied to select only the quadruplets that simultaneously satisfy the three conditions:

• V ∈ [3.52,3.60], where the broadened error bar allows for some uncertainty in the bolometric correction itself; – log g ∊ [3.6,4.0];

• a mass M such that |M – M_L\/M < 0.15, where M_L is the mass derived through a L(M) relation valid for main sequence stars with M > 2 M_⊙.

The L(M) relation (in solar units) is derived, for consistency, by fitting L as a function of M using the same tables:

Although the mean fractional error of the fitted points is about 4%, we opted for a less stringent 15% tolerance to account for possible deviations from the standard behavior of the underlying sample. We underline that this test does not allow to constrain the age of the star, since the tables are built by averaging over stars of different ages.

The posterior distributions of (M, R, T_eff, π) are shown in Figure A.1. While the mean values of effective temperature and distance do not change significantly, a strong preference for a large value of radius appears (R = 5.8 ± 0.3 R_⊙), while the mass distribution shifts to very high values .

Considering that older stars on the main sequence are brighter, that is, age and mass are anti-correlated for MS stars of the same magnitudes, we expect our mass estimate to be corrected toward higher or lower values depending on the adopted age. In order to introduce the age into the discussion, we ran a similar but independent test that only relies on PARSEC isochrones (Marigo et al. 2017) of varying age and metallicity; drawing 10⁷ masses and parallaxes as above, we selected the combinations that simultaneously satisfied:

• apparent Gaia magnitudes (G, G_bp, G_rp) within 0.1 mag from their observed values; absorption coefficients were taken from Wang & Chen (2019);

• a derived log(g) ∈ [3.6,4.0] as before;

• T_eff ∈ [20900,22900] K.

The test could never be passed using the nominal age (t = 16 Myr) and metallicity ([Fe/H]=0.00), but increasing at least one of the two constraints (t ∈ [16,25] Myr or [Fe/H] ∈ [0.10,0.20]) a family of solutions appeared (Figure A.2–A.3); the posterior mass distribution is consistent with the result of the first test (M = 8.8 ± 0.3 M_⊙ for the run at constant metallicity, M = 9.2 ± 0.2 M_⊙ for the run at constant age; we may combine the two to get M = 9.0 ± 0.3 M_⊙). We average the two mass determinations from the independent methods to a derive a final M = 9.1 ± 0.3 M_⊙. In the same way, we derive R = 5.6 ± 0.2 R_⊙; the median values of the posteriors of the first test will be used as best-fit estimates for T_eff and π. Finally, while we are not able to solve the age-metallicity degeneracy, we encapsulate the general tendency for a larger age in a revised age estimate t = 20 ± 4 Myr.

Fig. A.1 Results of the first optimization test of stellar parameters. Corner plot showing the posterior distribution of the quadruplets (M, R, Teff, π) consistent with the observational filter for the optimization test based on Pecaut & Mamajek (2013).

Looking at Figure B.1, we can reasonably exclude the presence of a close stellar companion with q > 0.08; additionally, we do not have any evidence to support the existence of a smaller stellar companion. We notice that the parallax error in the latest Gaia release is in line with that expected in Gaia EDR3 for a source with G = 3.5 (Lindegren et al. 2021): the high errors associated with Gaia’s astrometric solution are likely a mere consequence of the extreme brightness of the star.

Fig. A.2 Results of the second optimization test of stellar parameters (constant metallicity). Corner plot showing the posterior distribution of the doublets (M, π) consistent with the observational filter for the optimization test based on PARSEC isochrones at fixed [Fe/H]=0.00. The corresponding Teff, R and age distribution are shown too.

Fig. A.3 Results of the second optimization test of stellar parameters (constant age). Corner plot showing the posterior distribution of the doublets (Μ,π) consistent with the observational filter for the optimization test based on PARSEC isochrones at fixed t = 16 Myr. The corresponding Teff, R and [Fe/H] distribution are shown too.

Appendix B: Possible binarity of µ² Sco

The astrometric solution of µ² Sco in the literature appears highly problematic (Table B.1). As shown by Kervella et al. (2019), a significant deviation of proper motion components from the long-term motion reconstructed by cross-matching Hip-parcos and Gaia can be suggestive of the gravitational effects from an unseen stellar companion. Given the importance of the system architecture for a correct characterization and interpretation of our results, we tried to put constraints on possible unresolved stellar companions.

Fig. B.1 Limits on the mass of a possible unresolved stellar companion coming from several techniques. A M ≥ 0.3 - 0.5 M⊙ companion is excluded at sub-au separations by RVs (95% confidence interval), a M > 0.7 M⊙ is excluded at 1-5 au by the lack of X-ray detection. Astrometry (PMA) and coronagraphic observations place even stronger limits at d > 5 au. At ~ 290 au, μ2 Sco b lies outside the x range.

Although Kervella et al. (2019) reports a significant

PMA was not significant at Gaia DR2 era. Thus, we checked the consistency of the astrometric solution from Hipparcos, finding that the 2007 reduction yields a PMA along right ascension with an opposite sign and similar magnitude (about 0.7 mas yr⁻¹) compared to the 1997 reduction based on the same raw data (Perryman et al. 1997). Therefore, the tabulated PMA of Hipparcos cannot be trusted for this star. We use the nonsignificant Gaia DR2 PMA as a separation-dependent upper limit on the stellar companion mass.

µ² Sco is neither a visual binary (Gaia Collaboration 2021) nor an interferometric binary (Rizzuto et al. 2013); additionally, it does not appear as an X-ray source in ROSAT (Berghoefer et al. 1996). The upper limit on the X-ray luminosity, given in units of L_X = log f_x [erg s⁻¹], is L_X = 29.82. Since the estimate assumes a distance d ~ 206 pc, recalibration with the distance adopted in this work yields L_X = 29.65. Such an emission is below the predicted plateau of a K5 star (Poppenhaeger & Wolk 2014), posing an upper limit, independent of separation, of M₂ < 0.7 M_⊙.

We then analyzed HARPS (Mayor et al. 2003) and FEROS (Kaufer et al. 1999) radial velocity data, spanning about 10 yearsretrieved from the ESO Science Archive¹⁹: after extraction of RV from 37 He I and atomic lines, we verified that the rms of the observations (0.21 km s⁻¹) is smaller than individual uncertainties (~ 0.3 km s⁻¹). Spectra taken at short temporal separation can push the sensitivity further: the four HARPS observations provide a small scatter of 74 m/s, and the two FEROS observations acquired on JD=56523.66 only differ by 90 m/s. The overall shallow trend of ~ 50 m s⁻¹ yr⁻¹ can be used to derive an additional upper mass limit: we used the Exo-DMC code (Bonavita 2020) under default assumptions to derive a 95% confidence interval for this mass limit.

Appendix C: Stellar magnetic fields and rotation

As the presence of a strong magnetic field can hinder a reliable comparison with model isochrones (see Section 2), we analyzed three high-resolution spectra of µ² Sco obtained as part of the Magnetism in Massive Stars (MiMeS) project (Wade et al. 2016): two taken by the spectropolarimeter ESPaDOnS (Donati et al. 2006) in 2014 and 2015, and one taken by spectropolarimeter HARPSpol (Piskunov et al. 2011) in 2011. To begin with, comparison of the intensity spectrum (Stokes I) with a non-LTE TLUSTY synthetic model (Hubeny & Lanz 2011) yields values for the stellar parameters that are fully consistent with those derived through our optimization tests. We analyzed both Stokes I and Stokes V (circular polarization) spectra applying the Least-Squares Deconvolution technique (LSD) to perform a sort of weighted mean in all the spectral lines (Donati et al. 1997). This method provides mean photospheric Stokes I and V profiles with a S/N much better than in the individual lines, and allows therefore to put stringent constraints on the surface magnetic field.

The weights used in LSD are the predicted central depth of the intensity lines, the wavelength and the Landé factor. The mask, that is the list of predicted lines, has been obtained by using the VALD atomic line database20, removing only some usually very strong lines (Balmer and He I) affected by non-LTE effects.

The resulting LSD profiles are shown in Figure C.1. The bottom profile shows the mean Stokes I profile, typical of a photospheric profile of a massive star, with broad wings indicating a significant macroturbulent velocity. The profile is clearly disturbed by the β Cephei pulsations. The top profile is Stokes V, and the middle one is the null profile N, computed in a way that allows us to check that no spurious polarization is present in our data. All the curves are normalized to the mean continuum intensity, Ic.

The Stokes V profile is totally flat, indicating that there is no Zeeman detection in the data. All the measurements of the line-of-sight component of the magnetic field averaged over the surface of the star, called the longitudinal magnetic field Bl, are consistent with 0 G with an uncertainty of ~ 15G. Using a Monte Carlo simulation (Alecian et al. 2016), we estimate that if a dipole field at the pole with B ≥ 170 G exists, we would have had a 90% chance to have detected it. As the typical magnetic fields of massive stars are between a few hundreds of Gauss to few kG, μ2 Sco is not likely to host a strong fossil field.

Fig. C.1 Resulting LSD profiles for the three μ2 Sco spectra (black: ESPaDOnS 2014; red: ESPaDOnS 2015; green: HARPSpol 2011). Top profile: Stokes V. Middle profile: null profile. Bottom profile: mean Stokes I profile. For each observation, the typical errorbars of individual spectral points are shown on the right. A vertical offset was applied to the top and the middle curve for better visibility.

With respect to rotational velocity, the low observed v sin i = 52 km s-1 implies a geometric probability of having v > 100 km s-1 of about 13%. By comparison, the median v sin i for B stars belonging to ~ 20 Myr regions has been estimated as ~ 125 km s 1 (Wolff et al. 2007). μ2 Sco appears therefore to be a slow rotator.

Appendix D: Formation and dynamical analysis of the system under a GI scenario

Although the large eccentricities of both the robust µ² Sco b and the probable CC0 can be suggestive of a strong dynamical evolution after their birth, it is nevertheless instructive to see whether an in-situ formation at the position of CC0 is possible according to current formation models. With a best-fit semimajor axis of about 20 au, the mean irradiation level of CC0 is comparable to that of Mercury, preventing the presence of ice grains which are fundamental in a CA scenario (Mordasini et al. 2012). Whether the same argument applies to GI is not clear, crucially depending on the values of Toomre’s Q parameter and the cooling time, which determine if the protoplanetary disk can fragment or not at this position. Fragmentation occurs if Toomre’s Q parameter, defined as: (D.1)

is less than unity. Here Ω* is the epicyclic frequency, cs is the sound speed, G is the gravitational constant and Σ is the surface density. Since the radial dependencies of Ω* and Σ are similar (α r-3/2), the radial variation of Q depends on that of cs, and the treatment can be simplified. Recalling that a full modeling of GI around stars as massive as μ2 Sco is still lacking, we might start from simple scaling laws from the results obtained for solar hosts. The relevant proportionalities in Q are given by: . As regards the disk mass, we conservatively assume a linear proportionality between Mdisk and Mstar. For the gas temperature, we might start from our best-fit estimates for stellar radius and Teff to derive the effect due to stellar irradiation at 20 au (Eq. 4, 96) that is Tirr = 690 ± 50 K, supposing a flared disk in vertical hydrostatic equilibrium. The onset of fragmentation for 0.1 M⊙ disks around 1M© stars in theoretical simulations occurs for T ⊙ 50 K (Helled et al. 2014), so we may write: (D.2)

Since fragmentation happens if Qc < 1, a formation from GI would not be possible at the location of CC0. However, if Mdisk/Mstar increases by more than 20% than expected from the linear proportionality in the mass range of B stars, as assumed above, in situ fragmentation at the position of CC0 becomes possible. Indeed, indications exist that - at least for M < 2M⊙ – Mdisk α Mq with q = 1.3 – 2.0 (Pascucci et al. 2016): assuming that q > 1.1 for B stars, the value of Qc would become lower than 1. Also, a warmer disk reduces the cooling time, favoring the onset of fragmentation at fixed Q. On the other hand, given the much lower irradiation temperature at the distance of μ2 Sco b, μ2 Sco b could have easily formed in situ via GI: not only the Toomre parameter is likely much smaller than unity in the approximations considered above; but also, more importantly, the conditions in protoplanetary disks at R > 50 au are most favorable to fragmentation since gas, even in a massive disk, is expected to be optically thin, thus leading to short cooling timescales (Boley et al. 2010).

An interesting consequence of this scenario would be that disks around stars similar to μ2 Sco might be unstable at some phase of their evolution over a wide range of separations. However, the formation of planet-like objects such as the companions of b Cen and μ2 Sco rather than more massive (even stellar) objects would imply an extremely low efficiency of companion accretion, with the rest either accumulated on other objects (the star itself or a stellar companion) or lost from the system.

Appendix E: Companion candidates

List of Objects

‘mu2 Scorpii’ on page 1

‘b Centauri’ on page 1

Acknowledgments

Acknowledgements.

We are extremely grateful for the extremely helpful comments during the refereeing process, which significantly helped raise the quality of this paper. The observations have been acquired at the ESO VLT telescope (program 1101.C-0258). This work has made use of the SPHERE Data Center, jointly operated by OSUG/IPAG (Grenoble), PYTHEAS/LAM/CeSAM (Marseille), OCA/Lagrange (Nice), Observatoire de Lyon and Observatoire de Paris/LESIA (Paris). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. 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 work was supported by the Programme National de Planétologie (PNP) and ASHRA of CNRS/INSU, co-funded by CNES. V.S., R.G., D.M., S.D. and V.D. acknowledge the support of PRIN-INAF 2019 Planetary Systems At Early Ages (PLATEA). Part of this work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation. N.E., S.P.Q. and L.M. acknowledge the financial support of the SNSF. G.D.M. acknowledges the support of the Deutsche Forschungsgemeinschaft (DFG) priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (MA 9185/1-1). G.D.M. also acknowledges the support from the Swiss National Science Foundation under grant BSSGI0_155816 “PlanetsInTime”. A.M. acknowledges support from the ERC Consolidator Grant funding scheme (project ASTEROCHRONOMETRY, https://www.asterochronometry.eu, G.A. n. 772293). For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

References

All Tables

All Figures

	Fig. 1 Proper motion components of LS (blue) and ELS (orange) members. Both →1 Sco and →2 Sco are consistent with membership to ELS.
In the text

	Fig. 2 Sky coordinates of bona fide Sco-Cen members from Damiani et al. (2019), shown as black dots. Classical boundaries for UCL are indicated by the red box. A zoom of LS (blue) and ELS (orange) region is drawn in the upper-right corner; →2 Sco is indicated by the black star.
In the text

	Fig. 3 Parallax distribution of ELS stars. We adopt → = 5.9 ± 0.2 mas as parallax prior for →2 Sco.
In the text

	Fig. 4 IRDIS images of →2 Sco obtained with TLOCI (left panel: first epoch; right panel: second epoch). The star, artificially obscured by the coronagraphic mask, is at the center of the image. Several background sources can be easily seen as bright point sources. →2 Sco b is the source inside the white circle.
In the text

Fig. 5 IFS images of →2 Sco obtained with ASDI-PCA with 25 modes (left panel: first epoch; right panel: second epoch). The star is at the center of the image but not visible due to the presence of the coronagraphic mask – which artificially masks the innermost ~90 mas of the image – and the aggressive data analysis method used. The probable companion CC0 is the source inside the white circle. A background source is visible on the lower left of the first-epoch image.

In the text

Fig. 6 Position of CC0 at different wavelengths. (A) J-band image, first epoch; (B) H-band image, first epoch; (C) J-band image, second epoch; (B) H-band image, second epoch. The Y-band image is not shown due to the extreme faintness of the source. White circles indicate the best-fit positions from Table E.1; circle radii have been enlarged to ~10 times the uncertainties on best-fit positions to enhance visibility. Unlike a speckle (red circle), the separation of CC0 does not increase with wavelength.

In the text

	Fig. 7 (K1 → K2, K) color magnitude diagram of known substellar companions and field objects. IRDIS CCs only seen in one epoch are shown as red dots. The position of CC0 ( = →2 Sco b) is consistent with that expected for a young object.
In the text

Fig. 8 Astrometric motion of CCs present in both epochs. As the plot shows the difference in separation between the two epochs, a comoving source should be close to the origin (labeled by “2018”). A background source with null proper motion will move according to the dashed curve as a reflection of the target’s motion, ending in the blue circle labeled as ‘2021’ at the second epoch. The motion of CC0 (labeled by '0') and CC2 (labeled by ‘2’) is distinct from the cloud of background contaminants.

In the text

Fig. 9 Contrast spectrum for CC0 obtained with PACO-ASDI from 2018 IFS data (black dots) and from PACO-ASDI for IRDIS data acquired at the same epoch (5→ upper limits; red squares). The second-epoch ASDI-PCA spectrum (gray squares) is shown for reference. The solid line is the contrast spectrum expected for a 20 Myr old, 19.5 MJ brown dwarf using AMES-Dusty models. The dot-dashed line is a black body curve with the same temperature and radius as this model.

In the text

Fig. 10 Subset of possible orbits for the robust companion candidate →2 Sco b. (A) Sky-projected orbital fit (orange star icon) based on the 2018 and 2021 epochs which overlap at this scale, showing 100 randomly drawn orbits from the orbitize! MCMC chains. The position of →2 Sco is indicated by the yellow star. (B) Evolution over time of the planet-star separation (→); measured points and the corresponding error bars are shown in purple. (C) Evolution over time of the position angle (PA); again, measured points and the corresponding error bars are shown in purple.

In the text

Fig. 11 Posterior distributions of orbital parameters for the robust companion µ2 Sco b. Corner plot showing the posterior distributions from orbitize! for semi-major axis (a), eccentricity (e) and inclination (i). The 16th (left dashed line), 50th (middle dashed line) and 84th (right dashed line) percentiles of each parameter distribution are indicated, and the 50th percentile (median) value is listed with 1→ uncertainties derived from the lower and upper percentiles. One percent of all chains, representing long-tail values, have been excluded from the corner plots for clarity, but are still considered for the percentile calculations. No priors or constraints have been given for any of the parameters.

In the text

Fig. 12 Subset of possible orbits for the probable CC0. (A) Sky-projected orbital fit (orange star icons) based on the 2018 and 2021 epochs, showing 100 randomly drawn orbits from the orbitize! MCMC chains. The position of →2 Sco is indicated by the yellow star. (B) Evolution over time of the planet-star separation (→); measured points and the corresponding error bars are shown in purple. (C) Evolution over time of the position angle (PA); again, measured points and the corresponding error bars are shown in purple.

In the text

	Fig. 13 Posterior distributions of orbital parameters for the probable CC0. Corner plot showing the posterior distributions from orbitize! for semi-major axis (a), eccentricity (e) and inclination (i) derived in the same way as in Fig. 11.
In the text

Fig. 14 Mass ratio vs. irradiation for known exoplanets. Only exoplanets whose stellar host mass is known to a precision of at least 30% are shown. Each planet is labeled according to its detection method: transits in green, radial velocity in red, microlensing in orange and direct imaging in blue. BEAST discoveries are overplotted with larger blue symbols, and circular orbits with radius equal to the observed projected distance are assumed; Solar System planets (images from NASA) are also shown for reference. µ2 Sco b can be considered a Jupiter analog both in terms of irradiation and mass ratio (similar to the directly imaged 51 Eridani b; Macintosh et al. 2015, obscured here by the icon for Jupiter) while the irradiation received by probable CC0 is similar to Mercury’s. Sources: NASA Exoplanet Archive (https://exoplanetarchive.ipac.caltech.edu/), the Extrasolar Planets Encyclopaedia (http://exoplanet.eu/).

In the text

	Fig. A.1 Results of the first optimization test of stellar parameters. Corner plot showing the posterior distribution of the quadruplets (M, R, Teff, π) consistent with the observational filter for the optimization test based on Pecaut & Mamajek (2013).
In the text

	Fig. A.2 Results of the second optimization test of stellar parameters (constant metallicity). Corner plot showing the posterior distribution of the doublets (M, π) consistent with the observational filter for the optimization test based on PARSEC isochrones at fixed [Fe/H]=0.00. The corresponding Teff, R and age distribution are shown too.
In the text

	Fig. A.3 Results of the second optimization test of stellar parameters (constant age). Corner plot showing the posterior distribution of the doublets (Μ,π) consistent with the observational filter for the optimization test based on PARSEC isochrones at fixed t = 16 Myr. The corresponding Teff, R and [Fe/H] distribution are shown too.
In the text

Fig. B.1 Limits on the mass of a possible unresolved stellar companion coming from several techniques. A M ≥ 0.3 - 0.5 M⊙ companion is excluded at sub-au separations by RVs (95% confidence interval), a M > 0.7 M⊙ is excluded at 1-5 au by the lack of X-ray detection. Astrometry (PMA) and coronagraphic observations place even stronger limits at d > 5 au. At ~ 290 au, μ2 Sco b lies outside the x range.

In the text

Fig. C.1 Resulting LSD profiles for the three μ2 Sco spectra (black: ESPaDOnS 2014; red: ESPaDOnS 2015; green: HARPSpol 2011). Top profile: Stokes V. Middle profile: null profile. Bottom profile: mean Stokes I profile. For each observation, the typical errorbars of individual spectral points are shown on the right. A vertical offset was applied to the top and the middle curve for better visibility.

In the text