© E. Zari et. al 2021

1. Introduction

Luminous and hot stars are massive, and hence rare and short lived. However, they play decisive roles across different fields of astrophysics. They dominate the interaction between stars and the interstellar gas and dust (interstellar medium; ISM) in their host galaxies by heating or ionising those components; those stars also interact with the ISM through powerful stellar winds and eventually supernova explosions (Mac Low & Klessen 2004; Hopkins et al. 2014). Luminous and hot stars must indeed have played an important role in galaxy evolution throughout cosmic time, via their intense winds, ultraviolet radiation fields, chemical processing, and explosions (Haiman & Loeb 1997; Douglas et al. 2010; Bouwens et al. 2011).

Luminous and hot stars are inevitably young, and therefore they can serve as tracers of recent massive star formation. Although they make up an insignificant fraction of the overall stellar mass, they contribute a major portion of the light of the disc and can therefore probe the spiral structure and young disc kinematics of our Galaxy and others (Xu et al. 2018; Chen et al. 2019; Dobbs & Baba 2014; Kendall et al. 2011, 2015).

Also, these objects are – directly or through their remnants – decisive drivers of the chemical evolution for many elements in the periodic table and serve as crucial laboratories for stellar evolution in this arguably least well-understood regime of stellar physics.

Massive stars are born predominantly as members of binary and multiple systems (Sana et al. 2012, 2014; Kobulnicky et al. 2014; Moe & Di Stefano 2017). As a consequence, most of them are expected to undergo strong binary interaction, which drastically alters their evolution (Podsiadlowski et al. 1992; Van Bever & Vanbeveren 2000; O’Shaughnessy et al. 2008; de Mink et al. 2013; Langer et al. 2020).

Finally, massive (binary) stars are the only channel to yield binaries that involve black holes and neutron stars in the disc of the Milky Way. Therefore, understanding massive star binaries, including as they evolve away from their zero-age main sequence, is indispensable for understanding the distribution of gravitational wave events. For all of the above reasons, a multi-epoch, spectroscopic census of luminous and hot stars across the Galaxy is needed, providing spatial and dynamic information together with estimates for masses, ages, metallicity, multiplicity, and other spectroscopic information. The SDSS-V survey (Kollmeier et al. 2017, and in prep.; Bowen & Vaughan 1973; Gunn et al. 2006; Smee et al. 2013; Wilson et al. 2019) will provide such a comprehensive, multi-epoch spectroscopic program on hot and massive stars, but this program must be based on a clear and quantitative selection function to enable rigorous subsequent population studies.

There is no universal precise ‘definition’ across different sub-communities within astrophysics of when a star is ‘luminous’, ‘hot’, and ‘massive’. ‘Hot’ may either mean that its T_eff is sufficiently high that the spectrum is dominated by H and He lines, rather than metal lines (i.e. OBA stars); or it may mean that T_eff is sufficiently high that the star produces significant amounts of ionising radiation (O and early B stars). ‘Massive’ may either mean, massive enough to ‘go supernova’ or to have a convective core.

For the current context we choose an operative definition of luminous and hot as T_eff ≳ 8000 K and M_K ≲ 0 mag, or – loosely speaking – OBA stars. This is the set of stars that addresses the science sketched out above. Some aspects, such as ISM ionisation, are only affected by O and early B stars; but it is important to remember that even the most massive stars do not produce ionising radiation during all of their evolution. Therefore we aim in our sample definition for a sample that is hot enough to eliminate all the luminous red giant branch (RGB) and asymptotic giant branch (AGB) stars, but include most evolved phases of massive stars. Stars with T_eff ≥ 8000 K can be photometrically pre-selected across the entire sky; yet efficient discrimination between 10 000 K and 20 000 K by means of photometry only is almost impossible in the presence of dust.

In light of the science goals above, in this paper we present the selection of the target sample of massive stars of the SDSS-V survey. We describe an approach to use kinematics to improve distance estimates. Finally we use a ‘clean’ subset of the target sample to study the structure of the Milky Way disc as traced by young, massive stars.

The selection of our target sample is based on the combination of Gaia EDR3 (Gaia Collaboration 2016, 2021a) astrometry and photometry and 2MASS (Skrutskie et al. 2006) photometry. Similar samples of stars in the upper main sequence (UMS) were already selected by Poggio et al. (2018) and Romero-Gómez et al. (2019) to characterise the structure and the kinematic properties of the Galactic warp. Poggio et al. (2018) used a combination of Gaia DR2 and 2MASS colours to select stars of spectral type earlier than B3, and obtained a sample of 599 494 stars. Romero-Gómez et al. (2019) used only Gaia DR2 photometry to select stars brighter than M_G, 0 = 2 mag (thus likely including also A-type stars), and obtained a sample of 1 860 651 stars. Our sample definition (described in Sect. 2) lies somewhat in the middle between those by Poggio et al. (2018) and Romero-Gómez et al. (2019). In Sect. 3 we describe the method that we use to estimate distances. In Sect. 4 we clean the target sample of sources with spurious astrometry and we filter out older contaminants. In Sect. 5 we use this filtered sample to study the 3D space distribution of the sources and we compare the distribution of OBA stars with other tracers of massive star formation and spiral arm structure. In Sect. 6 we discuss our findings in the context of the structure and nature of the spiral arms of the Milky Way. Finally, we summarise our results and draw our conclusions in Sect. 7.

2. Target sample

In this section we present the criteria that we applied to select the target sample for the SDSS-V spectroscopic survey, and we describe its characteristics in terms of sky distribution, magnitude and luminosity distribution, variability, purity, and completeness.

2.1. Astro-photometric selection

To select our sample, we use Gaia EDR3 photometry (Riello et al. 2021) and astrometry (Lindegren et al. 2021a) combined with 2MASS photometry. The cross-match between Gaia and 2MASS is not yet provided in the Gaia archive (Marrese et al., in prep.)¹. For this reason we first performed a cross-match with Gaia DR2, followed by a cross-match with 2MASS. The ADQL query for this cross-match is provided in Appendix A.

We restrict our query to the sources with G < 16 mag. The motivation behind this condition is twofold. On the one hand, it is a SDSS-V technical requirement: at G = 16 mag we can obtain a S/N = 75 with 15 min exposures with the BOSS spectrograph (Smee et al. 2013). On the other hand, for G > 16 mag the Gaia EDR3 parallax precision rapidly deteriorates (see e.g. Table 3 in Gaia Collaboration 2021a), thus (a) estimating absolute magnitudes (which are crucial to our selection) becomes nontrivial and (b) inferring distances strongly depends on the choice of the prior (see Appendix B and Sect. 3).

We define a proxy for the absolute magnitude of a star in the K_s band, , and we require . This translates to the following parallax condition:

(1)

and aims at selecting a reasonably sized sample of bright stars. The value M_Ks = 0 mag corresponds roughly to a B7V-type star (Pecaut & Mamajek 2013). To include also later B-type stars we should require , as M_Ks = 1 mag roughly corresponds to the absolute magnitude of a A0V-type star. We chose however as a good compromise between completeness and contamination (see Sect. 2.3).

We stress that for the survey target list we do not require any condition on the photometric quality of 2MASS, or on the photometric errors, the parallax accuracy, or the goodness of fit of the astrometric solution. This is motivated by the fact that we are aiming to obtain a known, well-defined selection function: not only do measurement errors not describe intrinsic physical properties of the sources, and therefore should not enter the selection criteria, but they are also difficult to model.

At this point our sample, shown in the (J − K_s)₀ versus G_BP − G_RP colour–magnitude diagram in Fig. 2, consists of luminous stars, which could either be OBA, RGB, or AGB stars. OBA stars fall in the region highlighted by the orange ellipse. The stars in the horizontal stripes are stars with large 2MASS photometric errors or 2MASS photometric quality flags (ph_flag) different from “AAA”.

Fig. 1. Two-dimensional histogram of J − Ks vs. G − Ks (left), (J − Ks)0, and G − Ks (right) for the stars in the LAMOST sample by Liu et al. (2019); see main text for more details on the sample.

Fig. 2. (J − Ks)0 vs. GBP − GRP colour–magnitude diagram. OBA stars roughly fall in the region highlighted by the orange ellipse. This representation shows that hot stars stand out qualitatively in colour–colour space.

To select OBA stars we use simple photometric cuts, roughly following the procedure outlined in Poggio et al. (2018, see in particular their Fig. 1). The first step of our selection consists in selecting sources with:

(2)

The first condition excludes most of the red sources, and the second condition removes sources with unnaturally blue colours. We obtained the (J − K_s)₀ colour by applying the following Equation:

(3)

which we derived by noticing that, for the spectroscopically confirmed O- and B-type stars selected by Liu et al. (2019), the (J − K_s) and (G − K_s) colours are linearly correlated (see Fig. 1, left), and by assuming that the difference between the J and K_s magnitudes for such stars is close to zero (see e.g. Pecaut & Mamajek 2013).

Figure 3 shows the distribution of sources in the G − K_s versus J − H colour–colour diagram. O- and B-type stars lie on a sequence in the J − H versus G − K_s colour–colour diagram as a consequence of interstellar reddening and are clearly separated from redder turn-off stars and giants. We therefore select stars located in the region defined by the equations:

(4)

Fig. 3. J − H vs. G − Ks colour–colour diagram of the sources with (J − Ks)0 < 0.1 mag and (J − Ks)0 > −0.3 mag. The solid lines delineate the criteria for selecting OBA stars and are defined in the text.

(solid grey lines in Fig. 3). Finally, we select stars in the G versus G − K_s colour–magnitude diagram, where giants still contaminating our sample can easily be separated from OBA stars (see Fig. 4). In particular we select stars satisfying:

(5)

Fig. 4. G − Ks vs. G colour–magnitude diagram of the sources selected after applying the criteria in the (J − Ks)0 vs. G − Ks and G − Ks vs. J − H colour–magnitude diagrams (see text). The solid orange line shows the equation G = 2(G − Ks)+3, and is designed to exclude the few cool giants that may still contaminate the sample.

(solid orange line in Fig. 4). The catalogue consists of 988 202 entries.

2.2. Characteristics of the target sample

Figure 5 shows the sky distribution of the sample selected in the previous sections in Galactic coordinates (l, b).

Fig. 5. On-sky density distribution of the candidate OBA stars selected in Sect. 2, in Galactic coordinates. Most of the sources are located in the Galactic plane (b = 0°) and in the Large and Small Magellanic clouds.

The stars in our catalogue are distributed in our Galaxy (mainly on the Galactic plane) and in the Large (LMC, l, b ∼ 280° , − 34°) and the Small (SMC, l, b ∼ 303° , − 45°) Magellanic Clouds. We stress that massive stars in the Magellanic clouds are included in the target sample of the SDSS-V survey. However, in the present paper we focus on the Milky Way disc, thus in the following sections we restrict our sample to stars with |b|< 25°. Massive star-forming regions can be recognised as over-densities in the source distribution.

Dust features are visible as ‘gaps’ in the star distribution. For example, the Aquila rift, the Pipe Nebula, and the Ophiucus clouds can be easily identified towards the Galactic centre (l, b = 0, 0 deg). By visually inspecting the features of Fig. 5, it is possible to identify old stars contaminating the sample. For instance, at l, b ∼ 309° ,15° the ω Centaurus globular cluster is visible, while the stars towards the Galactic centre delineate the shape of the bulge. We further study the purity and completeness of our sample in Sect. 2.3.

Figure 6 shows the distribution in apparent (left) and absolute magnitude (right) of the sources. The orange histograms show the distribution of all the sources in our sample. The grey histograms show the distribution of the sources selected in Sect. 3. In the right hand panel, absolute magnitudes are computed using the inversion of parallaxes as the distance (orange histogram) or astro-kinematic distances (grey histogram; see Sect. 3). The grey distribution appears shifted towards larger M_G values. This is due to the fact that, for stars with small parallaxes (and large parallax errors), our astro-kinematic distances are shorter than those obtained by inverting the parallax.

Fig. 6. Left: G magnitude distribution of all the OBA stars selected in Sect. 2.1 (orange histogram) and of the sources selected in Sect. 3 (grey histogram), with astro-kinematic distances comparable with those estimated by Bailer-Jones et al. (2018). Right: absolute magnitude distribution of all the OBA stars selected in Sect. 2.1, computed using their parallax (orange histogram), and of the sources selected in Sect. 3, computed using their astro-kinematic distances (grey histogram).

Figure 7 (left) shows the variability distribution of the sources in the G band vs. the ratio of the variability in the G_BP and G_RP bands. Following roughly Belokurov et al. (2017), Deason et al. (2017), and Iorio et al. (2018), the variability in a given filter x is defined as:

(6)

Fig. 7. G-band variability vs. GRP/GBP variability of all the sources selected in Sect. 2.1 (left) and of the sources selected in Sect. 4 (right). Eclipsing binaries have the same level variability in GBP and GRP and G−variability lower than 0.3. RR Lyrae stars have GBP/GRP variability ∼1.6 and G-band variability lower than 0.4.

where phot_x_n_obs is the number of observations contributing to the x-band photometry, and

(7)

While most of our sources are not variable, some cluster in two regions: RR Lyrae stars have G-band variability ranging from 5 to 40% (Belokurov et al. 2017) and G_BP/G_RP variability around 1.6; eclipsing binaries have G variability up to 30% and roughly equal G_BP and G_RP variability. Figure 7 (right) shows the same as Fig. 7 (left) but for the sources selected in Sect. 3 (see below), that have astro-kinematic distances compatible with those estimated by Bailer-Jones et al. (2018). This selection removes RR Lyrae stars, while retaining eclipsing binaries: this is expected as RR Lyrae stars do not follow the same kinematics as the OBA stars that we are interested in.

2.3. Completeness and purity

To study the purity level of our sample we cross-matched it with LAMOST DR6², and find 36 617 sources in common. We derived effective temperatures (T_eff) and surface gravity (log g) from the LAMOST spectra by fitting ab initio model spectra generated with the 1D-LTE ATLAS12 model atmosphere (Xiang et al., in prep.). Figure 8 shows the distribution of sources in the log T_eff vs. M_Ks plane. Around 45% of the sources are hotter than 9700 K (orange solid line in Fig. 8, left, corresponding to the temperature of a A0V-type star) and around 20% are hotter than 14 000 K (roughly the temperature of a B7V-type star). We assume that we can extrapolate the same numbers to the entire sample. The PARSEC isochrones (A_V = 0 mag and solar metallicity) in Fig. 8 (Bressan et al. 2012; Tang et al. 2014) show that our main sources of contamination are evolved massive or intermediate mass stars (e.g. yellow and blue super-giants) and that the sample is substantially free from RGB and AGB stars, which are colder and more luminous than the stars in our sample. In the left-hand panel, the isochrones are colour-coded by their log age (from 1 Myr to 500 Myr), in the right-hand panel they are colour-coded by mass.

Fig. 8. Left: distribution of the sources with LAMOST DR6 spectra in the log Teff vs. MKs (grey). The orange solid line indicates Teff = 9700 K, which corresponds to the temperature of a A0V-type star (Pecaut & Mamajek 2013). The isochrones are colour-coded according to their age (in Myr). Right: same as left, with the isochrones colour-coded according to their mass.

To estimate the completeness of the sample we cross-matched with the O and B-type star catalogue by Liu et al. (2019) and the O-type star catalogue by Sota et al. (2014) and Maíz Apellániz et al. (2016). The Liu et al. (2019) sample consists of 16 032 stars that reduce to 15 206 after cross-matching with Gaia EDR3 and 2MASS and further checking for duplicated sources (1″ cross-match radius), and to 9083 stars by applying Eq. (1). The cross-match between the Liu et al. (2019) catalogue and ours gives 8212 sources. The stars that are not included in the cross-match are those whose infrared colours are not consistent with our selection. The sample of Sota et al. (2014) and Maíz Apellániz et al. (2016, the Galactic O-star spectroscopic survey, GOSSS) consists of 590 stars that reduce to 580 after cross-matching with Gaia EDR3 and 2MASS. Our catalogue contains 503 (∼86%) of the GOSSS stars. As with the LAMOST sample, the stars missing from our selection do not follow some of our photometric selection criteria. For example, many of the missing GOSSS sources (49%, 54 stars) have J − H > 0.15(G − K_s)+0.05 (see Eq. (4)), that is 2MASS-Gaia EDR3 colours consistent with being giants.

3. Distance estimates

To study the 3D space distribution of our sources, precise distance estimations are needed. We estimate distances using a model designed to reproduce the properties of our data-set in terms of spatial and luminosity distribution, and the additional information that stars belonging to our sample should follow Galactic rotation, with a small, typical velocity dispersion. By making such an assumption we neglect non-circular streaming motions. We can then predict the true proper motions () of the stars in our sample and compare them with the observed ones, thus adding a further constraint to the distance estimation and deriving ‘astro-kinematic’ distances. The probability density function (pdf) for a star in our sample to be at a certain distance d_kin is given by:

(8)

where o = (ϖ−ϖ₀,μ_l*,μ_b) is the array of the astrometric observables, i.e. parallax ϖ (with ϖ₀ the parallax zero-point) and proper motions components in l and b, μ_l*, and μ_b; m_Ks is the apparent magnitude in the K_s band; Θ_KM represents our kinematic model; Θ_SM represents our model for the distribution of stars in the Galaxy; and Θ_CMD accounts for the observational effects due to our selection function on the spatial distribution of our sample.

The method and all the terms in Eq. (8) are described in detail in Appendix B. We adopt the mode of the pdf of Eq. (8) as the astro-kinematic distance estimate d_kin, and we use the 16th and 84th percentiles to estimate distance errors. We also tried using the median of the pdf as a distance estimate, but this did not cause significant differences in the maps presented in Sect. 5. Figure 9 (left) shows the comparison between our astro-kinematic distances and the photo-geometric distances estimated by Bailer-Jones et al. (2020). More than 85% of the sources have astro-kinematic distances consistent within 1σ with the photo-geometric distances from Bailer-Jones et al. (2020). We comment on the sources with inconsistent distances in Sect. 4.2.

Fig. 9. Comparison between the photo-geometric distances estimated by Bailer-Jones et al. (2020) and our astro-kinematic distances (see Sect. 3), after removing spurious sources (left), and after removing sources with vertical velocity not consistent with our kinematic model and with MKs < 0 mag (right). The orange dashed line has equation y = x.

4. Filtered sample

In this section we define a subset of the target sample (Sect. 2), which we use to map the structure of the young Milky Way disc. Such a “filtered” sample is obtained by cleaning the target sample for sources with likely spurious astrometric solutions (Sect. 4.1) and by removing sources with kinematic properties not consistent with the model that we used to estimate astro-kinematic distances (Sect. 4.2).

4.1. Removal of spurious sources

To remove spurious sources, we use a method developed by Rybizki et al. (2021) based on Gaia Collaboration (2021b). Spurious sources have poor astrometric solutions that can be due to the inconsistent matching of the observations to different physical sources. This is more likely to occur in regions of high source surface density (e.g. in the Galactic plane) or for close binary systems (either real or due to perspective effects). To identify poor astrometric solutions in their solar neighbourhood catalogue, Gaia Collaboration (2021b) constructed a random-forest classifier that assigns a probability to each source of having a good (or bad) astrometric solution based on astrometric quantities and quality indicators. The classification probability is approximately zero for sources with poor astrometric solutions, and approximately one for sources with good astrometric solutions. Gaia Collaboration (2021b) applied this classifier to stars with distances less than 100 pc. Rybizki et al. (2021) used similar methods to develop a neural network classifier, which they applied to the entire Gaia EDR3 catalogue. The distribution (in logarithmic scale) of ‘astrometric fidelities’ for the sources in our target sample is shown in Fig. 10. We select sources with astrometric fidelities > 0.5 (around 75% of the target sample) following the definition of ‘good’ and ‘bad’ astrometric solutions in Rybizki et al. (2021). The distribution in the Galactic plane of sources with good and bad astrometric solutions is shown in Fig. 8 of Rybizki et al. (2021).

Fig. 10. Astrometric fidelity distribution for the sources in our target sample. Low astrometric fidelity values (< 0.5) correspond to bad astrometric solutions, and high astrometric fidelity values (> 0.5) correspond to good astrometric solutions.

While this procedure allows us to obtain a more reliable sample, it might also eliminate binaries. Binaries however do not constitute the focus of this paper and do not affect the space distribution of our filtered sample.

4.2. Kinematic cleaning

Our astro-kinematic distance estimates are not accurate for sources that do not follow the disc kinematics assumed by our model. Such sources can be grouped into two categories:

1. Stars whose kinematics is influenced by the bar. These are located mostly towards the Galactic centre, and have astro-kinematic distances larger than Bailer-Jones photo-geometric distances. These stars could in principle be young. Their astro-kinematic distance estimates could however be wrong due to the kinematic effects of the bar itself, which have not been included in our model (see Sect. 3 and Appendix B). In the following sections we focus on the in-plane distribution of stars within d_kin ≲ 5 kpc, and thus we do not consider sources in the bar/bulge regions.

2. Stars with heated vertical kinematics. These stars have Bailer-Jones photo-geometric distances d_BJ ≳ 3 kpc and astro-kinematic distances d_kin ∼ 2 kpc. They could correspond to old contaminants, such as the RR Lyrae-type variable stars discussed in Sect. 2.2. To remove them we proceed as follows. First, we compute their vertical velocity, v_z (see Appendix C). We fit the vertical velocity distribution with a Gaussian Mixture model with two components, the first (component A) with mean vertical velocity ⟨v_z⟩_A = −0.5 km s⁻¹ and velocity dispersion σ_vz, A = 7.0 km s⁻¹ which has properties comparable with the young population that we are interested in studying, and the second (component B), with ⟨v_z⟩_B = 1.6 km s⁻¹ and velocity dispersion σ_vz, B = 29. km s⁻¹, which instead better reproduces the properties of an old population. We estimate the probability of each star belonging to either component A or B, and we select those stars with a greater probability of belonging to component A.

Figure 9 (right) shows the comparison between the astro-kinematic distances and Bailer-Jones photo-geometric distances after cleaning the target sample as explained above. At this point, around 95% of the sources have astro-kinematic distances consistent with Bailer-Jones photo-geometric distances within 1σ.

Finally, using our astro-kinematic distances, we estimate the absolute magnitude of each star in the K_s band, M_Ks, and we select stars with M_Ks < 0 mag. This condition refines Eq. (1), and is designed to select stars of type earlier than B7V. After applying the quality cuts described in this section, our filtered sample reduces to 435 273 stars.

5. Three-dimensional space distribution

To create the 3D density maps, we follow the method outlined in Zari et al. (2018). We compute galactic Cartesian coordinates, x, y, and z, for all the sources using our newly estimated astro-kinematic distances and we define a box V = 12 × 12 × 1.5 kpc centred on the Sun. We divide the cube into volume elements v = 5 × 5 × 5 pc. After computing the number of stars in each volume, we estimate the star density D(x, y, z) by smoothing the distribution by means of a 3D Gaussian filter, using a technique similar to that used by Bouy & Alves (2015). The Gaussian width (equal on the three axes) is w = 5 pc and the Gaussian is truncated at 3σ. The choice of a certain w value is arbitrary. A high w value produces a smooth, less detailed map, while a low w value results in a noisy map.

Figure 11 shows the projection on the Galactic plane of the filtered sample selected in Sect. 4. A number of features are visible, showing different structures at different scales and distances. At the centre of the map (x, y ∼ 0, 0), a low-density gap is visible. This is due to a combination of reasons. First, in the solar neighbourhood the number of OBA stars is low. The Sun is indeed currently located in a cavity, partially filled with hot, low-density gas, which is usually called the Local Bubble (LB). The LB was likely created 10–20 Myr ago by a number of supernova explosions that likely halted any low- or high-mass star formation episode in the region. Second, our selection criteria exclude by construction stars of later spectral type than ∼B7V. Thus, even if some late-type B stars are present within 200−300 pc of the Sun, they are excluded from our selection. Third, there are selection effects due to the fact that (a) nearby bright sources are not included in Gaia EDR3 and that (b) some of those that are included might have poor 2MASS photometry; these latter are excluded by the colour selection in Eq. (2).

Fig. 11. Surface density of the stars selected in Sect. 3 projected on the Galactic plane. The Sun is in (0,0), the x-axis is directed towards the Galactic centre, and the y-axis towards Galactic rotation. The z-axis is perpendicular to the plane. The density is displayed in arbitrary units. The dashed circles have radii from 1 to 6 kpc, in steps of 1 kpc. The labels follow the nomenclature proposed by K. Jardine (http://gruze.org/galaxymap/map_2020/). The gaps in the density distribution mentioned in Sect. 5 are indicated by grey circles.

The small dense clumps can be associated to well-studied OB associations, for example Cygnus, Carina, Cassiopeia, and Vela are visible. The distribution of the lower density contours traces the Milky Way young disc structure, in particular the spiral arm location (see Sects. 5.2 and 6). The density distribution presents numerous low-density gaps. These could be due to the fact that our view might be obscured by interstellar dust towards certain lines of sight, however some gaps are located in regions of relatively low extinction (see also Sect. 5.3). For example, the Perseus gap (x, y ∼ −2, 1 kpc) has been detected in the distribution of other young tracers, which are traditionally associated with spiral arms, such as the HII maps (Hou & Han 2014) and the Cepheid distribution (Skowron et al. 2019, and see Sect. 5.4), and in the OB stars shown by Romero-Gómez et al. (2019) and Poggio et al. (2018, and see Sect. 5.2), and the high-mass star-forming regions of Reid et al. (2014) and Reid et al. (2019, and see Sect. 5.5). By tracing the open cluster distribution in the Galactic disc, Cantat-Gaudin et al. (2020) suggest that a similar gap might be observed around x, y ∼ 1, −1 kpc, where our map shows an under-density, although less pronounced than the Perseus gap. In general, the distribution of stars along the ‘arms’ is not homogeneous, but characterised by under- and over-densities. For distances larger than 3–4 kpc, the number of sources (and thus their spatial density) decreases due to our magnitude limit at G = 16 mag (see Sect. 2).

The radial features are caused by ‘shadow cones’ produced by foreground extinction, and by distance uncertainties. As shown in Fig. B.3, distance uncertainties are lower than 10% for stars within 5 kpc, and therefore are only partially responsible for the elongation.

5.1. The Galactic warp

As mentioned in Sect. 1, early-type star samples have been used to study the properties of the Galactic warp (see Poggio et al. 2018; Romero-Gómez et al. 2019, and references therein). Although in this paper we do not focus specifically on the warp, we studied the median height of our filtered sample with respect to the plane of the Galaxy and the median vertical velocity distribution, and compared our results to those by Poggio et al. (2018) and Romero-Gómez et al. (2019) to further validate our selection.

Figure 12 shows a map of the median height z of our sample. The radial features are likely an artefact due to uneven sampling above or below the Galactic plane due to foreground extinction. Such features make the interpretation of Fig. 12 uncertain and may prevent us from seeing the warp shape, especially for X > −2.5 kpc. For X < −2.5 kpc, the impact of extinction is less dramatic (see also Fig. 14), and therefore the displacement visible above and below the plane might indeed trace the warp (see also Fig. 5 in Romero-Gómez et al. 2019).

Fig. 12. Median height for the filtered sample projected on the Galactic plane. We divided the XY plane into bins of 200 pc width, and we only show the ones containing more than ten stars. The dashed circles have radii from 1 to 9 kpc, in steps of 1 kpc.

The kinematic signature of the warp is evident in Fig. 13 (see e.g. Poggio et al. 2018; Romero-Gómez et al. 2019), which shows the median vertical velocity (v_z) distribution of our filtered sample on the Galactic plane (see Appendix C for more details on how the vertical velocity was computed). For distances larger than 4 kpc from the Sun, the vertical velocities have positive values that seem to peak at Y ∼ 0 kpc and decrease towards both sides. At smaller distances from the Sun, some features also present positive vertical velocity. While some correspond to density enhancements in the 3D spatial distribution of our sources, this is not the case for all of them. Such variation in v_z is likely related to the perturbations experienced by the disc (e.g. interaction with satellite galaxies).

Fig. 13. Median vertical velocity vz distribution on the Galactic plane. We divided the XY plane into bins of 200 pc width, and we only show the ones containing more than ten stars. The dashed circles have radii from 1 to 9 kpc, in steps of 1 kpc. The Galactic warp is visible for R > 4 kpc towards the Galactic anti-centre.

5.2. Comparison with other O and B star maps

By combining Gaia EDR3 and VPHAS data with literature catalogues, Chen et al. (2019) obtained a sample of 14 880 OB stars, earlier than B3V, which they used to describe the morphology of the spiral structure of the Milky Way. The method used by these latter authors to select O and B stars substantially differs from ours, but a qualitative comparison between their Fig. 4 and our maps shows substantially the same gross structures. This increases our confidence that both maps are revealing the same features in the Galactic O and B star distribution. The larger number of sources in our catalogue enables us to reveal more details in the source distribution. Interestingly, neither map shows a clear spiral structure. This is further discussed in the Sect. 6.

Figure 3 of Romero-Gómez et al. (2019) shows the distribution on the Galactic plane of their upper main sequence sources. Figure 11 of the present paper shows essentially the same over-densities, albeit in greater detail. This greater detail is mainly attributable to two factors: (a) to create their maps, Romero-Gómez et al. (2019) used a larger bin size than in the present work, and (b) their sample likely includes more later type stars than ours because of their selection criteria: this causes the distribution of sources to be smoother.

Figure 3 (panel A) in Poggio et al. (2018) also shows the distribution of upper main sequence sources on the Galactic plane. Similarly to Romero-Gómez et al. (2019), the pixel size of Poggio et al. (2018) is higher than ours. As mentioned above, this creates a smoother map, where it is not possible to identify small-scale over-densities. However, the comparison between large-scale structures shows many similarities.

5.3. Comparison with dust distribution

Dust is also a tracer of Galactic structure. Figure 14 shows the projection of the 3D dust distribution from Lallement et al. (2019) on the Galactic plane, over-plotted on top of our density map (see Fig. 11). The units of the density distribution are arbitrary. The dust density distribution shows discrepancies with respect to the star distribution. For example, the two elongated structures at the centre of the map are not prominent in the OB star distribution, and the dust features in the first and fourth Galactic quadrants are off-set with respect to the star density distribution, and seem to have a different separation and relative inclination. Such offsets are expected, as newly formed stars will drift from their birthplaces as they follow galactic rotation while the spiral arms move with a given pattern speed.

Fig. 14. Same as Fig. 11, with dust density contours from Lallement et al. (2019) over-plotted. Both densities are displayed in arbitrary units. The dashed circles have radii from 1 to 5 kpc, in steps of 1 kpc.

5.4. Comparison with Cepheid distribution

Classical Cepheids are young (< 400 Myr) variable stars, whose distances can be estimated thanks to their period–luminosity relation. Figure 15 shows the Cepheids identified by Skowron et al. (2019) over-plotted on the density distribution of O and B stars (same as in Figs. 11 and 14). The colour-bar represents the ages (in Myr) determined by Skowron et al. (2019). The Cepheid distribution traces the density enhancements corresponding to the Sagittarius-Carina arm reasonably well (see Sect. 5.5), while the correspondence with the other density enhancements is not as tight. Although the age distribution of the Cepheids and our filtered sample is similar, the selection criteria are different: this makes a more direct comparison difficult.

Fig. 15. Cepheids from Skowron et al. (2019), colour-coded by age (Myr). The density map is the same as Fig. 11. The dashed circles have radii from 1 to 6 kpc, in steps of 1 kpc.

5.5. Comparison with maser distribution

The spiral structure of the Milky Way can be traced also by water and methanol masers associated with high-mass-star-forming regions (HMSFRs). By measuring parallaxes and proper motions of such HMSFRs, Reid et al. (2019) provided estimates of the fundamental Galactic parameters, among which the pitch angle and the arm width. Figure 16 shows the masers and the location of the arms from the fit from Reid et al. (2019) on top of the OB star density map shown in Fig. 11. Similarly to the Cepheids distribution, there is good agreement between the position of the masers and the location of the most prominent over-densities, which makes us confident of our distance estimates. Unfortunately there are no available maser data tracing the spiral structure in the fourth Galactic quadrant (X > 0, Y < 0), making a more detailed comparison in that region non-trivial. However, by tracing the spiral arms beyond the data while keeping the parameters fixed (dashed lines, and lighter contours), more discrepancies become visible, especially towards the inner Galaxy, towards the Sagittarius-Carina (purple) and Scutum-Centaurus arm (orange).

Fig. 16. Masers (dots) and the fit to the spiral arms (solid lines) from Reid et al. (2019) are over-plotted on top of our density map (same as Fig. 11). The maser distances were computed by naively inverting the parallax. The dashed circles have radii from 1 to 6 kpc, in steps of 1 kpc. The shaded regions correspond to the arm width. Different colours correspond to different spiral arms: Outer, green; Perseus, grey; Local, pink; Sagittarius-Carina, purple; Scutum-Centaurus, orange.

6. Discussion

We devised a large and systematically selected sample of massive, young stars and studied its properties focusing on the spatial distribution of our sources on the Galactic plane. In this section we discuss our findings in the context of the spiral arm structure of the Milky Way.

Traditionally, the four main spiral arms of the Milky Way are considered to be Perseus, Sagittarius, Scutum, and Norma, with an additional Outer arm, which might be the outer part of one of the inner arms, and the Local or Orion arm, which is shorter and may be a spur or a bridge between two arms (see for instance Churchwell et al. 2009). In the Milky Way, the term ‘spiral arm’ has been used very broadly to refer to regions encompassing HI, molecular gas (and masers), dust, and young stars, while it has been shown from external galaxies that these tracers have quite distinct morphology. This has led to some inconsistencies in defining the spiral arms in our Galaxy. For example, the structure parameters derived by Reid et al. (2019) (see Fig. 16 and Sect. 5.5 of the present paper) do not agree with those presented in Chen et al. (2019) (see Sect. 5.2). Further, certain features might be interpreted as spiral arms or spurs connecting different arms: this is for example the case for the dust complex labelled ‘Lower Sagittarius-Carina’ by Lallement et al. (2019) (see Sects. 5.3 and Fig. 14) and ‘Lower Sagittarius-Carina Spur’ by Chen et al. (2020). Further, by studying the correlation between the location of young clusters and molecular clouds in NGC 7793 and M51 (respectively a flocculent and a grand-design spiral galaxy), Grasha et al. (2018, 2019) found that the star clusters that are associated (i.e. located within the footprint of a giant molecular cloud) are young, with a median age of 2 Myr in NGC 7793 and of 4 Myr in M51. Older clusters are mostly unassociated with any molecular cloud. Therefore, equating the same spiral arm morphology with different tracers (such as early type stars and dust) might be inappropriate.

The over-density towards the inner Galaxy (visible for positive x values, i.e. in the first and fourth Galactic quadrants) in Fig. 11 is associated with the Sagittarius-Carina and Scutum-Centaurus arm, and is much more prominent than the others, containing numerous HMSFRs. On the contrary, the distribution of massive stars associated with the Perseus arm peaks mainly in the second Galactic quadrant (towards Cassiopeia). This is consistent with the findings of Reid et al. (2019), and would point to the conclusion that the Perseus arm as traced by O- and B-type stars is not a dominant arm, and might be dispersing in the field. This suggests that a recent, large-scale episode of (massive) star formation occurred in Sagittarius-Carina and Scutum-Centaurus, while in the other arms (Perseus and Local) stars formed earlier, except in a few isolated massive associations, such as Cygnus and Cassiopeia. To confirm this, we selected stars brighter than M_Ks = −1 mag and M_Ks = −2 mag (corresponding respectively to the K_s absolute magnitude of B2V-type and B1V-type stars), and we evaluated their density on the plane following the procedure described in Sect. 5. The density maps obtained with these samples are shown in Fig. 17, where also the same map of Fig. 11 is shown for clarity. The density contours in Fig. 17 (left) trace the same dense structures as in Fig. 11 (we note that the contour levels are different), while the low-density contours substantially disappear. This is even more evident in Fig. 17 (right). This is expected, as intrinsically brighter stars are on average younger and thus have had less time to disperse in the field.

Fig. 17. Left: same as Fig. 11. Centre: density of stars with MKs < −1 mag projected on the Galactic plane. Right: density of stars with MKs < −2 mag projected on the Galactic plane. We note that the density scales are different.

Drimmel (2000) and Drimmel & Spergel (2001) studied the spiral arm features of the Milky Way in the far-infrared and near-infrared and concluded that the near-infrared features were consistent with a two-armed spiral model, while the other two arms were traced only by the far-infrared emission and thus could perhaps be present only in gas and young stars. More recently, Xu et al. (2018) and Chen et al. (2019) noted that the spiral arm structure traced by O- and early B-type stars may have many substructures in addition to the four major arms. These latter authors concluded that the Milky Way might not be a pure grand-design spiral galaxy with two or four well-defined, dominant arms. This might suggest that the Milky Way could exhibit characteristics of both flocculent and grand-design arms, with grand-design features seen in the infrared (old stellar populations) and more flocculent (or multi- arm) features seen in the optical (young stars). As mentioned above, these characteristics have been observed in spiral galaxies for instance by Kendall et al. (2011, 2015, and references therein), who concluded that galaxies that exhibit a grand-design structure in the optical also exhibit such structure in the near infra-red (NIR), while there are optically flocculent galaxies which also exhibit NIR grand-design structure. As pointed out by Dobbs & Baba (2014), Elmegreen et al. (1999, 2011) noted nevertheless that most flocculent galaxies do not exhibit grand-design structure and those that do have very weak spiral arms. Elmegreen et al. (2003) further suggested that both grand-design and flocculent spirals (as seen in the old stars) exhibit a similar structure in the gas and young stars (independent of the underlying old stellar population) which is driven by turbulence in the disc. The co-existence of different spiral-arm patterns in the Milky Way disc might have different explanations. For example, Martos et al. (2004) concluded that secondary arms could be resonance-related features in the quasi-stationary density wave picture, while Steiman-Cameron et al. (2010) and Drimmel (2000) proposed that the traditional four spiral arms may represent the dynamical response of a gas disc to a two-arm spiral perturbation in the mass distribution.

As mentioned above, the density distribution presented in Fig. 11 exhibits a prominent density enhancement roughly corresponding to the Sagittarius-Carina and Scutum-Centaurus arms. This alone does not allow us to draw firm conclusions on the nature (or on the number) of the spiral arms of the Milky Way. On the one hand, star formation occurs in a clumpy and patchy fashion: this might explain the observed density distribution without necessarily assuming a flocculent structure. On the other hand, the map shown in Fig. 15 might suggest a different picture, in which the Milky Way, as traced by young stars (≲300 Myr), does not show a discrete set of spiral arms.

The spiral arm structure of the Milky Way can also be investigated using stellar kinematics (see e.g. Antoja et al. 2016). A spiral mass density perturbation can have a very different morphology from spiral-like density enhancements in molecular gas or young stars. Eilers et al. (2020) studied the kinematics of red giant stars and showed a spiral feature in their radial velocities that they interpret as arising from a perturbation to the potential of the Milky Way caused by a two-armed logarithmic spiral. Their model predicts the locations of the spiral perturbation in the Milky Way, with one roughly co-spatial to the Orion (or Local) arm, and the other roughly corresponding to the Outer arm. The location of the perturbation predicted by Eilers et al. (2020) does not closely follow the over-densities of the map in Fig. 11. The comparison between the radial velocities of the O- and B-type star sample and the giant sample will be crucial to make progress in our understanding of the mechanism giving rise to the spiral signature presented by Eilers et al. (2020). The data of the SDSS-V survey will make such a comparison possible, and will also likely allow the construction of a surface density map of the giant sample which could be directly compared with the map presented in this work. The complete kinematic information that we will obtain by combining Gaia proper motions with SDSS-V radial velocities will also allow us to better separate stars belonging to different spiral arms, and to study their internal motions. Indeed, as already mentioned in Sect. 5, vertical velocities are not related to the perturbation in the disc kinematics induced by the spiral arms.

7. Conclusions

In this study we analyse the 3D space distribution of a sample of hot and luminous OBA stars, focusing in particular on their configuration in the Galactic disc. Our target selection is based on the combination of Gaia EDR3 astrometry and photometry and 2MASS photometry, and is designed to provide a well-defined selection function. We describe the properties of the target sample in terms of distribution in the sky, brightness, and variability. We estimate the purity and contamination of the catalogue by comparing with existing O- and B-type star catalogues, such as those by Liu et al. (2019), and Sota et al. (2014) and Maíz Apellániz et al. (2016), and by deriving stellar parameters (T_eff and log g) for stars in common with LAMOST DR6.

By assuming that young massive stars are on near-circular orbits with a small velocity dispersion, we compute astro-kinematic distances, and compare them to those derived by Bailer-Jones et al. (2020). We use these distances to study the distribution of sources in the Galactic disc of a subset of the target sample, which we obtain by filtering out sources with spurious astrometric solutions and kinematic properties inconsistent with our model.

We find that the distribution of sources in the plane of the Milky Way is highly structured, characterised by over- and under-densities. Some of the density enhancements correspond to massive star forming regions, such as Carina, Cygnus, and Cassiopeia. With these associations, the inner arms (Sagittarius-Carina and Scutum-Centaurus) are strikingly more prominent in OBA stars than any of the outer arms, which we find little evidence for.

The distribution of O- and B-type stars from previous catalogues, classical Cepheids, dust, and HMSFRs shows similarities with our density distribution. However, the picture of the spiral arm structure of the Galaxy that we obtain in this study is complex and may suggest that young stars show little tendency to be neatly organised in distinct spiral features.

To assess different spiral arm models, it would be necessary to extend our map beyond the volume currently studied and, perhaps more importantly, to combine it with better Gaia data and spectroscopic information.

The target sample that we have devised is indeed optimised for spectroscopic follow-up with the SDSS-V survey. The information that we will be able to obtain with SDSS-V will be crucial for studying the kinematics and dynamics of the spiral arms, and thus the nature of the spiral arms themselves. Finally, this sample will further enable the study of the properties of massive stars, for example in terms of multiplicity, internal structure, and (binary) evolution, in a way that is statistically robust.

Acknowledgments

We thank the referee for their comments, which improved the quality of this manuscript. We would also like to thank: J. Rybizki and G. Green for making the catalogue of spurious sources available in advance of publication; A. Gould for discussions that prompted us to use kinematics for distance estimations; M. Sormani for discussions on the kinematics of the MW bar and bulge. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. Funding for the Sloan Digital Sky Survey V has been provided by the Alfred P. Sloan Foundation, the Heising-Simons Foundation, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS web site is www.sdss.org. SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration, including the Carnegie Institution for Science, the Chilean Participation Group, the Gotham Participation Group, Harvard University, The Johns Hopkins University, L’Ecole polytechnique fédérale de Lausanne (EPFL), Leibniz-Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Extraterrestrische Physik (MPE), Nanjing University, National Astronomical Observatories of China (NAOC), New Mexico State University, The Ohio State University, Pennsylvania State University, Space Telescope Science Institute (STScI), the Stellar Astrophysics Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Illinois at Urbana-Champaign, University of Toronto, University of Utah, University of Virginia, and Yale University. The research leading to these results has (partially) received funding from the KU Leuven Research Council (grant C16/18/005: PARADISE), from the Research Foundation Flanders (FWO) under grant agreement G0H5416N (ERC Runner Up Project), as well as from the BELgian federal Science Policy Office (BELSPO) through PRODEX grant PLATO. E.P. acknowledges support from the Centre national d’études spatiales (CNES). This research made use of TOPCAT (Taylor 2005), Astropy, (Astropy Collaboration 2013, 2018), matplotlib (Hunter 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), and scikit-learn (Pedregosa et al. 2011). This work would have not been possible without the countless hours put in by members of the open-source community all around the world.

References

Appendix A: Query of the Gaia archive

Here we provide an example query for cross-matching 2MASS and Gaia EDR3 in the Gaia archive using the cross-match with Gaia DR2.

        select edr3.* xdr2.*, tm.*
from gaiaedr3.gaia_source as edr3
inner join gaiaedr3.dr2_neighbourhood as xdr2
on edr3.source_id = xdr2.dr3_source_id
inner join gaiadr2.tmass_best_neighbour as xtm
on xdr2.dr2_source_id = xtm.source_id
inner join gaiadr1.tmass_original_valid AS tm
on tm.tmass_oid = xtm.tmass_oid
WHERE xtm.angular_distance < 1.
AND xdr2.angular_distance < 100.
AND edr3.phot_g_mean_mag < 16.
AND edr3.parallax < power(10,(10-tm.ks_m)/5)

Appendix B: Astro-kinematic distances

In this section we describe all the terms of Eq. (8), which we report here for convenience:

(B.1)

The prior p(d, m_Ks | l, b, Θ_SM, Θ_CMD) can be written as:

(B.2)

The term p(d|l, b, Θ_SM) represents the probability density function of observing a star in the direction (l, b) at distance d from the Sun according to our assumed model for the spatial distribution of stars in the Galaxy. The term f(d, m_Ks|Θ_CMD) specifies the fraction of stars that can be observed at a distance d and magnitude m_Ks given the distribution of stars in colour–magnitude space and our selection criteria. In the following sections, we describe the different components of Eq. (B.1).

Figure B.1 shows example pdfs for four randomly selected stars. Both the parallaxes of the stars in the top row have σ_ϖ/ϖ < 20%. In the left panel, the distance estimate is mostly constrained by the term P(d | μ_l * ,Θ_KM), which represent the probability of observing a star at distance d given the observed proper motion along Galactic latitude and the kinematic model Θ_KM. P(d | μ_l * ,Θ_KM) shows bi-modality, with a primary maximum at ∼6 kpc, and a secondary maximum at ∼1.8 kpc. The total distance pdf shows traces of such bi-modality. In the right panel of the top row on the other hand, the parallax component dominates the distance determination. The star in the left panel of the bottom row has σ_ϖ/ϖ < 20%, however both P(d | μ_l * ,Θ_KM and P(d | μ_b, Θ_KM peak at closer distances than the parallax term, and thus put a strong constraint on the final distance estimate. Finally, the star in the right panel of the bottom row has σ_ϖ/ϖ ≈ 70%. The final distance pdf is nevertheless still quite narrow, as a result of the combination of the different terms.

Fig. B.1. Posterior pdfs for four stars in our sample (their Gaia EDR3 source_id is displayed at top of each panel). The thick black solid line represents p(d | o, mKsΘKM, ΘSM, ΘCMD). The coloured thin lines illustrate the different components of the pdfs. The correlation terms in the covariance matrix of Eq. (B.12) are neglected here for illustrative purposes. The distributions are scaled so that they peak at unity.

B.1. Kinematic model

We assume that stars in our sample follow the rotation curve determined by Eilers et al. (2019),

(B.3)

with R the Galactocentric radius and R_⊙ = 8.122 kpc (Gravity Collaboration 2018). The intrinsic velocity dispersions at the Sun’s location in the components of the velocity in Galactic coordinates U, V, W are σ_U = 16.7 km s⁻¹, σ_V = 10.8 km s⁻¹, and σ_W = 6 km s⁻¹ (Robin et al. 2003). These values are assumed to be constant over the disc.

In Cartesian Galactic coordinates, the pdf of the space velocity for one star in our sample is:

(B.4)

the prime signifying the transpose of the vector. The quantity u takes into account Galactic rotation and solar motion, and can be written as:

(B.5)

where

(B.6)

is the circular velocity in Galactocentric Cartesian coordinates of a star at a Galactocentric distance r_star and x_star = [x_star, y_star, z_star]. The unit vectors and e_Z are defined as and e_Z = [0, 0, 1]. At the Sun’s position x_⊙,

(B.7)

where e_Y = [0, 1, 0] is the unit vector in the Y direction in Galactic coordinates. We assumed the Sun’s position with respect to the Galactic centre to be x_⊙ = [ − 8.122, 0.,0.025] kpc and the Sun’s peculiar motion with respect to the Local Standard of Rest to be v_⊙ = [11.1, 12.24, 7.25] km s⁻¹ (Schönrich et al. 2010). We write the velocity dispersion matrix S as:

(B.8)

For each star, the astrometric observables are the parallax ϖ, and the proper motions components in l and b, μ_l* = μ_l cos b, and μ_b respectively. These are collected in arrays:

(B.9)

where ϖ₀ is the systematic zero point offset of Gaia EDR3 parallaxes. There is not a unique ϖ₀ for all the Gaia EDR3 stars. Following Bailer-Jones et al. (2020), we applied the parallax zero-point correction derived by Lindegren et al. (2021b).

By assuming Gaussian errors, the pdf for the observables o given the true values , is then:

(B.10)

The true values are defined as:

(B.11)

where d is the true distance to the star and A = 4.74047 km yr s⁻¹. The two vectors p and q are the two components of the normal triad in longitude and latitude, and, as above, the prime signify their transpose. The elements of the astrometric covariance matrix C are:

(B.12)

and are provided in the Gaia archive.

The joint pdf of the observables with the velocity is therefore:

(B.13)

The pdf for the observables is obtained by marginalising over the velocity:

(B.14)

The integral can be resolved analytically. Since the product of two normal pdf is normal and the marginal density of a normal pdf is also normal, we can write:

(B.15)

where

(B.16)

and

(B.17)

B.2. Structural model

The probability of a star to be at a true distance d is proportional to the stellar density ρ predicted by our spatial model Θ_SM, so that we can write:

(B.18)

where the Jacobian term d² takes volume effects into account. For convenience we write our model in Galactocentric coordinates (R, ϕ, z), where the stellar density is modelled as an exponential disc:

(B.19)

Here, R and z are the Galactocentric radius and height above the plane. The quantity h_z = 0.15 kpc is the disc scale height (cf. Poggio et al. 2020), and L = 3.5 kpc is the best fit scale-length of the young (< 1 Gyr) disc, where we have used the model presented in Frankel et al. (2018) with the data and method in Frankel et al. (2019).

B.3. Luminosity function

We assume that the Milky Way has a universal luminosity function ϕ(M_Ks) = ϕ(M_Ks(m_Ks, d)), where M_Ks is the absolute magnitude in the 2MASS K_s band, and M_Ks = m_Ks − 5log₁₀(d)+5. Contrary to Poggio et al. (2020) or Astraatmadja & Bailer-Jones (2016), we use the near-infrared K_s band to be able to neglect extinction, at least as a first approximation. Following the procedure outlined by Poggio et al. (2020), we construct the distribution of stars in colour–magnitude space ϕ(m_Ks, d, c) by assuming a constant star formation rate, the two-part power law Kroupa initial mass function corrected for unresolved binaries, and solar metallicity. The term c indicates a colour (for instance, G − K_s). The distribution of stars in colour–magnitude space is then modified by setting to zero the parts of the distribution where T_eff < 8000 K. This roughly corresponds to applying the colour cuts aiming at selecting hot stars as described in Sect. 2.1, and effectively changes the term S(c) of Eq. (B.20) into a selection approximately based on effective temperature, S(T_eff). We can thus write the effective luminosity function, shown in Fig. B.2, f(d, m_Ks | Θ_CMD) as:

(B.20)

Fig. B.2. Effective luminosity function f(d, mKs|ΘCMD) used to determine our prior, normalised to unity.

where S(c) = 1 if T_eff > 8000 K and S(c) = 0 otherwise.

B.4. Distance uncertainties

Figure B.3 shows the fractional distance uncertainty for stars in the filtered sample (see Sect. 4). For distances smaller than ∼4 kpc, parallax uncertainties dominate the fractional error distribution. This is a regime where σ_ϖ/ϖ < 0.1, and the inverse parallax is a good distance estimate (see also Bailer-Jones et al. 2020). For distances larger than ∼4 kpc, the kinematic model and the other priors described above keep the distance uncertainties from exploding. The fractional errors are below the 10% for d_kin < 5 kpc, and remain well below the 20% for d_kin < 10 kpc.

Fig. B.3. Fractional distance uncertainty dkin, 84 − dkin, 16/2 dkin for stars in the clean sample. The dashed light grey lines correspond to the 16th and 84th percentiles, the solid grey line corresponds to the median fractional error.

B.5. Distance validation

To validate our distances, we identified the stars in our sample that belong to the three clusters IC 4756, NGC 2112, and King 21, respectively at 469 pc, 1090 pc, and 3005 pc, as determined by Cantat-Gaudin et al. (2020). Figure B.4 shows the distance pdfs of such stars. As expected, the pdfs for stars at larger distances are broader than for those at closer distances. The distance estimates for all the cluster members are compatible within 2σ with the average distance value. The dispersion around the median distance is between 10 and 20 pc for IC 4756, around 65 pc for NGC 2112, and around 300 pc for King 21. This indicates that distance uncertainties are within the 10% for d ∼ 3 kpc (where most of the over-densities are located), in agreement with Appendix B.4.

Fig. B.4. Single star distance pdfs for the members of the three clusters IC 4756 (black), NGC 2112 (blue), and King 21 (orange). The grey vertical lines represent the cluster distances (respectively 469 pc, 1090 pc, and 3005 pc) derived by Cantat-Gaudin et al. (2020).

Appendix C: Vertical velocity

As described in Drimmel et al. (2000), the vertical velocity v_z of a star is computed by:

(C.1)

where b is Galactic longitude, μ_b is the proper motion along b, d is the distance estimate, W_⊙ is the Sun’s velocity along Z, and v_r is the line-of-sight velocity. The majority of stars in our sample lack line-of-sight velocities, and therefore it is not possible to directly calculate their vertical velocity. However, for stars at low Galactic latitudes (|b|< 20°), the term v_rsin(b) can be approximately computed by assigning to each star the velocity it would have if it followed the Galactic rotation curve exactly:

(C.2)

where: S_⊙ = U_⊙ cos l + V_⊙sin(l), with l the galactic longitude of a star, and U_⊙ and V_⊙ the components of the Sun’s motion along X and Y respectively; S = v_ϕ R_⊙/R − v_LSRsin(l), with R_⊙ the distance of the Sun from the Galactic Centre, R a star’s Galactocentric radius, v_ϕ the azymutal velocity, and v_LSR the standard of rest velocity. We assumed (U, V, W)_⊙ = (11.1, 12.24, 7.25) km s⁻¹ from Schönrich et al. (2010), and the rotation curve derived by Eilers et al. (2019), reported in Eq. (B.3).

All Figures

	Fig. 1. Two-dimensional histogram of J − Ks vs. G − Ks (left), (J − Ks)0, and G − Ks (right) for the stars in the LAMOST sample by Liu et al. (2019); see main text for more details on the sample.
In the text

	Fig. 2. (J − Ks)0 vs. GBP − GRP colour–magnitude diagram. OBA stars roughly fall in the region highlighted by the orange ellipse. This representation shows that hot stars stand out qualitatively in colour–colour space.
In the text

	Fig. 3. J − H vs. G − Ks colour–colour diagram of the sources with (J − Ks)0 < 0.1 mag and (J − Ks)0 > −0.3 mag. The solid lines delineate the criteria for selecting OBA stars and are defined in the text.
In the text

	Fig. 4. G − Ks vs. G colour–magnitude diagram of the sources selected after applying the criteria in the (J − Ks)0 vs. G − Ks and G − Ks vs. J − H colour–magnitude diagrams (see text). The solid orange line shows the equation G = 2(G − Ks)+3, and is designed to exclude the few cool giants that may still contaminate the sample.
In the text

	Fig. 5. On-sky density distribution of the candidate OBA stars selected in Sect. 2, in Galactic coordinates. Most of the sources are located in the Galactic plane (b = 0°) and in the Large and Small Magellanic clouds.
In the text

Fig. 6. Left: G magnitude distribution of all the OBA stars selected in Sect. 2.1 (orange histogram) and of the sources selected in Sect. 3 (grey histogram), with astro-kinematic distances comparable with those estimated by Bailer-Jones et al. (2018). Right: absolute magnitude distribution of all the OBA stars selected in Sect. 2.1, computed using their parallax (orange histogram), and of the sources selected in Sect. 3, computed using their astro-kinematic distances (grey histogram).

In the text

	Fig. 7. G-band variability vs. GRP/GBP variability of all the sources selected in Sect. 2.1 (left) and of the sources selected in Sect. 4 (right). Eclipsing binaries have the same level variability in GBP and GRP and G−variability lower than 0.3. RR Lyrae stars have GBP/GRP variability ∼1.6 and G-band variability lower than 0.4.
In the text

	Fig. 8. Left: distribution of the sources with LAMOST DR6 spectra in the log Teff vs. MKs (grey). The orange solid line indicates Teff = 9700 K, which corresponds to the temperature of a A0V-type star (Pecaut & Mamajek 2013). The isochrones are colour-coded according to their age (in Myr). Right: same as left, with the isochrones colour-coded according to their mass.
In the text

	Fig. 9. Comparison between the photo-geometric distances estimated by Bailer-Jones et al. (2020) and our astro-kinematic distances (see Sect. 3), after removing spurious sources (left), and after removing sources with vertical velocity not consistent with our kinematic model and with MKs < 0 mag (right). The orange dashed line has equation y = x.
In the text

	Fig. 10. Astrometric fidelity distribution for the sources in our target sample. Low astrometric fidelity values (< 0.5) correspond to bad astrometric solutions, and high astrometric fidelity values (> 0.5) correspond to good astrometric solutions.
In the text

Fig. 11. Surface density of the stars selected in Sect. 3 projected on the Galactic plane. The Sun is in (0,0), the x-axis is directed towards the Galactic centre, and the y-axis towards Galactic rotation. The z-axis is perpendicular to the plane. The density is displayed in arbitrary units. The dashed circles have radii from 1 to 6 kpc, in steps of 1 kpc. The labels follow the nomenclature proposed by K. Jardine (http://gruze.org/galaxymap/map_2020/). The gaps in the density distribution mentioned in Sect. 5 are indicated by grey circles.

In the text

	Fig. 12. Median height for the filtered sample projected on the Galactic plane. We divided the XY plane into bins of 200 pc width, and we only show the ones containing more than ten stars. The dashed circles have radii from 1 to 9 kpc, in steps of 1 kpc.
In the text

	Fig. 13. Median vertical velocity vz distribution on the Galactic plane. We divided the XY plane into bins of 200 pc width, and we only show the ones containing more than ten stars. The dashed circles have radii from 1 to 9 kpc, in steps of 1 kpc. The Galactic warp is visible for R > 4 kpc towards the Galactic anti-centre.
In the text

	Fig. 14. Same as Fig. 11, with dust density contours from Lallement et al. (2019) over-plotted. Both densities are displayed in arbitrary units. The dashed circles have radii from 1 to 5 kpc, in steps of 1 kpc.
In the text

	Fig. 15. Cepheids from Skowron et al. (2019), colour-coded by age (Myr). The density map is the same as Fig. 11. The dashed circles have radii from 1 to 6 kpc, in steps of 1 kpc.
In the text

Fig. 16. Masers (dots) and the fit to the spiral arms (solid lines) from Reid et al. (2019) are over-plotted on top of our density map (same as Fig. 11). The maser distances were computed by naively inverting the parallax. The dashed circles have radii from 1 to 6 kpc, in steps of 1 kpc. The shaded regions correspond to the arm width. Different colours correspond to different spiral arms: Outer, green; Perseus, grey; Local, pink; Sagittarius-Carina, purple; Scutum-Centaurus, orange.

In the text

	Fig. 17. Left: same as Fig. 11. Centre: density of stars with MKs < −1 mag projected on the Galactic plane. Right: density of stars with MKs < −2 mag projected on the Galactic plane. We note that the density scales are different.
In the text

Fig. B.1. Posterior pdfs for four stars in our sample (their Gaia EDR3 source_id is displayed at top of each panel). The thick black solid line represents p(d | o, mKsΘKM, ΘSM, ΘCMD). The coloured thin lines illustrate the different components of the pdfs. The correlation terms in the covariance matrix of Eq. (B.12) are neglected here for illustrative purposes. The distributions are scaled so that they peak at unity.

In the text

	Fig. B.2. Effective luminosity function f(d, mKs|ΘCMD) used to determine our prior, normalised to unity.
In the text

	Fig. B.3. Fractional distance uncertainty dkin, 84 − dkin, 16/2 dkin for stars in the clean sample. The dashed light grey lines correspond to the 16th and 84th percentiles, the solid grey line corresponds to the median fractional error.
In the text

	Fig. B.4. Single star distance pdfs for the members of the three clusters IC 4756 (black), NGC 2112 (blue), and King 21 (orange). The grey vertical lines represent the cluster distances (respectively 469 pc, 1090 pc, and 3005 pc) derived by Cantat-Gaudin et al. (2020).
In the text