© ESO 2022

1. Introduction

Westerlund 1 (Wd 1) is believed to be the most massive young cluster in the Milky Way, with mass estimates ranging from about 5 × 10⁴ to ≳​10⁵ M_⊙ (Clark et al. 2005; Brandner et al. 2008; Gennaro et al. 2011). Beyond the anecdotal interest of being a top contender for the title of most massive young cluster in the Local Group – in all likelihood beaten by R136 in the Large Magellanic Cloud – the huge mass of Wd 1 turns it into a prime laboratory for massive star evolution. There are so many massive stellar systems in the cluster that rare evolutionary phases are sampled, and different avenues for binary interaction are probed, resulting in a variety of outcomes (see, e.g., Clark et al. 2011).

The population associated with the cluster encompasses all sorts of massive stars, from O-type giants and supergiants to red supergiants (Clark et al. 2020), including a number of extremely luminous B-type hypergiants (Negueruela et al. 2010), an LBV (Ritchie et al. 2009b), several yellow hypergiants (Clark et al. 2010), some transitional B-type supergiants with emission lines (Ritchie et al. 2010; Clark et al. 2014), an extremely luminous B[e]-like object (Clark et al. 2013), and a large population of Wolf-Rayet stars (Crowther et al. 2006). The vast majority of these objects lie within less than 1 arcmin from the nominal cluster centre, although some outliers, such as the Wolf-Rayet stars known as N (=WR 77n) or X (=WR 77sd), are more than 3 arcmin away.

All these varieties of massive objects are generally assumed to correspond to stars with initial masses ≳30 M_⊙ and imply ages ≲5 Ma for stars evolving in isolation (if binary evolution is considered, the age is somewhat older). Nevertheless, an empirical determination of their masses will depend, to a large degree, on the parameters derived for the cluster. Observed stellar magnitudes are transformed into a luminosity and hence a mass, impacting heavily on the age adopted for Wd 1. The two main unknowns are the true distance to Wd 1 and the effect of extinction on the determination of luminosities. The cluster is affected by heavy reddening. Different estimates of the average extinction towards Wd 1 encompass the A_V ≈ 10–12 range, although there are substantial variations across the face of the cluster (Negueruela et al. 2010; Damineli et al. 2016).

With such high extinction, derivation of absolute magnitudes is uncertain, not only because of the possibility of a non-standard extinction law, but also because of the extreme colour terms involved in the colour transformations (Clark et al. 2005), which may lead to systematic errors. Nevertheless, Negueruela et al. (2010) found no strong evidence for large deviations from a standard extinction law from the analysis of VRI photometry. Contrarily, the analysis of near-IR data suggests a much steeper extinction law (Lim et al. 2013; Damineli et al. 2016).

Under the assumption of a standard extinction law, different estimates based on the observed massive star population favour distances around 5 kpc for Wd 1 (Clark et al. 2005; Crowther et al. 2006; Negueruela et al. 2010). Conversely, the analysis of the near-IR colour magnitude diagram with pre-main sequence isochrones by Brandner et al. (2008) resulted in a lower distance of 3.6 kpc. A close value of 3.8 kpc was also obtained from NIR photometry by Lim et al. (2013). Kothes & Dougherty (2007) obtained an indirect estimate of the cluster distance by resorting to the gas clouds in its immediate neighbourhood. By assuming that the clouds were physically connected to Wd 1 and that they followed the Galactic rotation curve, they determined d ≈ 3.9 ± 0.7 kpc. Although these are sensible assumptions, the Galactic spiral pattern delineated by maser radial velocities and trigonometric parallaxes (Reid et al. 2014, 2019) shows a very complex situation along this sightline (l ≈ 340°), where the Scutum–Centaurus and Norma–Outer arms are not easily distinguishable, and tracers of a given arm show a large dispersion in velocity.

Gaia DR2 data for the region of Wd 1 were affected by very strong systematic effects, which rendered any distance estimation very unreliable (Clark et al. 2020). Aghakhanloo et al. (2020) claimed to be able to measure an accurate parallax to the cluster from a Bayesian analysis of parallaxes to stars along the line of sight, obtaining a value of $\pi =0.{35}_{-0.06}^{+0.07}$ mas. Such a result is surprising, given that the median parallax for known members is 0.19 mas (Clark et al. 2020), and implies a distance of only $2.6_{-0.4}^{+0.6}$kpc. This value is substantially shorter than any previous determination and would, according to these authors, imply a total mass not much higher than 2 × 10⁴ M_⊙. Davies & Beasor (2019) carried out a more careful analysis, by selecting only stars whose proper motions were similar to those of known cluster members, coming to a distance of ${3.9}_{-0.64}^{+1.0}$kpc, more in line with previous authors.

More recently, Aghakhanloo et al. (2021) have insisted on a distance of ${2.8}_{-0.6}^{+0.7}$kpc derived from an analysis of the distribution of parallaxes in Gaia EDR3 for all the stars in the area, without any reference to membership. Meanwhile, Beasor et al. (2021) have used EDR3 data on OB stars with a coherent proper motion distribution to derive a substantially longer distance of 4.1${}_{-0.36}^{+0.66}$ kpc. Moreover, based on an estimation of the luminosity of the cool supergiants in Wd 1, Beasor et al. (2021) have argued that all these objects have ages around 10 Ma, implying masses of only ∼17 M_⊙, while the presence of some more massive stars and the younger age derived from pre-main-sequence isochrone fits can be explained by non-coeval or extended star formation.

The aim of this paper is to provide a well-founded estimation of the extent of Westerlund 1 and its possible connection to other populations in its surroundings, with the purpose of evaluating the likelihood of a complex (multi-age or multi-cluster) population. For this, we use two complementary tools: Gaia EDR3 astrometric data and spectroscopy of a large sample of luminous stars in a field of radius one degree surrounding the cluster. After introducing the data used in Sect. 2, we carry out membership analysis on the Gaia EDR3 data in Sect. 3 and then proceed to calculate an accurate distance to the cluster in Sect. 4. We then describe the population surrounding the cluster in Sect. 5 and explore the consequences of our findings in the discussion. We close the paper with our conclusions.

2. Data collection

2.1. Spectroscopy data

Observations of the field surrounding Westerlund 1 were obtained with the fibre-fed dual-beam AAΩ spectrograph mounted on the 3.9 m Anglo Australian Telescope (Siding Springs, Australia) on the nights of 2011, July 20 and 21. The Two Degree Field (2dF) multi-object system was utilised to position fibres. The instrument allows simultaneous observations in two different arms by using a dichroic beam-splitter with crossover at 5700 Å. Each arm of the AAOmega system was equipped with a 2k × 4k E2V CCD detector and an AAO2 CCD controller. In the blue arm, we used grating 580V, which gives a resolving power R = 1300 over ∼2100 Å (with central wavelength at 4500 Å). Given the extremely high extinction to most of our targets, useful blue spectra were obtained for only a handful of objects. In the red arm, the 1700D grating was used, providing R ≈ 11 000 in a 500 Å wide spectral window containing the near-infrared Ca II triplet (CaT). This grating must be centred on 8700 Å, but the actual central position of each spectrum is determined by the position of the target on the sky, with a maximum shift of 20 Å.

The targets observed were selected by means of a combination of photometric catalogues, following the procedures outlined in Negueruela & Schurch (2007) and Negueruela et al. (2012), by means of the Q_IR index, defined as Q_IR = (J − H)−1.8 ⋅ (H − K_S). A sample of “early-type” stars was created by selecting objects from the 2MASS (Skrutskie et al. 2006) catalogue with good flags (i.e. “A” or “E”) in all filters whose Q_IR index falls in the range between −0.5 and +0.1, which is occupied by emission-line stars and early (OBA) stars. Some F stars may also be included, as room is made for photometric errors, but the main contaminant population are late-M AGB stars, whose SEDs are similar to those of reddened OB stars (e.g., Maíz Apellániz et al. 2020). This initial selection was then cross-matched with the DENIS i band (or USNO-B1.0 I band, if DENIS data did not exist) to guarantee that the target would be observable (i < 15). The list was completed with the very few catalogued early-type stars in the field and known members of Wd 1 (which, in most cases, have no reliable 2MASS magnitudes due to crowding and saturation).

Candidate cool luminous stars were selected with the same criteria that had worked successfully in other highly extincted areas (Negueruela et al. 2012). We picked very bright stars in the infrared (K_S < 7), with E(J − K_S) > 1.3 (i.e. redder than the intrinsic colour of any red supergiant) and Q_IR < 0.4 (a cut that leaves out most red giants). Again, a cross-match with DENIS and USNO-B1.0 was conducted to make sure that the star would be observable. Known cool supergiant members of Wd 1 were added by hand.

On the night of July 20, we observed a bright field, aimed at stars with I < 11.5, with three exposures of 300 s. The following night, we took four exposures of 1200 s, targeting mostly fainter targets. We used the standard reduction pipeline 2dfdr as provided by the AAT at the time. Wavelength calibration was attained by observing arc lamps before each target exposure. Details of the reduction can be found in González-Fernández et al. (2015).

2.2. Gaia astrometric and photometric data

Gaia EDR3 has provided unique quality photometric and astrometric data for more than 1.5 billion celestial objects (Gaia Collaboration 2021). The precision of the astrometric data (parallax and proper motions) is over 20% better, on average, than in the previous release, Gaia DR2. We have made use of the version of the EDR3 catalogue available from Vizier (I/350/gaiaedr3) to obtain astrometry and photometry for the stellar sample analysed in this work. We downloaded data for all stars within a circle of radius 12 arcmin around the nominal centre of Wd 1 given by the SIMBAD database (ICRS; Epoch 2000; RA = 251.76667°; Dec = −45.85136°; Brandner et al. 2008). The large radius chosen, more than six times larger than the 1.7 arcmin catalogued by Sampedro et al. (2017), is recommended by our objective of analysing the cluster and its environment in a highly crowded region of the sky whose line of sight could be crossing several spiral arms, and, therefore, a variety of stellar populations with different kinematic properties (e.g., Reid et al. 2014, 2019).

The selection has only been filtered by the condition that the re-normalized unit weight error (RUWE) is < 1.4. In addition, data are only taken from a proper motion box defined by μ_α [−10; 7] and μ_δ [−11; 5] in mas/a units¹. RUWE is a quality-fit parameter of Gaia astrometric solutions (Lindegren et al. 2021b). It is directly linked to the χ² statistic. Several authors suggest a RUWE threshold at 1.4 for the selection of large data samples, as this upper limit guarantees a well-behaved astrometric solution. This would be equivalent to choosing objects with χ²  ≤  2 (e.g., Stassun & Torres 2021, among others). Here, we adopt this criterion. The proper motion box is defined to remove the most obvious outliers in the proper motion space to prevent them from biasing the selection algorithm. The number of stars in the sample selected is 19 432. Several authors (e.g., Lindegren et al. 2021a; Maíz Apellániz 2022; Riello et al. 2021; Huang et al. 2021; Yang et al. 2021; Niu et al. 2021; Cantat-Gaudin & Brandt 2021), have discussed the possible systematic errors that could still be present in astrometry and photometry published in Gaia EDR3. We return to these corrections throughout the work, when we consider that their application could affect the physical variables estimated, and/or modify some selection criteria or some of the conclusions reached.

3. Membership analysis

3.1. Gaia EDR3: Membership exploratory analysis

As mentioned in the Introduction, the stellar population associated with Wd 1 contains highly reddened, intrinsically very luminous, massive stars. With this prior knowledge, we explore how Gaia EDR3 photometry informs us about the stellar population in the Wd 1 field. Figure 1 shows the G vs. (BP − RP) diagram for the full sample. We have marked in red the 134 stars for which spectral types have been provided by Clark et al. (2020) and Gaia EDR3 photometry exists in our sample. From all previous analyses (e.g., Clark et al. 2020, and references therein), all these objects are very likely cluster members. More than 90% of these stars have a (BP − RP) colour redder than 4 mag, with the remainder having very slightly bluer values. For an initial exploration of cluster parameters, we chose stars with (BP − RP) > 4, according to their Gaia EDR3 photometry. We call these objects the Red Branch (RB) stars. It is noteworthy that any correction to the Gaia EDR3 photometry, as suggested by some authors (Yang et al. 2021; Niu et al. 2021, and references therein), is too small to alter in any significant way the sample of RB stars selected.

Fig. 1. G vs. (BP − RP) colour-magnitude diagram for the total sample with RUWE ≤ 1.4 (grey small dots). Stars with spectroscopic information in Clark et al. (2020) and Gaia EDR3 photometry are shown in red. This sample is clearly seen to consist mostly of stars redder than (BP − RP) > 4.

As a first approach, we explored how this data set behaves in different subspaces of the phase space formed by the astrometric variables (α, δ, μ_α, μ_δ, ϖ). Initially, we did not use any parallax correction, as we were simply trying to determine the shape and structure of the distributions. For a more quantitative analysis, we later apply the necessary systematic corrections to parallax at the cost of substantially reducing the sample. In Fig. 2, we show the spatial distribution of the 975 RB stars in the field of Wd 1. The coordinates are angular distances to the centre of the cluster in RA (X) and Dec (Y), both measured in arcminutes. The stars distribute forming an elliptical and very dense central core whose major axis is ≈7′ long. Surrounding the core, a halo can be seen, mainly in the north-west quadrant, with a density of objects much lower (a factor 50) than the core density maximum.

Fig. 2. Red Branch stars, defined as those with (BP − RP) > 4, plotted on the plane of the sky (dots). Coordinates are angular distance to the cluster centre along RA (X axis) and Dec (Y axis). Red squares identify stars with spectra reported in the catalogue of Clark et al. (2020). The vast majority of red squares are within a radius of $\approx 3\stackrel{\prime }{.}5$.

The Vector Point Diagram (VPD) in Fig. 3 displays a similar pattern. Here, we also include the complete EDR3 sample to show how the distribution of the RB stars differs from that of the bulk of data. The distribution of the whole sample is more elongated and centred on a nearby, but significantly different, position to the centre of the RB stars.

Fig. 3. Isodensity contours (black line) of the total sample (grey small dots). Red dots are RB stars, which are in an eccentric position with respect to the whole sample distribution. They appear centred on one focus of the internal dispersion ellipse.

The histogram of the parallaxes for the total sample and for the Red Branch (Fig. 4) clearly indicates that the stars with redder (BP − RP) colours present a mode well separated from that of the total sample (by ≈0.20 mas, attending to the position of the Kernel density maxima for both populations). Moreover, the RB FWHM is twice as narrow as that of stars of bluer colour.

Fig. 4. Histogram showing the parallax distribution for the total sample (grey) and the RB stars (red). Frequency is normalised to the maximum value for better visualising the difference in modes. Attending to the histogram, this difference is 0.25 mas (bin size). Kernel density estimators lead to a slightly shorter separation of ≈0.20 mas.

Our exploratory analysis thus suggests that the RB stellar population shows distributions in the different planes that are compatible with, and representative of, a young and rich stellar cluster, from whose population we are just detecting the tip of the iceberg, i.e. the most luminous stars. We therefore conclude that the (BP − RP) colour is an excellent discriminator between cluster members and field stars for this stellar system.

3.2. Quantitative membership analysis: SALSON

As a more sophisticated approach, we then proceeded to obtain individual membership probabilities by applying a non-parametric method based on the direct estimation of the probability density functions (pdf) by Kernel functions (SALSON; Cabrera-Cano & Alfaro 1990). This method is based on the assumption that cluster members are more densely concentrated than field stars in the space of variables. Eq. (1) defines the pdf of the class k for the location of the object i, in a variable space S, using a Kernel, K, where S_i represents the vectorial coordinates of the object i in the space S, as

$$\begin{array}{c}\hfill D_{\mathbf{Si}}(k)={\displaystyle \sum _{i\ne j}^{\forall j\in k\mathrm{class}}{K}_{\mathrm{S}}({\mathbf{S}}_i,{\mathbf{S}}_j,h_{\mathrm{S}}).}\end{array}$$(1)

We assume that there are only two classes (k): cluster (c) and field (f). We work with a Gaussian Kernel that only depends on the smoothing parameter h_S. The space S is formed by the four variables (X,  Y,  μ_α,  μ_δ), but considering that the subspaces P= (X,  Y) and V= (μ_α,  μ_δ) are statistically independent, in such a way that D_PVi(k) can be estimated as D_PVi(k)  =  D_Pi(k)×D_Vi(k). The separation of the space of variables allows us to determine three membership probabilities, one each corresponding to the subspaces P and V, both in 2D, and the third associated with the 4D space. These probabilities are defined by the expression

$$\begin{array}{c}\hfill {\displaystyle P{r}_{\mathbf{Si}}(k)=\frac{Pr(k)\times D_{\mathbf{Si}}(k)}{{\sum }_{\forall k}Pr(k)\times D_{\mathbf{Si}}(k)},}\end{array}$$(2)

where Pr(k) are the a priori probabilities that are estimated, in a frequentist inference approach, as the relative frequency of elements belonging to each class (k). Obtaining the pdfs, as well as the a posteriori probabilities, follows this procedure:

1. Removing the most evident outliers in both subspaces.

2. Estimating the Kernel smoothing parameters (h_P and h_V) as those that maximise the likelihood in both subspaces (P and V), respectively (Silverman 1986). This requires performing the previous step; otherwise, the smoothing parameter estimates could be biased, and so could be the final probability assignment (Cabrera-Cano & Alfaro 1985).

3. Estimating the pdf in subspace V for each object in the whole sample, and sort them from highest to lowest densities.

4. Making the first classification between both classes, by assigning as cluster members the higher 0.4 percentile of the pdf in V (D_Vi).

5. Calculating D_Pi(c), D_Pi(f), D_Vi(c), and D_Vi(f), by using Eq. (1). With these values, the probabilities Pr_Pi(c), Pr_Vi(c), and Pr_PVi(c) are then estimated via Eq. (2).

6. Classifying as cluster members those objects with Pr_Vi(c)≥0.5 and Pr_PVi(c)≥ 0.5. This criterion takes into account both subspaces, giving more weight to the kinematic variables.

7. Comparing this classification with the previous one, and, if they do not coincide, returning to step 5. If they match, the run is stopped and an output file is written, containing the D_Si(k) for the two subspaces and both classes, as well as the three estimated probabilities Pr_Pi, Pr_Vi, and Pr_PVi.

From the 19 432 objects in our initial sample, SALSON identifies 8010 as cluster members. In Figs. 5 and 6 we show the pdfs of the cluster class in subspaces P and V. Stars identified as members distribute in positional space (P; see Fig. 5) with two different shapes, a central concentration of stars, within a raw external radius of around $3\stackrel{\prime }{.}5$, and an extended halo, where the pdfs for the cluster (D_P(c); blue) and field (D_P(f); green) classes are much smaller and very close to each other, so that the probabilities given by Eq. (2) can be very uncertain for these objects.

Fig. 5. DP pdfs on the X–Y plane for the probable clusters members selected by SALSON. DP(c) in blue and DP(f) in green. DP(c) and DP(f) show similar values for objects well outside the central core, which leads to larger uncertainties in the estimated cluster membership probabilities.

Fig. 6. DV pdfs in the proper motion space for the probable clusters members selected by SALSON. DV(c) is in blue and DV(f), in green. DV(c) shows a non-single peak distribution with a rounded bump at lower proper motions. This is likely indicative of a field contamination not removed by SALSON. Possible causes are discussed in the text.

In the proper motion subspace (Fig. 6), there is a clear separation between both pdfs, defined by the plane D_V(c) = 0.02. The pdf in V shows a rounded peak, but there is also a bump of stars adjacent to the main maximum, which suggests that we could have a non-negligible contamination from other stellar populations with different, but close, angular velocities.

In any event, the number of objects selected is much higher than the number of objects with red colours that we identified in the previous section. Moreover, the stars selected as cluster members are distributed all over the photometric CMD. Therefore, we proceed to prune the dataset of those objects, which, despite their SALSON cluster member classification, show lower pdf values in both subspaces, as well as very similar values for both classes. To this aim, we plot the radial spatial distribution of the cluster by drawing the D_P(c) (blue) and D_P(f) (green) densities versus the radial angular distance to cluster centre for the objects selected as members (Fig. 7).

Fig. 7. Wd 1 radial density profile. The DP pdfs versus angular distance (arcmin) to cluster centre are shown for the members selected by SALSON. DP(c) is shown in blue and DP(f), in green. We identify as most likely cluster members those stars within the upper-left rectangle (red lines): DP(c) > 0.002 and r < 3.5 arcmin.

This plot suggests that most conspicuous members are within a radius of $\approx 3\stackrel{\prime }{.}5$ and show D_P(c) larger than 0.002. Thus, stars enclosed in the upper-left rectangle of Fig. 7 are our best candidates according to their spatial distribution. By imposing this criterion plus that of the kinematic subspace (D_V(c) > 0.02), we extract 1154 high probability members supported by, and estimated from, astrometric Gaia EDR3 data. Still, these objects display a variety of photometric colours and magnitudes that are incompatible with a single population. If we now apply the colour criterion that was identified in the previous section, ((BP − RP) > 4), we are left with a set of 401 bona fide cluster members that we later use as probes for obtaining accurate values of the physical parameters of Wd 1. We have to note that the selected cluster radius is larger than any other previous estimation, but we prefer to be conservative in this initial selection and come back later to a more refined calculation of the cluster radius.

3.3. Wd 1 and its surroundings: Is there anything else?

In the previous section, we have selected a list of bona fide cluster members by applying astrometric and photometric criteria. However, the simple application of a non-parametric method to the astrometric variables separates a much larger population, prompting two questions: Is there a second kinematic population hidden in this selection? Is our application of a colour cut fully justified?

To answer these questions, we proceed to analyse the spatial distribution of the (BP − RP) colour for the 8010 probable cluster members selected by SALSON. The goal is to find whether there are narrow ranges of colour that are spatially concentrated in the RA-Dec subspace. For this, we use the Spectrum of Kinematic Groupings algorithm (SKG; Alfaro & González 2016; Alfaro & Román-Zúñiga 2018) based on the Minimum Spanning Tree (MST) which was first proposed for the analysis of mass segregation (Allison et al. 2009; Maschberger & Clarke 2011), but has been easily extended to other physical variables that may also be spatially segregated (González & Alfaro 2017; Costado et al. 2017; Costado & Alfaro 2018; Carballo-Bello et al. 2021).

The basic foundations of the method rely on sorting the variable to be analysed (ascending or descending order does not matter), selecting the number of points per bin, which is tentatively taken as the integer closer to the square root of the number of total points, and choosing another integer that stands for the number of overlapping points between two consecutive intervals. The algorithm estimates the MST of each bin in the position space and the median of the distances between adjacent points of the MST. This value is compared with that obtained for a sample of the same size extracted from the original sample (full range of the variable). If the quotient between this last median and that of the corresponding bin is significantly greater than 1, we say that this interval of the variable is spatially segregated, or, in other words, that condition

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_i-(2\times {\sigma }_{{\mathrm{\Lambda }}_i})>1\end{array}$$(3)

is fulfilled. For a more detailed explanation on how Λ_i and σ_Λi are estimated, we refer the reader to Alfaro & González (2016). A more conservative criterion can be chosen by simply increasing the multiplicative factor of σ_Λi in Eq. (3). The latest implementation of the algorithm (Alfaro & Román-Zúñiga 2018) also calculates the value of Q (Cartwright & Whitworth 2004) as a quantitative description of the internal spatial structure of the objects in each bin. Values above 0.8 indicate central distributions. Clumpy patterns present values below 0.8 and a homogeneous distribution would show values around 0.8 (Cartwright & Whitworth 2004; Sánchez & Alfaro 2009). The precision in the estimation of Q is obtained by bootstrapping (Alfaro & Román-Zúñiga 2018).

By applying this technique to the 8010 probable members outlined by SALSON, we obtain the SKG shown in Fig. 8. Along the x axis we plot the central (BP − RP) colour of each bin, while on the vertical axis we display Λ_i − (2  ×  σ_Λi). The plot clearly highlights the increase of Λ_i with growing colour (starting around (BP − RP)≈ 4), which indicates a sharp spatial segregation for the redder stars, peaking at (BP − RP)≈5². Figure 9 shows the plot of Q vs. (BP − RP) colour. Again, the reddest stars show a significant increasing trend of Q with growing colour that indicates that they are more centrally concentrated the redder they are, in agreement with the results observed in the SKG (Fig. 8). On the blue side of the colour interval ((BP − RP) < 1.5), a spatial segregation with a central concentration (Q > 0.8) is also observed, although with Λ and Q values much lower than those shown by the red tip of the colour interval. Based on this result, we then select two data sets of spatially segregated stars: (a) the red group (RG) formed by stars with (BP − RP) > 4, and (b) the blue group (BG) formed by objects with (BP − RP) < 1.5. Both groups present well differentiated distributions of the astrometric variables in Gaia EDR3.

Fig. 8. Spectrum of the spatial segregation by (BP − RP) colour. Error bars are σ values. Λ–(2 × σ) values higher than 1 are indicative of colour bins spatially grouped. A strong spatial segregation is observed for stars redder than (BP − RP) > 4. A blue grouping is also detected with (BP − RP) < 1.5.

Fig. 9. Q spatial concentration parameter as a function of colour. Error bars have been obtained by bootstrapping. Q ≈ 0.8 separates the clumpy and central spatial distributions. The higher Q is, the more spatially concentrated the distributions. This plot reinforces the conclusions derived from Fig. 8.

To visualise this, we show the distribution of both groups in equatorial coordinates (Fig. 10), proper motions (Fig. 11) and parallaxes (Fig. 12). The isodensity contours for both groups have also been drawn for descriptive purposes in the position and proper motion diagrams. The RG is mainly associated with the known population of Wd 1, as already discussed in previous sections. The BG appears to be distributed over the area located to the south-west of the cluster. The density local maximum is around (l, b) = ($339\stackrel{°}{.}41$, $-0\stackrel{°}{.}38$), but the population is much less concentrated. Similarly, the RG is densely concentrated in the Vector Point Diagram, while the BG seems to be associated with a second bump of lower height, well separated from the central value of cluster proper motion. However, the most conclusive evidence that we are observing two different stellar systems is provided by the histogram of the parallaxes (Fig. 12) which shows two clearly different distributions. We performed a two-sample Kolmogorov-Smirnov test for the parallaxes of the groups RG (771) and BG (222). The statistic KS is 8.198, which leads to a probability p < 0.00001, whose value is much lower than 0.05, therefore rejecting the hypothesis that the two samples come from the same distribution. It is clear then that the BG represents a grouping located closer to the Sun, with a central parallax (mode) implying a raw distance of around 2 kpc.

Fig. 10. Stars (dots) and isodensity contours (lines) for the two groups spatially segregated by (BP − RP). BG (in blue) contains the objects in segregated bins with colour bluer than 1.5. RG (in red) is the same for colours redder than 4. The RG is mainly associated with the cluster centre, while the BG shows a stellar population with its maximum density well separated from the core of Wd 1. The axes show angular distances in equatorial coordinates.

Fig. 11. Stars and isodensity contours in proper motion space for the two groups spatially segregated by (BP − RP). Symbols and colours as in Fig. 10. The spatial segregation by colour is mimicked in the VPD. Although both distributions show a non-null overlap, their density maxima are well separated in proper motion.

Fig. 12. Histogram of corrected (see Sect. 4.1 for parallax bias corrections) Gaia EDR3 parallaxes for the RG (red) and BG (blue) groups, showing two different statistical distributions. RG shows smaller parallaxes than BG. This result, together with those shown in Figs. 8–12, strongly suggests that both populations are well separated by colour, position, kinematics, and distance. Gaussian Kernel density estimators are also drawn.

The distribution of the RG and BG groups on the plane of the sky appears to be primarily drawn by the interstellar extinction pattern in the field. Both groups extend mainly over the western region of the area analysed (Fig. 10), while to the east of the cluster the surface density of objects decreases significantly. Going into more detail, while the RG is distributed over most of the western semicircle, the BG is concentrated towards the south-west of the cluster.

The nature of the BG group is unclear. We discuss it in Appendix A. Independently of this, the analysis presented in this section leaves two conclusions: (1) the simple application of a non-parametric method to the astrometric variables (positions and proper motions) does not separate effectively the cluster population from its surroundings, and (2) the population with redder colours appears as clearly distinct in both physical space and astrometric parameter space, fully justifying our choice of colour as a discriminant.

4. Distance estimation and kinematic parameters

Once we have selected a safe set of cluster members, we can use their Gaia EDR3 parallaxes to estimate the distance to the cluster. Unfortunately, the quality of astrometric data for bona fide cluster members is low, because of crowding and faintness, and we cannot calculate bias corrections for a majority.

4.1. Parallax bias correction

Parallaxes catalogued in Gaia EDR3 present systematic biases (or zero points) that depend on several factors, such as the photometry of the star, its location on the sky, ecliptic latitude and number of orbits, as well as the set of free parameters used in fitting the astrometric solution (Lindegren et al. 2021a; Maíz Apellániz 2022). Those two references provide independent estimates of the systematic biases that are quite similar for faint stars (G > 13) but are somewhat different for brighter objects. The biases are of the order of a few tens of mas and can be corrected by employing the Python algorithm provided by the first set of authors³ and by the IDL algorithm provided by the second author as an appendix to his paper (which can be used for both zero points). From now on, we refer to the first systematic bias as the Lindegren zero point (Lztp) and to the second one as the Maíz Apellániz zero point (MAztp).

We applied both algorithms to the whole sample of 19 432 objects. A total of 12 070 stars have zero-point estimates, which represents 62% of the initial sample. This drastic reduction in the number of objects is mainly due to the position of Wd 1 in a very crowded region, showing a wide range of magnitudes and colours, but also because we have only recovered interpolated solutions (see Lindegren et al. 2021a, and the corresponding tutorial, for a more detailed explanation of these issues). In other words, we do not work with extrapolated zero points for the sake of accuracy, at the cost of reducing the final number of cluster members. The effect on the bona fide members is much worse than on the whole sample, as we are dealing with some of the most crowded and reddest stars. From 401 members selected in Sect. 3.2, only 142 stars have interpolated solutions. Of those, 138 have five-parameter astrometric solutions and four have six-parameter astrometric solutions.

4.2. Cluster parallax

Once we have corrected the individual parallaxes from the zero point biases estimated in the previous section, we proceed to determine the cluster parallax and its uncertainty. Until now, we have not filtered the data by astrometric errors, either in proper motions or parallaxes; we only applied a RUWE threshold. Now, we do introduce the errors catalogued in Gaia EDR3 to estimate the total uncertainty of the parallax and use them as a weight parameter (alone or together with other parameters connected to the cluster membership) in determining the mean parallax of the cluster and its error. To this goal, we follow the procedure proposed by Maíz Apellániz (2022, and references therein).

Parallax errors catalogued in Gaia EDR3 data release do not include the actual uncertainty in the measurement of that variable. However, the true value can be approximated (see discussion in Maíz Apellániz et al. 2021) by the expression:

$$\begin{array}{c}\hfill {\sigma }_{{\varpi }_{\mathrm{e}}}=\sqrt{k^2{\sigma }_{\varpi ,i}^2+{\sigma }_{\varpi ,\mathrm{s}}^2}.\end{array}$$(4)

The values involved in this prescription can be different for different releases and the subscripts i, s, and e indicate formal, systematic, and total errors, respectively. The multiplicative constant k depends on several factors, and there are different empirical approaches for its estimation. We consider two possibilities: (a) the estimate from Maíz Apellániz et al. (2021) when applying Lzpt and (b) the estimate from Maíz Apellániz (2022) when applying MAzpt.

Both approaches (and others, such as Vasiliev & Baumgardt 2021, as well) yield similar values for k, which do not differ significantly from 1.3 for most of the stars in our sample with high weights. For σ_ϖ, s, we use the value of 10.3 μas from Maíz Apellániz et al. (2021).

In Fig. 13 we show the histograms of the Lzpt correction for the entire sample (red) and for the bona fide members (blue). Both histograms have been normalized to their respective maximum to better distinguish distribution modes. While the majority of stars have corrections of ≈−0.03 mas, close to the average value for Gaia EDR3 given by Lindegren et al. (2021a), cluster members typically have corrections between −0.05 and −0.085, with the mode at −0.08 mas. The photometric values (G and (BP − RP)) are the main variables that influence the high absolute value of this correction for cluster stars, as can be seen in Fig. 14.

Fig. 13. Histogram of parallax zero-point corrections according to Lindegren et al. (2021a). Bona fide members are in blue, and the global sample is in red. The histograms have been normalized to their maxima for better visualizing the different distribution modes.

Fig. 14. Cluster CMD showing the value of the zero point applied to each one of the bona fide members of Wd 1 for which such correction has been estimated. Left panel: Lindegren. Right panel: Maíz Apellániz. The scale is in μas and larger symbols have larger full weights.

The value of the cluster unbiased parallax can be estimated by the expression

$$\begin{array}{c}\hfill {\varpi }_c={\displaystyle \sum _{i=1}^n{w}_i{\varpi }_i,}\end{array}$$(5)

where the weights w_i can be defined in two different ways by assigning different weights: (a) by taking into account only the true uncertainty of the parallaxes, as in Maíz Apellániz et al. (2021), or (b) by also including information related to the membership analysis that we have carried out in Sect. 3.2. In this way, two different aspects are involved in the weighted average, the parallax error, which is intrinsic to the star and its measurement, and the congruence of its positional and kinematic values within the cluster, which somehow measures its connection to the stellar system itself. Both conditions constrain and refine the final cluster parallax. The expressions that we use are:

$$\begin{array}{c}\hfill {\displaystyle w_i=\frac{1/{\sigma }_{\mathrm{e},i}^2}{{\sum }_{i=1}^{n}1/{\sigma }_{\mathrm{e},i}^2}}\end{array}$$(6)

and

$$\begin{array}{c}\hfill {\displaystyle w_{t,i}=\frac{(1/{\sigma }_{\mathrm{e},i}^2)\times D_{P,i}(c)\times D_{V,i}(c)}{{\sum }_{i=1}^n(1/{\sigma }_{\mathrm{e},i}^2)\times D_{P,i}(c)\times D_{V,i}(c)}}\end{array}$$(7)

for the case when only parallax error is taken into account (a; Eq. (6)), and for the case when we also include the pdfs for the cluster class in both subspaces (positional and kinematic) in the computation (b; Eq. (7)). We refer to Eq. (6) as the parallax weights and to Eq. (7) as the full weights.

The error of the weighted mean parallax is obtained with the formula

$$\begin{array}{c}\hfill {\sigma }_{{\varpi }_c}^2={\displaystyle \sum _{i=1}^n{w}_i^2+2\sum _{i=1}^{n-1}\sum _{j=i+1}^n{w}_i{w}_j{V}_{\varpi }({\theta }_{\mathit{ij}}).}\end{array}$$(8)

Here, the second term on the right side of the equation represents the angular covariance, where V_ϖ(θ_ij) is given by Eq. (8) in Maíz Apellániz et al. (2021). This analytic function depends on the covariance angular value, V_ϖ(0), at the limit θ → 0, which is taken as 100 μas in this work.

We present the results for the parallax to Westerlund 1 in Table 1. We show eight different combinations, by varying the type of solution used (all or only interpolated), the zero point applied (MAzpt or Lzpt) the number of astrometric parameters for the stars (five and six or just five), and the weights used (parallax or full). In each case, the number of stars in each sample, the cluster parallax, and the cluster distance are given. For the distances we use the prior of Maíz Apellániz (2001), Maíz Apellániz et al. (2005), which is specific for OB stars, with the parameters of Maíz Apellániz et al. (2008).

The results for the eight cases in Table 1 are very similar, with the two members of any given pair within one sigma or less of each other. This indicates that our distances are robust and independent of the selection we choose from the sample of probable members, or which zero point is used. From now on, we adopt the case C distance of $4.{23}_{-0.21}^{+0.23}$ kpc.

4.3. Cluster kinematics

Once the sample of bona fide members of Wd 1 has been selected and its distance estimated with the procedure detailed in the previous section, we now proceed to determine the main kinematic variables of the cluster. For this, we now use the whole sample of 401 cluster members.

As with the parallaxes, the proper motions in Gaia EDR3 show systematic errors that, although negligible in most cases, need to be corrected if we want to determine the motion of the cluster in the Galaxy with adequate precision to analyse its current and future dynamic state. The zero points for proper motion corrections have been determined following the procedure described in Cantat-Gaudin & Brandt (2021). The estimates depend mainly on the celestial positions and on the G magnitude of the stars. In the case of the 401 stars in our sample, this correction is limited to an average zero point of −0.02 mas/a in μ_δ. Given the cluster coordinates, and the stellar brightness range, the zero point in μ_α is negligible. The mean values and variances of the μ_α and μ_δ components have been estimated following the same procedure as for the determination of parallax. A value of μ = [ − 2.231 ±  0.008, −3.697 ±  0.008], in Equatorial coordinates, has been obtained in mas/a units by using the full weights from the previous subsection.

5. Spectroscopic data

5.1. Description of the sample

In total, our AAOmega observations include 14 early-type members of Westerlund 1, the LBV Wd1-W243, two yellow hypergiants, Wd1-W4 and Wd1-W265, and the four red supergiants. Spectra of the most luminous stars are displayed in Fig. 15 (left panel). The stars that were observed are indicated in the penultimate column of Table B.1. Five of the early-type stars are located outside the well studied core region (of radius around 2′). Four of them were subsequently observed with FLAMES and are thus included in the catalogue of spectroscopic members in Clark et al. (2020), as W1049, W1053, W1067, and W1069. The fifth is more than 5′ away from the centre, but its proper motions are within 1-σ of the cluster average and its location in the CMD confirms it as the most distant known member. Following the naming convention in Clark et al. (2020), we identify this object (Gaia EDR3 5 940 104 594 038 054 144) as Wd1-1070. It is further discussed in Appendix B.

Fig. 15. AAOmega spectra of luminous supergiants. All the spectra have been normalised and the continuum has been shifted for clarity. Left panel: members of Wd 1. The (green) dashes above the spectrum of Wd1-W57a indicate the positions of Paschen lines. The three thin (orange) vertical lines show the position of the Ca II triplet. The short (red) dashes below the spectrum of Wd1-W265 identify some of the strongest N I lines. The dotted line marks the position of the 8621 Å DIB, whose strength correlates well with extinction up to a saturation value. Right panel: other stars in the field, not connected to the cluster. Symbols as in the left panel.

Outside this central concentration, which corresponds to the cluster, there are 29 other “blue” targets, none of them closer than 20′ from the cluster centre. Of these, four objects are catalogued early-type stars with usable blue spectra. All other sources, with only I-band spectra, are new identifications, except for star #119 = [GKF2010] MN48, which has subsequently been signalled as an LBV (Kniazev et al. 2016), and is further discussed in Appendix C. The list of targets observed is shown in Table C.1. Spectral types have been estimated following the criteria discussed in Negueruela et al. (2010). See Appendix C for further details.

In addition, we have spectra for 321 “red” targets. Of them, only two are within 5′ of the cluster centre. As is typical of the criteria used, a large fraction of the stars observed turn out to be late-M giants, in all likelihood AGB stars, whose colours are indistinguishable from those of red supergiants. The list of targets observed and their characteristics can be found in Table C.2. Spectral types have been estimated following the criteria discussed in Negueruela et al. (2012). We find around 100 stars whose spectral features allow a supergiant classification (luminosity class Ib-II or above). Of these, close to 40 can be unambiguously identified as luminous red supergiants (RSGs), i.e. massive stars.

Despite a moderate success rate among RSG candidates (fully within the expectations of the method), between the blue and red targets, we identify at least 70 previously unknown massive stars within 1° of Wd 1. A representative sample of good quality spectra is shown in Fig. 15 (right panel), where the whole range of spectral types is represented.

5.2. Radial velocities

For the majority of the stars in the “red” sample (spectral types G and later), it was possible to calculate radial velocities by means of the method described in detail in González-Fernández et al. (2015). In short, we cross-correlated each observed spectrum against the most similar MARCS synthetic spectrum in the grid available from the POLLUX database (Palacios et al. 2010). In González-Fernández et al. (2015), this method resulted in a typical uncertainty (velocity dispersion for stars observed more than once during that run) around 1.0 km s⁻¹. As that sample was observed with the same instrumental configuration and on the same dates as our data, we can assume with confidence the same uncertainty for our measurements. Finally, we calculated the heliocentric and Local Standard of Rest corrections through the spectroscopic-analysis software ispec (Blanco-Cuaresma et al. 2014; Blanco-Cuaresma 2019).

For 115 objects, we can compare the RVs measured from our spectra to those measured by Gaia DR2. The average difference (our heliocentric velocities minus those from DR2) is −0.7 km s⁻¹ (below our uncertainty of 1.0 km s⁻¹), with a standard deviation of 3.0 km s⁻¹. There are only five stars whose differences are > 6 km s⁻¹ in modulus. These values suggest that our RVs are in exactly the same reference system as those from Gaia DR2, and we can safely use our much larger sample.

Unfortunately, the vast majority of our targets are faint and have very poor astrometric solutions, and any attempt to derive individual distances results in very large uncertainties. For completeness, in Table C.2, we give the distances derived by Bailer-Jones et al. (2021). If we select only objects that we have classified as RSGs and impose a quality criterion ϖ/σ_ϖ > 6, there might be a trend to more negative radial velocities with increasing distance, as expected from the radial velocity curve in this direction. However, no object with a Bailer-Jones et al. (2021) distance > 3 kpc passes the quality criterion, and thus our objects cannot be used to explore the radial velocity distribution up to the distance of Westerlund 1 with the current EDR3 dataset.

6. Discussion

We have investigated the population of massive stars in the field surrounding Westerlund 1 and the kinematic properties of the cluster and its surroundings. The first conclusion that may be drawn from the previous analysis is that this is a very complex field. A direct application of non-parametric Bayesian inference (SALSON) to separate the cluster from the field results in a sample with an important field contamination. The stars identified as cluster members, based only on sky positions and proper motions, represent 40% of the initial sample, but are spread over the parameter space both in parallax and in the CMD. This field contamination is mainly caused by four reasons: (1) the cluster proper motions are quite similar to those of many foreground stars along the line of sight (see Fig. 3), (2) the actual cluster members represent a very small fraction of the total sample, (3) almost all actual members have large errors in their astrometric parameters (see Table B.1), precluding their identification as a distinct population, and (4) the large sampling radius, much larger than the cluster size, results in a higher pollution by field stars (Sánchez et al. 2010, 2020), which we finally remove by restricting the sample to the central $3\stackrel{\prime }{.}5$.

In spite of this, the spatial distribution of members selected by SALSON is strongly concentrated towards the position of Westerlund 1. There are two main reasons for this: (1) actual cluster members are very heavily concentrated, driving the behaviour of the whole subsample, and (2) Westerlund 1 is seen through a hole in the extinction, so that the density of foreground stars is also much higher in its vicinity. In fact, the density of stars in the SALSON sample is high over the whole south-western quadrant, including the cluster itself, and much lower in the other three quadrants.

Given this high contamination, some additional information is needed to allow a better selection of cluster members. For this, we resort to the (BP − RP) colour, after our analysis shows that (1) redder stars are much more strongly concentrated on the sky and in the vector point diagram than the rest of the sample and (2) all known members are among the reddest stars in the field. We note, however, that – unlike Beasor et al. (2021) – we do not start from the known members and we do not limit our analysis to stars with similar proper motions. As our primary aim is determining the distance to the cluster, we do not use parallax information in the selection of cluster members. We allow the non-parametric analysis of proper motion and spatial distribution, together with a colour cut, to select the cluster members. Nevertheless, our sample includes most of the members listed in Clark et al. (2020). As can be seen in Table B.1, most of those not included have poor astrometric parameters, thus confirming the validity of our approach. Without the additional information provided by the spatial distribution and colour cuts, the sample would be dominated by Scutum Arm stars whose proper motions are moderately similar to those of cluster members, while their astrometric errors are much smaller. These circumstances very likely explain the results of Aghakhanloo et al. (2021).

To evaluate how effective our member selection is, we consider the catalogued members brighter than RP = 13.5, a range for which Clark et al. (2020) are likely complete and where we can expect astrometric parameters to present reasonable errors. There are 39 such objects in Clark et al. (2020). Our algorithm identifies 27 as cluster members. Of the 12 objects not identified as bona fide members, one (Wd1-W9) has no astrometric solution in EDR3, one (WR T) is outside the $r=3\stackrel{\prime }{.}5$ threshold and eight have RUWE > 1.4. The only two objects with good astrometry that are not considered bona fide members by the algorithm are Wd1-1049 and Wd1-1067, two halo members whose proper motions are compatible with the cluster average, but that just fail to pass the cut in positional probability.

Conversely, essentially all the astrometric bona fide members are part of the cluster population. We find only two interlopers among the > 40 astrometric members brighter than G = 15.5 (all the other ones have counterparts in the catalogue of Clark et al. 2020). The first one is Gaia EDR3 5 940 106 208 947 388 416, which has colours similar to the RSGs in the optical, but is much fainter in the near-infrared. We have a spectrum of this object (our target #238) which shows it to be a very late red giant. The second one is Gaia EDR3 5 940 106 625 576 004 992, which we identify with object F1 in Ritchie et al. (2009a), a field red giant that happens to have proper motions compatible with membership. Therefore, we can safely assume that most of our 401 astrometric members are OB cluster members.

Nevertheless, our sample of bona fide members is very far from complete. There are some bright members with poor astrometry, among which we can cite the RSG Wd1-W75, or the luminous supergiants Wd1-W28, Wd1-41 and Wd1-43b. But the effects of crowding become much more important for fainter stars. To illustrate this point, in Fig. 16, we plot the spatial distribution of our members, divided into a bright (G < 17) and a faint (G ≥ 17) sample. Almost all the stars in the faint sample come from the periphery of the cluster. There is not a single faint member selected in the densest regions, demonstrating how complicated this field is for Gaia. This can be easily visualised by looking at Table B.1, which contains the EDR3 parameters for the cluster members listed in Clark et al. (2020). Stars with G ≈ 17 that lie in crowded regions are lacking parameters (e.g., Wd1-W6b, Wd1-W14c, WR J, WR K) or have large errors in their proper motions that prevent their identification as members (e.g., WR G, WR H), while others of similar brightness in less crowded regions have smaller errors and are picked up by SALSON (e.g., WR Q, or even Wd1-1064, with G = 17.6, and WR W with G = 18.2).

Fig. 16. Spatial distribution of Gaia EDR3 bona fide (circles) and likely (defined as objects picked by SALSON as possible members and having 3.5 ≤ BP − RP < 4.0; triangles) members of Westerlund 1, divided by brightness. Blue symbols represent stars brighter than G = 17, while red symbols are fainter objects. There are only 5 bright likely members. See text for details.

It could be argued that our colour cut at (BP − RP) = 4.0 is to some degree arbitrary, and a slightly bluer or redder threshold could have been chosen. The analysis presented in Sect. 3.3 suggests that cluster members start to be numerous around (BP − RP) = 3.5. There are a few catalogued members with (BP − RP) just below 4 mag (we can cite Wd1-1010, Wd1-1013, Wd1-1020 and Wd1-W228b). They are all located in the south-western tip of the cluster, clearly the area of the lowest extinction. Moving our threshold to (BP − RP) > 3.5 adds about one hundred extra likely members. Given their spatial distribution, heavily concentrated towards the cluster, their membership is very likely. To illustrate this, they have been added to Fig. 16, where they appear as triangles. Interestingly, they are all faint. Except for the four objects mentioned above and one other uncatalogued source, they are all fainter than G = 17. The absence of any bright stars at lower reddening indicates that stars with (BP − RP) < 4 are a minor component of the cluster population. Like their more reddened siblings, these faint objects form a halo around the cluster core. Again, no faint astrometric members are detected in the densest regions, strongly hinting that crowding is preventing the detection of faint stars in the cluster core. The degree of incompleteness in the magnitude range G = 17–19 is very difficult to evaluate, as the quality of Gaia data (whether a star has an astrometric solution, useful uncertainties or even BP and RP magnitudes) depends strongly on the object’s local environment.

In any event, the main source of bias in our sample is the faintness of stars in the BP band, due to the enormous reddening. The CMD for our bona fide members (Fig. 17) clearly shows that we are missing a very large number of stars because they are too faint in BP. This is the only sensible explanation for the diagonal edge to the cluster sequence on its red side, while the (J − K_S) vs. K_S CMD is essentially a vertical strip extending over six magnitudes (e.g., Gennaro et al. 2017). As a rough guide to the level of incompleteness, we cross-matched the list of bona fide members displaying G < 17 with the photometry of Gennaro et al. (2017). We then counted the number of objects lying in the same region of the (J − K_S) vs. K_S CMD as these stars. Taking into account that the photometry of Gennaro et al. (2017) is also incomplete (due to both crowding and saturation in the vicinity of the brightest members), we estimate that the Gaia sample is > 50% incomplete in this luminosity range.

Fig. 17. Cluster CMD for all the Gaia EDR3 bona fide members of Westerlund 1 and spectroscopic members from Clark et al. (2020) not selected by the algorithm. Diamonds represent the most luminous hyper- or supergiants (according to colour, blue, yellow – including the LBV Wd1-W243 – and red). Triangles are luminous B-type supergiants (those classified as Ia). Green triangles have mid-B types, while turquoise triangles have early-B types. The bulk of OB members is represented by circles. The larger blue circles represent the astrometric members, while the smaller brown circles are stars in Clark et al. (2020) not selected by the algorithm. The orange circles represent the candidate high-reddening sample. The red squares are objects that present much redder infrared colour than typical, (J − KS) > 2.2.

6.1. The size of Wd 1

In optical images, Wd 1 has a distinctive shape, resembling a crescent, with two well separate groups, the main one to the north and a smaller aggregate to the south. Contrarily, our selection of members is distributed in an elliptical shape, with a rather smooth distribution, as also found by Gennaro et al. (2017), based on near-infrared star counts. These two disparate observations need to be reconciled.

When we look at the CMD for bona fide members (Fig. 17), we find a very large spread of colours. Leaving aside the strip of cool hypergiants at the top, we find a significant number of objects whose (BP − RP) colour is much higher than the average or mean⁴ for the members (4.7 and 4.8 mag, respectively). To understand this distribution, we selected a sample of candidate OB stars with high reddening, taking the red edge of the cluster sequence, at (BP − RP)≳5.4. We note that this edge may not truly represent the reddest cluster members, as it is likely determined by the observability of objects in the BP band. When we look at the spatial location of these objects (Fig. 18), we find that they are very tightly concentrated: most of them are distributed in a narrow strip between the two main concentrations mentioned in the previous paragraph or immediately to the east of this strip. Such concentration immediately suggests that these objects represent a population with higher extinction.

Fig. 18. Spatial distribution for the objects shown in Fig. 17 (symbols as in that figure). The spatial concentration of the orange circles, presumably high-reddening OB stars, stands out. A few objects of interest discussed in the text are labelled. Wd1-1039 is the only object with two labels, as it is both a Ia supergiant and a highly reddened star (Wd1-1069 is much more heavily reddened, but is not marked with an orange circle, as it lies outside the main distribution). A handful of the members in Clark et al. (2020) are outside the area shown.

To confirm this point, we cross-matched our selection of bona fide members with the near-infrared photometry of Gennaro et al. (2017). When this photometry is plotted, these candidate high-reddening stars display values of (J − K_S)≈1.8–2.0, as opposed to the 1.5–1.6 shown by the bulk of the cluster members. This difference of ≈0.3 mag in E(J − K_S) corresponds, by using the extinction law of Damineli et al. (2016), to an extra extinction of A_V ≳ 2 mag. Given their position in Fig. 17, this additional reddening implies that several of these stars have intrinsic magnitudes comparable to those of the luminous OB supergiants at the top of the cluster sequence. This is confirmed by their K_S < 10 mag. Despite this, most of them were not identified as candidate members by Clark et al. (2005) or subsequent works, presumably because they are not detectable in their UBV photometry. They have never been observed spectroscopically, and thus they were not used by Damineli et al. (2016) to calculate the average extinction to the cluster. As an example, in Fig. 17, we identify the position of Wd1-1039, which is classified as B1 Ia. Unfortunately, this object has no infrared photometry in Gennaro et al. (2017).

Interestingly, there is a significant number of stars that show even redder infrared colours, with (J − K_S)≈2.3. They are also concentrated to the east of the cluster, although not so strongly. Unexpectedly, some of them have Gaia counterparts, which are faint objects of moderately red colours (marked as stripped squares in Figs. 17 and 18). Their nature is unclear. They may be objects with intrinsic infrared excess, such as Herbig Be stars⁵. Alternatively, given their very faint BP ≈ 22 magnitudes and large photometric errors, their optical colours might be unreliable. Spectroscopic observations will be needed to explore these possibilities.

In any event, there is a substantial population of Gaia members that have redder colours than the bulk of the population, indicative of higher extinction. Their spatial distribution leads us to conclude that the appearance of Wd 1 in optical images is driven by differential extinction. The range of extinction values present in the cluster is broader than previously assumed, with many infrared-bright sources simply not visible in the blue. Moreover, since Damineli et al. (2016) used only objects with spectral classification from Negueruela et al. (2010) to calculate individual reddenings and extinctions, their sample is strongly biased to low values. The true average extinction across the face of the cluster must be higher than their estimate. In this respect, it is worth noting that Andersen et al. (2017) found a significantly higher value for the average cluster extinction and a higher dispersion than Damineli et al. (2016) when using a technique based on individual dereddening of a sample of PMS stars observed with HST⁶.

This realisation has strong implications for our knowledge of the cluster basic properties. On the one hand, the combination of Gaia and near-IR colours suggest that about 15 intrinsically luminous OB supergiants had not been previously identified. This compares to the approximately 55 OB supergiants of luminosity class Ia or Iab listed in Clark et al. (2020), resulting in an increase of over 25% in the cluster post-MS population.

Moreover, the assumption used by Beasor et al. (2021) to calculate bolometric luminosities for red and yellow hypergiants of an almost homogenous extinction for the cluster is shown not to hold. In fact, the RSGs Wd1-W20 and Wd1-W75, which these authors find to be of lower luminosity (implying an older population), are very evidently much more reddened than most of the cluster members, as can be seen in Fig. 17. Their near-IR colours, (J − K_S)≈3.8 and 3.6, respectively, are about a whole magnitude redder than those of the two other cluster RSGs⁷. Wd1-W20 is surrounded by some of the most heavily reddened OB stars in the new sample (see Fig. 18). Likewise, from their position in Fig. 17, the YHGs W12a and W265, which Beasor et al. (2021) find to be of low luminosity, are obviously more reddened than other stars of similar types. In contrast, Wd1-W4 or Wd1-W26, which according to Beasor et al. (2021) are much more luminous, happen to be rather less extinguished.

In fact, given the bias against stars with faint BP magnitudes discussed above, we might be missing cluster members with very high reddening. In support of this, we mention the case of the BHG Wd1-1069. This star, which was first located with the observations reported here (our target #355), and later repeatedly targeted during FLAMES monitoring (see Clark et al. 2020; their Fig. 3), lies about 4′ to the south-east of the cluster centre (see Fig. 18 for reference). Since it is outside the radius of $3\stackrel{\prime }{.}5$ imposed to select bona fide members, it is not in our list. Its proper motions are consistent with those of the cluster within their errors. As these uncertainties are large, we cannot decide if it is a halo member or a slow ejection from the cluster. Wd1-1069 is an extremely bright blue hypergiant, which displays (J − K_S) = 2.2 and an extreme (BP − RP) = 7.4 (see its position in Fig. 17), which implies an absolute magnitude around M_KS = −8 at the cluster distance, similarly to its spectral twins Wd1-W7 and Wd1-W33. With BP = 21.1, this object is at the limit of detectability with Gaia, suggesting that any objects fainter in the blue (essentially any, except for other blue hypergiants) and affected by the same amount of reddening will not be detected by Gaia in the blue band.

6.2. The sightline to Wd 1

The STILISM tool⁸ (Capitanio et al. 2017) indicates that in this general direction ($l=339\stackrel{°}{.}5$, $b=-0\stackrel{°}{.}4$) extinction is low at short distances, reaching on average only E(B − V)≈0.3 at 1 kpc and ≈0.5 ± 0.1 at 2 kpc. Extinction is, however, highly patchy, as the bright O9 Ib supergiant HD 151 018, which is projected exactly on top of the cluster and has a Gaia EDR3 distance of $2.{05}_{-0.06}^{+0.08}$kpc (Bailer-Jones et al. 2021), suffers A_V = 2.9 (Maíz Apellániz & Barbá 2018).

To investigate how extinction grows in this direction, we implement the following procedure: from the 2MASS catalogue, we take only sources with high precision photometry, by requiring all three magnitudes to have photometric errors < 0.03 mag. We then select objects with −0.1 < Q_IR < +0.1, corresponding to early-type stars⁹. By imposing high photometric quality, we are more likely to select objects that are in fact of early type (although, as discussed earlier, there can be some degree of contamination by dusty AGBs, especially among the reddest objects), but we are restricting our sample to intrinsically bright stars as we move to higher distances. We obtain a list of 1972 sources, reaching K_S = 13.6, which is the magnitude of a late-B main sequence star in Wd 1. This list is then cross-matched with Gaia EDR3 and the Bailer-Jones et al. (2021) list of inferred distances¹⁰. In Fig. 19, we plot the observed (BP − RP) against the distances from Bailer-Jones et al. (2021). Only objects with BP < 21 are included, as fainter sources lack reliable colours, resulting in a total of 1721 sources.

Fig. 19. Observed colour of stars in the field likely to be of early type (as estimated from their 2MASS colours) against their Gaia EDR3 distances (Bailer-Jones et al. 2021). The thick (red) line is a 200-point running average. There is an obvious decrease in the number of objects around d ≈ 2.8 kpc, which coincides with a steep increase in the average colour. The shaded area represents the running standard deviation. See the text for the many biases present.

Since most stars selected will be of early type, the observed colour should be an approximate measurement of the extinction (bar the intrinsic colour of the star, which ranges from −0.2 to −0.5 mag from late-B to early-O stars). For distances lower than ∼1 kpc, all stars have low extinction. Around d ∼ 1 kpc, the extinction suddenly grows, with no values lower than (BP − RP)≈ 0.5 found any more, and values above (BP − RP)≈ 1.5 appearing for the first time. There is an abrupt change in the distribution at d ≈ 2.8 kpc. The number of sources in the graph abruptly decreases, while stars with (BP − RP) < 2 almost disappear. Given the selection procedure, in particular the high photometric quality required, the only sensible interpretation is the presence of an extinction wall at this distance. This is corroborated by the abrupt increase in the average colour between 2.6 and 3.0 kpc as seen in the running average displayed in Fig. 19. The behaviour of the (J − K_S) colour is identical, although the range of values is much smaller, as expected. Comparison to recent models of galactic structure (see, e.g., Reid et al. 2019; Hou 2021) suggests that this extinction wall at ∼2.8 kpc is associated with molecular clouds located in the Scutum-Centaurus arm, while the smaller rise in extinction around 1 kpc is due to the Sagittarius arm, which is not very prominent along this sightline. Interestingly, the distance to the cluster obtained by Aghakhanloo et al. (2021) coincides with the distance to the extinction wall.

The distribution of extinction may be further probed with our observations and the characteristics of known clusters in the region. Figure 20 shows the classification spectra of a few objects that could be observed with the blue arm of AAOmega. We classify HDE 328 811, whose EDR3 distance is 1.6 kpc (Bailer-Jones et al. 2021), as B2 V. It has UBV photometry from Whiteoak (1963), who estimates E(B − V) = 0.5. HD 150 533 is listed as a Be star in SIMBAD, despite being classified as an O9 Ib supergiant by Garrison et al. (1977). We find a slightly lower luminosity for this object, whose EDR3 distance is 2.5 kpc (Bailer-Jones et al. 2021). UBV photometry is available from a number of sources, starting with Whiteoak (1963), and indicates E(B − V) = 1.0. Cantat-Gaudin et al. (2020) identify a number of moderately distant clusters around our field: NGC 6216 (d = 2.6 kpc, A_V = 2.15), Ruprecht 121 (d = 2.0 kpc, A_V = 2.54), UBC 548 (d = 2.4 kpc, A_V = 2.27) and UBC 669 (d = 2.4 kpc, A_V = 1.9). All these observations confirm the presence of patchy and relatively low extinction out to around 2.5 kpc.

Fig. 20. Spectra of stars in the foreground to Wd 1 whose extinction is sufficiently low to allow blue spectroscopy.

Contrarily, the distances to the high-reddening blue sources identified with our AAOmega observations (see Table C.1 and Fig. 21) are above the discontinuity observed. Unfortunately, all these sources are very faint, and their parallax measurements are subject to large uncertainties. However, only one object, #365, has a nominal distance below 2.8 kpc, namely $2.0_{-0.3}^{+0.7}$kpc (Bailer-Jones et al. 2021), although, with π/e_π = 4.6, this value is highly uncertain. Interestingly, many sources with 4 ≲ (BP − RP)≲5 have nominal distances in the vicinity of 3 kpc. Of particular interest is a close group of three early-B supergiants, #236, #255 and #263, within less than 1′ of each other and having essentially the same proper motions (within their errors), which are probably signposting a young open cluster or association, about 25′ north of Wd 1.

Fig. 21. As in Fig. 19, but for the sample of “blue” stars observed spectroscopically (cluster members are not included). The thick discontinuous lines represent the running standard deviation from Fig. 19, which suggests that the distances in Bailer-Jones et al. (2021) for our sample may be underestimated.

6.3. Cluster radial velocity

The brightest members of Wd 1 are extremely luminous supergiants with complex atmospheric structures. Many of them show substantial spectroscopic variability, ranging from line-profile changes – indicative of pulsation – to changes in spectral types by several subtypes (Ritchie et al. 2009a; Clark et al. 2010). This results in rather significant variations in the RV measured at different epochs (Ritchie et al. 2009a). Furthermore, some of them can even show substantial differences in the RV measured from lines corresponding to different ions, as is the case of the yellow hypergiant Wd1-W265 (Ritchie et al. 2009a). Because of this, defining an average cluster RV value is very difficult, unless a large sample is observed at a given time or – even preferably – over a large number of epochs.

Gaia DR2 provides radial velocities for five bright cluster members, the YHGs Wd1-W4, Wd1-W8a, and Wd1-W265, and the RSGs Wd1-W26, and Wd1-W75. Except for Wd1-W4, all the RVs have large dispersions, reflecting the variability. They are, however, long-term averages and therefore quite reliable. The RSG Wd1-W75 is not part of our Gaia-based sample, as its RUWE > 1.4 removed it from the selection process. The weighted (full weights) average heliocentric radial velocity for the four objects in the bona fide sample is −47.4 ± 0.4 km s⁻¹. After the solar motion correction (Elias et al. 2006), the average radial velocity with respect to the Local Standard of Rest is ν_LSR = −43.0 ± 0.4 km s⁻¹. If weight is not taken into account, the heliocentric mean value is −45.4 ± 1.6 km s⁻¹. Including Wd1-W75 in the calculation of the cluster radial velocity changes the unweighted mean by about −1 km s⁻¹, bringing it closer to the weighted mean of −47.4 km s⁻¹. Although the Gaia RV values are almost certainly very good estimates of the stellar average radial velocities, the formal errors likely underestimate our uncertainty in the definition of a radial velocity for the cluster. Based on the dispersion of individual RVs, we take ±4.0 km s⁻¹ as an upper limit on this uncertainty.

From our AAOmega data, we obtain single-epoch RVs for six cluster members, the four RSGs plus Wd1-W4 and Wd1-W265. Their average is −48.7 km s⁻¹, in excellent agreement with the Gaia DR2 value. Both values are also in good agreement with the RVs measured by Ritchie et al. (2022) from a sample of putatively non-binary early-type stars (which average −43 km s⁻¹), although we cannot make a direct comparison to ensure that their velocities are in the same system as those from Gaia. Given the correction to the Local Standard of Rest in this direction (+4.1 km s⁻¹), we can conclude that the velocity of the cluster is not very far away from v_LSR = −43 km s⁻¹. This value is substantially higher than the −55 km s⁻¹ assumed by Kothes & Dougherty (2007), and, according to the Galactic rotation curve in this direction (e.g., Reid et al. 2019), would place the cluster at a very low distance (≲3 kpc).

However, the gas emission maps displayed by Reid et al. (2019, e.g., their Fig. 3) very strongly suggest that clouds in this direction do not follow closely the rotation curve. On the other hand, a combination of recent studies of various spiral tracers, including ethanol and water masers, draws a l − ν_LSR map (Hou 2021) where the radial velocity obtained for Wd 1 would place it very close to the Norma-Outer arm¹¹. Moreover, the distance that we derive for the cluster also agrees with estimates to the arm in this direction. Given that the extinction wall around 2.8–3.0 kpc is in all likelihood associated with the Scutum-Crux arm, Wd 1 most likely lies in the Norma arm, probably on its near side.

6.4. The field around Westerlund 1: a hive of RSGs

The RSGs in Wd 1 are extremely bright near-IR sources, with Wd1-W26 and Wd1-W237 presenting K ≲ 2.5, Wd1-W20 at K ≈ 3.1, and Wd1-W75 about 0.3 mag fainter (Borgman et al. 1970). Recently, Beasor et al. (2021) have argued that these objects represent an older burst of star formation than the bulk of the early-type stars. Although the assumption of a pretty homogenous extinction towards the cluster may play a major role in this conclusion, the possibility that the cluster is part of a larger association formed along a moderately long timespan is thus worth exploring. To this aim, we looked for stars of similar infrared brightness that could be RSGs related to the cluster in the field surrounding it.

There are 24 stars brighter than K_S = 4 in the (2 degree diameter) field that we observed. Of them, only two are brighter than Wd1-W26 and Wd1-W237, the catalogued Mira variable V823 Ara and the OH/IR star HDE 328 913, which we observed (our target #3), finding that it is a very luminous star, which we classify as M4.5 I. This is the only one among these 24 objects whose proper motions suggest it could have been ejected from the cluster (although it is almost one degree away and a more careful analysis would be necessary to ascertain this possibility). However, its EDR3 parallax, ϖ = 0.58 ± 0.04 mas seems to completely rule out this possibility. In fact, only 5 of the 24 sources with K_S < 4 have a Bailer-Jones et al. (2021) distance ≳3 kpc. These four objects have proper motions roughly in agreement with the cluster averages, but they are all located > 20′ away from the cluster centre. We have a spectrum for one of them, #390, an M3.5 Ib supergiant with RV = −50 km s⁻¹. Another object of interest for which we have a spectrum is #152, with K_S = 3.2 ± 0.3. This star has essentially the same ppms as the cluster, but its parallax ϖ = 0.54 ± 0.14 places it more than 2 σ away from the cluster average. However, given that ϖ/e_ϖ is only 4.0, its high luminosity (M3 Ia), and a radial velocity similar to those of cluster members (RV = −44 km s⁻¹), we cannot fully rule out a connection with the cluster until a more accurate value of the parallax is available. In any event, this object is > 20′ away from the cluster.

As we move towards fainter near-IR magnitudes, we fail to find the combination of high luminosity, RV and proper motion that would allow us to claim a connection to the cluster. In Fig. 22, we carry out a simple experiment, by subtracting from the proper motions of all the RSGs that we have observed the cluster average values found above. As can be seen, many of the RSGs have proper motions differing significantly from those of the cluster (although we must note that errors are in some cases very large). The very few objects for which this difference in proper motions might be compatible with an ejection from the cluster are very faint in K_S (e.g., the K1 Ia star #214, with K_S = 6.82 ± 0.02) or not very luminous (e.g., the M1 Ib object #248 or the M2 Ib #334). Therefore, we find little to no evidence of RSGs ejected from the cluster either in our spectroscopic sample or in the sample of stars brighter than K_S = 4. Likewise, we find no evidence for an extended population surrounding the cluster that could be indicative of a large OB association centred on the cluster à la Per OB1 (see, e.g., de Burgos et al. 2020). A few of the RSGs and some of the blue stars have proper motions and parallaxes that agree very well with those of Westerlund 1, and may be indicating the presence of a diffuse OB association with a range of ages, although they could also be part of the general population of massive stars associated with the Norma arm. In any event, their presence very strongly suggests that the proper motions of Westerlund 1 are indistinguishable from those of stars in its neighbourhood. Since they are also quite similar to those of the much more numerous population associated with the foreground Scutum-Crux arm, we are tempted to conclude that proper motions in this direction mainly reflect the relative motion between stars following their Galactic orbits and the Sun.

Fig. 22. Relative proper motions of massive stars observed with AAOmega. The arrows are proportional to the difference between the proper motions of the star and the average proper motions of the cluster. Blue dots represent our “blue” targets. Red dots are objects that we classify as red supergiants.

Despite this lack of connection to the cluster, the large population of RSGs found is remarkable in itself. Among the hundreds of red luminous stars observed, close to 40 have spectral classifications compatible with being evolved massive stars. As a comparison, a circle of radius 1° centred on χ Per (i.e. the densest region of Per OB1) contains nine RSGs. Such a bountiful harvest of RSGs is testimony to the richness of the Scutum-Crux arm in this direction. Within the small field covered by our AAOmega observations, we find a very substantial population of massive stars, evidence for ongoing star formation and thick dark clouds that obscure our view to Westerlund 1.

7. Conclusions

Our analysis of Gaia EDR3 data in a circle of radius 12′ surrounding Westerlund 1 reveals the extraordinary complexity of this field. We find extreme crowding, huge variations in foreground extinction and a very complex kinematic configuration. As a consequence, the direct application of a non-parametric method to the positional and kinematic data of stars does not provide us with a clean sample of cluster members, but needs further pruning based on the distance to the cluster centre and the (BP − RP) colour of the stars. Selecting the reddest stars ((BP − RP) > 4) and restricting ourselves to the central $3\stackrel{\prime }{.}5$ of the circle, we find that:

1. The population associated with the cluster is well separated. The Bayesian membership analysis leads to 401 bona fide cluster members, 142 of them with parallaxes corrected from EDR3 systematic zero-points. This sample includes the majority of spectroscopically catalogued members in Clark et al. (2020).

2. The parallax of the cluster is estimated to be ϖ_c, t = 0.238 ± 0.012 mas. This value is statistically indistinguishable from that obtained by Beasor et al. (2021) for a selection, based on proper motion analysis, of catalogued OB members with data in Gaia EDR3 (ϖ_c, t = 0.243 ± 0.027 mas). Different selection methods and samples of Wd 1 members lead to the same results for the cluster parallax, using the best available astrometric data. Nevertheless, we note that the parallax corrections for cluster members are generally quite extreme (typically, −0.08 mas), and therefore systematic effects cannot be completely ruled out. The analysis presented by Maíz Apellániz (2022) suggests that a large systematic is unlikely for the range of colours and magnitudes found in the cluster, but this needs to be confirmed. A distance of 4.2 kpc likely places Wd 1 in the Norma arm.

3. The cluster distributes on the plane of the sky (RA, Dec) with an elliptical shape, where the principal diameter is ≈7 arcmin. This value must be considered an upper limit on the cluster size, which would lead to a major semiaxis of 3.5 arcmin, well above the radius determined by other authors (see, for example, Dias et al. 2002, who give a value $R=1\stackrel{\prime }{.}2$). The radius of a cluster is a fuzzy concept that admits different definitions and approaches (i.e. Sánchez et al. 2020). If we assume that the minor axis of the ellipse is $2\stackrel{\prime }{.}5$, the effective radius would be close to 3′, which at the distance estimated for Wd 1 corresponds to 3.7 pc.

4. The proper motion of the cluster derived from the 401 bona fide cluster members is μ_α(c) = − 2.231 and μ_δ(c) = − 3.697 mas/a, which at the distance of 4.2 kpc, yields V_α(c) and V_δ(c) space velocities of 44.5 ± 3.2 and 73.8 ± 5.3 km s⁻¹, respectively.

5. There are indications of the presence of a low-density halo surrounding the cluster, as many of the very red stars plotted in Fig. 2 have proper motions compatible with those of the cluster (see Fig. 3). However, we also find stars with fully compatible astrometric parameters at much higher angular distances (for instance, stars #246, #359, #389, #395, and #421 in Table C.1 have parallax and proper motions fully compatible with membership within their errors, but are all located > 20′ away from the cluster). This suggests that the proper motions of Wd 1 are not distinguishable from those of the field population in the Norma arm with the current level of accuracy. More accurate proper motions, such as those in DR4 (perhaps helped by spectroscopy), will be needed to prove the existence of a halo and search for runaway stars ejected from the cluster.

6. Gaia DR2 provides radial velocity values for four stars considered as members of the cluster. The weighted average of the radial velocities provides different mean values depending on the chosen weights, in such a way that we adopt a heliocentric radial velocity of −45.4 km s⁻¹ with an internal precision of 1.4 km s⁻¹, corresponding to the unweighted average estimation. This may be an underestimation of the true uncertainty in the RV value, despite the quality of Gaia data, given the small number of objects for averaging, and the large number of RV variable stars. This value is fully consistent with the average RV of 58 OB stars without convincing signs of binarity (−43.1 km s⁻¹; Ritchie et al. 2022).

7. Many of the astrometric members detected have colours indicative of higher reddening than most of the members listed by Clark et al. (2020). The majority of these objects are located very close together in the sky, in a region to the east of the cluster centre and a strip that extends between the two main groupings seen in optical images. About 15 of these very red objects have infrared magnitudes and colours that identify them as luminous OB supergiants. This finding implies that extinction across the face of Wd 1 is much more variable than previously assumed. As a consequence, the average extinction and its dispersion are likely higher than found by Damineli et al. (2016), in agreement with Andersen et al. (2017).

8. The sightline to Wd 1 goes through an extremely rich portion of the Scutum-Crux arm, which contains a huge population of massive stars (see Appendix C). The distribution of sources in this direction suggests that this population is mostly located at a distance ∼3 kpc, behind a wall of dark clouds that contribute substantially to the extinction to the cluster.

In all, Gaia EDR3 data provide us with a much improved view of the cluster properties. Future data releases will likely sharpen this view.

Acknowledgments

This paper is dedicated to the memory of Simon Clark, who tragically passed away as this work was beginning. For twenty years he and the first author, together with many others, explored Westerlund 1. The expedition is far from over, but the captain is no longer at the helm. All new crews that take up this task will be sailing on his trail. We thank the anonymous referee for helpful suggestions. This research is partially supported by the Spanish Government Ministerio de Ciencia, Innovación y Universidades under grants PGC2018-093 741-B-C21/C22 and PGC2018-095 049-B-C21/C22 (MICIU/AEI/FEDER, UE) and the Generalitat Valenciana under grant PROMETEO/2019/041. EJA acknowledges financial support from the State Agency for Research of the Spanish MCIU through the “Center of Excellence Severo Ochoa” award to the Instituto de Astrofísica de Andalucía (SEV-2017-0709)”. RD acknowledges further support from the Spanish Government Ministerio de Ciencia e Innovación through grant SEV 2015-0548, and from the Canarian Agency for Research, Innovation and Information Society (ACIISI), of the Canary Islands Government, and the European Regional Development Fund (ERDF), under grant with reference ProID2 017 010 115. The AAOmega observations have been supported by the OPTICON project (observing proposal 2011A/014), which was funded by the European Commission under the Seventh Framework Programme (FP7). This research has made use of the Simbad, Vizier and Aladin services developed at the Centre de Données Astronomiques de Strasbourg, France. 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. The Gaia data are processed with the computer resources at Mare Nostrum and the technical support provided by BSC-CNS. In this work, we have made extensive use of TOPCAT (Taylor et al. 2005). We thank its author, and subsequent contributors, for the creation and development of this tool.

References

Appendix A: Is there a cluster in front of Westerlund 1?

The nature of the BG is an open question. The open cluster vandenBergh-Hagen 197 (MWSC 2458; BH 197; van den Bergh & Hagen 1975) is catalogued to lie in the Northern part (in Galactic coordinates) of the area occupied by the BG. Catalogue data for this object show disparate values for the basic physical properties. While the two studies that have recently estimated the distance and age of this cluster coincide (Kharchenko et al. 2013; Dias et al. 2014) in placing it between 1.6 and 2 kpc and assign a similar extinction, E(B − V)∼1, they differ widely in the estimation of age: Kharchenko et al. (2013) dates it close to 1 Ga, while Dias et al. (2002) give 30 Ma, which seems more compatible with the colour observed for the BG. With this information, we can speculate that the area west of Wd 1, at least within the radius of 12′ chosen for this study, may contain a star formation region with a minimum size of 9 pc, which includes the cluster BH 197 surrounded by a halo of low-density blue stars at an average distance of 2 kpc.

The membership analysis performed by Sampedro et al. (2017), based on UCAC4 (Zacharias et al. 2013) data, and on the previous work by Dias et al. (2014), selects a set of 274 members (considered as such by the four methods used in that work) within a radius of 4′ around the coordinates of the cluster centre given by Kharchenko et al. (2013).

Of the 274 BH 197 members catalogued by Sampedro et al. (2017), 253 are in our initial catalogue. The G vs (BP − RP) CMD and proper-motion VPD drawn with Gaia EDR3 data for this sample show distributions far away from those expected for a stellar cluster. The CMD presents different sub-groupings at different colour ranges that cannot be explained by a stellar system with a single stellar population, even advocating a very complicated extinction pattern. The same can be said for its distribution in proper motion space. Members span more than 10 mas/a along the main axis, and are far from resembling what is expected for a star cluster at the distance of 2 kpc. Thus, we conclude that BH 197 cannot be considered an open cluster in the classical sense of the term, but rather, in the best event, a subregion of a larger star-forming region that could extend beyond 9 pc, this lower limit being defined by the sampling radius of our data.

The distribution of YSOs along the line of sight, which crosses the Galactic plane in the direction l ≈ 340°, shows a clumpy pattern (Kuhn et al. 2021), with groupings of YSOs (identified from Spitzer’s catalogue) associated with the Carina-Saggitarius and Scutum-Centaurus arms, but with some clustering in the inter-arm region. Our BG appears to be located in the vicinity of these concentrations and could represent the more evolved optical counterpart of the ongoing star-formation process. There is an additional argument that supports the hypothesis that the BG is a relatively extensive, active star-forming region. Dense-gas hub-filament systems (HFSs) appear to be good tracers of high-mass star-formation activity (Schneider et al. 2012; Peretto et al. 2013; Kumar et al. 2020). One HFS (HFS14082) detected by Kumar et al. (2020) is located within the area covered by the BG. This hub is at a distance compatible with that estimated for our Blue Group, especially if we take into account the different distance estimation methods and their associated uncertainties.

Appendix B: Cluster members

In Table B.1 (available at the CDS), we list the members of Wd 1 with spectroscopic observations in the catalogue of Clark et al. (2020) and their Gaia EDR3 basic information. The coordinates in Clark et al. (2020) were collected from a variety of sources, and had different accuracies. A majority (114) have accurate coordinates, and the nearest EDR3 source is within $0\stackrel{″}{.}5$. A few (8) are > 1″ away from the nearest EDR3, but all are within $1\stackrel{″}{.}5$. The only exception is the RSG Wd1-W20 for which Clark et al. (2020) erroneously list the same coordinates as for Wd1-W237.

In addition to the members already known, we list the newly confirmed member W1070. This object is outside the $3\stackrel{\prime }{.}5$ region where we have selected bona fide members, but its astrometric values are typical of cluster members. Although, its spectrum is very poor, the narrow Paschen lines and presence of a very prominent C III 8502 Å line identify it as an O9 I star.

Appendix C: Stars observed with AAOmega

In Tables C.1 and C.2 (available at the CDS), we list all the stars that were observed with AAOmega, divided according to their nature. Table C.1 contains blue stars and a handful of yellow supergiants. Many of these objects are quite faint, and their spectra are quite noisy (sky subtraction is worse in spectra with low signal-to-noise). For some, the only features seen are broad and shallow Paschen lines, which identify them as OB stars (mid and late B stars have deeper Paschen lines; see Negueruela et al. 2010). In a few, a moderately strong C III 8502 Å line is seen, signalling a luminous star with spectral type close to O9.

Our targets include 9 B-type supergiants, 3 A-type luminous supergiants and five F-type supergiants. This is a surprisingly high number, considering the small field surveyed. The high number of supergiants observed indicates that we are only observing the intrinsically most luminous early-type stars. Most of them have a 2MASS (J − K_S) colour > 1.5, indicative of very high extinction – the intrinsic colour of all these objects is close to (J − K_S)≈0, with only the F-type supergiants approaching 0.3. A few of the most luminous supergiants have colours (J − K_S) > 2.0, suggesting, and parallaxes compatible with being, background objects behind the Norma arm.

Table C.2 contains cool stars (with spectral types G and later, corresponding to our "red" targets). Targets were ordered by their brightness in the i or I band. This is reflected in the quite strong correlation between ID number and G band brightness. Spectral types have been estimated following the criteria discussed in Negueruela et al. (2012). The quality of the spectra differs greatly, depending on many factors. For average count rates below about 1000, sky subtraction becomes substantially worse, and many of the fainter targets only have approximate classifications, or simply an M classification that indicates that TiO bands are visible, but the quality of the spectrum prevents further classification.

C.1 Some interesting early-type stars

Star #412 is coincident with the mid-IR source IRAS 16475−4609. Based on the morphology of the WISE images and presence of the nearby candidate young stellar object IRAS 16474−4610, #412 could be the ionising star of a small H II region at a very high distance.

Star #168 is a luminous B5 Ia supergiant (see Fig. 15). It is located approximately $2\stackrel{\prime }{.}5$ to the north of the compact group formed by #236, #255 and #263, and its proper motions, though similar to those of the other three, are compatible with having travelled from their vicinity. All these stars have parallaxes compatible with a distance ∼3 kpc.

Target #119 is the central star of the ring nebula [GKF2010] MN48 (Gvaramadze et al. 2010). These authors proposed it as a candidate LBV star. Wachter et al. (2010) classified it as an unusual emission-line star, based on H- and K-band spectra. Kniazev et al. (2016) presented comprehensive photometry and spectroscopy of the source, confirming its LBV nature. Based on the evolution of brightness and colours, Kniazev et al. (2016) argued that the star become cooler and more luminous between 2009 and 2011, and then faded as its temperature increased again. If this interpretation is correct, our spectrum was taken close to maximum light. The spectrum is shown in Fig. C.1, compared to that of the known cluster LBV Wd1-W243. The I-band spectrum of Wd1-W243 is very similar to those reported for 2005 – 2009 by Ritchie et al. (2009b, their fig. 11), although the emission cores of the Paschen lines seem to have weakened, an effect that is more clearly seen in Pa 11, which is essentially a pure absorption profile now. The absorption spectrum is dominated by N I lines, which correspond well with a 8 500 K hypergiant model (Ritchie et al. 2009b, their fig. 14).

Fig. C.1. AAOmega spectrum of the LBV MN48, compared to that of Wd1-W243. In both stars, the strongest absorption lines correspond to N I. The Paschen lines in Wd1-W243 are almost entirely in absorption, while those in MN48 display P-Cygni profiles. See text for a description.

The spectrum of MN48 is dominated by much stronger emission lines. The Paschen lines show weak P-Cygni profiles. The Ca II lines have much stronger P-Cygni profiles, although at this resolution it cannot be discarded to attribute the absorption profiles to the broader Paschen lines. The N I spectrum in absorption is somewhat weaker than in Wd1-W243, but still compatible with an early-A super/hypergiant. A low-resolution optical spectrum taken two months before our observations is dominated by emission lines corresponding to singly ionised metals, mainly Fe II (Kniazev et al. 2016). According to these authors, the star became noticeably hotter over the following four years. A spectrum taken in 2015 shows the upper Paschen lines essentially in absorption, while the O I 8446 Å line, powered by Lyman β fluorescence, remains strong (peaking around 1.4 times the continuum level, as in our spectrum). Our spectrum also shows a strong Fe II 8490 Å line, which is pumped by Lyman α.

The proper motions of MN48 are almost identical to the cluster averages, but – given the moderate errors – may indicate a small difference around 0.3 mas in pmRA. As the star is located ∼25′ to the east and ∼9′ to the north, travelling this distance would require around 5 Ma, which makes association with Wd 1 possible. Given the range of ages currently being considered for Wd 1, the star should have been ejected early in the history of the cluster. Kniazev et al. (2016) measure an RV =  − 37 km s⁻¹, which is also compatible with a slow ejection from Wd 1. It is worthwhile commenting that Gvaramadze (2018) identifies another LBV (MN44) that may also have been ejected from Wd 1 (in approximately the opposite direction) and that also requires a long (∼4 Ma) travel time. If these associations are confirmed by more accurate proper motions, they would imply that both stars were ejected early in their lives, when they were O-type main sequence stars, and evolved into LBVs close to their current locations.

All Tables

All Figures

	Fig. 1. G vs. (BP − RP) colour-magnitude diagram for the total sample with RUWE ≤ 1.4 (grey small dots). Stars with spectroscopic information in Clark et al. (2020) and Gaia EDR3 photometry are shown in red. This sample is clearly seen to consist mostly of stars redder than (BP − RP) > 4.
In the text

	Fig. 2. Red Branch stars, defined as those with (BP − RP) > 4, plotted on the plane of the sky (dots). Coordinates are angular distance to the cluster centre along RA (X axis) and Dec (Y axis). Red squares identify stars with spectra reported in the catalogue of Clark et al. (2020). The vast majority of red squares are within a radius of $\approx 3\stackrel{\prime }{.}5$.
In the text

	Fig. 3. Isodensity contours (black line) of the total sample (grey small dots). Red dots are RB stars, which are in an eccentric position with respect to the whole sample distribution. They appear centred on one focus of the internal dispersion ellipse.
In the text

	Fig. 4. Histogram showing the parallax distribution for the total sample (grey) and the RB stars (red). Frequency is normalised to the maximum value for better visualising the difference in modes. Attending to the histogram, this difference is 0.25 mas (bin size). Kernel density estimators lead to a slightly shorter separation of ≈0.20 mas.
In the text

	Fig. 5. DP pdfs on the X–Y plane for the probable clusters members selected by SALSON. DP(c) in blue and DP(f) in green. DP(c) and DP(f) show similar values for objects well outside the central core, which leads to larger uncertainties in the estimated cluster membership probabilities.
In the text

	Fig. 6. DV pdfs in the proper motion space for the probable clusters members selected by SALSON. DV(c) is in blue and DV(f), in green. DV(c) shows a non-single peak distribution with a rounded bump at lower proper motions. This is likely indicative of a field contamination not removed by SALSON. Possible causes are discussed in the text.
In the text

	Fig. 7. Wd 1 radial density profile. The DP pdfs versus angular distance (arcmin) to cluster centre are shown for the members selected by SALSON. DP(c) is shown in blue and DP(f), in green. We identify as most likely cluster members those stars within the upper-left rectangle (red lines): DP(c) > 0.002 and r < 3.5 arcmin.
In the text

	Fig. 8. Spectrum of the spatial segregation by (BP − RP) colour. Error bars are σ values. Λ–(2 × σ) values higher than 1 are indicative of colour bins spatially grouped. A strong spatial segregation is observed for stars redder than (BP − RP) > 4. A blue grouping is also detected with (BP − RP) < 1.5.
In the text

	Fig. 9. Q spatial concentration parameter as a function of colour. Error bars have been obtained by bootstrapping. Q ≈ 0.8 separates the clumpy and central spatial distributions. The higher Q is, the more spatially concentrated the distributions. This plot reinforces the conclusions derived from Fig. 8.
In the text

Fig. 10. Stars (dots) and isodensity contours (lines) for the two groups spatially segregated by (BP − RP). BG (in blue) contains the objects in segregated bins with colour bluer than 1.5. RG (in red) is the same for colours redder than 4. The RG is mainly associated with the cluster centre, while the BG shows a stellar population with its maximum density well separated from the core of Wd 1. The axes show angular distances in equatorial coordinates.

In the text

	Fig. 11. Stars and isodensity contours in proper motion space for the two groups spatially segregated by (BP − RP). Symbols and colours as in Fig. 10. The spatial segregation by colour is mimicked in the VPD. Although both distributions show a non-null overlap, their density maxima are well separated in proper motion.
In the text

Fig. 12. Histogram of corrected (see Sect. 4.1 for parallax bias corrections) Gaia EDR3 parallaxes for the RG (red) and BG (blue) groups, showing two different statistical distributions. RG shows smaller parallaxes than BG. This result, together with those shown in Figs. 8–12, strongly suggests that both populations are well separated by colour, position, kinematics, and distance. Gaussian Kernel density estimators are also drawn.

In the text

	Fig. 13. Histogram of parallax zero-point corrections according to Lindegren et al. (2021a). Bona fide members are in blue, and the global sample is in red. The histograms have been normalized to their maxima for better visualizing the different distribution modes.
In the text

	Fig. 14. Cluster CMD showing the value of the zero point applied to each one of the bona fide members of Wd 1 for which such correction has been estimated. Left panel: Lindegren. Right panel: Maíz Apellániz. The scale is in μas and larger symbols have larger full weights.
In the text

Fig. 15. AAOmega spectra of luminous supergiants. All the spectra have been normalised and the continuum has been shifted for clarity. Left panel: members of Wd 1. The (green) dashes above the spectrum of Wd1-W57a indicate the positions of Paschen lines. The three thin (orange) vertical lines show the position of the Ca II triplet. The short (red) dashes below the spectrum of Wd1-W265 identify some of the strongest N I lines. The dotted line marks the position of the 8621 Å DIB, whose strength correlates well with extinction up to a saturation value. Right panel: other stars in the field, not connected to the cluster. Symbols as in the left panel.

In the text

	Fig. 16. Spatial distribution of Gaia EDR3 bona fide (circles) and likely (defined as objects picked by SALSON as possible members and having 3.5 ≤ BP − RP < 4.0; triangles) members of Westerlund 1, divided by brightness. Blue symbols represent stars brighter than G = 17, while red symbols are fainter objects. There are only 5 bright likely members. See text for details.
In the text

Fig. 17. Cluster CMD for all the Gaia EDR3 bona fide members of Westerlund 1 and spectroscopic members from Clark et al. (2020) not selected by the algorithm. Diamonds represent the most luminous hyper- or supergiants (according to colour, blue, yellow – including the LBV Wd1-W243 – and red). Triangles are luminous B-type supergiants (those classified as Ia). Green triangles have mid-B types, while turquoise triangles have early-B types. The bulk of OB members is represented by circles. The larger blue circles represent the astrometric members, while the smaller brown circles are stars in Clark et al. (2020) not selected by the algorithm. The orange circles represent the candidate high-reddening sample. The red squares are objects that present much redder infrared colour than typical, (J − KS) > 2.2.

In the text

Fig. 18. Spatial distribution for the objects shown in Fig. 17 (symbols as in that figure). The spatial concentration of the orange circles, presumably high-reddening OB stars, stands out. A few objects of interest discussed in the text are labelled. Wd1-1039 is the only object with two labels, as it is both a Ia supergiant and a highly reddened star (Wd1-1069 is much more heavily reddened, but is not marked with an orange circle, as it lies outside the main distribution). A handful of the members in Clark et al. (2020) are outside the area shown.

In the text

Fig. 19. Observed colour of stars in the field likely to be of early type (as estimated from their 2MASS colours) against their Gaia EDR3 distances (Bailer-Jones et al. 2021). The thick (red) line is a 200-point running average. There is an obvious decrease in the number of objects around d ≈ 2.8 kpc, which coincides with a steep increase in the average colour. The shaded area represents the running standard deviation. See the text for the many biases present.

In the text

	Fig. 20. Spectra of stars in the foreground to Wd 1 whose extinction is sufficiently low to allow blue spectroscopy.
In the text

	Fig. 21. As in Fig. 19, but for the sample of “blue” stars observed spectroscopically (cluster members are not included). The thick discontinuous lines represent the running standard deviation from Fig. 19, which suggests that the distances in Bailer-Jones et al. (2021) for our sample may be underestimated.
In the text

	Fig. 22. Relative proper motions of massive stars observed with AAOmega. The arrows are proportional to the difference between the proper motions of the star and the average proper motions of the cluster. Blue dots represent our “blue” targets. Red dots are objects that we classify as red supergiants.
In the text

	Fig. C.1. AAOmega spectrum of the LBV MN48, compared to that of Wd1-W243. In both stars, the strongest absorption lines correspond to N I. The Paschen lines in Wd1-W243 are almost entirely in absorption, while those in MN48 display P-Cygni profiles. See text for a description.
In the text