© The Authors 2024

1 Introduction

Young exoplanets with a favourable brightness and separation to their host star can be directly imaged. However, such exoplanets can be confused with more distant background stars that happen to lie in the line of sight. A common way to distinguish these scenarios is to observe both the host star and candidate over time to look for a common proper motion and/or parallax. An example of doing this is shown in Fig. 1 (data from Janson et al. 2021b and published in Squicciarini et al. 2022). The positions of the host star and a number of exoplanet candidates were measured at two epochs, 2018 and 2021. The orange crosses show the measured change in position (relative to the host star) of candidates between these two epochs. The host star has some proper motion and parallax between 2018 and 2021: the black dashed line shows how objects with zero proper motion and parallax would move over this time period (as our view is centred on the host star). Those orange points clustered near the black point labelled 2021 are therefore consistent with being distant background stars. The candidate labelled ‘b’, on the other hand, has a motion more consistent with the host star and so is more likely to be a true companion.

To make this procedure quantitative, we must take into account the measurement uncertainties and any covariance between them. Background stars do not have zero parallax and proper motion, so we need a proper model for their motions too. There may also be more than two epochs, so we want to take into account all of the data within a single assessment. In the literature these confounding factors are usually not considered (Lagrange et al. 2010; Lafrenière et al. 2011; Carson et al. 2013; Kuzuhara et al. 2013; De Rosa et al. 2015; Konopacky et al. 2016; Chauvin et al. 2017; Keppler et al. 2018; Bohn et al. 2021; Janson et al. 2021a; Squicciarini et al. 2022; Franson et al. 2023; Mesa et al. 2023; Chomez et al. 2023). Such studies often compare the second epoch position with that expected for a stationary background star. If they deviate significantly, the null hypothesis of being a background star is rejected and some alternative model, that is, a companion, is implicitly accepted. This classical hypothesis-testing approach does not, however, assess whether the data might be even more unlikely under the companion model.

The goal of this paper is to put this astrometric confirmation approach on a solid statistical footing. We developed a model to evaluate whether the multi-epoch motion of the candidate is more likely to be a co-moving exoplanet or a coincidental field star. Our model is based on the proper motion and parallax distributions of field stars in the same area of the sky as the candidate, and with similar magnitudes to the candidate being tested. Using the odds ratio, we compare this background model with a model in which the candidate is co-moving with the host star. We applied our method to multi-epoch measurements of candidates observed in the B-Star Exoplanet Abundance Study (BEAST) survey (Janson et al. 2021b), which has an abundance of candidate objects due to observing close to the Galactic plane. A python package implementing our method is accessible via GitHub¹.

Fig. 1 Change in position of exoplanet candidates (orange crosses) relative to the star µ2 Sco between two measurement epochs. A co-moving source should be close to the origin (labelled ‘2018’). A background source with zero proper motion will move according to the dashed curve (a reflection of the host star’s parallax and proper motion) ending in the black star labelled ‘2021’. The motion of µ2 Sco b is distinct from the cloud of background stars in the field that are (through this plot) deemed not to be exoplanets. This figure has been adapted from Squicciarini et al. (2022).

2 Methods

We developed a probabilistic method that compares the likelihoods of position measurements of an exoplanet candidate under two models: the first assumes the candidate is a co-moving companion, and the second assumes it is a field star. The former uses the proper motion and parallax of the exoplanet candidate’s host star, and the latter additionally uses a magnitude-dependent fit to the parallax and proper motion distributions of field stars close to the host star’s line of sight. Our method can use an arbitrary number of epochs of astrometric observations to derive the likelihoods. We also implemented a special case of the general method that neglects the parallax. This is less realistic but is easier to visualize. The field star astrometric model is currently built on data from the Gaia (Gaia Collaboration 2023) and 2MASS (Two Micron All Sky Survey; Skrutskie et al. 2006) surveys, but could use other present or future surveys.

2.1 Positional model

We consider the position of a candidate object relative to the host star as a function of time. Assuming that the candidate does not move relative to the star, then in a geocentric coordinate system this can be described as a linear motion with a superimposed par-allactic motion. In a Cartesian plane projection, we write this as $$\begin{array}{l}\text{Δ}{x^{\prime }}_i=\left({x^{\prime }}_{a,0}-{x^{\prime }}_{\star ,0}\right)+\left({\mu }_{a,x}-{\mu }_{\star ,x}\right)t_i+\left({\varpi }_a-{\varpi }_{\star }\right)s_x\left(t_i\right)\\ \text{Δ}{y^{\prime }}_i=\left({y^{\prime }}_{a,0}-{y^{\prime }}_{\star ,0}\right)+\left({\mu }_{a,y}-{\mu }_{\star ,y}\right)t_i+\left({\varpi }_a-{\varpi }_{\star }\right)s_y\left(t_i\right),\end{array}$$(1)

where µ is the proper motion and ϖ is the parallax. Subscript a refers to the candidate object, subscript ⋆ to the host star, subscript 0 to the true position of the candidate at a reference epoch, and subscript i to the ith epoch relative to the reference epoch. The functions s_x and s_y are periodic phase factors for the parallax motion given by the orbital motion of Earth around the Sun at epoch i. The primed variables denote the true position of the candidate, while the unprimed variables denote the measured position of the candidate. In this paper ‘true positions’ denote those we would obtain in the absence of noise in our positional measurements at the epochs, but they are still based on some externally provided values of parallax and proper motion. Those external values are considered to be noisy, and their covariances are taken into account.

Equation (1) can be generalized to N measurements of the candidate’s position over time (∆x, ∆y) = {∆x_i, ∆y_i} at times t = {t_i} with i ∈ [1, N], where the bold font indicates a vector of measurements. The first position measurement defines the true position of the candidate in the first epoch (∆x₁, ∆y₁) = $\left(\text{Δ}{x^{\prime }}_1,\text{Δ}{y^{\prime }}_1\right)$ without loss of generality.

2.2 Overall probabilistic model

Given the measured positions of the exoplanet candidate over time, we computed the likelihood of the data under two models. The first model, denoted M_c, assumes the exoplanet candidate is a co-moving companion, sharing the same proper motion and parallax as the host star. We made the simplifying assumption that the orbital period is long compared to the observational baseline and so neglected the candidate’s orbital motion. This may not be justified in all cases, which we elaborate on in Sect. 4. The second model, denoted M_b, assumes the candidate to be a background object² with a proper motion and parallax distribution constructed from a set of background stars in a narrow field of view around the host star (described in Sect. 2.3). We then compared these two models via the odds (likelihood) ratio $$r_{c,b}=\frac{P\left(\text{Δ}x,\text{Δ}y\mid M_c\right)}{P\left(\text{Δ}x,\text{Δ}y\mid M_b\right)},$$(2)

which indicates which model favours the data more. This does not give the posterior probability of the model given the data, however, for which we would first need to establish model prior probabilities.

Let us first consider the denominator in Eq. (2), the likelihood under the background model. This can be computed as a combination of two probability density functions (PDFs). The first of these is the probability of the noisy position measurements given the true positions, P(∆x, ∆y | ∆x′, ∆y′), which reflects the noise in the determination of the centroid of a point spread function on the detector. We assumed this to be a Gaussian with mean (∆x′, ∆y′) and a covariance matrix Λ that reflects the accuracy of, and correlations in, the measurements. The second term is P(∆x′, ∆y′ | t, M_b). This represents the spread in possible true positions of a background star at specific times arising from the spread in the proper motion and parallax of the background star population. We assumed this term to be a Gaussian distribution with a mean given by the right side of Eq. (1). The terms in that equation, as well as the covariance matrix of the Gaussian (which we denote A′), come from our fit to the background star population (Sect. 2.3). These two PDFs we then combined via a marginalization to give the required likelihood $$\begin{array}{l}P\left(\text{Δ}x,\text{Δ}y\mid t,M_b\right)={\displaystyle {\int }_{\text{Δ}{x}^{\prime },\text{Δ}{y}^{\prime }}P}\left(\text{Δ}{x}^{\prime },\text{Δ}y\mid \text{Δ}{x}^{\prime },\text{Δ}{y}^{\prime }\right)\\ \times P\left(\text{Δ}{x}^{\prime },\text{Δ}{y}^{\prime }\mid t,M_b\right)\text{dΔ}{x}^{\prime }\text{\hspace{0.17em}dΔ}{y}^{\prime }.\end{array}$$(3)

This shows that the background model likelihood is a convolution of the two PDFs. Given that these are both Gaussian, the result of the convolution is also a Gaussian, with mean (Ax′, Ay′) and covariance matrix (Λ + Λ′). More details are provided in Appendix A.

We took the same approach to computing the likelihood of the data under the co-moving companion model M_c, the numerator in Eq. (2). Here, however, the relative position of the exoplanet candidate was assumed to be constant, meaning that the relative position measurements have zero variance and so P(Ax′, Ay′ | t, M_c) is simply a delta function. The convolution in Eq. (3) is then trivial, resulting in P(∆x, ∆y | t, M_c) being a Gaussian with mean (∆x′, ∆y′) and covariance matrix Λ.

Having now expressions for both likelihoods in terms of the measured relative positions and their covariances, and in terms of the parallax and proper motion distribution of the background model, we could compute the odds ratio in Eq. (2) and thus decide which model better explains the data.

We refer to the above general method as the proper motion and parallax covariance method. If the parallax is negligible, or if the observations are separated by almost exactly a year, then we can set the parallax in Eq. (1) to zero. If we also have only two epochs, then instead of using the two position measurements directly, we could convert them into a single proper motion. We call this the ‘proper motion-only’ method. In this special case, the ‘true’ proper motion $\text{Δ}{\mu }^{\prime }=\left(\text{Δ}{{\mu }^{\prime }}_x,\text{Δ}{{\mu }^{\prime }}_y\right)$ relative to the host star can be written $$\text{Δ}{{\mu }^{\prime }}_x=\frac{\text{Δ}{x^{\prime }}_2-\text{Δ}{x^{\prime }}_1}{t_2-t_1}$$(4)

which we compared to the measured relative proper motions from the data (same equation without the accents). The likelihoods for the two models are then $$P\left(\text{Δ}\mu \mid M_c\right)=\mathcal{N}\left(\left[\begin{array}{l}0\hfill \\ 0\hfill \end{array}\right],\left[\begin{array}{cc}\mathrm{Var}\left(\text{Δ}{\mu }_x\right)& \mathrm{Cov}\left(\text{Δ}{\mu }_x,\text{Δ}{\mu }_y\right)\\ \mathrm{Cov}\left(\text{Δ}{\mu }_x,\text{Δ}{\mu }_y\right)& \mathrm{Var}\left(\text{Δ}{\mu }_y\right)\end{array}\right]\right)$$(5)

and $$P\left(\text{Δ}\mu \mid M_b\right)={\displaystyle \int P}\left(\text{Δ}\mu \mid \text{Δ}{\mu }^{\prime }\right)P\left(\text{Δ}{\mu }^{\prime }\mid M_b\right)\text{dΔ}{\mu }^{\prime }.$$(6)

Both equations are Gaussian distributions. In the co-moving companion model the likelihood is centred at zero with a covariance matrix given by the uncertainties in position measurements. In the background model the likelihood is centred at the proper motion ∆µ′ of the background star distribution relative to the host star, with a covariance matrix given by the convolution of the proper motions of this background star population with the measurement uncertainties.

A more complete derivation of both the full model and this simpler (but less representative) proper motion-only model is given in Appendix A³.

Fig. 2 KS -band magnitude distribution of stellar objects in Gaia DR3 in the area of the sky within 0.3° of the 23 BEAST stars. Those in red have measured 2MASS KS-band photometry. Those in black stripes have KS predicted by a colour transformation.

2.3 Proper motion and parallax distribution model for field stars

In order to assess whether a candidate is astrometrically consistent with a given field star population (that is, to evaluate the likelihood P (∆x, ∆y | t, M_b)), we have to create a model for the parallax and proper motion distributions of the population. Ideally this model would be a function of the true distance and velocity of the star, but we of course do not have this information for an arbitrary candidate. However, we know that distance and velocity – and therefore parallax and proper motion – depend on the measurable direction in the Galaxy (see, for example, Bailer-Jones 2023). We can also condition our model on other relevant measurements, the most obvious being magnitude, as this contains some information about distance and stellar population. Ideally we would also use colour (as a crude proxy of mass and age and thus having velocity dependence), but this is often not available for many exoplanet surveys.

We built our model empirically using Gaia Data Release 3 (DR3; Gaia Collaboration 2023). We selected stars that are near to the candidate in Galactic coordinates and that have similar brightness. The latter can be hard to achieve because Gaia observes in the optical and is not as deep as infrared exoplanet surveys. To address this we used a positional cross-match of Gaia with 2MASS (Skrutskie et al. 2006) to assign infrared magnitudes to Gaia stars where possible. For those Gaia objects that did not have a 2MASS counterpart, we assigned synthetic KS magnitudes using the Gaia colour-colour transformation based on the G-band magnitude and G_BP – G_RP colour from Riello et al. (2021). While this transformation is not perfect, it significantly increased our field star sample size. Figure 2 shows in red all bright field stars covered by 2MASS within 0.3° of the 23 target stars of the BEAST. The distribution of stars for which we calculated synthetic K_S -band magnitudes via the colour transformations are shown with black stripes.

Using our sample, we then build a smooth model of the astro-metric distributions as a function of magnitude. This allows us to evaluate (if necessary by extrapolation) the astrometric model at the candidate’s brightness. To build the model we first group the data into magnitude bins, then fit a two-dimensional Gaussian distribution to the proper motions and parallaxes in each bin. An example of one magnitude bin is shown in Fig. 3. We then fit a linear function to the mean values of the two-dimensional (2D) Gaussian as a function of magnitude (Fig. 4), excluding points that lie outside the 10th to 90th percentile range of the magnitudes to achieve a more robust fit. This fit smooths out the variations and allowed us to extrapolate the model to fainter candidates than are in Gaia. The standard deviations as a function of magnitude were fit as an exponential $${\sigma }_{\text{fit }}(m)={\sigma }_{\min }+a\exp \left(-b\left(m-m_0\right)\right)$$(7)

with a > 0. This fit has a minimum value of σ_min, a fixed constant. The magnitude m₀ is the mean of the observed magnitudes. An example of this fit is shown in Fig. 5. A linear fit was used instead if the exponential fit residuals were larger than those of a linear fit. In making this linear fit we prevented σ_fit (m) from becoming smaller than 1 mas yr⁻¹.

We did not identify any definitive and consistent trend in the correlation between parallax and proper motion. We therefore opted to fit a constant value for the correlation instead.

Fig. 3 Proper motion distribution of stellar objects with a KS -magnitude between 18 and 19 in a 0.3° sky area around HIP 82545 (µ2 Sco). The elliptical contours are the boundaries that encompass 50, 90, and 99% of the stellar objects.

Fig. 4 Variation of the mean of the proper motion distribution in field stars as a function on stellar magnitude. Each bin includes 200 stellar objects in a 0.3° sky area around HIP 82545 (µ2 Sco). Each point is the mean value of the 2D Gaussian fit in each magnitude bin. Only those points lying within the 10th–90th percentile range of magnitudes were used for the linear fit.

Fig. 5 Variation of the standard deviation and correlation of the proper motion distributions for the field stars around HIP 82545 (µ2 Sco) as a function of magnitude. Each point comes from a 2D Gaussian fit over a narrow magnitude bin (shown in Fig. 3). Each bin includes 200 stars in a 0.3° sky area around the target star. Crosses denote fits using stars with 2MASS magnitudes, and circles denote those with magnitudes computed from the Gaia colour transformation.

2.4 Verification of the method using simulations

A verification of our model would ideally be based on a large set of real data in which the true nature of each candidate is known, but this was not available. We therefore used simulations to assess the reliability of the odds ratio. We simulated 1000 trajectories of a co-moving companion and 1000 trajectories of background model objects.

Starting with the real star µ² Sco, we constructed a background model for proper motion and parallax. The initial position of the candidate relative to µ² Sco was randomly generated. The position at subsequent epochs (one year apart) is given by $$x_{i+1}=x_i+\left({\mu }_x+\mathcal{N}\left(0,\sigma =3{\mathrm{mas }}_{{\text{yr}}^{-1}}\right)\right)\text{Δ}t$$(8)

and similarly for y. This equation represents either a co-moving companion with a zero (relative) proper motion (µ_x = µ_y = 0), or a relative proper motion of the background model to µ² Sco (${\mu }_x={\mu }_{\text{RA}}(m)-{\mu }_{x,{\mu }^2}\text{Sco}$, where µ_RA(m) is a fixed value defined by the brightness of the candidate), both with a normal distribution 𝒩 to simulate the propagation of the proper motion uncertainty. We chose a standard deviation of 3 mas yr⁻¹. We simulated 1000 paths the candidate could take under the companion model M_c using Eq. (8) with µ = 0, and another 1000 paths under the background object model M_b using a single fixed proper motion for a background object with a similar brightness to the candidate. The spread seen in Fig. 6 thus arises from the noise we add at each epoch. Our method distinguishes between those paths via a logarithmic odds ratios greater than zero (companion model favoured, show in blue) and smaller than zero (background model favoured, shown in black). We found that all 2000 trajectories were correctly identified as coming from the model from which they were generated.

Fig. 6 Simulated test of our method by propagating the positions of 2000 simulated candidates to µ2 Sco over four epochs. Half are propagated according to the co-moving model, which just adds zero mean noise at each epoch. The other half are propagated according to the proper motion of background stars, plus noise. The colours denote the logarithmic odds ratio that our method computes for each trajectory. All 2000 are correctly identified.

3 Results

3.1 General results

BEAST (Janson et al. 2021b) employs the VLT (Very Large Telescope)/SPHERE (Spectro-Polarimetric High-contrast Exo-planet REsearch) Beuzit et al. (2019) extreme adaptive optics instrument to conduct a direct imaging survey of 85 B-type stars in the young (5–20 Myr) and relatively nearby (120–150 pc) Scorpius-Centaurus (Sco-Cen) region (Fig. 7). The targets have similar distance and G-band magnitude, but vary in colour. By targeting young B-type stars, their exoplanets should be relatively bright in the near-infrared. Before this survey, B-stars had not been systematically surveyed for exoplanets. Radial velocity surveys, which are more sensitive to close-in planets, reveal comparatively few planets around massive stars (Reffert et al. 2015). Thus BEAST addresses the question of whether massive stars can form massive planets at larger separations. Certain formation scenarios, such as disk instability, might primarily take place in the outermost areas of extensive disks surrounding massive stars (Helled et al. 2014). Sco-Cen is located close to the Galactic plane, so observations often include many spurious candidates in the field-of-view. This makes the survey a prime target for testing our method of distinguishing bound companions from faint field stars.

Of the 85 stars in BEAST, 23 have candidates with at least two epochs of observation (as of 2022-04-01). These 23 stars have in total 263 candidates. The projected separations between the candidates and their host stars range from 49 to 1457 au (Fig. 8).

We applied both our proper motion-only model and our proper motion and parallax covariance model to all 263 candidates. The results are shown in Table 1, sorted by descending logarithmic odds ratio log₁₀ r_cb(µ, ϖ), that is, the most likely real companions appear first. In what follows we focus on the results using the parallax and proper motion model in the final column. Visualizations of all 20 candidates are shown later. The complete table of candidate astrometry and photometry from the BEAST survey will be published in Delorme et al. (in prep.).

Most interesting are the seven candidates that have an odds ratio greater than ten (log₁₀ r_cb(µ, ϖ) > 1). Two of these, b Cen b and µ² Sco b, are well-established as exoplanets (Janson et al. 2021b; Squicciarini et al. 2022). We further infer that the candidates of HIP 81208 (two candidates), HIP 61257, HIP 52742, and HIP 76048 are likely true companions (but not necessarily in the exoplanet mass regime). The remaining candidates in Table 1 are still formally favoured by the co-moving companion model, but not by much. These will require more epochs or other observations in order to determine their nature. We now examine the results on some of the individual, high odds ratio candidates, which also serve to illustrate how our method works.

Fig. 7 Colour-magnitude diagram of all BEAST stars, except for η Cen at (BP – RP, G) = (4.75, 2.25) mag.

Fig. 8 Distribution of the separations of all candidate from all targets identified in BEAST. This assumes the candidates are at the same distance as their respective target stars.

3.2 Results on individual candidates

µ² Sco b. The astrometric field star model for the target µ² Sco is shown in Fig. 9. The black line is our fit to the background stars (described in Sect. 2.3), which can be compared to the various BEAST candidates in this field in orange. Nearly all of the candidates agree with the field star model, showing the large degree of field star contamination that can be present in exoplanet searches. Just one does not agree; this is the exoplanet µ² Sco b. The PDF over the proper motion (the likelihood) for this candidate under the background model is shown by the black curves and contours in Fig. 10. The likelihood for the companion model is shown by the blue curves and contours, and the measured proper motion is shown as the orange point and line. In this example we see that the measurement is far more consistent with the companion model than with the background model. This is a rather clear-cut case even by visual inspection; many other cases in Appendix B are more ambiguous.

b Cen b. This candidate, which was identified as an exo-planet by Janson et al. (2021a), has three epochs in our analysis.

Figure 11 shows the change in positions of the candidate under the background model M_b (black), the (unchanged) position of the candidate under the co-moving companion model M_c (blue), and the measured positions of the candidate (orange). The first observation of this candidate was at the epoch 2000.4 (circle) based on archival images by Shatsky & Tokovinin (2002). In the relative reference frame to the host star, a co-moving companion would not move from its 2000.4 position in our model, and so is not shown for other epochs. The next observation of the candidate took place in 2019.2 with BEAST. Based on the background proper motion and parallax distribution, as well as the host star’s proper motion, a background object in the same area of the sky as the host star would have moved about 850 mas since 2000.4. However, the figure shows that by 2019.2 the candidate has moved much less (orange point). Another epoch at 2021.3 confirms that the motion of the candidate is much closer to the host star’s motion than it is to that of the background population. The background star hypothesis is clearly ruled out in this case. But because the measured positions lie near the 99% contour of our co-moving companion model, some residual motion compared to the co-moving companion model likely remains. This can be explained by orbital motion - which is not taken into account – over the 21-yr observational baseline (Janson et al. 2021a).

HIP 61257 ‘B’. This highly probable candidate of HIP 61257 can be further examined by including non-BEAST data. Kouwenhoven et al. (2005) discuss a potential companion around HIP 61257, but they ultimately classify it as a background object. This potential companion has a K_S -band magnitude of 12.43 mag and a separation of 5540 mas, which coincides with the magnitude and separation of one of our candidates with logarithmic odds ratio greater than zero in both modelling frameworks. The astrometric motion plot (Fig. 12) shows that this candidate is co-moving with the host star. As the astrometric data from Kouwenhoven et al. (2005) are reported without separation and position angle uncertainties, we adopted an uncertainty of 10 mas in the relative right ascension and declination. This companion was also discussed by Gratton et al. (2023), who identified it to be a low-mass star (0.083 ± 0.01 M_⊙) based on its Gaia and K-band magnitudes. The logarithmic odds ratios without the archival data are log r_cb(µ, ϖ) = 12.91 and log r_cb(µ) = 2.45. We therefore confirm this candidate as a genuine binary companion based on its astrometry.

HIP 81208 B and C. These two candidates are identified as co-moving companions by Viswanath et al. (2023). They identify HIP 81208 B as a 67 MJup object, mostly likely a brown dwarf, and HIP 81208 C as a 0.135 M_⊙ low-mass star. Our results for these are shown in Figs. 13a,b (using astrometry and photometry from Table C.1 in Viswanath et al. 2023). Both candidates are favoured by the co-moving companion model with our more sophisticated parallax and proper motion method. The proper motion-only model rejects HIP 81208 B but accepts HIP 81208 C. This system was recently analysed further and found to be a gravitationally bound hierarchical quadruple system comprized of low mass objects with a newly discovered companion to the C component (Cb; Chomez et al. 2023).

HIP 52742 x9ld1uh0 and mm9uw3ha. Two candidates with magnitudes of K_S = 12.03 mag and K_S = 19.49 mag are formally favoured by our companion model. The brighter one we identify as a companion with high significance (log odds ratio 12.6). Gratton et al. (2023) also identified this as a companion: Adopting an age of 82.5 Myr (obtained from assuming membership of a comoving group outside Sco-Cen; Janson et al. 2021b) they estimated a mass of 0.51 M_⊙ at a projected separation of 176 au. In a search for astrometric acceleration from a comparison of Hipparcos and Gaia DR3 proper motions, Brandt (2021) found marginally significant evidence for an acceleration of this star, which might be evidence of this companion. The fainter candidate companion to this star that we identify has a log odds ratio of just 0.16, which is not significant.

HIP 76048 vslvj1zp. We identify one potential companion to this star with an odds ratio of 100 in our analysis, albeit with a very short baseline of just 0.2 yr. Brandt (2021) found no significant evidence of acceleration from their HIPPARCOS-Gaia DR3 proper motion study.

Fig. 9 Proper motion model of the background stars for µ2 Sco for the right ascension direction (top) and the declination direction (bottom) in the International Celestial Reference System (ICRS). This is similar to Fig. 4, but now includes the BEAST candidates shown in orange and transformed to ICRS.

Fig. 10 Demonstration of the proper motion-only method for the known exoplanet µ2 Sco b. The likelihood under the background model is shown in black and the likelihood under the co-moving companion model is shown in blue. The bottom left panel shows these as two-dimensional (Gaussian) distributions. The other two panels show the one-dimensional marginal distributions. The measured relative proper motion is shown in orange: the uncertainties in this is not shown because it is included in the two likelihoods (see Eqs. (5) and (6)).

Fig. 11 Visualization of the predicted positions of the candidate companion b of the star HIP 71865 (b Cen) under the proper motion and parallax model. A co-moving companion would remain at the position of the first epoch (blue circle) because orbital motion is not included. A field star with the modelled proper motion and parallax of nearby (mostly background) stars would be measured at the two later epochs at the two positions shown by the black triangles. The actual measured change in positions of the candidate are shown as orange triangles. The fact that these are much nearer to the blue distribution means this is likely to be a true companion, something that is properly quantified by our method. The contour lines show 50, 90, and 99% of the enclosed probability, reflecting the propagated uncertainty in the parallaxes, proper motions, and BEAST position measurements. The marginal likelihoods are shown on both axes. This visualization does not show the covariances between the measurements at different epochs, which are nonetheless taken into account by our method (see Eq. (A.4)).

Fig. 12 Results of our model applied to the position measurements of HIP 61257 ‘B’ from Kouwenhoven et al. (2005) and Janson et al. (2021b) over a 17-yr baseline. See the caption to Fig. 11 for a description.

4 Discussion

4.1 Interpretation

Our analysis favours the co-moving companion model for five candidates from the BEAST data set accessible for this study, as well as two candidates in the HIP 81208 system that have a second epoch from Viswanath et al. (2023). Of these seven candidates, two are confirmed exoplanets: b Cen AB b (Janson et al. 2021a) and µ² Sco b (Squicciarini et al. 2022). HIP 61257 ‘B’ is very likely a stellar-mass companion (Gratton et al. 2023). The two candidates of HIP 81208 are discussed as being a brown dwarf and a low-mass stellar companion by Viswanath et al. (2023). The remaining two companions of the host stars HIP 52742 and HIP 76048 - with a measurement baseline of only 0.2 yr – are unconfirmed at the time of writing.

Most of the targets in our analysis have only two observation epochs with small temporal baselines (of order one year), so many candidates show little motion relative to the host star between the epochs. It is precisely these cases where ad hoc approaches to assessing companionship are inconclusive, and our statistical framework is most useful. For many of these short baselines, the co-moving model is not favoured. Longer temporal baselines make it easier to distinguish between the models, especially as the targets of direct imaging campaigns tend to be nearby stars with large proper motion.

4.2 Model assumptions

Our model does not currently take into account orbital motion, so the co-moving companion model may not be favoured even if the candidate is a true companion with significant orbital motion over the observational baseline. Orbital periods for directly imaged planets tend to be long, of the order of centuries or even millennia. In those cases, neglecting the orbital motion does not affect our analysis. Planets with relatively short orbits, in contrast, such as β Pic b and c with periods of 23.6 and 3.3 yr respectively (Lacour et al. 2021), can show significant orbital motion.

If orbital motion can be directly observed then the candidate is very likely to be confirmed as a companion (for example Marois et al. 2008). This could be included as a relatively straight forward extension to our modelling approach by including a model with path curvature. Plausible paths (priors) could be generated from a Keplerian model (using the model of Blunt et al. 2020, for example) and then marginalizing over them to determine the model likelihood. An alternative approach is to show that a candidate’s motion is consistent with that of field stars while also showing that the range of orbits that can explain the data are beyond the escape velocity of the system (see for example Nielsen et al. 2017; Wagner et al. 2022).

We built our field star astrometric model using Gaia parallaxes and proper motions. On account of the limited depth of Gaia, we then had to extrapolate our model to the fainter magnitudes of our candidates. This can introduce biases, and thus incorrectly favour or disfavour the background model.

Another drawback of Gaia is that it observes in the optical, whereas direct imaging surveys for planets are currently done mostly in the near-infrared. To obtain the necessary infrared magnitudes of the Gaia sources we had to crossmatch to an infrared survey. We chose 2MASS because it is all-sky. However, it is not very deep, so for many field stars we instead had to predict the K_S -band magnitude using colour transformations from Gaia. We used the transformations of Riello et al. (2021) for objects in the colour range −0.5 < G_BP − G_RP < 2.5. The accuracy of this for brighter stars around µ² Sco is demonstrated in Fig. 14. The 2MASS objects span a K_S -band range of only 10 to 18 mag while the transformed photometry extends to 21 mag. Some exoplanet surveys are also conducted at shorter wavelengths, such as J and H, for which colour transformations are more reliable as they are closer to the observed Gaia bands. Ideal for our purposes, of course, would be an infrared version of Gaia.

As explained in Sect. 2.3, in our background model the parallax and proper motion distributions depend only on direction and magnitude. A further improvement would be to add an additional dependence on the colour, if that is provided by the imaging survey. We could then also use the measured colour to infer something about the intrinsic properties of the object, such as its spectral type (Parviainen et al. 2019). Although that may help determine whether or not it is a low-mass object, it would not tell us whether it is gravitationally bound.

In computing the odds ratio we only considered kinematics. We do not take into account the number density of background stars or the angular separation between the candidate and target star. Yet for a given separation, the less dense the background, the more likely the candidate is a genuine companion (Tamura 2016; Squicciarini et al. 2022). This information could be incorporated as an additional multiplicative odds ratio, although it requires a model (or measurement) of the stellar density at the faint magnitudes at which we conduct our exoplanet survey.

Finally, we emphasize that we report our results as a ratio of likelihoods of two models, where each likelihood is the probability of the data given the model. To convert each likelihood into a posterior probability of the model given the data, we would need to adopt a prior probability for the model, and know that our models are exhaustive. The latter is not yet the case, as we have neglected orbital motion, for example. The prior could incorporate the direction-dependent variation of the background star number density.

Fig. 13 Astrometric motion of the two candidates of HIP 81208 with an observation baseline of three years. Viswanath et al. (2023) report both objects as co-moving companions based on their proper motion analysis. Our analysis supports this claim.

Fig. 14 Difference between the KS magnitude calculated from Gaia via colour transformations and the and measured 2MASS KS magnitudes for 934 023 objects from Gaia with 2MASS counterparts in the vicinity of µ2 Sco.

5 Summary

This work has introduced a statistical method that uses multi-epoch astrometry of an imaged exoplanet candidate to compare a co-moving companion model with a chance-aligned field star model. It puts what is commonly referred to as the ‘common proper motion test’ on a probabilistic footing.

Our statistical model enables a quantitative analysis of an arbitrary number of epochs, a task that cannot be achieved effectively through visual inspection. We consider the proper motion and parallax of the host star and the candidate and evaluate the likelihoods under two different models for the candidate: one in which it is a co-moving companion with negligible orbital motion, the other in which it is a member of the field star population. For the latter we built a probabilistic model of the distribution of the proper motions and parallaxes of field stars as a function of magnitude, using a fit to Gaia data in the field of each target star.

We applied our method to a sample of 263 candidates around 23 stars from the B-Star Exoplanet Abundance Study. We first developed a purely proper motion based method, which we then extended to take into account the parallax. This model accommodates the covariance in the astrometry both between the measurements and across multiple epochs, for both Gaia astrom-etry and the direct measurements. We identify seven candidates as co-moving companions. Five of these have been identified as real companions in the literature, including the two exoplanets µ² Sco b and b Cen(AB) b. The remaining two candidates are priority targets for further investigation.

Our modelling approach is publicly available as an open-source Python package on GitHub⁴, allowing for easy evaluation and visualization of existing and new data. While this work presents an improvement over current practices, there is scope for further improvement. Most useful would be the inclusion of exoplanet orbital motion in the companion model, the incorporation of stellar number densities, and discriminating field stars from exoplanets based on their spectral information.

Acknowledgements

We are grateful for the helpful comments from the referee that improved the paper. This study is based on observations acquired at the European Southern Observatory (ESO) VLT telescope (program 1101.C-025). This publication makes use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://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 publication also makes use of data 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.

Appendix A Mathematical description of the multi-epoch probabilistic model

We describe here in more detail the probabilistic model introduced in Sect. 2.2. In high-contrast imaging, the set of positions of candidates are generally determined relative to the host star, namely $$\begin{array}{l}\text{Δ}{x}^{\prime }=x_a^{\prime }-x_{\star }^{\prime }\\ \text{Δ}{y}^{\prime }=y_a^{\prime }-y_{\star }^{\prime }\end{array}$$(A.1)

where the subscript a refers to the candidate and ⋆ to the host star. The accents on the coordinates ${x^{\prime }}_a$ and ${x^{\prime }}_{\star }$ indicate that they are ‘true’, that is, without positional measurement uncertainties. The noisy, observed versions are x_a and x_⋆ at the (noise-free) times t. The times we use in our models are measured relative to the Gaia epoch of observation of the host star, because we use Gaia data to build our proper motion and parallax distribution models for the field stars.

As the true positions are nonetheless predictions based on noisy parallaxes and proper motions (and propagated in time using Eq. 1), we write the probability density functions of these true positions as $$P\left(\text{Δ}{x}^{\prime },\text{Δ}{y}^{\prime }\mid t,M(\varpi ,\mu )\right)$$(A.2)

where M can be M_c or M_b and the PDF can be written as Gaussians with the variances and covariances between the two positions ∆x′ and ∆y′. For the companion model, these (co)variances are derived from the Gaia data for the host star. For the background model, these (co)variances come from our fit to the distribution of the Gaia parallaxes and proper motions of a set of field stars, as described in Sect. 2.3. The likelihoods for both of these model (which appear in the odds ratio of Eq. 2) are obtained by marginalizing over the unknown true positions (∆x′, ∆y′) predicted by the model, as follows $$\begin{array}{l}P(\text{Δ}x,\text{Δ}y\mid t,M(\varpi ,\mu ))\\ ={\displaystyle \int P}\left(\text{Δ}x,\text{Δ}y,\text{Δ}{x}^{\prime },\text{Δ}{y}^{\prime }\mid t,M(\varpi ,\mu )\right)d\text{Δ}{x}^{\prime }d\text{Δ}{y}^{\prime }\\ ={\displaystyle \int P}\left(\text{Δ}x,\text{Δ}y\mid \text{Δ}{x}^{\prime },\text{Δ}{y}^{\prime }\right)\times \\ P\left(\text{Δ}{x}^{\prime },\text{Δ}{y}^{\prime }\mid t,M(\varpi ,\mu )\right)d\text{Δ}{x}^{\prime }d\text{Δ}{y}^{\prime }.\end{array}$$(A.3)

P(∆x, ∆y | ∆x′, ∆y′) is the likelihood of the measured positions given the true positions, and so reflects the noise in the determination of the image centroid. We take this to be a Gaussian with mean (∆x′, ∆y′) and 2N-dimensional Gaussian Λ, where N is the number of measurement epochs. Equation A.3 is a convolution of two Gaussians, which results in another Gaussian.

In the background model M_b, P (∆x′, ∆y′ | t_i, M_b) is a Gaussian with mean (∆x′, ∆y′) and covariance matrix $${\Lambda }^{\prime }=\left[\begin{array}{ccccc}\mathrm{Var}\left(\text{Δ}{x}_1^{\prime }\right)& \mathrm{Cov}\left(\text{Δ}{x}_1^{\prime },\text{Δ}{y}_1^{\prime }\right)& \mathrm{Cov}\left(\text{Δ}{x}_1^{\prime },\text{Δ}{x}_2^{\prime }\right)& \mathrm{Cov}\left(\text{Δ}{x}_1^{\prime },\text{Δ}{y}_2^{\prime }\right)& \cdots \\ & \mathrm{Var}\left(\text{Δ}{y}_1^{\prime }\right)& \mathrm{Cov}\left(\text{Δ}{y}_1^{\prime },\text{Δ}{x}_2^{\prime }\right)& \mathrm{Cov}\left(\text{Δ}{y}_1^{\prime },\text{Δ}{y}_2^{\prime }\right)& \dots \\ & & \mathrm{Var}\left(\text{Δ}{x}_2^{\prime }\right)& \mathrm{Cov}\left(\text{Δ}{x}_2^{\prime },\text{Δ}{y}_2^{\prime }\right)& \cdots \\ & & & \mathrm{Var}\left(\text{Δ}{y}_2^{\prime }\right)& \dots \\ ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]$$(A.4)

where terms are given in the order $\left(\text{Δ}{x^{\prime }}_1,\text{Δ}{y^{\prime }}_1,\text{Δ}{x^{\prime }}_2,\text{Δ}{y^{\prime }}_2,\dots \right)$. The symmetry of the covariance matrix fills the lower triangle of Eq. A.4. The individual terms of Λ′ are computed from Eq. 1 to be $$\begin{array}{l}\mathrm{Var}\left(\text{Δ}{x}_i^{\prime }\right)=t_i^2\left[\mathrm{Var}\left({\mu }_{b,x}\right)+\mathrm{Var}\left({\mu }_{\star ,x}\right)\right]\\ +s_x{\left(t_i\right)}^2\left[\mathrm{Var}\left({\varpi }_b\right)+\mathrm{Var}\left({\varpi }_{\star }\right)\right]\\ +2t_i{s}_x\left(t_i\right)\left[\mathrm{Cov}\left({\varpi }_b,{\mu }_{b,x}\right)+\mathrm{Cov}\left({\varpi }_{\star },{\mu }_{\star ,x}\right)\right]\\ \mathrm{Var}\left(\text{Δ}{y}_i^{\prime }\right)=t_i^2\left[\mathrm{Var}\left({\mu }_{b,y}\right)+\mathrm{Var}\left({\mu }_{\star ,y}\right)\right]\\ +s_y{\left(t_i\right)}^2\left[\mathrm{Var}\left({\varpi }_b\right)+\mathrm{Var}\left({\varpi }_{\star }\right)\right]\\ +2t_i{s}_y\left(t_i\right)\left[\mathrm{Cov}\left({\varpi }_b,{\mu }_{b,y}\right)+\mathrm{Cov}\left({\varpi }_{\star },{\mu }_{\star ,y}\right)\right]\\ \mathrm{Cov}\left(\text{Δ}{x}_i^{\prime },\text{Δ}{y}_i^{\prime }\right)=t_i^2\left[\mathrm{Cov}\left({\mu }_{b,x},{\mu }_{b,y}\right)+\mathrm{Cov}\left({\mu }_{\star ,x},{\mu }_{\star ,y}\right)\right]\\ +s_x\left(t_i\right)s_y\left(t_i\right)\left[\mathrm{Var}\left({\varpi }_b\right)+\mathrm{Var}\left({\varpi }_{\star }\right)\right]\\ +t_i{s}_y\left(t_i\right)\left[\mathrm{Cov}\left({\varpi }_b,{\mu }_{b,x}\right)+\mathrm{Cov}\left({\varpi }_{\star },{\mu }_{\star ,x}\right)\right]\\ +t_i{s}_x\left(t_i\right)\left[\mathrm{Cov}\left({\varpi }_b,{\mu }_{b,y}\right)+\mathrm{Cov}\left({\varpi }_{\star },{\mu }_{\star ,y}\right)\right]\\ \mathrm{Cov}\left(\text{Δ}{x}_i^{\prime },\text{Δ}{x}_j^{\prime }\right)=t_i{t}_j\left[\mathrm{Var}\left({\mu }_{b,x}\right)+\mathrm{Var}\left({\mu }_{\star ,x}\right)\right]\\ +t_i{s}_x\left(t_j\right)\left[\mathrm{Cov}\left({\varpi }_b,{\mu }_{b,x}\right)+\mathrm{Cov}\left({\varpi }_{\star },{\mu }_{\star ,x}\right)\right]\\ +t_j{s}_x\left(t_i\right)\left[\mathrm{Cov}\left({\varpi }_b,{\mu }_{b,x}\right)+\mathrm{Cov}\left({\varpi }_{\star },{\mu }_{\star ,x}\right)\right]\\ +s_x\left(t_i\right)s_x\left(t_j\right)\left[\mathrm{Var}\left({\varpi }_b\right)+\mathrm{Var}\left({\varpi }_{\star }\right)\right]\\ \mathrm{Cov}\left(\text{Δ}{y}_i^{\prime },\text{Δ}{y}_j^{\prime }\right)=t_i{t}_j\left[\mathrm{Var}\left({\mu }_{b,y}\right)+\mathrm{Var}\left({\mu }_{\star ,y}\right)\right]\\ +t_i{s}_y\left(t_j\right)\left[\mathrm{Cov}\left({\varpi }_b,{\mu }_{b,y}\right)+\mathrm{Cov}\left({\varpi }_{\star },{\mu }_{\star ,y}\right)\right]\\ +t_j{s}_y\left(t_i\right)\left[\mathrm{Cov}\left({\varpi }_b,{\mu }_{b,y}\right)+\mathrm{Cov}\left({\varpi }_{\star },{\mu }_{\star ,y}\right)\right]\\ +s_y\left(t_i\right)s_y\left(t_j\right)\left[\mathrm{Var}\left({\varpi }_b\right)+\mathrm{Var}\left({\varpi }_{\star }\right)\right].\\ \mathrm{Cov}\left(\text{Δ}{x}_i^{\prime },\text{Δ}{y}_j^{\prime }\right)=t_i{t}_j\left[\mathrm{Cov}\left({\mu }_{b,x},{\mu }_{b,y}\right)+\mathrm{Cov}\left({\mu }_{\star ,x},{\mu }_{\star ,y}\right)\right]\\ +t_i{s}_y\left(t_j\right)\left[\mathrm{Cov}\left({\varpi }_b,{\mu }_{b,x}\right)+\mathrm{Cov}\left({\varpi }_{\star },{\mu }_{\star ,x}\right)\right]\\ +t_j{s}_x\left(t_i\right)\left[\mathrm{Cov}\left({\varpi }_b,{\mu }_{b,y}\right)+\mathrm{Cov}\left({\varpi }_{\star },{\mu }_{\star ,y}\right)\right]\\ +s_x\left(t_i\right)s_y\left(t_j\right)\left[\mathrm{Var}\left({\varpi }_b\right)+\mathrm{Var}\left({\varpi }_{\star }\right)\right].\end{array}$$(A.5)

From the convolution in Eq. A.3, we see that P (∆x, ∆y | t, M_b) is a 2N-dimensional Gaussian with mean (∆x′, ∆y′) and covariance matrix Λ + Λ′. Equations A.5 assume that there is no covariance between the field star proper motion and parallax and the host star proper motion and parallax. This is not strictly true, because both the background model and the host star data are drawn from Gaia, and Gaia shows correlations between sources separated by small angles on the sky. But in practice our background model describes the distribution of an ensemble of stars, and so is less affected by Gaia’s correlations between individual sources.

In the co-moving companion model M_c, the position of the companion relative to the star is constant, leading to a zero covariance for the second term under the integral in Eq. A.3 making P (∆x′, ∆y′ | t, M_c) a delta function. P (∆x, ∆y | t, M_c) is therefore a 2N-dimensional Gaussian with mean (∆x′, ∆y′) and covariance Λ.

Appendix B Candidates with positive log odds ratios

Figures B.1, B.2, and B.3 show the modelling results for the 18 BEAST-only candidates that have logarithmic odds ratio log r_tcb > 0. The results for the two targets with additional data are shown in Fig. 13.

Fig. B.1 Candidates from BEAST with a logarithmic odds ratio log rtcb > 0. The schematics of the plots are explained in Fig. 11.

Fig. B.2 Candidates from BEAST with a logarithmic odds ratio log rtcb > 0. The schematics of the plots are explained in Fig. 11.

Fig. B.3 Candidates from BEAST with a logarithmic odds ratio log rtcb > 0. The schematics of the plots are explained in Fig. 11.

References

All Tables

All Figures

Fig. 1 Change in position of exoplanet candidates (orange crosses) relative to the star µ2 Sco between two measurement epochs. A co-moving source should be close to the origin (labelled ‘2018’). A background source with zero proper motion will move according to the dashed curve (a reflection of the host star’s parallax and proper motion) ending in the black star labelled ‘2021’. The motion of µ2 Sco b is distinct from the cloud of background stars in the field that are (through this plot) deemed not to be exoplanets. This figure has been adapted from Squicciarini et al. (2022).

In the text

	Fig. 2 KS -band magnitude distribution of stellar objects in Gaia DR3 in the area of the sky within 0.3° of the 23 BEAST stars. Those in red have measured 2MASS KS-band photometry. Those in black stripes have KS predicted by a colour transformation.
In the text

	Fig. 3 Proper motion distribution of stellar objects with a KS -magnitude between 18 and 19 in a 0.3° sky area around HIP 82545 (µ2 Sco). The elliptical contours are the boundaries that encompass 50, 90, and 99% of the stellar objects.
In the text

	Fig. 4 Variation of the mean of the proper motion distribution in field stars as a function on stellar magnitude. Each bin includes 200 stellar objects in a 0.3° sky area around HIP 82545 (µ2 Sco). Each point is the mean value of the 2D Gaussian fit in each magnitude bin. Only those points lying within the 10th–90th percentile range of magnitudes were used for the linear fit.
In the text

Fig. 5 Variation of the standard deviation and correlation of the proper motion distributions for the field stars around HIP 82545 (µ2 Sco) as a function of magnitude. Each point comes from a 2D Gaussian fit over a narrow magnitude bin (shown in Fig. 3). Each bin includes 200 stars in a 0.3° sky area around the target star. Crosses denote fits using stars with 2MASS magnitudes, and circles denote those with magnitudes computed from the Gaia colour transformation.

In the text

Fig. 6 Simulated test of our method by propagating the positions of 2000 simulated candidates to µ2 Sco over four epochs. Half are propagated according to the co-moving model, which just adds zero mean noise at each epoch. The other half are propagated according to the proper motion of background stars, plus noise. The colours denote the logarithmic odds ratio that our method computes for each trajectory. All 2000 are correctly identified.

In the text

	Fig. 7 Colour-magnitude diagram of all BEAST stars, except for η Cen at (BP – RP, G) = (4.75, 2.25) mag.
In the text

	Fig. 8 Distribution of the separations of all candidate from all targets identified in BEAST. This assumes the candidates are at the same distance as their respective target stars.
In the text

	Fig. 9 Proper motion model of the background stars for µ2 Sco for the right ascension direction (top) and the declination direction (bottom) in the International Celestial Reference System (ICRS). This is similar to Fig. 4, but now includes the BEAST candidates shown in orange and transformed to ICRS.
In the text

Fig. 10 Demonstration of the proper motion-only method for the known exoplanet µ2 Sco b. The likelihood under the background model is shown in black and the likelihood under the co-moving companion model is shown in blue. The bottom left panel shows these as two-dimensional (Gaussian) distributions. The other two panels show the one-dimensional marginal distributions. The measured relative proper motion is shown in orange: the uncertainties in this is not shown because it is included in the two likelihoods (see Eqs. (5) and (6)).

In the text

Fig. 11 Visualization of the predicted positions of the candidate companion b of the star HIP 71865 (b Cen) under the proper motion and parallax model. A co-moving companion would remain at the position of the first epoch (blue circle) because orbital motion is not included. A field star with the modelled proper motion and parallax of nearby (mostly background) stars would be measured at the two later epochs at the two positions shown by the black triangles. The actual measured change in positions of the candidate are shown as orange triangles. The fact that these are much nearer to the blue distribution means this is likely to be a true companion, something that is properly quantified by our method. The contour lines show 50, 90, and 99% of the enclosed probability, reflecting the propagated uncertainty in the parallaxes, proper motions, and BEAST position measurements. The marginal likelihoods are shown on both axes. This visualization does not show the covariances between the measurements at different epochs, which are nonetheless taken into account by our method (see Eq. (A.4)).

In the text

	Fig. 12 Results of our model applied to the position measurements of HIP 61257 ‘B’ from Kouwenhoven et al. (2005) and Janson et al. (2021b) over a 17-yr baseline. See the caption to Fig. 11 for a description.
In the text

	Fig. 13 Astrometric motion of the two candidates of HIP 81208 with an observation baseline of three years. Viswanath et al. (2023) report both objects as co-moving companions based on their proper motion analysis. Our analysis supports this claim.
In the text

	Fig. 14 Difference between the KS magnitude calculated from Gaia via colour transformations and the and measured 2MASS KS magnitudes for 934 023 objects from Gaia with 2MASS counterparts in the vicinity of µ2 Sco.
In the text

	Fig. B.1 Candidates from BEAST with a logarithmic odds ratio log rtcb > 0. The schematics of the plots are explained in Fig. 11.
In the text

	Fig. B.2 Candidates from BEAST with a logarithmic odds ratio log rtcb > 0. The schematics of the plots are explained in Fig. 11.
In the text

	Fig. B.3 Candidates from BEAST with a logarithmic odds ratio log rtcb > 0. The schematics of the plots are explained in Fig. 11.
In the text