11institutetext: Department of Physics and Astronomy ‘Galileo Galilei’,
University of Padova; Vicolo dell’Osservatorio 3, I-35122 Padova, Italy
11email<EMAIL_ADDRESS>22institutetext: INAF – Osservatorio
Astronomico di Padova; Vicolo dell’Osservatorio 5, I-35122 Padova, Italy
33institutetext: Institutionen för astronomi, Stockholms universitet; AlbaNova
universitetscentrum, SE-106 91 Stockholm, Sweden 44institutetext: Jet
Propulsion Laboratory, California Institute of Technology; 4800 Oak Grove
Drive, Pasadena CA 91109, USA. 55institutetext: Univ. Grenoble Alpes, CNRS,
IPAG; F-38000 Grenoble, France 66institutetext: Unidad Mixta Internacional
Franco-Chilena de Astronomía, CNRS/INSU UMI 3386 and Departamento de
Astronomía, Universidad de Chile; Casilla 36-D, Santiago, Chile
77institutetext: CRAL, UMR 5574, CNRS, Université Lyon 1; 9 avenue Charles
André, F-69561 Saint-Genis-Laval Cedex, France 88institutetext: Aix-Marseille
Université, CNRS, CNES, LAM (Laboratoire d’Astrophysique de Marseille); UMR
7326, F-13388 Marseille, France. 99institutetext: Departamento de Física
Teórica, Universidad Autónoma de Madrid; S-28049 Madrid, Spain
1010institutetext: Centro de Investigación Avanzada en Física Fundamental
(CIAFF), Facultad de Ciencias, UAM; S-28049 Madrid, Spain 1111institutetext:
Landessternwarte, Zentrum für Astronomie der Universität Heidelberg;
Königstuhl 12, D-69117 Heidelberg, Germany 1212institutetext: Leiden
Observatory, Leiden University; Postbus 9513, NL-2300 RA Leiden, The
Netherlands 1313institutetext: Department of Astronomy, University of
Michigan; 1085 S. University Ave, Ann Arbor MI 48109, USA 1414institutetext:
School of Physical Sciences, Faculty of Science, Technology, Engineering and
Mathematics, The Open University; Walton Hall, Milton Keynes, MK7 6AA, UK
1515institutetext: ETH Zürich, Institute for Particle Physics and
Astrophysics; Wolfgang-Pauli-Strasse 27, CH-8093 Zürich, Switzerland
1616institutetext: School of Physics and Astronomy, University of Birmingham;
Edgbaston, B15 2TT, UK 1717institutetext: Max-Planck-Institut für Astronomie,
Königstuhl 17, D-69117 Heidelberg, Germany 1818institutetext: Center for
Theoretical Physics and Cosmology, Institute for Computational Science,
Universität Zürich; Winterthurerstrasse 190, CH-8056 Zürich, Switzerland
1919institutetext: LESIA, Observatoire de Paris, Université PSL, CNRS,
Sorbonne Université, Univ. Paris Diderot; Sorbonne Paris Cité, 5 Place Jules
Janssen, F-92195 Meudon, France 2020institutetext: Institut für Astronomie und
Astrophysik, Universität Tübingen; Auf der Morgenstelle 10, D-72076 Tübingen,
Germany 2121institutetext: Physikalisches Institut, Universität Bern;
Gesellschaftsstr. 6, CH-3012 Bern, Switzerland

# A scaled-up planetary system around a supernova progenitor††thanks: Based on
observation made with European Southern Observatory (ESO) telescopes at
Paranal Observatory in Chile, under program 1101.C0258

V. Squicciarini 1122 R. Gratton 22 M. Janson 33 E. E. Mamajek 44 G. Chauvin
5566 P. Delorme 55 M. Langlois 77 A. Vigan 88 S. C. Ringqvist 33 G. Meeus
991010 S. Reffert 1111 M. Kenworthy 1212 M. R. Meyer 1313 M. Bonnefoy 77 M.
Bonavita 1414 D. Mesa 22 M. Samland 33 S. Desidera 22 V. D’Orazi 22 N. Engler
1515 E. Alecian 55 A. Miglio 1616 T. Henning 1717 S. P. Quanz 1515 L. Mayer
1818 O. Flasseur 1919 G.-D. Marleau 202021211717

(Received –; accepted –)

###### Abstract

Context.  Virtually all known exoplanets reside around stars with
$M<2.3\leavevmode\nobreak\ M_{\odot}$ either due to the rapid evaporation of
the protostellar disks or to selection effects impeding detections around more
massive stellar hosts.

Aims. To clarify if this dearth of planets is real or a selection effect, we
launched the planet-hunting B-star Exoplanet Abundance STudy (BEAST) survey
targeting B stars ($M>2.4\leavevmode\nobreak\ M_{\odot}$) in the young (5-20
Myr) Scorpius-Centaurus association by means of the high-contrast spectro-
imager SPHERE at the Very Large Telescope.

Methods. In this paper we present the analysis of high-contrast images of the
massive ($M\sim 9\leavevmode\nobreak\ M_{\odot}$) star $\mu^{2}$ Sco obtained
within BEAST. We carefully examined the properties of this star, combining
data from Gaia and from the literature, and used state-of-the-art algorithms
for the reduction and analysis of our observations.

Results. Based on kinematic information, we found that $\mu^{2}$ Sco is a
member of a small group which we label Eastern Lower Scorpius within the
Scorpius-Centaurus association. We were thus able to constrain its distance,
refining in turn the precision on stellar parameters. Around this star we
identify a robustly detected substellar companion ($14.4\pm 0.8$
$M_{\mathrm{J}}$) at a projected separation of $290\pm 10$ au, and a probable
second similar object ($18.5\pm 1.5$ $M_{\mathrm{J}}$) at $21\pm 1$ au. The
planet-to-star mass ratios of these objects are similar to that of Jupiter to
the Sun, and the flux they receive from the star is similar to those of
Jupiter and Mercury, respectively.

Conclusions. The robust and the probable companions of $\mu^{2}$ Sco are
naturally added to the giant $10.9\leavevmode\nobreak\ M_{\mathrm{J}}$ planet
recently discovered by BEAST around the binary b Cen system. While these
objects are slightly more massive than the deuterium burning limit, their
properties are similar to those of giant planets around less massive stars and
they are better reproduced by assuming that they formed under a planet-like,
rather than a star-like scenario. Irrespective of the (needed) confirmation of
the inner companion, $\mu^{2}$ Sco is the first star that would end its life
as a supernova that hosts such a system. The tentative high frequency of BEAST
discoveries is unexpected, and it shows that systems with giant planets or
small-mass brown dwarfs can form around B stars. When putting this finding in
the context of core accretion and gravitational instability formation
scenarios, we conclude that the current modeling of both mechanisms is not
able to produce this kind of companion. The completion of BEAST will pave the
way for the first time to an extension of these models to intermediate and
massive stars.

###### Key Words.:

planetary systems – Stars: early-type – Stars: individual: mu2 Sco – Stars:
individual: b Cen – Techniques: high angular resolution –

##  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_{\mathrm{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 years, 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 space111Which comprises, but is not limited to,
planetary (mass, radius), stellar (mass, radius, age, activity), and orbital
(orbital distance, period, eccentricity, inclination) properties. (Gaudi et
al., 2021). Most previous and ongoing surveys have focused on stars as massive
as (or less massive than) the Sun222More massive stars than the Sun have
mostly been targeted in later evolutionary phases, once they have left the
main sequence to become giant stars., 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 ($\lesssim 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_{\odot}$ and goes down to zero at $M>3M_{\odot}$ (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
($\sim 10^{4}$ 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\lesssim 2.4M_{\odot}$); 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\pm 1.6M_{\mathrm{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
$\mu^{2}$ Scorpii ($\mu^{2}$ Sco), composed of a robustly-detected $14.4\pm
0.8\leavevmode\nobreak\ M_{\mathrm{J}}$ companion and a probable closer
$18.5\pm 1.5\leavevmode\nobreak\ M_{\mathrm{J}}$ companion, and discuss its
implication for the current scenarios of planetary formation. We present an
extensive discussion of stellar properties in Section 2. Section 3 details the
observations and Section 4 discusses the data reduction methods. In Section 5
we identify the companion candidates, treating with special care the innermost
one. In Section 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 Section
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

$\mu^{2}$ Sco (also known as HR 6252, HD 151985, HIP 82545,
Pipirima333Pipirima is the Tahitian name for the pair $\mu^{2}$ and $\mu^{1}$
Sco, referring to mythological twin siblings. The IAU Working Group on Star
Names (WGSN) adopted Pipirima as the proper name for $\mu^{2}$ Sco in 2017
(https://www.iau.org/public/themes/naming_stars/).) is a naked-eye star
belonging to the young Scorpius-Centaurus association (de Zeeuw et al., 1999).
The main astrometric, 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 $\mu^{2}$ 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 Section 2.1-2.6, and the final
derived parameters are reported in Section 2.7. Complete details on the
derivation are provided in Appendix A.

Table 1: Main astrometric, kinematic and photometric properties of $\mu^{2}$ Sco collected from the literature. References: (1): Gaia Collaboration et al. (2021); (2): Kervella et al. (2019); (3): Gontcharov (2006); (4): Brown & Verschueren (1997); (5) Mermilliod (2006); (6): Cutri et al. (2003); (7) Hiltner et al. (1969). Name | Value | Reference
---|---|---
$\alpha$ (∘, J2016.0) | 253.083869820(64) | (1)
$\delta$ (∘, J2016.0) | -38.017636551(38) | (1)
$\mu_{\alpha}^{*}$ (mas yr-1) | $-11.772\pm 0.022$ | (2)
$\mu_{\delta}$ (mas yr-1) | $-23.105\pm 0.021$ | (2)
RV (km s-1) | $1.3\pm 0.8$ | (3)
$v\sin{i}$ (km s-1) | $52\pm 5$ | (4)
V (mag) | $3.565\pm 0.005$ | (5)
G (mag) | $3.543\pm 0.003$ | (1)
J (mag) | $4.15\pm 0.28$ | (6)
H (mag) | $4.159\pm 0.25$ | (6)
K (mag) | $4.292\pm 0.31$ | (6)
spectral type | B2IV | (7)

###  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. $\mu^{2}$ Sco, in
particular, appears to lie within a small clump of stars (a region centered on
galactic coordinates $(l,b)=(343.1^{\circ},4.7^{\circ})$ with a radius of
about $2^{\circ}$ 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 et al., 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 8). 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 (Figure 1).

Figure 1: Proper motion components of LS (blue) and ELS (orange) members. Both
$\mu^{1}$ Sco and $\mu^{2}$ Sco are consistent with membership to ELS.

We notice that our target (together with its sibling $\mu^{1}$ Sco) actually
lies in a peripheral area of LS, $\sim 2^{\circ}$ 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, Figure 2). ELS appears
somewhat closer to the Sun than the whole LS: we adopt its parallax
distribution (modeled as a normal, $\pi=5.9\pm 0.2$ mas) as our parallax prior
(Figure 3).

Figure 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; $\mu^{2}$ Sco is indicated by the black star. Figure 3: Parallax distribution of ELS stars. We adopt $\pi=5.9\pm 0.2$ mas as parallax prior for $\mu^{2}$ Sco. Table 2: Defining criteria for the UCL, LS and ELS samples used throughout this paper. The criteria should be interpreted as an addition to the cuts already applied in the definition of the initial catalog. Equatorial coordinates and proper motions refer to the ICRS at epoch J2016.0. Criterion | UCL | LS | ELS
---|---|---|---
initial catalog |  Damiani et al. (2019), |  Damiani et al. (2019), | Gaia EDR3
bona fide sources | bona fide sources
galactic longitude $(^{\circ})$ | $313<l<343$ | $339<l<350$ | $345<l<347.5$
galactic latitude $(^{\circ})$ | — | $1.5<b<7.7$ | $3.3<b<5.0$
right ascension $(^{\circ})$ | — | — | $251<\alpha<255$
declination $(^{\circ})$ | — | — | $-40<\alpha<-36$
parallax (mas) | — | $5.3<\pi<6.2$ | $5.3<\pi<6.5$
proper motion along $\alpha$ (mas yr-1) | — | — | $-18.7<\mu_{\alpha}^{*}<-10.2$
proper motion along $\delta$ (mas yr-1) | — | — | $-25.0<\mu_{\delta}<-16.1$
velocity along $\alpha$ (km s-1) | — | $-11.2<v_{\alpha}<-8.5$ | $-10<v_{\alpha}<-7$
velocity along $\delta$ (km s-1) | — | $-19.3<v_{\delta}<-16.5$ | $-21<v_{\alpha}<-16$
no. of sources | 3842 | 575 | 73

###  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\in[5,12]M_{\odot}$. Likewise, we do not favor a particular
initial value for the radius: $R\in[3.5,10.0]\leavevmode\nobreak\ R_{\odot}$.
For a B2V star on the main sequence we expect $R\sim 4$ R⊙ (Pecaut & Mamajek,
2013), but the classification of the star as a B2IV (Hiltner et al., 1969)
argues for a larger radius.

###  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
$\sim 10^{5}$ 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
preparation), we collected photometry from Gaia 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 ($\sim 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
($\alpha=0.05$)444The higher ($\sim 5$ and $\sim 4$ times, respectively)
density of LS and ELS with respect to the coeval UCL is likely due to
significant self-gravity, slowing down the expansion caused by random
turbulent motion within the natal molecular cloud. Indeed, the dispersion of
proper motions of ELS (0.71 mas yr-1 $\approx 0.55$ km s-1 is lower than the
escape velocity of $\sim 0.72$ km s-1 obtained by summing the masses of all
observed members ($\sim 50\leavevmode\nobreak\ M_{\odot}$) and considering its
projected size ($\rho\sim 17^{\prime}\sim 0.85$ pc). The lifetime of a small
cluster such as ELS against tidal disruption is about 35 Myr (Lamers et al.,
2005), about twice its estimated age..

In the absence of tighter constraints, we adopt the age of UCL ($t=16\pm 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 $\beta$ Cephei type variability and in
particular for that of a slowly pulsating B-star (SPB) $\beta$ 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.

###  2.4 Reddening

To estimate the reddening $E(B-V)$ in the direction of $\mu^{2}$ Sco in a
robust way, we considered six different determinations:

  * •

starting from $UBV$ photometry from Mermilliod (2006) and $uvby\beta$
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 $uvby\beta$ photometry through the calibration of Castelli
(1991): $E(B-V)=0.019$ mag;

  * •

dereddening the $uvby\beta$ photometry with the older, independent calibration
of Shobbrook (1983): $E(B-V)=0.020$ mag;

  * •

from Lyman $\alpha$ observations, Savage & Panek (1974) estimate the
interstellar column density of HI toward $\mu^{2}$ Sco to be $N(H)=2.5\cdot
10^{20}\text{ cm}^{-2}$. 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 $\mu^{2}$ 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\pm 0.024$ mag.

By averaging these estimates, we obtain $E(B-V)=0.022\pm 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\pm 0.015$ mag
and $A(K)=0.0062\pm 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, $\mu^{2}$ 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 $\delta$ Cet
($B-V=-0.219$, $U-B=-0.850$, $E(B-V)=0.018$) which has a median
$T_{\mathrm{eff}}$ in the recent literature of $T_{\mathrm{eff}}\simeq 21600$
K (Cardiel et al., 2021). Hence we do not expect the effective temperature of
$\mu^{2}$ Sco to differ too much from that of $\delta$ Cet.

Table 3: New effective temperature estimates for $\mu^{2}$ Sco. Details on the derivation of each estimate are provided in the text. # |  $T_{\mathrm{eff}}$ (K)
---|---
1 | 20913
2 | 21083
3 | 21655
4 | 21989
5 | 22063
6 | $22700\pm 270$
7 | 25978
adopted | $21900\pm 1000$

In Table 3 we list several new $T_{\mathrm{eff}}$ estimates for $\mu^{2}$ Sco
based on photometric data from the literature. The estimates are derived in
the following way:

  1. 1.

photometry from Hauck & Mermilliod (1997), dereddened through Castelli (1991)
using Balona (1994) $T_{\mathrm{eff}}$ calibration;

  2. 2.

photometry from Hauck & Mermilliod (1997), $T_{\mathrm{eff}}/H_{\beta}$
relation by Balona (1984);

  3. 3.

employing the $(U-B)_{o}$ vs $T_{\mathrm{eff}}$ trend based on B2IV standard
stars Pecaut & Mamajek (2013);

  4. 4.

photometry from Hauck & Mermilliod (1997), dereddened through Castelli (1991)
using Napiwotzki et al. (1993) $T_{\mathrm{eff}}$ calibration;

  5. 5.

photometry from Hauck & Mermilliod (1997), using [c1] index adopting Nieva
(2013) $T_{\mathrm{eff}}$ scale;

  6. 6.

comparing IUE spectrophotometry555Taken from https://archive.stsci.edu/iue/.
in the wavelength range [110-195] nm with the grid of model atmospheres by
Castelli & Kurucz (2003);

  7. 7.

photometry from Hauck & Mermilliod (1997), using [u-b] index adopting de Geus
et al. (1989) $T_{\mathrm{eff}}$ scale.

We adopted the averaged value ($T_{\mathrm{eff}}=21900\pm 1000$ K) as our
$T_{\mathrm{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}\in[3.6,4.0]$.

Table 4: Literature surface gravity estimates for $\mu^{2}$ Sco. References for each estimate are provided in the first column. References: (1): de Geus et al. (1989); (2): Wolff (1990); (3): Grigsby et al. (1992); (4): Nieva (2013). Source | $\log{g}$
---|---
(1) | $3.9\pm 0.2$
(2) | 3.9
(3) | $3.6\pm 0.2$
(4) | 3.916
adopted | $3.8\pm 0.2$

###  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
$\mu^{2}$ Sco.

Table 5: Stellar parameters of $\mu^{2}$ Sco derived in this work. Name | Value
---|---
parallax (mas) | $5.9\pm 0.2$
$T_{\mathrm{eff}}$ (K) | $21700\pm 900$
$E(B-V)$ (mag) | $0.022\pm 0.006$
$A(V)$ (mag) | $0.068\pm 0.015$
age (Myr) | $20\pm 4$
mass ($M_{\odot}$) | $9.1\pm 0.3$
radius ($R_{\odot}$) | $5.6\pm 0.2$

With a spectral type of B2IV, corresponding to a best-fit effective
temperature of $21700\pm 900$ K, $\mu^{2}$ Sco is one of the brightest stars
in the BEAST sample (Janson et al., 2021b). Indeed, its fitted mass of
$M=9.1\pm 0.3\leavevmode\nobreak\ M_{\odot}$ qualifies the star, from a
stellar evolution standpoint, as a ”massive star” ($M>8M_{\odot}$): a star
that could explode as a supernova in the next 10–20 million years (Nomoto &
Leung, 2017), in particular as an electron-capture supernova (Nomoto, 1984).
As a $20\pm 4$ Myr star, $\mu^{2}$ Sco is currently evolving off the main
sequence, consistently with its spectral classification.

We note that in the specific case of $\mu^{2}$ 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 $\mu^{2}$ 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
($\lambda_{K_{1}}=2.1025\pm 0.1020\leavevmode\nobreak\ \mu$m,
$\lambda_{K_{2}}=2.2550\pm\leavevmode\nobreak\ 0.1090\leavevmode\nobreak\
\mu$m; see Vigan et al. 2010) over a field of view (FoV) of $11"\times 11"$,
while IFS collected images – which can be combined to build low-resolution
spectra ($R=30$) – in the $\lambda\in[0.96-1.64]\leavevmode\nobreak\ \mu$m
wavelength range for each position over a field of about $1.7"\times 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.

Table 6: Details of the observations of $\mu^{2}$ Sco. DIT is the integration time, $N_{exp}$ the number of frames after selection (out of 192 frames), Rot is the total FoV rotation. The seeing and the Strehl ratio are averages over the duration of the exposures. UT Date | Instr. | Filter | DIT | $N_{exp}$ | Rot | Seeing | Strehl | Airmass | True north corr. | Platescale
---|---|---|---|---|---|---|---|---|---|---
|  |  | s |  | ∘ | $"$ |  $H$ band |  | ∘ | mas pixel-1
2018-04-24 | IFS | YJH | 16 | 174 | 53.0 | 0.35 | 0.91 | 1.03 | $-1.761\pm 0.06$ | $7.46\pm 0.02$
2018-04-24 | IRDIS | K12 | 16 | 104 | 43.8 | 0.35 | 0.91 | 1.03 | $-1.761\pm 0.06$ | $12.256\pm 0.016$
2021-06-04 | IFS | YJH | 16 | 186 | 53.2 | 0.33 | 0.90 | 1.03 | $-1.780\pm 0.07$ | $7.46\pm 0.02$
2021-06-04 | IRDIS | K12 | 16 | 161 | 51.3 | 0.33 | 0.90 | 1.03 | $-1.780\pm 0.07$ | $12.261\pm 0.005$

##  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 (Section 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 Figure 4 and Figure 5,
respectively.

Figure 4: IRDIS images of $\mu^{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. $\mu^{2}$ Sco b is the source
inside the white circle. Figure 5: IFS images of $\mu^{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 $\sim 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.

##  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 $\sim 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 9.

###  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\times 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 $\sim$ 8.0,
$\text{separation}=128\pm 2$ mas and $\text{PA}=99.5\pm 1.1^{\circ}$. Although
not significant (S/N=2.5), a source at a similar location
($\text{separation}=123\pm 7$ mas, $\text{PA}=114\pm 4^{\circ}$) 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,
$\text{separation}\sim 126\pm 2$, $\text{PA}=99.9\pm 0.8^{\circ}$ for the
first epoch, S/N=5.0, $\text{separation}=125\pm 3$, $\text{PA}=111.0\pm
1.6^{\circ}$ 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 ($\sim 2\times 10^{4}$), the S/N
obtained with PCA-ASDI and PACO-ASDI corresponds to values of the False Alarm
Probability (FAP) equal to $\sim 10^{-11}$, 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 Gaussian666We ran the Shapiro-Wilk test
for Gaussianity (using the on-line calculator available at
https://www.statskingdom.com/shapiro-wilk-test-calculator.html) for the inner
50$\times$50 pixels (about 0.37$\times$0.37 arcsec) region of all these S/N
map, eliminating the area within 0.5 $\lambda$/D of the coronagraphic edge and
that around CC0 and using pixels having even indices in both coordinates, to
avoid the concern related to covariance of adjacent pixels. Depending on the
case, this test found or not some small deviations from Gaussianity. However,
the kurtosis is almost always higher than expected for a Gaussian
distribution, suggesting that outliers are more common than expected with a
Gaussian distribution. (likely because of the edge effects introduced by the
coronagraphic 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 $\sim 11^{\circ}$ from the first-epoch PA, the separation being
the same ($d_{CC0}\approx 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 $\Delta r$ not to be larger than
the fraction $\Delta r_{max}$ of the face-on circular orbit with $r=d_{CC0}$
covered over the timespan $\Delta t\approx 3.1$ yr separating the two epochs:

$\Delta r\lesssim\Delta r_{max}=\frac{2\pi d_{CC0}}{T_{orb}}\cdot\Delta
t\approx 78\text{ mas},$ (1)

where $T_{orb}=1\text{ yr}\cdot(s$
[au]$)^{3/2}\cdot(M_{*}[M_{\odot}])^{-1/2})=1\text{
yr}\cdot\sqrt{21^{3}/9.1}\approx 32$ yr is the period of the circular orbit
with $s=1\text{ au}\cdot d_{CC0}/\text{parallax}\approx 21$ au around
$\mu^{2}$ Sco ($M_{*}=9.1M_{\odot}$). As a comparison, the observed shift of
CC0 is $\Delta r_{obs}=26\pm 4$ mas.

Due to the presence of the coronagraph at $\sim 100$ mas, we conveniently
define as ”interesting area” the annular sector, centered on CC0’s first-epoch
PA, with separation $100\text{ mas}<d<200\text{ mas}$ and semiangular width
$\Delta\text{PA}=\frac{\Delta r_{max}}{2\pi d_{CC0}}\cdot 360^{\circ}\approx
35^{\circ}$.

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 $p_{RPA}\sim 19/99\cdot
35^{\circ}/360^{\circ}\sim 4\%$.

Finally, the separation of CC0 is constant with wavelength at both epochs
(Figure 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 $K_{1}$ and $K_{2}$ data at first epoch – although at a
very low S/N$\sim$2– However, given the very low S/N of these IRDIS
detections, we consider them as upper limits in the following discussion.

Figure 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 9; circle radii have been enlarged to $\sim 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.

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.

##  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^{\circ}$) of the target,
the final images are abundant in detected sources: two companion candidates
are seen inside the IFS FoV Figure 5, while in the wider IRDIS FoV there are
46 objects (Figure( 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.

####  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
(Figure 7). However, owing to the large color uncertainty, their position in
the ($K_{1}-K_{2}$, $K_{1}$) color-magnitude diagram (CMD) does not provide us
with a definitive answer. The same applies for CC29 and CC30, only seen in
$K_{1}$ in one epoch. Given that the geometric probability of a random
alignment between a background source and a target star with separation $d$ is
$\propto d^{2}$, we are only able to argue that the presence of several secure
background sources that have a smaller separation to $\mu^{2}$ Sco than each
of these 11 sources7775 sources for CC7 and CC8, 6 for CC10 and CC11, 14 for
CC20, 18 for CC25, 21 for CC29 and CC30, 24 for CC34 and CC35, 30 for CC42.,
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.

Figure 7: ($K_{1}-K_{2}$, $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 (=$\mu^{2}$ Sco b) is consistent with that
expected for a young object.

####  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 Figure 8).

Figure 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.

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\pm 1.2$
mas (rms: 6.7 mas) along right ascension, $-15.2\pm 1.8$ mas (rms: 10.2 mas)
along declination. The astrometric shift of CC0 is $-44.4\pm 3.5$ mas along
RA, $-96.1\pm 4.0$ mas along declination; that of CC2 is $-40.7\pm 2.1$ mas
along RA, $-77.9\pm 2.1$ mas along declination. Hence, the observed shifts of
CC0 and CC2 are $\sim 8\sigma$ and $\sim 9\sigma$ 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 Section 6.2.

####  6.1.3 Motion of background sources

As shown by Figure 8, the CC cloud is shifted on average by $\sim-16$ mas with
respect to the position expected on the basis of the motion of $\mu^{2}$ 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. 1.

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

  2. 2.

the photocenter of $\mu^{2}$ 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. 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 $\mu^{2}$ 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 $\mu^{2}$ Sco is
very similar to that of its neighbor $\mu^{1}$ Sco ($\mu_{\alpha}=-11.867\pm
0.043$ mas yr-1 in RA and $\mu_{\delta}=-22.511\pm 0.035$ mas yr-1 in Dec) and
to the average value for ELS ($\mu_{\alpha}=-11.57\pm 0.48$ mas yr-1 and
$\mu_{\delta}=-22.12\pm 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 $\mu^{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 $\sim$ 10 years
(see Appendix B).

To investigate scenario 3, we looked for Gaia EDR3 data of field stars near
$\mu^{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:

$\displaystyle<G>$ $\displaystyle=$ $\displaystyle 20.09\leavevmode\nobreak\
\text{mag,}$ $\displaystyle<\pi>$ $\displaystyle=$ $\displaystyle
0.30\leavevmode\nobreak\ \text{mas,}$ $\displaystyle<\mu_{\alpha}^{*}>$
$\displaystyle=$ $\displaystyle-2.681\pm 0.042\leavevmode\nobreak\ \text{ mas
yr}^{-1}\leavevmode\nobreak\ (rms:3.367\leavevmode\nobreak\ \text{ mas
yr}^{-1}),$ $\displaystyle<\mu_{\delta}>$ $\displaystyle=$
$\displaystyle-3.848\pm 0.043\leavevmode\nobreak\ \text{ mas
yr}^{-1}\leavevmode\nobreak\ (rms:3.431\leavevmode\nobreak\ \text{ mas
yr}^{-1}).$

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\pm 0.37$ mas
yr-1 (rms: 2.17 mas yr-1) along RA and $-4.87\pm 0.56$ mas yr-1 (rms: 3.29 mas
yr-1) along dec. We anticipate that a small systematic offset, of about
$-0.6\pm 0.4$ mas along RA and $+1.7\pm 0.4$ mas along dec, is likely to
affect the proper motions of our CCs (see Section 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 $\sim$3 kpc, and absolute magnitudes $M_{G}\sim 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
$M_{K}\sim 4.8$ according to the tables by Pecaut & Mamajek (2013). The
corresponding apparent magnitude $K\sim 17.4$ gives a contrast with respect to
$\mu^{2}$ Sco of d$K=13.1$ mag; these stars are roughly two magnitudes
brighter than the CC cloud (median contrast: d$K\sim 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 $M_{K}\sim 6.7$; according to the tables by
Pecaut & Mamajek (2013) this corresponds to M3V stars that have an absolute
$M_{G}\sim 10.0$, and to an apparent magnitude $G\sim 22.6$. IRDIS CCs are too
faint to be detected by Gaia (that has a limiting magnitude of $G\sim 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 $\mu^{2}$ Sco

In this subsection, we examine the possibility that CC2 is not a bound
companion of $\mu^{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 $\mu^{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 $\mu^{2}$ Sco because of a ”remarkable”
astrometric shift relative to $\mu^{2}$ Sco. Two factors must be taken into
account: the probability of finding an object within certain boundaries of
$\mu_{\alpha}^{*}$ and $\mu_{\delta}$; the fact that we have fully888That is,
at least twice. As already mentioned, two observations are needed to compute
relative proper motions. observed 25 stars in the survey.

Starting from the Gaia bg sample, we identify as interesting the stars with
$\mu_{\alpha}^{*}<-9.8$ mas yr-1 and $\mu_{\delta}<-21.1$ mas yr-1 999This is
a conservative approach that considers interesting all the objects having a
larger relative proper motion – with respect to the CC cloud – than $\mu^{2}$
Sco itself; see below.. After excluding a few objects that, despite passing
the test, have parallax $>5$ mas (which would correspond to $4\lesssim
dK\lesssim 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 $f_{G}=1.67\cdot 10^{-5}$. The FAP of finding one star
with these features in our observations of $\mu^{2}$ Sco (that is, within the
IRDIS field of view), given that we observe 36 CCs in both epochs, is then
$6\cdot 10^{-4}$.

A similar argument can be made by creating a synthetic sample of background
stars by means of the Besançon Galaxy Model interface101010Available at
https://model.obs-besancon.fr/modele_home.php. (Czekaj et al., 2014). We
selected a sample of stars with distance $d\in[0,50]$ kpc, apparent $K$
magnitude similar to that observed for CC0 ($K\in[15,17]$ mag), radius
$\rho=5$ arcmin around the position of $\mu^{2}$ Sco. Out of 6595 stars, just
4 (corresponding to $f_{B}=1.34\cdot 10^{-5}$) pass the proper motion test
described above; this yields a FAP -given 36 CCs- of $4.8\cdot 10^{-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
$\mu^{2}$ Sco, but we assume for simplicity a strict equality; if we further
assume that the fractions $f_{B}$ and $f_{G}$ 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 $\sim
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\cdot 10^{-3}$ (Gaia-based test) and $3.7\cdot
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 $\mu^{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\leavevmode\nobreak\ M_{\mathrm{J}}$) 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 association111111Given
that any Gaia-based census of association members is not complete at such low
masses..

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\leavevmode\nobreak\ M_{\odot}$. 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
$[5M_{\mathrm{J}},75M_{\mathrm{J}}]$. We do not consider objects below
$5M_{\mathrm{J}}$ 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 $\mu^{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\leavevmode\nobreak\ M_{\odot}$ (see “Stellar system analysis”). Since
Gaia is complete within this magnitude range (its limiting magnitude is $G\sim
21$), we assume a sharp transition between 100% completeness above 15
$M_{\mathrm{J}}$ to 0% completeness below 15 $M_{\mathrm{J}}$. 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_{\odot}$, obtaining a complete census
of approximately $\sim 4021$ sources for UCL, $\sim 602$ for LS, $\sim 76$ for
ELS.

Multiplying these numbers by the integral of the PDF from $5M_{\mathrm{J}}$ to
$75M_{\mathrm{J}}$, we obtain the expected number of interesting objects
$\tilde{N}$:

$\displaystyle\tilde{N}_{LS}$ $\displaystyle\sim$ $\displaystyle 123,$
$\displaystyle\tilde{N}_{ELS}$ $\displaystyle\sim$ $\displaystyle 16,$
$\displaystyle\tilde{N}_{UCL}$ $\displaystyle\sim$ $\displaystyle 823.$

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^{\circ},343^{\circ}]\times[2^{\circ},28^{\circ}]$ as in Pecaut &
Mamajek (2016) were used. For LS and ELS we computed the mean ¡l¿, ¡b¿ and the
related standard deviations $\sigma_{l}$ and $\sigma_{b}$, and defined the
areas as
$[<l>-2\sigma_{l},<l>+2\sigma_{l}]\times[<b>-2\sigma_{b},<b>+2\sigma_{b}]$. We
get:

$\displaystyle A_{LS}$ $\displaystyle=$ $\displaystyle
23.2\leavevmode\nobreak\ \text{deg}^{2}=3.0\cdot 10^{8}\leavevmode\nobreak\
\text{arcsec}^{2},$ $\displaystyle A_{ELS}$ $\displaystyle=$ $\displaystyle
3.8\leavevmode\nobreak\ \text{deg}^{2}=4.9\cdot 10^{7}\leavevmode\nobreak\
\text{arcsec}^{2},$ $\displaystyle A_{UCL}$ $\displaystyle=$ $\displaystyle
747.0\leavevmode\nobreak\ \text{deg}^{2}=9.7\cdot 10^{9}\leavevmode\nobreak\
\text{arcsec}^{2},$

so that the projected densities are:

$\displaystyle\Sigma_{LS}$ $\displaystyle=$
$\displaystyle\tilde{N}_{LS}/A_{LS}=4\cdot 10^{-7}\leavevmode\nobreak\
\text{arcsec}^{2},$ $\displaystyle\Sigma_{ELS}$ $\displaystyle=$
$\displaystyle\tilde{N}_{ELS}/A_{ELS}=3.2\cdot 10^{-7}\leavevmode\nobreak\
\text{arcsec}^{2},$ $\displaystyle\Sigma_{UCL}$ $\displaystyle=$
$\displaystyle\tilde{N}_{UCL}/A_{UCL}=8.5\cdot 10^{-8}\leavevmode\nobreak\
\text{arcsec}^{2}.$

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"\times 11"$):

$\displaystyle N_{BD,UCL}$ $\displaystyle=$
$\displaystyle\Sigma_{UCL}\cdot(11")^{2}=1.0\cdot 10^{-5},$ $\displaystyle
N_{BD,LS}$ $\displaystyle=$ $\displaystyle\Sigma_{LS}\cdot(11")^{2}=5.0\cdot
10^{-5},$

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:

$\displaystyle p_{BD,UCL}$ $\displaystyle=$ $\displaystyle 2.6\cdot 10^{-4},$
$\displaystyle p_{BD,LS}$ $\displaystyle=$ $\displaystyle 1.0\cdot 10^{-3}.$

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 $\mu^{2}$ Sco ($\mu^{2}$
Sco b). As regards CC0, based on the reliability arguments presented in the
Section 5.2, we consider it to be a probable detection that would make it a
second physical companion to $\mu^{2}$ Sco.

###  6.3 Companion characterization

The star $\mu^{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 $\mu^{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 one121212As a consistency check, we stress that in the H-band region,
where the signal is highest, the two spectra are consistent with one another..
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
$K_{1}$ and $K_{2}$ band provided by the same algorithm applied to the IRDIS
first-epoch dataset (see Figure 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 $\mu^{2}$
Sco. The best-fit model has a mass of $M_{c}=19.5\pm 0.9$ $M_{\mathrm{J}}$, a
temperature of $T_{\mathrm{eff}}$=$2262\pm 28$ K and
$\log{L/L_{\odot}}=-3.08\pm 0.03$. The comparison is generally fairly good
(reduced $\chi^{2}=2.55$); most of the contribution to $\chi^{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\leavevmode\nobreak\ \mu$m according to any realistic atmospheric model
(compare, e.g., Marley et al., 2021), we attribute this discrepancy to
residual speckle noise.

Figure 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\sigma$ 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.5M_{\mathrm{J}}$ brown dwarf using AMES-Dusty
models. The dot-dashed line is a black body curve with the same temperature
and radius as this model.

A similar fit with BHAC15 models yields $M_{c}=17.4\pm 0.9M_{\mathrm{J}}$,
$T_{\mathrm{eff}}=2274\pm 28$ K and $\log{L/L_{\odot}}=-3.04\pm 0.03$. We
account for theoretical uncertainties by averaging the two estimates to obtain
a final mass of $M_{c}=18.5\pm 1.5M_{\mathrm{J}}$.

We highlight that, despite the high level of irradiation at its expected
orbital distance ($21\pm 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 $\mu^{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 $K_{1}$ and $K_{2}$ band; combining them with
the $K$ magnitude of the primary, we derive absolute ($K_{1}$, $K_{2}$)
magnitudes. No significant photometric variation is observed between the two
epochs. The position of the companion in the ($K_{1}-K_{2}$, $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 Figure 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 $M_{b}=14.4\pm 0.8$ $M_{\mathrm{J}}$.

As regards CC0, the reported magnitudes are computed by collapsing IFS
spectral channels 12-21 (1.159-1.333 $\mu$m) and 30-38 (1.504-1.636 $\mu$m)
from the first epoch to build the two bands $J_{IFS}=1.246\leavevmode\nobreak\
\mu$m (band width=0.174 $\mu$m) and $H_{IFS}=1.570\leavevmode\nobreak\ \mu$m
(band width=0.132 $\mu$m); the contrasts in the two bands are translated into
absolute magnitudes by adding 2MASS $J$ and $H$ magnitudes of the primary.

####  6.3.2 Astrometry and orbits

The position of $\mu^{2}$ Sco b and CC0 relative to the star is not exactly
the same in the two epochs, showing a significant displacement of $\sim 10$
mas and $\sim 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 Section 6.1.3 (6286 sources), we
select the 2741 sources with parallax 0.1 mas$<\pi<$5 mas (i.e.,
$0.2\leavevmode\nobreak\ \text{kpc}<d<10\leavevmode\nobreak\ \text{kpc}$). The
mean proper motion of this ”restricted Gaia bg sample” is $-2.51\pm 0.07$ mas
yr-1 (rms=3.53 mas yr-1) along right ascension, and $-3.82\pm 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
$\mu_{\alpha}^{*}$, one for $\mu_{\delta}$) at level $\alpha=0.05$. While the
null hypothesis could not be rejected for $\mu_{\alpha}^{*}$
($p=0.13>\alpha$), a significant difference exists with respect to
$\mu_{\delta}$ ($p=0.008<\alpha$). In particular, the CC cloud appears to have
a somewhat higher $\mu_{\alpha}^{*}$ and lower $\mu_{\delta}$ than the
restricted Gaia bg sample. To quantify this, we identified the range of
$\mu_{\alpha}^{*}$ and $\mu_{\delta}$ shifts to be solidly applied to the
whole all CC cloud so that the KS test is passed at level $\alpha=0.05$.

The range of shifts needed for $\mu_{\alpha}^{*}$ is
$\Delta\mu_{\alpha}^{*}\in[-1.44,+0.22]$ mas yr-1, while for $\mu_{\delta}$ it
is $\Delta\mu_{\delta}\in[0.55,2.81]$ 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 $\Delta\alpha=-1.34$ mas, $\Delta\delta=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 Section 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 $\mu^{2}$ Sco b and
CC0, we conservatively opted for treating $\Delta\alpha$ and $\Delta\delta$ as
an additional source of random uncertainty on their relative astrometry, and
propagate it to derive final uncertainties $\sigma_{d}$ and
$\sigma_{\text{PA}}$ on separation and PA, respectively, that are somewhat
broadened with respect to those in Table 9: for CC0,
$(\sigma_{d},\sigma_{\text{PA}})=(2.5\text{ mas},2.0^{\circ})$ and $(3.3\text{
mas},2.2^{\circ})$ for the first and second epoch, respectively; for CC2 =
$\mu^{2}$ Sco b, $(\sigma_{d},\sigma_{\text{PA}})=(3.1\text{
mas},0.12^{\circ})$ and $(2.8\text{ mas},0.11^{\circ})$ 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 $\mu^{2}$ Sco b using the orbitize!
Python package131313Available at https://github.com/sblunt/orbitize ., run
with the recommended parameters for reliable convergence, for a total of
$2\cdot 10^{6}$ orbits (Blunt et al., 2020). The priors used for the MCMC were
stellar mass and parallax, taken from Table 5, as well as the star-planet
separations and position angles, taken from Table 9. A subsample of suitable
orbits, as well as the posterior distributions of semimajor axis,
eccentricity, and inclination, are shown in Figure 10-13.

Figure 10: Subset of possible orbits for the robust companion candidate
$\mu^{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 $\mu^{2}$ Sco is
indicated by the yellow star. (B) Evolution over time of the planet-star
separation ($\rho$); 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. Figure
11: Posterior distributions of orbital parameters for the robust companion
$\mu^{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\sigma$ 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. Figure 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 $\mu^{2}$ Sco is indicated by the
yellow star. (B) Evolution over time of the planet-star separation ($\rho$);
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. Figure 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 Figure 11.

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, ($m_{1}$, $a_{1}$, $e_{1}$) and ($m_{2}$, $a_{2}$, $e_{2}$)
the mass, semimajor axis and eccentricity of the inner and outer companion,
respectively. Under the hypothesis that $m_{1}\ll M_{*}$ and $m_{2}\ll M_{*}$,
a system is Hill stable, meaning that the two companions will avoid close
approaches at all times, if:

$\alpha^{-3}\left(\mu_{1}+\frac{\mu_{2}}{\delta^{2}}\right)(\mu_{1}\gamma_{1}+\mu_{2}\gamma_{2}\delta)^{2}\geq
1+3^{4/3}\frac{\mu_{1}\mu_{2}}{\alpha^{4/3}}$ (2)

Gladman (1993). Here $\mu_{1}=m_{1}/M_{*}$, $\mu_{2}=m_{2}/M_{*}$,
$\alpha=\mu_{1}+\mu_{2}$, $\Delta=a_{2}-a_{1}$,
$\delta=\sqrt{1+\Delta/a_{1}}$, $\gamma_{1}=\sqrt{1-e_{1}^{2}}$,
$\gamma_{2}=\sqrt{1-e_{2}^{2}}$.

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 $\sim 5$ au
up to $\sim 800$ au), though this cannot be excluded at present due to the
large uncertainties still existing on the orbital parameters.

The final parameters of $\mu^{2}$ Sco b and the candidate CC0 derived
throughout this Section are reported in Table 7.

Table 7: Absolute magnitudes, mass, projected separation, semimajor axis, eccentricity and inclination for the robust $\mu^{2}$ Sco b and the probable CC0. Here ($F_{1}$, $F_{2}$) = ($K_{1}$, $K_{2}$) for $\mu^{2}$ Sco b, ($J$, $H_{2}$) for CC0. | $F_{1}$ | $F_{2}$ | Mass | $q$ | Proj. sep | $a$ | $e$ | $i$
---|---|---|---|---|---|---|---|---
| mag | mag | $M_{\mathrm{J}}$ |  | au | au |  | ∘
$\mu^{2}$ Sco b | $9.86\pm 0.33$ | $9.48\pm 0.32$ | $14.4\pm 0.8$ | 0.0015(1) | $290\pm 10$ | $242.4^{+114.5}_{-52.1}$ | $0.56^{+0.27}_{-0.26}$ | $96.6^{+21.5}_{-20.5}$
CC0 | $11.51\pm 1.10$ | $9.69\pm 0.38$ | $18.5\pm 1.5$ | 0.0019(2) | $21\pm 1$ | $18.9^{+11.7}_{-5.0}$ | $0.61^{+0.19}_{-0.32}$ | $62.8^{+9.9}_{-16.6}$

##  7 Discussion

The two objects, having $M_{b}=14.4\pm 0.8M_{\mathrm{J}}$ and $M_{c}=18.5\pm
1.5M_{\mathrm{J}}$, are just above the deuterium-burning limit ($M\sim
13M_{\mathrm{J}}$) 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 ($q_{b}=0.0015$, $q_{c}=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 $\mu^{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 (Figure
14). The $\mu^{2}$ Sco system appears to be in many regards a scaled-up
version of the Solar System.

Figure 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. $\mu^{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/).

Leaving aside the case to be confirmed of CC0, the robust discovery of a
physical companion to $\mu^{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 ($M=10.9\pm 1.6M_{\mathrm{J}}$,
$q=0.0013_{-0.0003}^{+0.0004}$) around the binary b Cen system. In that case
the primary star has a mass of $\sim 6$ M⊙ and the overall mass is similar to
that of $\mu^{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 system141414Even in the worst case of the relatively compact LS, the
timescale for a close encounter between $\mu^{2}$ Sco and another star within
1000 au is $\sim 200$ Myr, using Eq. 1 from Malmberg et al. (2007) with
$n=1pc^{-3}$, $M_{tot}=10M_{\odot}$ and $v=0.7$ km s-1.; 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.

$\mu^{2}$ 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 Figure 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 $1M_{\mathrm{J}}\lesssim M\lesssim 40M_{\mathrm{J}}$ 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 surveys151515Mass and separation
ranges are highly varying across different direct imaging studies. See Table 3
from Vigan et al. (2021) for a comprehensive comparison of the cited studies.
(Bowler, 2016; Nielsen et al., 2019; Vigan et al., 2021). As mentioned in
Section 1, radial velocity studies indicate a turnover in this trend at about
2 $M_{\odot}$; 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 $\mu^{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 $M_{\mathrm{J}}$is determined by the possibility to ignite hydrogen
burning in their core, the transition from giant planets to brown dwarfs is
often set to $\sim 13\leavevmode\nobreak\ M_{\mathrm{J}}$, the so-called
deuterium-burning limit (DBL). According to this definition, the $\mu^{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<4M_{\mathrm{J}}$ and stellar metallicity, hint of formation within a
protoplanetary disk; more massive ($M\sim 10$ $M_{\mathrm{J}}$) 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\leavevmode\nobreak\ M_{\mathrm{J}}$ and $25\leavevmode\nobreak\
M_{\mathrm{J}}$) to the evolved 2.7 $M_{\odot}$ star $\nu$ 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<13M_{\mathrm{J}}$. 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 $\mu^{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 (Figure
14); the mass ratios are smaller than those in the $\nu$ 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),
$\mu^{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,
$\mu^{2}$ Sco b and b Cen b should both be considered as planets.

A third way to assess whether $\mu^{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 $\mu^{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 disk161616The latter mechanism, according to the
second distinction, would instead fall under the ”planet” definition.. One
might wonder if objects like $\mu^{2}$ Sco b and b Cen b constitute the high-
mass 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):

$\frac{\partial^{2}N}{\partial\log_{10}(a)\partial
q}\propto\frac{(\log_{10}(a[au])-2.59)^{2}}{2\cdot 0.79^{2}}\cdot q^{0.25},$
(3)

normalizing it to the median frequency of $1-75\leavevmode\nobreak\
M_{\mathrm{J}}$ 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:

$\frac{\partial^{2}N}{\partial\log_{10}(a)\partial
q}\propto\frac{(\log_{10}(a[au])-0.79)^{2}}{2\cdot 0.77^{2}}\cdot q^{-1.31}.$
(4)

The expected number of companions around one star with $a\in[100,1000]$ au,
$q\in[0.0005,0.0030]$ is $1.4\cdot 10^{-4}$ for the SP scenario, $9.0\cdot
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\cdot 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.5M_{\odot}$ (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 $\mu^{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 $\sim 9\leavevmode\nobreak\
M_{\odot}$ 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
($\sim 10\leavevmode\nobreak\ M_{\oplus}$), 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 $\mu^{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\leavevmode\nobreak\ {\rm
eV}<h\nu<13.6\leavevmode\nobreak\ {\rm eV}$; EUV: $13.6\leavevmode\nobreak\
{\rm eV}<h\nu<0.1\leavevmode\nobreak\ {\rm keV}$) and X-ray
($0.1\leavevmode\nobreak\ {\rm keV}<h\nu<2\leavevmode\nobreak\ {\rm 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 $\sim 10^{5}$ years) 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 ($M_{A}\sim 5.5\leavevmode\nobreak\ M_{\odot}$,
$M_{B}\sim 3.5\leavevmode\nobreak\ M_{\odot}$) and the single $\mu^{2}$ Sco
($M\sim 9.1\leavevmode\nobreak\ M_{\odot}$). 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 $\mu^{2}$ Sco would yield a three
times shorter disk lifetime around $\mu^{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 $\mu^{2}$ Sco. Combining Eq. (11) from Cesaroni et al. (2007) for the
photoevaporation outflow rate with the expected ionizing photon flux $\Phi$
for the three stars, the ratio171717We interpolate between $\log{\Phi_{\rm
EUV}}$ values for $M=3\leavevmode\nobreak\ M_{\odot}$ and
$M=7\leavevmode\nobreak\ M_{\odot}$ (see Table 2 by Gorti & Hollenbach, 2009)
and $\phi_{i}$ values for $M=16.8\leavevmode\nobreak\ M_{\odot}$,
$M=25.6\leavevmode\nobreak\ M_{\odot}$ and $M=65\leavevmode\nobreak\
M_{\odot}$ (Hollenbach et al., 1994), deriving the empirical relation
$\log{\Phi}=-5.9705\log{M}+21.553\log{M}+30.1$. The b Cen photoevaporation
outflow rate is the sum of the individual contributions of b Cen A and b Cen
B. between the two disk survival timescales should be around $\sim 20$. While
the impact of this on planet formation is difficult to be properly assessed,
the presence of (at least) a companion with $q\sim 0.002$ around $\mu^{2}$ 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)181818Alternatively, the problem might be alleviated if
$\mu^{2}$ Sco were actually formed by the merging of two nearly equal mass
stars, which is possibly not an exotic case (de Mink et al., 2014). On this
respect, it is notable that $\mu^{2}$ Sco appears to be a slow rotator (see
Appendix C); in fact it has been argued recently that post-mergers should
appear as slow rotators Schneider et al. (2019); Wang et al. (2022). On the
other hand, the same studies suggest that mergers might have strong magnetic
fields, but there is no evidence for this in $\mu^{2}$ Sco (see Appendix C).

On the other hand, giant planets may possibly form inside massive
protoplanetary disks on a much shorter timescale of $\sim 10^{4}$ 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\leavevmode\nobreak\ M_{\odot}$ 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 $\mu^{2}$ 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\leavevmode\nobreak\
M_{\odot}$) has hinted at an interesting eccentricity dichotomy between the
two populations: while the former usually has low eccentricity values
($e=0.13^{+0.12}_{-0.08}$), 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 $\mu^{2}$ 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 $\mu^{2}$ 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.2M_{\odot}$
(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 $\mu^{2}$ 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 $\mu^{2}$ 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.

##  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\gtrsim
2.4\leavevmode\nobreak\ M_{\odot}$) 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 $\mu^{2}$ 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 $\beta$ 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\pm 4$ Myr and a
mass of $9.1\pm 0.3$ M⊙, $\mu^{2}$ 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 $\Delta K_{1}\approx
11.7$ mag and $\Delta K_{2}\approx 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 $\mu^{2}$ Sco and found both hypothesis to be highly
improbable. We have therefore concluded that this object – $\mu^{2}$ Sco b– is
a real physical substellar companion to the star. With a best-fit $M=14.4\pm
0.8\leavevmode\nobreak\ M_{\mathrm{J}}$, $\mu^{2}$ Sco b is the first ”planet-
like”, in the sense defined in Section 7, companion to a supernova progenitor.

A second $18.5\pm 1.5\leavevmode\nobreak\ M_{\mathrm{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 $\mu^{2}$ Sco b have angular separations of $0.128\pm 0.002$ arcsec
and $1.709\pm 0.005$ arcsec from the star that translate, at a distance
$d=169\pm 6$ pc, into projected separations of $21\pm 1$ au and $290\pm 10$
au, respectively. The small astrometric displacements over the three years
separating the two epochs (Figure 8) are a consequence of their orbital motion
around $\mu^{2}$ 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
(Figure 10-13).

$\mu^{2}$ Sco b (and possibly CC0) are naturally added to b Cen b, a $10.9\pm
1.6M_{\mathrm{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 $\mu^{2}$ 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 $\mu^{2}$ 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 $\mu^{2}$ 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.

###### Acknowledgements.

We are extremely grateful to the anonymous referee for the constructive
comments, 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. VS, RG, DM, SD and VD 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. NE, SPQ
and LM acknowledge the financial support of the SNSF. GDM acknowledges the
support of the Deutsche Forschungsgemeinschaft (DFG) priority program SPP 1992
“Exploring the Diversity of Extrasolar Planets” (MA 9185/1-1). GDM also
acknowledges the support from the Swiss National Science Foundation under
grant BSSGI0_155816 “PlanetsInTime”. 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

  * Alecian et al. (2016) Alecian, E., Tkachenko, A., Neiner, C., Folsom, C. P., & Leroy, B. 2016, A&A, 589, A47
  * Allard et al. (2001) Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357
  * Allard et al. (2011) Allard, F., Homeier, D., & Freytag, B. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 448, 16th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. C. Johns-Krull, M. K. Browning, & A. A. West, 91
  * Balona (1984) Balona, L. A. 1984, MNRAS, 211, 973
  * Balona (1994) Balona, L. A. 1994, MNRAS, 268, 119
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42
  * Barrado (2016) Barrado, D. 2016, in EAS Publications Series, Vol. 80-81, EAS Publications Series, 115–175
  * Beltrán & de Wit (2016) Beltrán, M. T. & de Wit, W. J. 2016, A&A Rev., 24, 6
  * Berghoefer et al. (1996) Berghoefer, T. W., Schmitt, J. H. M. M., & Cassinelli, J. P. 1996, A&AS, 118, 481
  * Beuzit et al. (2019) Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155
  * Blunt et al. (2020) Blunt, S., Wang, J. J., Angelo, I., et al. 2020, AJ, 159, 89
  * Boley et al. (2010) Boley, A. C., Hayfield, T., Mayer, L., & Durisen, R. H. 2010, Icarus, 207, 509
  * Bonavita (2020) Bonavita, M. 2020, Exo-DMC: Exoplanet Detection Map Calculator
  * Bonomo et al. (2017) Bonomo, A. S., Desidera, S., Benatti, S., et al. 2017, A&A, 602, A107
  * Boss (2003) Boss, A. P. 2003, ApJ, 599, 577
  * Bowler (2016) Bowler, B. P. 2016, PASP, 128, 102001
  * Bowler et al. (2020) Bowler, B. P., Blunt, S. C., & Nielsen, E. L. 2020, AJ, 159, 63
  * Brown & Verschueren (1997) Brown, A. G. A. & Verschueren, W. 1997, A&A, 319, 811
  * Cadman et al. (2020) Cadman, J., Rice, K., Hall, C., Haworth, T. J., & Biller, B. 2020, MNRAS, 492, 5041
  * Carbillet et al. (2011) Carbillet, M., Bendjoya, P., Abe, L., et al. 2011, Experimental Astronomy, 30, 39
  * Cardiel et al. (2021) Cardiel, N., Zamorano, J., Bará, S., et al. 2021, MNRAS, 504, 3730
  * Castelli (1991) Castelli, F. 1991, A&A, 251, 106
  * Castelli & Kurucz (2003) Castelli, F. & Kurucz, R. L. 2003, in Modelling of Stellar Atmospheres, ed. N. Piskunov, W. W. Weiss, & D. F. Gray, Vol. 210, A20
  * Cesaroni et al. (2007) Cesaroni, R., Galli, D., Lodato, G., Walmsley, C. M., & Zhang, Q. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 197
  * Chauvin et al. (2017) Chauvin, G., Desidera, S., Lagrange, A. M., et al. 2017, in SF2A-2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. C. Reylé, P. Di Matteo, F. Herpin, E. Lagadec, A. Lançon, Z. Meliani, & F. Royer, Di
  * Claudi et al. (2008) Claudi, R. U., Turatto, M., Gratton, R. G., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II, ed. I. S. McLean & M. M. Casali, 70143E
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, II/246
  * Czekaj et al. (2014) Czekaj, M. A., Robin, A. C., Figueras, F., Luri, X., & Haywood, M. 2014, A&A, 564, A102
  * Damiani et al. (2019) Damiani, F., Prisinzano, L., Pillitteri, I., Micela, G., & Sciortino, S. 2019, A&A, 623, A112
  * de Geus et al. (1989) de Geus, E. J., de Zeeuw, P. T., & Lub, J. 1989, A&A, 216, 44
  * de Mink et al. (2014) de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7
  * De Rosa et al. (2014) De Rosa, R. J., Patience, J., Wilson, P. A., et al. 2014, MNRAS, 437, 1216
  * de Zeeuw et al. (1999) de Zeeuw, P. T., Hoogerwerf, R., de Bruijne, J. H. J., Brown, A. G. A., & Blaauw, A. 1999, AJ, 117, 354
  * Delorme et al. (2017) Delorme, P., Meunier, N., Albert, D., et al. 2017, in SF2A-2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. C. Reylé, P. Di Matteo, F. Herpin, E. Lagadec, A. Lançon, Z. Meliani, & F. Royer, Di
  * Deng et al. (2021) Deng, H., Mayer, L., & Helled, R. 2021, Nature Astronomy, 5, 440
  * Dohlen et al. (2008) Dohlen, K., Langlois, M., Saisse, M., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II, ed. I. S. McLean & M. M. Casali, 70143L
  * Donati et al. (2006) Donati, J. F., Catala, C., Landstreet, J. D., & Petit, P. 2006, in Astronomical Society of the Pacific Conference Series, Vol. 358, Solar Polarization 4, ed. R. Casini & B. W. Lites, 362
  * Donati et al. (1997) Donati, J. F., Semel, M., Carter, B. D., Rees, D. E., & Collier Cameron, A. 1997, MNRAS, 291, 658
  * Flasseur et al. (2018) Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2018, A&A, 618, A138
  * Flasseur et al. (2020) Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2020, A&A, 637, A9
  * Fontanive et al. (2020) Fontanive, C., Allers, K. N., Pantoja, B., et al. 2020, ApJ, 905, L14
  * Forgan & Rice (2013) Forgan, D. & Rice, K. 2013, MNRAS, 432, 3168
  * Forgan et al. (2018) Forgan, D. H., Hall, C., Meru, F., & Rice, W. K. M. 2018, MNRAS, 474, 5036
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Galicher et al. (2018) Galicher, R., Boccaletti, A., Mesa, D., et al. 2018, A&A, 615, A92
  * Gaudi et al. (2021) Gaudi, B. S., Meyer, M., & Christiansen, J. 2021, in ExoFrontiers; Big Questions in Exoplanetary Science, ed. N. Madhusudhan, 2–1
  * Gladman (1993) Gladman, B. 1993, Icarus, 106, 247
  * Gontcharov (2006) Gontcharov, G. A. 2006, Astronomy Letters, 32, 759
  * Gorti et al. (2009) Gorti, U., Dullemond, C. P., & Hollenbach, D. 2009, ApJ, 705, 1237
  * Gorti & Hollenbach (2009) Gorti, U. & Hollenbach, D. 2009, ApJ, 690, 1539
  * Grigsby et al. (1992) Grigsby, J. A., Morrison, N. D., & Anderson, L. S. 1992, ApJS, 78, 205
  * Hauck & Mermilliod (1997) Hauck, B. & Mermilliod, M. 1997, VizieR Online Data Catalog, II/215
  * Hayashi et al. (1985) Hayashi, C., Nakazawa, K., & Nakagawa, Y. 1985, in Protostars and Planets II, ed. D. C. Black & M. S. Matthews, 1100–1153
  * Helled et al. (2014) Helled, R., Bodenheimer, P., Podolak, M., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 643
  * Hiltner et al. (1969) Hiltner, W. A., Garrison, R. F., & Schild, R. E. 1969, ApJ, 157, 313
  * Hobbs (1978) Hobbs, L. M. 1978, ApJ, 222, 491
  * Hollenbach et al. (1994) Hollenbach, D., Johnstone, D., Lizano, S., & Shu, F. 1994, ApJ, 428, 654
  * Hubeny & Lanz (2011) Hubeny, I. & Lanz, T. 2011, TLUSTY: Stellar Atmospheres, Accretion Disks, and Spectroscopic Diagnostics
  * Janson et al. (2021a) Janson, M., Gratton, R., Rodet, L., et al. 2021a, Nature, 600, 231
  * Janson et al. (2021b) Janson, M., Squicciarini, V., Delorme, P., et al. 2021b, A&A, 646, A164
  * Johnson & Morgan (1953) Johnson, H. L. & Morgan, W. W. 1953, ApJ, 117, 313
  * Johnson et al. (2010) Johnson, J. A., Aller, K. M., Howard, A. W., & Crepp, J. R. 2010, PASP, 122, 905
  * Kaufer et al. (1999) Kaufer, A., Stahl, O., Tubbesing, S., et al. 1999, The Messenger, 95, 8
  * Kennedy & Kenyon (2008) Kennedy, G. M. & Kenyon, S. J. 2008, ApJ, 673, 502
  * Kervella et al. (2019) Kervella, P., Arenou, F., Mignard, F., & Thévenin, F. 2019, A&A, 623, A72
  * Kratter et al. (2010) Kratter, K. M., Murray-Clay, R. A., & Youdin, A. N. 2010, ApJ, 710, 1375
  * Lafrenière et al. (2007) Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É. 2007, ApJ, 660, 770
  * Lallement et al. (2018) Lallement, R., Capitanio, L., Ruiz-Dern, L., et al. 2018, A&A, 616, A132
  * Lambrechts & Johansen (2012) Lambrechts, M. & Johansen, A. 2012, A&A, 544, A32
  * Lamers et al. (2005) Lamers, H. J. G. L. M., Gieles, M., Bastian, N., et al. 2005, A&A, 441, 117
  * Lecavelier des Etangs & Lissauer (2022) Lecavelier des Etangs, A. & Lissauer, J. J. 2022, arXiv e-prints, arXiv:2203.09520
  * Lenz et al. (2017) Lenz, D., Hensley, B. S., & Doré, O. 2017, ApJ, 846, 38
  * Lindegren et al. (2021) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2
  * Macintosh et al. (2015) Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Science, 350, 64
  * Maire et al. (2016) Maire, A.-L., Langlois, M., Dohlen, K., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI, ed. C. J. Evans, L. Simard, & H. Takami, 990834
  * Malik et al. (2015) Malik, M., Meru, F., Mayer, L., & Meyer, M. 2015, ApJ, 802, 56
  * Malmberg et al. (2007) Malmberg, D., de Angeli, F., Davies, M. B., et al. 2007, MNRAS, 378, 1207
  * Marigo et al. (2017) Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77
  * Marleau et al. (2019) Marleau, G.-D., Coleman, G. A. L., Leleu, A., & Mordasini, C. 2019, A&A, 624, A20
  * Marley et al. (2021) Marley, M. S., Saumon, D., Visscher, C., et al. 2021, ApJ, 920, 85
  * Mawet et al. (2014) Mawet, D., Milli, J., Wahhaj, Z., et al. 2014, ApJ, 792, 97
  * Mayor et al. (2003) Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20
  * Mayor & Queloz (1995) Mayor, M. & Queloz, D. 1995, Nature, 378, 355
  * McCall (2004) McCall, M. L. 2004, AJ, 128, 2144
  * Mermilliod (2006) Mermilliod, J. C. 2006, VizieR Online Data Catalog, II/168
  * Mesa et al. (2015) Mesa, D., Gratton, R., Zurlo, A., et al. 2015, A&A, 576, A121
  * Meyer et al. (2021) Meyer, M., Calissendorff, P., & Amara, A. 2021, submitted on A&A
  * Miret-Roig et al. (2022) Miret-Roig, N., Bouy, H., Raymond, S. N., et al. 2022, Nature Astronomy, 6, 89
  * Moe & Di Stefano (2017) Moe, M. & Di Stefano, R. 2017, ApJS, 230, 15
  * Mordasini et al. (2012) Mordasini, C., Alibert, Y., Benz, W., Klahr, H., & Henning, T. 2012, A&A, 541, A97
  * Napiwotzki et al. (1993) Napiwotzki, R., Schoenberner, D., & Wenske, V. 1993, A&A, 268, 653
  * Nielsen et al. (2019) Nielsen, E., De Rosa, R., Macintosh, B., et al. 2019, in AAS/Division for Extreme Solar Systems Abstracts, Vol. 51, AAS/Division for Extreme Solar Systems Abstracts, 100.02
  * Nieva (2013) Nieva, M. F. 2013, A&A, 550, A26
  * Nomoto (1984) Nomoto, K. 1984, ApJ, 277, 791
  * Nomoto & Leung (2017) Nomoto, K. & Leung, S.-C. 2017, in Handbook of Supernovae, ed. A. W. Alsabti & P. Murdin, 483
  * Pascucci et al. (2016) Pascucci, I., Testi, L., Herczeg, G. J., et al. 2016, ApJ, 831, 125
  * Pavlov et al. (2008) Pavlov, A., Feldt, M., & Henning, T. 2008, in Astronomical Society of the Pacific Conference Series, Vol. 394, Astronomical Data Analysis Software and Systems XVII, ed. R. W. Argyle, P. S. Bunclark, & J. R. Lewis, 581
  * Pecaut & Mamajek (2013) Pecaut, M. J. & Mamajek, E. E. 2013, ApJS, 208, 9
  * Pecaut & Mamajek (2016) Pecaut, M. J. & Mamajek, E. E. 2016, MNRAS, 461, 794
  * Perryman et al. (1997) Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49
  * Piskunov et al. (2011) Piskunov, N., Snik, F., Dolgopolov, A., et al. 2011, The Messenger, 143, 7
  * Pollack et al. (1996) Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62
  * Poppenhaeger & Wolk (2014) Poppenhaeger, K. & Wolk, S. J. 2014, A&A, 565, L1
  * Poznanski et al. (2012) Poznanski, D., Prochaska, J. X., & Bloom, J. S. 2012, MNRAS, 426, 1465
  * Quirrenbach et al. (2019) Quirrenbach, A., Trifonov, T., Lee, M. H., & Reffert, S. 2019, A&A, 624, A18
  * Rameau et al. (2013) Rameau, J., Chauvin, G., Lagrange, A. M., et al. 2013, A&A, 553, A60
  * Reffert et al. (2015) Reffert, S., Bergmann, C., Quirrenbach, A., Trifonov, T., & Künstler, A. 2015, A&A, 574, A116
  * Reggiani & Meyer (2013) Reggiani, M. & Meyer, M. R. 2013, A&A, 553, A124
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003
  * Rizzuto et al. (2013) Rizzuto, A. C., Ireland, M. J., Robertson, J. G., et al. 2013, MNRAS, 436, 1694
  * Röser et al. (2018) Röser, S., Schilbach, E., Goldman, B., et al. 2018, A&A, 614, A81
  * Samland et al. (2021) Samland, M., Bouwman, J., Hogg, D. W., et al. 2021, A&A, 646, A24
  * Savage & Panek (1974) Savage, B. D. & Panek, R. J. 1974, ApJ, 191, 659
  * Schlaufman (2018) Schlaufman, K. C. 2018, ApJ, 853, 37
  * Schneider et al. (2019) Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2019, Nature, 574, 211
  * Shobbrook (1983) Shobbrook, R. R. 1983, MNRAS, 205, 1215
  * Squicciarini et al. (2021) Squicciarini, V., Gratton, R., Bonavita, M., & Mesa, D. 2021, MNRAS, 507, 1381
  * Suri et al. (2021) Suri, S., Beuther, H., Gieser, C., et al. 2021, A&A, 655, A84
  * Vigan et al. (2021) Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72
  * Vigan et al. (2010) Vigan, A., Moutou, C., Langlois, M., et al. 2010, MNRAS, 407, 71
  * Wade et al. (2016) Wade, G. A., Neiner, C., Alecian, E., et al. 2016, MNRAS, 456, 2
  * Wagner et al. (2022) Wagner, K., Apai, D., Kasper, M., McClure, M., & Robberto, M. 2022, AJ, 163, 80
  * Wang et al. (2022) Wang, C., Langer, N., Schootemeijer, A., et al. 2022, Nature Astronomy [arXiv:2202.05552]
  * Wang & Chen (2019) Wang, S. & Chen, X. 2019, ApJ, 877, 116
  * Winn & Fabrycky (2015) Winn, J. N. & Fabrycky, D. C. 2015, ARA&A, 53, 409
  * Wolff (1990) Wolff, S. C. 1990, AJ, 100, 1994
  * Wolff et al. (2007) Wolff, S. C., Strom, S. E., Dror, D., & Venn, K. 2007, AJ, 133, 1092
  * Wolthoff et al. (2022) Wolthoff, V., Reffert, S., Quirrenbach, A., et al. 2022, arXiv e-prints, arXiv:2202.12800

##  Appendix A Optimization tests for stellar parameters

Our optimization tests for the parameters of $\mu^{2}$ Sco start from the
prior distributions of $(M,R,T_{\mathrm{eff}},\pi)$ described in Section 2. We
create a synthetic sample of $10^{7}$ stars, each one described by a
quadruplet $(M,R,T_{\mathrm{eff}},\pi)$, 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:

$\frac{L}{L_{\odot}}=\left(\frac{R}{R_{\odot}}\right)^{2}\cdot\left(\frac{T}{T_{\odot}}\right)^{4},$
(5)

we apply the interpolated bolometric correction as a function of
$T_{\mathrm{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\in[3.52,3.60]$, where the broadened error bar allows for some uncertainty
in the bolometric correction itself;

  * •

$\log{g}\in[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\leavevmode\nobreak\ M_{\odot}$.

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

$\log{L}=0.47+3.36\log{M}.$ (6)

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_{\mathrm{eff}},\pi)$ are shown in
Figure 15. 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\pm 0.3\leavevmode\nobreak\ R_{\odot}$), while the mass distribution
shifts to very high values ($M=10.0_{-0.9}^{+1.0}\leavevmode\nobreak\
M_{\odot}$).

Figure 15: Results of the first optimization test of stellar parameters.
Corner plot showing the posterior distribution of the quadruplets
$(M,R,T_{\mathrm{eff}},\pi)$ consistent with the observational filter for the
optimization test based on Pecaut & Mamajek (2013).

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^{7}$ 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)\in[3.6,4.0]$ as before;

  * •

$T_{\mathrm{eff}}\in[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\in[16,25]$ Myr or [Fe/H] $\in[0.10,0.20]$) a family of solutions appeared
(Figure 16-17); the posterior mass distribution is consistent with the result
of the first test ($M=8.8\pm 0.3\leavevmode\nobreak\ M_{\odot}$ for the run at
constant metallicity, $M=9.2\pm 0.2\leavevmode\nobreak\ M_{\odot}$ for the run
at constant age; we may combine the two to get $M=9.0\pm
0.3\leavevmode\nobreak\ M_{\odot}$). We average the two mass determinations
from the independent methods to a derive a final $M=9.1\pm
0.3\leavevmode\nobreak\ M_{\odot}$. In the same way, we derive $R=5.6\pm
0.2\leavevmode\nobreak\ R_{\odot}$; the median values of the posteriors of the
first test will be used as best-fit estimates for $T_{\mathrm{eff}}$ and
$\pi$. 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\pm 4$ Myr.

Figure 16: Results of the second optimization test of stellar parameters
(constant metallicity). Corner plot showing the posterior distribution of the
doublets $(M,\pi)$ consistent with the observational filter for the
optimization test based on PARSEC isochrones at fixed [Fe/H]=0.00. The
corresponding $T_{\mathrm{eff}}$, $R$ and age distribution are shown too.
Figure 17: Results of the second optimization test of stellar parameters
(constant age). Corner plot showing the posterior distribution of the doublets
$(M,\pi)$ consistent with the observational filter for the optimization test
based on PARSEC isochrones at fixed $t=16$ Myr. The corresponding
$T_{\mathrm{eff}}$, $R$ and [Fe/H] distribution are shown too.

##  Appendix B Possible binarity of $\mu^{2}$ Sco

The astrometric solution of $\mu^{2}$ Sco in the literature appears highly
problematic (Table 8). As shown by Kervella et al. (2019), a significant
deviation of proper motion components from the long-term motion reconstructed
by cross-matching Hipparcos 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.

Table 8: Astrometric solutions for $\mu^{2}$ Sco in the recent literature. The values adopted throughout this paper combine proper motions by Kervella et al. (2019) and the kinematic parallax based on membership to ELS. References: (1): Hipparcos (1997 reduction); (2): Hipparcos (2007 reduction); (3) Gaia DR2; (4) Gaia EDR3; (5) Kervella et al. (2019). Source | $\pi$ | $\mu_{\alpha}^{*}$ | $\mu_{\delta}$
---|---|---|---
| mas | mas yr-1 | mas yr-1
(1) | $6.31\pm 0.86$ | $-12.92\pm 0.66$ | $-23.80\pm 0.61$
(2) | $6.88\pm 0.12$ | $-11.09\pm 0.13$ | $-23.32\pm 0.11$
(3) | $7.92\pm 0.55$ | $-9.98\pm 0.96$ | $-19.87\pm 0.78$
(4) | $5.66\pm 0.28$ | $-12.11\pm 0.30$ | $-22.57\pm 0.27$
(5) | — | $-11.77\pm 0.02$ | $-23.11\pm 0.02$
adopted | $5.9\pm 0.2$ | $-11.77\pm 0.02$ | $-23.11\pm 0.02$

Although Kervella et al. (2019) reports a significant (S/N=4.1) proper motion
anomaly (PMA) at Hipparcos era, the 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-1) 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.

$\mu^{2}$ Sco is neither a visual binary (Gaia Collaboration et al. 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-1], is
$L_{X}=29.82$. Since the estimate assumes a distance $d\sim 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_{2}<0.7$
M⊙.

We then analyzed HARPS (Mayor et al. 2003) and FEROS (Kaufer et al. 1999)
radial velocity data, spanning about 10 years retrieved from the ESO Science
Archive191919Data available at http://archive.eso.org/scienceportal/home,
Programs 69.D-0677, 091.C-0713, 187.D-0917.: after extraction of RV from 37 He
I and atomic lines, we verified that the rms of the observations (0.21 km s-1)
is smaller than individual uncertainties ($\sim 0.3$ km s-1). 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
$\sim 50$ m s-1 yr-1 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.

Looking at Figure 18, 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.

Figure 18: Limits on the mass of a possible unresolved stellar companion
coming from several techniques. A $M>0.3-0.5\leavevmode\nobreak\ M_{\odot}$
companion is excluded at sub-au separations by RVs (95% confidence interval),
a $M>0.7\leavevmode\nobreak\ M_{\odot}$ 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 $\sim 290$ au, $\mu^{2}$ Sco b lies outside
the x range.

##  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 $\mu^{2}$ 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