© The Authors 2024

1. Introduction

The spatial and chemical distribution of stars in the stellar disc of the Milky Way results from the complex and intricate mechanisms that formed and shaped the evolution of our Galaxy. By exploring the chemical diversity across different stellar populations, we can glean valuable insights into the formation, enrichment, and evolution of galaxies, shedding light on fundamental astrophysical processes (Gilmore et al. 1989; Bovy et al. 2012, 2016; Rix & Bovy 2013; Hayden et al. 2015).

Investigation of the detailed chemical composition of stars throughout the Milky Way enables us to discern signatures of subsequent stellar generations and their respective nucleosynthetic contributions, and it provides insight into the formation history and enrichment mechanisms of various stellar populations within our Galaxy (e.g., Freeman & Bland-Hawthorn 2002). Ultimately, this knowledge allows us to build a timeline of galactic evolution, and to uncover the role played by different progenitors, such as massive stars, novae, and supernovae, in shaping the chemical landscape. By comparing observed chemical abundances with predictions from stellar evolution and galactic chemical evolution models, we can refine our understanding of the physical processes involved (e.g., Bensby et al. 2014; Martig et al. 2016). This iterative process allows us to improve our models, leading to more realistic simulations of stellar populations, galactic structures, and the chemical evolution of galaxies.

Many of these insights arise from asking how many stars of which abundance (and possibly which age) are observed in different parts of the Milky Way. This information is provided by large spectroscopic surveys. To analyse these surveys, their survey selection function is pivotal to understand which stars at which position in the Galaxy were able to be included in any survey under consideration (e.g., Rix et al. 2021). These functions consider various observational constraints, such as magnitude limits, sky coverage, signal-to-noise ratios, and target selection algorithms. Neglecting or misinterpreting selection functions can lead to distorted conclusions about the chemical distribution of stars and their implications for galactic evolution. The present paper is part of a series in the GaiaUnlimited project¹, which aims to provide the astronomical community with tools and methods to account for the selection function of various subsets of the Gaia data (see also Rix et al. 2021; Cantat-Gaudin et al. 2023; Castro-Ginard et al. 2023).

In this study, we focus on red clump stars, which are low-mass stars in the core helium-burning stage that are located in a distinct region of the Hertzsprung-Russell diagram and exhibit remarkable homogeneity. They are therefore excellent standard candles (Girardi 2016; Hawkins et al. 2017; Hall et al. 2019; Chan & Bovy 2020), and their chemical abundances allow us to investigate the radial abundance gradients in the Milky Way disc. In particular, examining the spatial variations in their chemical composition in different galactic regions provides insights into the kinematic properties of stellar populations, the role of migration, and the formation history of the galactic disc. We combine precise chemical abundances and parallaxes from Gaia DR3 (Gaia Collaboration 2016, 2021) with abundances from the Sloan Digital Sky Survey (SDSS) Apache Point Observatory Galactic Evolution Experiment (APOGEE) seventeenth data release (DR17; Majewski et al. 2017; Abdurro’uf et al. 2022). The combination of these two surveys offers a wide sampling of the Milky Way disc, with Gaia providing a dense coverage of the solar neighbourhood, while APOGEE offers a sparser but deeper coverage of the inner and outer disc. The aim of our paper is two-fold: We provide an example of how very different survey selection functions of surveys can be combined so that they can be modelled jointly. We also refine our understanding of the radial structure of mono-abundance populations (Bovy et al. 2012) in the Galactic disc.

This paper is structured as follows. Section 2.1 introduces the sample of APOGEE red clump stars and its selection function. Section 2.2 builds a sample of Gaia DR3 red clump stars with precise elemental abundances from the Gaia GSP-Spec pipeline, and Sect. 2.3 details the construction of the corresponding selection function. In Sect. 3, we then explore the density distribution of mono-abundance populations, and we discuss our results in Sect. 4. Section 5 closes with concluding remarks.

2. Data sets and their selection functions

In this section we introduce the two data sets that we used to determine the radial disc profiles for stellar mono-abundance populations in the Galactic disc. We start with the SDSS/APOGEE sample, and then we turn to the Gaia sample, for which the selection function has not been worked out before. For Gaia, determining the selection function, that is, the probability of any actual or counterfactual star to be in the sample, is a multistep process.

2.1. The SDSS/APOGEE Data

We employed the data from the 17th data release (Abdurro’uf et al. 2022, hereafter DR17) of the SDSS survey, focussing on the APOGEE subset (Majewski et al. 2017), which collected high-resolution near-infrared spectroscopic data for cool stars (in particular, red giants). All data used in this study were accessed through the SDSS portal². The identification of red clump stars in SDSS/APOGEE is described in Bovy et al. (2014). Figure 1 shows their spatial and radial distributions.

Fig. 1. Spatial distribution of the SDSS/APOGEE red clump stars and Gaia RC-GSP-Spec stars (see Sect. 2.2.1). Left: location projected on the Galactic plane. Middle: Galactic altitude Z vs. galactocentric distance R. Right: distribution of galactocentric distance R.

2.1.1. APOGEE Main Red Star Sample selection function

The complex selection function for the APOGEE targets in SDSS is described in Zasowski et al. (2013, 2017). The selection function of the entire survey is not tractable in principle, but we take advantage of the fact that the APOGEE Main Red Star Sample was selected by randomly choosing targets from the 2MASS catalogue (Two Micron All Sky Survey; Skrutskie et al. 2006) in specific regions of the ((J − K_s)₀, H) colour-magnitude diagram, corresponding to the photometric space of the reddening-corrected (J − K_s) colour and the observed (uncorrected) H magnitude. This enables us to construct its selection function.

The Main Red Sample stars can be identified in the APOGEE DR17 table with the EXTRATARG flag set to 0, PROGRAMNAME to bulge, disk, disk1, disk2 or apogee, and a FIELD name following the pattern L+B or L-B (e.g., 120+12), which indicates that it is part of the semi-regular grid footprint. About half of the 730 000 DR17 sources are part of pointings like this. To construct the selection function of this sample, we followed the general approach introduced in Sect. 4 of Bovy et al. (2014; see also Bovy et al. 2016; Frankel et al. 2019), which we summarise here.

In each field observed by APOGEE, the survey selected targets from three different magnitude ranges, which are referred to as short cohort, medium cohort, and long cohort. These cohorts typically correspond to H-magnitude ranges [7, 12.2] ([10, 12.2] in APOGEE-2), [12.2, 12.8], and [12.8, 13.3], respectively. The details of the cohort definition of each field are listed in the survey allDesign table³. The design table also specifies the colour range of each cohort. For APOGEE-1 observations, it is typically a single (J − K_s)₀ > 0.5 colour limit. For APOGEE-2, each cohort exists in two different colour ranges, typically, 0.5 < (J − K_s)₀ < 0.8 and (J − K_s)₀ > 0.8, which are referred to as blue and red (see Fig. 2).

Fig. 2. Binning in intrinsic colour-magnitude space used by APOGEE (left) and APOGEE-2 (right) to select their Main Red Star Sample targets from the 2MASS catalogue. The exact magnitude limits between the three cohorts (short, medium, and long) depend on the field and are provided by the APOGEE collaboration in auxiliary files. The selection is computed as the ratio of successfully observed targets to available sources in each of the three (in APOGEE) or six (in APOGEE-2) regions.

For each design, the APOGEE survey selected a random subset of ∼150 targets from the 2MASS data. The exact number depends on whether other science programmes have targets of interest in a given field of view, but it is ultimately limited by the number of available spectrograph fibers (230 in APOGEE, and 250 in APOGEE-2). In order to collect larger numbers of spectra, many fields have been observed multiple times with different designs (and a different corresponding target list). The selection function for these fields can only be constructed when every associated design is constructed with the same H and (J − K_s)₀ limits. In some rare cases, the radius of the field or the central coordinates of the pointing were also changed between designs, which means that the selection function is intractable for these fields as well.

The APOGEE survey conveniently shares⁴ the 2MASS data used for the target selection in each field, along with the extinction factor⁵A_K used to obtain (J − K_s)₀ = J − K − 1.5 ⋅ A_K (Indebetouw et al. 2005). To obtain the selection fraction in each cohort, we simply computed the ratio of successfully observed stars to those available in 2MASS in the corresponding region of the intrinsic colour-magnitude space. We also accounted for the fact that in some cases (especially in the long cohorts of bulge fields), a fraction of the 2MASS stars has no associated value of A_K and therefore could not be targeted by APOGEE.

We find a very good agreement between our results and those published by Frankel et al. (2019) for the short cohorts of APOGEE DR14 (except that DR17 covers more fields, and some of them are more complete). We make the final table available through the gaiaunlimited⁶ Python package, along with a module to query it at any location, magnitude, and colour.

2.1.2. Effective selection function of the APOGEE red clump sample

Section 2.1.1 expresses the APOGEE selection function in terms of survey observables (ℓ,b, H, (J − K)₀). For a stellar population of a given absolute magnitude, we can construct the effective selection function in physical coordinates (ℓ,b, d) where d is the distance from the Sun.

We assumed that red clump stars have an absolute magnitude H₀ = −1.5 mag and intrinsic colour (J − K)₀ = 0.60 mag, placing them in the bluer bin of the APOGEE-2 cohorts. We computed the expected apparent magnitude H for a red clump star at any location (ℓ,b) and distance using the Lallement et al. (2022) extinction map and assuming A_H = 0.1727 × A₀ + 0.0201⁷. Figure 3 shows a top-down view of the complex, discontinuous pattern of the effective selection function at b = 0. The difference in bright magnitude limit between APOGEE-1 (H = 7 mag) and APOGEE-2 (H = 10 mag) is especially visible in the fourth quadrant (lower right).

Fig. 3. Effective selection function for the APOGEE DR17 red clump stars shown in the Galactic (X, Y) plane at b = 0. Each pencil beam corresponds to a different APOGEE pointing. The difference in bright magnitude limit between APOGEE-1 (H = 7 mag) and APOGEE-2 (H = 10 mag) is visible in the fourth quadrant (lower left corner).

2.2. Gaia red clump stars and chemical abundances

In this section, we first build a sample of bona fide red clump stars, and we then select those with the most reliable metallicities and α abundances and calibrate these abundances to the APOGEE abundance scale.

2.2.1. Gaia red clump stars

We selected red clump stars from the Gaia DR3 catalogue from their location in a colour-absolute magnitude diagram, that is, the colours and magnitude were not corrected for extinction but corrected for distance according to their parallaxes (using the prescriptions of Lindegren et al. 2021a). This study later focuses on stars with Gaia DR3 chemical abundances (see Sect. 2.2.2), which are all brighter than G = 14 mag and closer than 3 kpc (80% are closer than 1.5 kpc). The mean parallax uncertainty at this magnitude is 0.016 mas, making the difference between the true and the estimated absolute magnitude negligible. This allows us to neglect the impact of parallax uncertainties on the selection of red clump stars.

To account for the effect of extinction on the magnitude and (G − G_RP) colour of the RC stars, we first visualised the colour-magnitude distribution of stars in common between Gaia and the APOGEE DR17 catalogue of red clump stars (introduced in Sect. 2.1). We fit a fifth-order polynomial to the locus traced by the APOGEE RC stars and considered stellar sources within ±0.3 mag⁸ of this trace that satisfied 0.6 < (G − G_RP) < 1.2 to be bona fide RC stars (Fig. 4). Assuming Gaussian parallax errors, our selection retains over 99.9% of the red clump stars with G < 14.

Fig. 4. Definition of the Gaia red clump star sample (see Sect. 2.2.1). The shaded area between the dashed lines represents our selection from the Gaia DR3 photometry described in Sect. 2.2.1, based upon the APOGEE DR17 red clump candidates indicated by the black dots. The grey background indicates the density of all sources in the “gradient sample” of Gaia Collaboration (2023) for reference.

2.2.2. The RC from the GSP-Spec abundance sample

The gaiadr3.astrophysical_parameters table from the Gaia archive⁹ provides stellar parameters and elemental abundances from the General Stellar Parametriser from spectroscopy (hereafter GSP-Spec; Recio-Blanco et al. 2023) analysis of five million sources observed with the Radial Velocity Spectrometer (RVS; Cropper et al. 2018). For Gaia DR3, the GSP-Spec pipeline characterised sources whose combined spectrum has a signal-to-noise ratio greater than 20, effectively favouring brighter sources and sources that received more visits. From this sample, we further selected the gradient analysis sample, which is defined in the chemical cartography study of Gaia Collaboration (2023). By construction, this sub-sample supports the analysis of the relations between chemistry and Milky Way structure, stellar kinematics, and orbital parameters. In particular, this sample corresponds to a subset of stars with good GSP-Spec quality flags (flags_gspspec¹⁰; Recio-Blanco et al. 2023) that satisfy the astrometric quality criterion ruwe< 1.4 (Lindegren et al. 2021b). Listing 3 in Appendix B of Gaia Collaboration (2023) gives the ADQL query to retrieve this sample.

We further selected the stars with total abundance uncertainties in metallicity and an α-element abundance smaller than 0.1 dex to be able to map the spatial distribution of mono-abundance populations,

$$\begin{array}{c}\hfill \sqrt{e_{[\mathrm{Fe}/\mathrm{H}]}^2+e_{[\alpha /\mathrm{Fe}]}^2}<0.1\mathrm{dex},\end{array}$$(1)

where

$$\begin{array}{c}\hfill e_{\text{[Fe/H]}}=\mathtt{mh\_gspspec\_upper}-\mathtt{mh\_gspspec\_lower}\end{array}$$(2)

and

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill e_{[\alpha /\mathrm{Fe}]}=& \mathtt{alphafe\_gspspec\_upper}\hfill \\ \hfill & -\mathtt{alphafe\_gspspec\_lower}\hfill \end{array}\end{array}$$(3)

are the widths of the 68% confidence intervals on iron abundance and α abundance listed in the Gaia DR3 astrophysical parameter table, respectively.

In the rest of this paper, we refer to this sample of 356 755 stars as the RC-GSP-Spec sample. Figure 1 shows the spatial distribution of these stars in blue.

The selection based on the abundance precision indirectly corresponds to a selection in signal-to-noise ratio in the RVS spectra, which to first order depends on the apparent G magnitude and on the location on the sky (the RVS spectra of stars in regions that are visited more often by Gaia are combined from a larger number of exposures). Because the spectra of metal-rich stars have more prominent absorption lines, their abundances tend to be more precise, while spectra with weaker lines require a high signal-to-noise ratio to reach the same abundance precision. Section 2.3.3 shows that the selection has a more stringent effect on metal-poor stars, and we account for the varying depth of the selection function with metallicity.

2.2.3. Calibrating the GSP-Spec abundances on the APOGEE scale

We find 2450 red clump stars in common in our sample and APOGEE DR17. We compared the metallicities reported by both catalogues, and we find that the GSP-Spec metallicities are 0.08 dex lower than the APOGEE values (with a dispersion of 0.07 dex) with a strong temperature dependence (left panel of Fig. 5), which is consistent with the biases reported by Recio-Blanco et al. (2023). We therefore calibrated the Gaia metallicity scale on APOGEE by fitting a linear relation between their offset and the GSP-Spec effective temperature,

$$\begin{array}{c}\hfill \mathrm{\Delta }[\mathrm{M}/\mathrm{H}]=2.988\times {10}^{-4}\times (T_{\mathrm{eff}}-4700)-7.264\times {10}^{-2}\end{array}$$(4)

Fig. 5. Abundance-scaling relations of Gaia DR3 GSP-Spec and APOGEE DR17 ASCAP abundances for a sample of 2450 red clump stars in common. Left: offset between Gaia DR3 mh_gspspec metallicities and APOGEE DR17 ASCAP [M/H]. Right: offset between [α/Fe] (from GSP-Spec) and [Ca/Fe] (from APOGEE). In both panels, the dashed line indicates the correction we applied to calibrate Gaia on the APOGEE scale, given in Eqs. (4) and (5).

where Δ[M/H] = [M/H]_{GSP − Spec} − [M/H]_APOGEE.

We also calibrated the GSP-Spec α abundances on the APOGEE scale. For stars of spectral type F or cooler, the spectral range covered by the Gaia Radial Velocity Spectrometer (845–872 nm) is dominated by calcium triplet lines (see e.g., Fig. 17 Cropper et al. 2018), which means that the alphafe_gspspec values listed in Gaia DR3 are mostly derived from Ca lines. We followed the approach of Recio-Blanco et al. (2023), and calibrated the GSP-Spec α abundances on the APOGEE DR17 [Ca/Fe] abundances rather than [α/Fe]. We find that the alphafe_gspspec values are lower by 0.07 dex on average (with a dispersion of 0.04 dex), and that the bias does not correlate with the effective temperature or surface gravity, but depends on metallicity (right panel of Fig. 5). We modelled this bias as a linear function of GSP-Spec metallicity,

$$\begin{array}{c}\hfill \mathrm{\Delta }[\alpha /\mathrm{H}]=-4.481\times {10}^{-2}\times {[\mathrm{M}/\mathrm{H}]}_{\mathrm{GSP}-\mathrm{Spec}}-4.414\times {10}^{-2}\end{array}$$(5)

where Δ[α/H]=[α/H]_{GSP − Spec} − [α/H]_APOGEE. We then applied the corrections in Eqs. (4) and (5) to the entire RC-GSP-Spec sample.

Figure 6 shows the distribution of APOGEE RC stars in the [α/Fe]-metallicity and in the [Ca/Fe]-metallicity planes. The APOGEE [α/Fe] (derived as the mean of [O/Fe], [Mg/Fe], [Si/Fe], [S/Fe], and [Ca/Fe]) estimates have a greater discriminating power than [Ca/Fe] alone¹¹. The latter barely allows us to distinguish stars between the high- and the low-α disc at low metallicity ([M/H] < −0.3 dex), and it does this not at all in the metal-rich regime.

Fig. 6. Red clump stars in chemical abundance space. Left: distribution of the APOGEE RC stars in the [α/Fe]-metallicity plane. The α-rich stars are highlighted manually. Middle: same stars in the [Ca/Fe]-metallicity plane. The [Ca/Fe] abundances do not have the same discriminating power as [α/Fe], which accounts for O, Mg, Si, S, and Ca for APOGEE. Right: distribution of the Gaia RC-GSP-Spec stars in the [α/Fe]-metallicity plane. The contours enclose 99%, 95%, 90%, 80%, and 50% of the sample. The Gaiaα-abundances are dominated by Ca lines and do not allow us to reliably distinguish the α-rich from the α-poor disc stars.

For the rest of this study, we focus on the α-poor disc, which we define as the stars with an APOGEE [Ca/Fe] or a calibrated GSP-Spec [α/Fe] lower than 0.1 dex. We discuss the consequences of possible contamination by the thin disc in Sect. 4.

2.3. Selection function of the Gaia RC-GSP-Spec sample

Building the selection function S for a given sample means identifying the different cuts and selections that cause objects of interest (in the present case, any red clump star in the Milky Way) to be absent from the sample, in order to link the observed number of objects N_obs to the true number N_true, where N_obs = N_true S.

If an RC star is not present in the RC-GSP-Spec sample, it can for one of these three reasons: (i) it is not in the Gaia DR3 catalogue; (ii) it is in this catalogue, but the selection performed in Sect. 2.2.1 failed to identify it as an RC star; and (iii) the spectroscopic measurements associated with this star are not of sufficient quality to pass the selections applied in Sect. 2.2.2. These are independent conditions that result in a product of selection function terms. When we wish to know how many RC stars lie in direction (ℓ,b) and at distance d, N_RC,obs, in the sample, this can be written as

$$\begin{array}{cc}\hfill N_{\text{RC,obs}}(\ell ,b,d)=& N_{\text{RC,MW}}\times \frac{N_{\text{RC,Gaia}}}{N_{\text{RC,MW}}}\times \frac{N_{\text{RC,selected}}}{N_{\text{RC,Gaia}}}\hfill \\ \hfill & \times \frac{N_{\text{flags}}}{N_{\text{RC,selected}}},\hfill \end{array}$$(6)

where N_RC,MW is the true number of RC stars at that location in the Milky Way (which we ultimately wish to infer), N_RC,Gaia is the number of RC stars at that location that is listed in the Gaia catalogue, N_RC,selected is the number of stars we were able to successfully identify as RC stars, and N_flags is the number of RC stars that also satisfy our chosen quality flags. (We omitted the location dependence on the right-hand side of the equal sign in Eq. (6) to simplify the notation.)

We need to express the fraction in terms of direct observables, but the distance d is not an observable. Because RC stars are standard candles with the same intrinsic magnitude M_G = 0.5 mag, however, we can rewrite Eq. (6) as a function of (ℓ,b, G). We estimate G at any (ℓ,b, d) through the distance d and a 3D extinction map. The completeness is then the product of the above selections,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{N_{\text{RC,obs}}}{N_{\text{RC,MW}}}(\ell ,b,G(d,A_G))=& \underset{=1\text{if}G<14}{\underbrace{S_{\mathit{Gaia}}}}\underset{=1\text{if}G<14}{\underbrace{S_{\text{RC selection}}}}\underset{=0\text{if}G>14}{\underbrace{{\overline{S}}_{\text{RC-GSP-Spec}}}},\hfill \end{array}\end{array}$$(7)

where we have used S as a shorthand for the ratios in Eq. (6), and all the terms on the right-hand side are functions of (ℓ,b, G); and G(d, A_G) is shorthand for G(d, A_G(ℓ, b, d)). We elaborate on the symbol choice of ${\overline{S}}_{\text{RC-GSP-Spec}}$ in Sect. 2.3.3.

In Eq. (7), the term S_Gaia is 1 if G < 14 because the overall Gaia catalogue is essentially complete at these bright magnitudes (Cantat-Gaudin et al. 2023). S_{RC selection} is also 1 in this magnitude range because the parallax errors are sufficiently small to enable the selection of virtually all red clump stars (see Sect. 2.2.1 and Fig. 4). The last term in this equation describes the fraction of identified RC stars that also satisfy the quality flags we imposed on the RC-GSP-Spec sample. We remark that none of our final RC-GSP-Spec sample stars is fainter than G = 14 mag. Since the completeness of the sample is 0% for G > 14 mag (setting the selection function to zero), we only need to describe the behaviour of the other terms in the brighter regime. We construct S_RC-GSP-Spec(ℓ, b, G) in Sect. 2.3.1.

2.3.1. The Gaia selection function dependence on magnitude and sky position

We investigated the completeness S_RC-GSP-Spec(ℓ, b, G) of the RC-GSP-Spec sample (constructed in Sect. 2.2.2) with respect to all 1.8 M Gaia DR3 red clump stars, which we assumed to be fully complete in the magnitude range in which precise GSP-Spec abundances are available (G < 14). Computing the selection function as ratios with respect to a more complete sample is a common approach in astronomy (see e.g., Rybizki et al. 2021; Castro-Ginard et al. 2023, for applications to Gaia data).

In order to understand the general behaviour of the selection function with G, we first computed the count ratios in bins of 0.5 mag computed over the entire sky and then in HEALPix regions at level 0, which correspond to 12 areas of about 3400 deg² each¹². We show both cases in Fig. 7. The completeness is close to 80% for sources with 8 < G < 10.5 mag, and it then decreases to zero near 13.5 mag. The completeness also decreases slightly for G < 8 mag. Because the parallax correction of Lindegren et al. (2021a) we employed in Sect. 2.2.1 is not available for stars brighter than G = 6 mag, our completeness is zero at the bright end (but only 0.05% of the sources with available GSP-Spec abundances are brighter than G = 6).

Fig. 7. Completeness of our RC-GSP-Spec sample with respect to all Gaia DR3 red clump stars computed in bins of 0.5 mag (black line). The grey lines, computed independently in the 12 HEALPixels at level 0, show the same quantity.

We aimed to construct a selection function that reflects both the spatial and the magnitude dependence of the completeness. We explored different binning schemes to find a compromise between spatial resolution, magnitude resolution, and reasonable noise levels. Because the magnitude range 6 < G < 8 mag only contains 2% of the RC-GSP-Spec sample, we grouped sources with 6 < G < 10.5 mag in a single magnitude bin and split the rest of the sample into bins of 0.5 mag. We tessellated the sky into HEALPix regions at level 3 (about 7.4 deg²). Although the fine details of the Gaia scanning pattern are better visualised with a finer spatial resolution, increasing the HEALPix level would produce noisier maps or require the use of broader magnitude bins (degrading the distance resolution of the selection function).

Figure 8 shows all-sky completeness maps in four magnitude intervals. Even at the bright end, Gaia is not complete. It is especially visible around (ℓ,b) = (5, −5) in the top left panel of Fig. 8, which contains a darker area of a few degrees below and to the east of the Galactic centre (around (ℓ,b) = (5, −5)). This region (which includes Baade’s Window) corresponds to a direction with particularly low extinction, which leads to an effectively high Gaia source density near the crowding limit of the survey. The source density is so high that it affects the completeness of the Gaia RVS catalogue itself (see Fig. 2 in Cantat-Gaudin et al. 2023). Patterns appear in our map for the magnitude range 11 − 11.5 mag. These correspond to the imprint of the Gaia scanning law and appear as stripes of lower completeness, mostly caused by the ruwe < 1.4 filter applied to the Gaia Collaboration (2023) GSP-Spec sample. We discuss this bias in Appendix A. At G > 11.5 mag, the more complete regions correspond to areas that the Gaia satellite observed more frequently (more visits; bottom panels of Fig. 8).

Fig. 8. Selection function of the RC-GSP-Spec sample, constructed from count ratios in spatial (healpix level 3) and magnitude bins. The maps are in Galactic coordinates, centred on (ℓ,b) = (0, 0), with ℓ increasing to the left.

Another possible approach that we did not follow here would have been to construct an analytic function of basic observables (Rix et al. 2021, for instance, only use the G magnitude) to predict the probability that a given source passes the selection filters. This approach is simple and powerful when the effect of cuts on, measurement errors or parallax signal-to-noise ratio is to be modelled, for example, which are mostly smooth functions of G in any given direction. The approach would also allow us to construct a selection function whose arguments are observable quantities even when they are not always available for all sources. The drawback of this approach is that it requires the user to accurately describe every single cut to the data with an appropriate model. In the present case, this would mean describing the mentioned ruwe quirk near G ∼ 11.5 (described in Appendix A) as a function of position and magnitude, as well as the various quality cuts applied to the GSP-Spec sample by Gaia Collaboration (2023). Our chosen approach allows us to compute S_RC-GSP-Spec in a single ratio step.

2.3.2. Effective red clump star selection function in Galactic coordinates

The selection function S_RC-GSP-Spec we constructed in Sect. 2.3 is a function of sky coordinates and apparent G magnitude. In the context of mapping the Milky Way disc, we are generally interested in knowing the completeness of a given sample as a function of spatial Galactic coordinates.

The apparent magnitude G of a star at any location in the Milky Way can be expressed as

$$\begin{array}{c}\hfill G(\ell ,b,d)=M_G+5{log}_{10}d-5+A_G(\ell ,b,d,T_{\mathrm{eff}}),\end{array}$$(8)

where M_G is the absolute magnitude of the star, d is its distance (in parsecs), and A_G is the total amount of interstellar extinction in the G band to the star, which we could estimate from a 3D extinction map and the effective temperature T_eff of that star¹³. For any particular choice of a star (defined by a value of reddening-free magnitude G₀ and temperature T_eff), we can use Eq. (8) to turn a function of (ℓ,b, G) into the effective selection function S_RC-GSP-Spec(ℓ, b, d). We refer to Sect. 3 of Bovy (2017) for an illustration of how a given selection function translates into a different effective selection function for different stellar populations in the Tycho-Gaia Astrometric Solution catalogue (Michalik et al. 2015).

In this paper, we study the spatial distribution of red clump stars, for which we assumed M_G = 0.5 mag (Ruiz-Dern et al. 2018; Hall et al. 2019) and T_eff = 4700 K (Girardi 2016). We used the all-sky dust maps of Lallement et al. (2022), and converted their extinction A₀ into A_G using the precomputed values provided by the Python package dustapprox¹⁴ for this effective temperature. Figure 9 shows the expected G magnitude and the corresponding effective selection function for red clump stars in the Galactic mid-plane.

Fig. 9. Effect of distance and extinction on the sample completeness. Top left: cumulative extinction at Z = 0, computed from the Lallement et al. (2022) dust map. Top right: predicted apparent G magnitude for a red clump star, assuming an intrinsic MG = 0.5 mag. Bottom left: selection function SGaia chem(X, Y, Z) plotted for b = 0. Bottom right: XY distribution of red clump stars with |Z|< 15 pc in our RC-GSP-Spec sample (Sect. 2.2.2). The Galactic centre is to the right in every panel.

We find no overall systematics between the observed G magnitudes and those predicted by the dust extinction model, except in a few distant regions (> 2 kpc from us) where the Lallement et al. (2022) map allows for zero extinction, and the observed red clump stars can be fainter than predicted by ∼0.1 mag. In high-extinction regions, some stars are slightly brighter than predicted by the model. This effect may be caused by the spatial resolution of the map, which smooths out the fine structure of the dust clouds at scales smaller than 5 pc. Because the selection function drops from 80% to 0% over the course of two magnitudes (Fig. 7) and is constructed in ranges of 0.5 in magnitudes, we did not consider these 0.1 mag biases (which only affect sparsely populated regions) to have a significant impact.

2.3.3. Metallicity dependence of the selection function

One of our quality flags for inclusion in the RC-GSP-Spec sample is that the combined metallicity and alpha-abundance uncertainty of a star be smaller than 0.1 dex (Eq. (1)). This depends foremost on the stellar brightness, but at a given brightness also on the stellar metallicity. Hence, our selection not only favours brighter sources and regions of the sky that were more frequently visited by Gaia (Fig. 8), but also favours high-metallicity stars whose spectral lines are inherently more prominent and yield smaller uncertainties. As a result, metal-poor stars included in the RC-GSP-Spec sample tend to be slightly brighter and more nearby than metal-rich stars, as Fig. 10 illustrates.

Fig. 10. Effect of metallicity on the selection function. Top left: distribution of metal-poor RC-GSP-Spec stars projected the Galactic (X, Y) plane for stars with |Z|< 200 pc alone. The Sun is at (X, Y) = (0, 0), and the Galactic centre is to the right. Top right: same for the most metal-rich stars in the sample. Metal-rich stars contained in the sample reach larger distances from the Sun. Bottom: G magnitude distribution in the entire sample (black) and in four ranges of metallicity. The offset ΔG used in Eq. (9) is indicated for the most metal-poor metallicity group.

The selection function ${\overline{S}}_{\text{RC-GSP-Spec}}$ shown in Figs. 8 and 9 only represents the metallicity-averaged selection function for the whole RC-GSP-Spec sample (black line in the lower panel of Fig. 10); hence the choice of $\overline{S}$. In practice, the apparent G-magnitude distribution of stars in the RC-GSP-Spec sample depends quite strongly on metallicity, as the coloured lines for different metallicity subsets in Fig. 10 show. These differences are most likely mostly due to differences in S_RC-GSP-Spec(G), not any intrinsic differences in the spatial distribution. We took an empirical approach to account for this metallicity dependence of the selection function. We did this by assuming that for any metallicity, the selection function is well described by the metallicity-averaged function (black line in Fig. 10), just shifted by the amount by which the peaks in the metallicity sub-bins differ from the peak in the metallicity-averaged distribution. For instance, for the metal-poor sub-sample ([M/H] < −0.3), we shifted the mean selection function so that it was 0.75 mag brighter.

Expressed mathematically, we did this via

$$\begin{array}{c}\hfill S_{}\text{RC-GSP-Spec}(\ell ,b,G,[\mathrm{M}/\mathrm{H}])\equiv \\ \hfill {\overline{S}}_{\text{RC-GSP-Spec}}(\ell ,b,G-\mathrm{\Delta }G),\end{array}$$(9)

where ΔG = ΔG([M/H]) is the metallicity-dependent shift of the peaks in the distributions in the bottom panel of Fig. 10. We resorted to this empirical approach because metallicities are not available for all stars in the Gaia catalogue, which would be needed to express the selection function as a ratio in the style of Eq. (6).

This approach could be refined even further by accounting for the fact that RC stars are not perfect standard candles and that their absolute magnitudes vary with age and metallicity (from M_G ∼ 0.35 for our most metal-poor stars to ∼0.6 according to PARSEC isochrones; Bressan et al. 2012). The mean absolute magnitude also depends on the age distribution within each metallicity group, with older RC stars being slightly fainter. Second-order corrections beyond the approach used in this study would therefore require deeper assumptions on the age-metallicity distribution at any point in the Milky Way. The overall effect of this finer correction would be to reduce the estimated effective volume of the metal-poor stars by ∼0.1 in distance modulus (because they are brighter, a given observed G corresponds to a shorter distance), and to increase the effective volume of the metal-rich group by a similar amount¹⁵, which we expect would have a negligible impact on the overall results of this study. We point out that this effect is present in the construction of the APOGEE selection function as well because the absolute H-band magnitude of RC stars also depends on their age and metallicity.

3. Modelling the Galactic disc profile for mono-abundance populations

In this section, we first qualitatively visualise the density profile of the disc (Sect. 3.1). This simple approach is statistically unsound, however. In Sect. 3.2 we build a forward model of the spatial distribution that accounts for the selection functions of both Gaia and APOGEE, and we provide quantitative results.

3.1. Qualitative Gaia density profile in the solar neighbourhood

We focused on the α-poor stars ([α/Fe]< 0.1 dex) in the RC-GSP-Spec sample. We binned them into six bins of [M/H], and a 200 pc grid in (X, Y, Z) coordinates and computed the mean selection function in each bin (by sampling it on a much finer grid and then averaging) to estimate the completeness. We obtained corrected counts in each bin by dividing the observed counts by the completeness.

To visualise the radial density profile in each metallicity bin, we plot the resulting density as a function of galactocentric distance R for the bins closest to the Galactic mid-plane ( − 200 < Z < 200 pc), assuming R_⊙ = 8.15 kpc (Reid et al. 2019).

Figure 11 shows the resulting completeness-corrected galactocentric distance distributions per [M/H] bin for the stars closest to the Galactic mid-plane ( − 200 < Z < 200 pc). The density increases towards the Galactic centre in all [M/H] bins, except for the most metal-poor ([M/H] < −0.3 dex) population, and the slope steepens with metallicity. The corrected density is noisier in regions of low completeness (darker points in Fig. 11) because they include fewer stars and because of the uncertainties on the selection function itself.

Fig. 11. Completeness-corrected number of α-poor RC-GSP-Spec stars in bins of Galactic (X, Y) bins of 200 × 200 × 200 pc near the Galactic mid-plane (−200 < Z < 200 pc) as a function of galactocentric distance in parsecs. We selected α-poor stars as [α/Fe]< 0.1 dex. A dot represents the corrected density for regions in which the expected completeness is greater than 20%. The top and bottom panels correspond to the data displayed in the top left and top right panels of Fig. 10, respectively. We stress that the observed counts should only be divided by completeness for qualitative visualisation and data exploration purposes. Section 3.3 implements a forward-modelling approach to apply the selection function to a density model.

An example of the vertical density profile traced by corrected counts is shown in Appendix C for the metal-rich ([M/H] > 0.3 dex) population. In the rest of this section, we fit a density model to the full 3D distribution of stars.

3.2. Gaia and APOGEE: Tracing the large-scale structure of the disc

In this section we use a forward-modelling approach to fit an analytical density model to the observed stellar distribution.

3.2.1. Density model

We describe the stellar density distribution ν_⋆ in the Milky Way disc with the model used by Bovy et al. (2016) and Mackereth et al. (2017), which assumes that the density is axisymmetric and that the radial and vertical density profiles are separable,

$$\begin{array}{cc}& {\nu }_{\star }(R,\varphi ,Z)={\nu }_{\star }(R,Z)=\mathrm{\Sigma }(R)\zeta (Z|R),\hfill \\ \hfill & \text{where}{\displaystyle \int }\zeta (Z|R)\mathrm{d}Z=1.\hfill \end{array}$$(10)

In this model, we assume that the radial density profile is a radially broken exponential,

$$\begin{array}{c}\hfill ln\mathrm{\Sigma }(R)\propto \{\begin{array}{cc}-h_{R,\text{in}}^{-1}(R-R_{\odot }),\hfill & \text{where}R\le R_{\text{peak}}\hfill \\ -h_{R,\text{out}}^{-1}(R-R_{\odot }),\hfill & \text{where}R>R_{\text{peak}}\hfill \end{array},\end{array}$$(11)

where R_⊙ = 8.15 kpc is the solar galactocentric radius (Reid et al. 2019), R_peak is the radius at which the change of slope occurs, and h_R, in and h_R, out are the scale length inside and outside the break radius, respectively. The density continuity at R = R_⊙ sets the relative normalisation of the two parts.

We assumed that the vertical profile is single exponential with a flare, making the scale height itself an exponential function of the galactocentric distance R,

$$\begin{array}{c}\hfill h_Z^{-1}(R)=h_{Z,\odot }^{-1}exp\left[R_{\mathrm{flare}}^{-1}(R-R_{\odot })\right],\end{array}$$(12)

with h_Z, ⊙ the scale height at the location of the Sun, and R_flare the flaring scale length. A negative value of R_flare corresponds to a $h_Z^{-1}$ decreasing with R, that is, a scale height h_Z increasing at larger radii. We illustrate the impact of flaring on the vertical structure in Appendix C. The (normalised) vertical density profile ζ at any R is given by

$$\begin{array}{c}\hfill ln\zeta (Z|R)=-h_Z^{-1}|Z|-ln2h_Z(R).\end{array}$$(13)

3.2.2. APOGEE and Gaia joint selection function

The APOGEE targets are chosen by randomly picking sources from the 2MASS catalogue, while the RC-GSP-Spec (Gaia) data do not rely on another input catalogue for its target selection or data processing. We can therefore consider the two samples as independent. According to set theory, the probability that a given RC star ends up in either or both of the catalogues is thus obtained as the sum, minus the probability that it is present in both,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill S_{\mathrm{combined}}=& S_{\text{RC-GSP-Spec}}+S_{\mathrm{APOGEE}}\hfill \\ \hfill & -S_{\text{RC-GSP-Spec}}\times S_{\mathrm{APOGEE}},\hfill \end{array}\end{array}$$(14)

where all terms are functions of (ℓ,b, d) (or, equivalently, X, Y, and Z), expressing the selection probability at any given location in the Milky Way, and S_RC-GSP-Spec also depends on metallicity (Eq. (9); Fig. 12).

Fig. 12. Effective selection function (shown at b = 0°) for the combination of our Gaia and APOGEE samples obtained from the effective selection functions displayed in Figs. 3 and 9 and Eq. (14). The completeness reaches over 80% in the solar neighbourhood dominated by the Gaia sample.

We stress that if the APOGEE survey assigned higher observing priorities to targets with available Gaia data, or if Gaia relied on APOGEE data to refine its GSP-Spec values, the assumption of independence would not hold, and it would be more complicated to calculate the joint selection function. For instance, Rybizki et al. (2021) have shown that in Gaia DR2, the selection function of the Radial Velocity Sample was strongly shaped by the Initial Gaia Source List (IGSL; Smart & Nicastro 2014) and is therefore correlated with the selection function of SDSS DR9 (York et al. 2000; Ahn et al. 2012) and GSC2.3 (Lasker et al. 2008) from which the IGSL was built. This issue was solved in the Gaia DR3 catalogue because the processing of the Radial Velocity Sample no longer relies on literature catalogues.

As Gaia has quickly become the reference catalogue for identifying objects of interest (from Milky Way stars to distant quasars), the next generation of large-scale spectroscopic surveys such as WEAVE (Dalton et al. 2014; Jin et al. 2023) or 4MOST (de Jong et al. 2012, 2014) will produce catalogues whose selection function will correlate strongly with that of Gaia, and in certain cases, they will even simply constitute subsamples of the Gaia data.

3.3. Model fitting

For a given choice of model parameters

$$\begin{array}{c}\hfill \theta =(h_{R,\text{in}},h_{R,\text{out}},R_{\text{peak}},h_{Z,\odot },R_{\mathrm{flare}}^{-1}),\end{array}$$

the model density at location (X, Y, Z) is ν_⋆(X, Y, Z|θ), and thus, the predicted observed number of stars is

$$\begin{array}{c}\hfill N_{\mathrm{pred}}(X,Y,Z|\theta )=N_0{\nu }_{\star }(X,Y,Z|\theta )S(X,Y,Z),\end{array}$$(15)

where S represents the selection function, which is the APOGEE selection function (see Sect. 2.1.1), the RC-GSP-Spec metallicity-dependent selection function (see Sect. 2.3.3), or the combined one (see Eq. (14)). We fixed N₀ for each model so that the total number of predicted counts matched the total number of observed stars in the bins where S > 0.

We fit the parameter set θ by maximising the Poisson likelihoods of the observed stellar counts with those predicted by our model per spatial bins (X, Y, Z),

$$\begin{array}{c}\hfill log\mathcal{L}(\theta )={\displaystyle \sum _{X,Y,Z}logP(N_{\mathrm{obs}}(X,Y,Z),N_{\mathrm{pred}}(X,Y,Z|\theta )),}\end{array}$$(16)

where the (X, Y, Z) bins are 200 × 200 × 200 pc cubes. To avoid biases due to a possibly inaccurate extrapolation of the Lallement et al. (2022) extinction map towards the inner disc, we also limited the fitting to bins with R > 4 kpc (this only discarded 2% of the APOGEE stars and none of the Gaia stars).

We sampled the log-posterior with the Markov chain Monte Carlo sampler emcee (Foreman-Mackey et al. 2013), exploring the parameter space for 1500 steps with 16 walkers and imposing very broad flat priors.

We performed the fitting procedure three times: (i) Using only the APOGEE RC stars, (ii) using only the Gaia RC-GSP-Spec stars, and (iii) using the combined sample. Figure 13 shows the resulting density profiles. In all cases, we discarded stars whose APOGEE [Ca/Fe] or Gaia GSP-Spec [α/Fe] estimates were greater than 0.1 dex to focus on the α-poor disc. We discuss these results in Sect. 4 and compare them to previous studies of the density profile of mono-abundance populations in the Milky Way.

Fig. 13. Density profiles of red clump stars from APOGEE, Gaia RC-GSP-Spec, and their combined data sets (from top to bottom). The left panels show the radial density profiles Σ(R) obtained in six ranges of metallicity in the α-poor disc ([α/Fe]< 0.1 dex). The profiles are shifted by an arbitrary offset to be visually separated. The panels on the right represent the flaring scale length and vertical scale height predicted by our model for each metallicity bin. The shaded areas (left) and error bars (right) represent the 95% confidence interval. The Gaia RC-GSP-Spec star sample covers a smaller range in R, especially at low metallicity. We list the best-fit parameters in Table D.1.

We also performed the analysis for the APOGEE stars by selecting the α-poor stars based on their APOGEE [α/Fe]. The results were very similar, suggesting that we are not strongly biased by our inability to separate the α-rich disc (whose relative weight is lower) from the α-poor disc.

To assess the impact of the uncertainty on the selection function, we performed the same fitting procedure on the Gaia RC-GSP-Spec sample restricted to regions in which the sample is at least 80% complete. The results are shown in Fig. B.1. This restriction reduces the radial coverage of the data greatly (to the solar neighbourhood), whicgh increases the parameter uncertainties. However, the radial density profiles and vertical scale lengths are statistically identical to those obtained with the entire RC-GSP-Spec sample.

4. Discussion

The results we obtained from APOGEE DR17 data (top row of Fig. 13) are overall consistent with those of Bovy et al. (2016) and Mackereth et al. (2017), who both used APOGEE DR12 data. In this section, we discuss the difference between the radial density profiles obtained with APOGEE and Gaia, which we attribute to their very different sampling of the disc. We also discuss the improved description of the vertical density profile obtained by combining the two surveys. Finally, we discuss some remaining sources of uncertainty that are not currently accounted for in our error budget.

4.1. Radial density profile

The radial profiles traced by APOGEE data peak at a galactocentric distance R_peak that continuously increases with decreasing metallicity. The radial density distribution traced by the Gaia RC-GSP-Spec sample (middle row of Fig. 13) appears to be broadly consistent, but with some major differences. The radial density profiles of the two populations with [M/H] < −0.15 appear mostly flat, which we attribute to the much shorter coverage of this sample, especially for low metallicities. All other populations peak at R_peak ∼ 7 kpc, which is at odds with the APOGEE results, where the near-solar metallicity samples peak between 8 and 10 kpc. The density profiles shown in Fig. 11 indicate that in this metallicity range, the density clearly increases with decreasing R. We attribute the difference between the profiles obtained from the APOGEE and Gaia data to their very different spatial coverage. Figure 1 shows that APOGEE mostly covers the disc outside the solar neighbourhood, while the Gaia RC-GSP-Spec sample provides a dense coverage for 7 < R < 9 kpc.

It is also likely that the true density profile of the α-poor population is not perfectly smooth (owing, e.g., to perturbations by the spiral arms and other non-axisymmetric features). The fit to the APOGEE data therefore describes the overall density distribution over a large scale, with a focus on the disc outside the solar radius, while Gaia provides stronger constraints on the local density. Mackereth et al. (2017) reported that their best-fit model to the APOGEE data underpredicts the stellar density in the solar neighbourhood compared to McKee et al. (2015), for example. The difference between APOGEE predictions for the solar neighbourhood and other data sets could be resolved with a more flexible (but more complex) model, which would allow us to capture the detailed small-scale variations revealed by Gaia in addition to the general trends traced by APOGEE.

4.2. Vertical density profile

The shallower depth of the Gaia RC-GSP-Spec sample is counterbalanced by its nearly isotropic coverage, while the deeper APOGEE footprint is mostly dominated by low-latitude fields. The Gaia RC-GSP-Spec provides better leverage than APOGEE for a constraint of the local vertical scale height of the disc (h_Z, ⊙). The uncertainties on the scale heights obtained from the Gaia sample are smaller by two to five times than those obtained with APOGEE, and they confirm previous reports that the scale height in the solar neighbourhood decreases with increasing metallicity (e.g., Bovy et al. 2016; Yu et al. 2021; Lian et al. 2022). The only population for which Gaia and APOGEE return discrepant values is the most metal-poor group, for which APOGEE suggests h_Z, ⊙ = 340 ± 70 pc, while Gaia indicates 507 ± 30 pc. This difference is likely due to contamination by the α-rich disc, whose scale height is significantly larger (∼1 kpc according to Bovy et al. 2016), but whose presence in the outer disc (and therefore in the APOGEE sample) is less significant due to its shorter scale length.

An unexpected result driven by the Gaia data is that the flaring of the α-poor (young) disc appears to increase with metallicity, with $R_{\mathrm{flare}}^{-1}$ ranging from ∼ − 0.05 to −0.20 kpc⁻¹ within the metallicity range of our sample. In other words, the scale height h_Z of the more metal-rich stellar populations increases more strongly with R (Eq. (12)). This effect has never been reported in the literature to our knowledge. The large uncertainties on $R_{\mathrm{flare}}^{-1}$ obtained from APOGEE (in this study and in Bovy et al. 2016) do not allow us to draw conclusions on the relation between metallicity and flaring. However, splitting the disc into age groups, Mackereth et al. (2017) and Anders et al. (2023) reported that the youngest populations, which exhibit the smallest scale heights, are also more strongly flared. This is expected if the flaring is due to external perturbations because in this scenario, the scale height of the perturbed populations become asymptotically similar in the outer disc, so that populations that are geometrically thinner in the inner disc are more strongly flared by definition. In the model we employed, the variation in the vertical scale height across the entire range of galactocentric distances covered by the data is described with a single number ($R_{\mathrm{flare}}^{-1}$), which makes it difficult to distinguish between specific scenarios. Lian et al. (2022) remarked that using a model in which the scale height is a linear function of R (rather than exponential) yields negligible flaring, and they proposed that a piece-wise model could highlight that flaring mostly affects the outer disc (beyond 10 kpc). Using this approach (but not resolving the Milky Way into mono-abundance or mono-age populations), Mateu & Vivas (2018) proposed that flaring sets on beyond 11 kpc. Numerical simulations have shown that flaring can be caused by multiple mechanisms (see e.g., Garcíadela Cruz et al. 2021, and references therein), such as external perturbations from satellite galaxies, stellar migration, or a tilt in the mass distribution of the Galactic halo (Han et al. 2023). Sotillo-Ramos et al. (2023) reported “a large variety of flaring flavours and amounts” in the TNG50 simulations (Pillepich et al. 2019), which indicates that the flaring of the Milky Way disc might be difficult to capture with a simple, smooth density model.

We show in Sect. 2.2.2 that the Gaia GSP-Spec abundances cannot distinguish between the high-α and low-α disc components. However, the strong flaring observed in the metal-rich population is probably not due to contamination by (geometrically thicker) α-rich stars, which are rare in the outer disc due to their short scale length (Bovy et al. 2016). At high metallicity ([M/H] > 0.2 dex), it is in fact difficult even for the APOGEE Mg-informed α abundances to distinguish the two disc components (e.g., Leung et al. 2023). The evolution of the flaring scale length with metallicity might then be investigated by simultaneously fitting two density models, providing a morphological definition of the α-rich and α-poor component of the disc at that metallicity (resembling the traditional thin/thick distinction, see e.g., Gilmore & Reid 1983; Jurić et al. 2008). A clearer separation of the discs might be obtained with a sample using Gaia XP abundances (Andrae et al. 2023) or the upcoming ground-based surveys WEAVE and 4MOST. This task is beyond the scope of this paper, however, which focuses on defining a selection function for the combination of two independent astronomical surveys. The α-rich contaminants could also be separated through a kinematic selection (e.g., Bensby et al. 2003), whose selection function would have to be constructed as well.

4.3. Uncertainty on the selection function

We did not consider the uncertainty of the selection function itself. A possible extension of this work could involve parametrising the selection function and simultaneously fitting for the selection function and the density distribution in a hierarchical fit. The model could include a correlation between metallicity and absolute magnitude of the RC, with priors describing some assumptions on the age-metallicity relation across the Galaxy (effectively acting as a nuisance parameter). For a more rigorous and quantitative comparison of the model uncertainties, we should also take the impact of the choice of extinction map into account: While we used the all-sky Lallement et al. (2022) map, Bovy et al. (2016) and Mackereth et al. (2017) used a combination of maps and did not cover the entire sky. We therefore stress that our reported best-fit values and corresponding uncertainties are ultimately tied to a given choice of extinction model.

5. Conclusion

We have presented the construction of a selection function for a sample of Gaia DR3 red clump stars with precise chemical abundances, for which we included a 3D extinction. We proposed a simple approximation to account for the metallicity dependence of the selection function. We provided a framework to define the combined selection function of Gaia and APOGEE red clump stars that accounts for the different abundance scales of the two surveys. The α abundances obtained from Gaia spectra are not sufficient to clearly separate the α-rich disc, and therefore, we focused on the α-poor disc.

The APOGEE survey provides a sparser but much deeper coverage of the disc, while Gaia brings better constraints on the vertical structure due to its dense all-sky coverage. We performed a model fitting using them separately abd then using both simultaneously to account for their respective selection functions.

We find a good overall agreement between the profiles traced by both surveys, and with previous studies. We also noted interesting differences stemming from the different spatial coverage of the surveys, notably some difference between the APOGEE and Gaia results in the solar neighbourhood. The explanation for this may be that the density model we employed was too smooth to simultaneously account for the details of the density distribution on a small and large scale.

The Gaia data allowed us to study the density profile of α-poor mono-abundance populations within 2 kpc. We confirmed the steepening of the density gradient with metallicity and the continuous decrease in scale height reported in studies based on APOGEE data.

We make the APOGEE selection function used in this study available through the GaiaUnlimited Python package¹⁶. The procedure we employed to construct the selection function of the Gaia subsample of red clump stars is also explained in detail in a Jupyter notebook¹⁷, which contains the code that produced Figs. 7–10.

Acknowledgments

T.C.G. thanks David Hogg and Adrian Price-Whelan for their feedback on this manuscript. 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, particularly those participating in the Gaia Multilateral Agreement. This work results from the Gaia Unlimited project, which has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 101004110. The Gaia Unlimited project was started at the 2019 Santa Barbara Gaia Sprint, hosted by the Kavli Institute for Theoretical Physics at the University of California, Santa Barbara. A. R. C. is partly supported by the Australian Research Council through a Discovery Early Career Researcher Award (DE190100656) and Discovery Project DP210100018. This work made use of the https://www.overleaf.com/ and the Python packages astropy (Astropy Collaboration 2018), scipy (Virtanen et al. 2020), numpy (Harris et al. 2020), healpy (http://healpix.sourceforge.net) (Górski et al. 2005; Zonca et al. 2019), and matplotlib (Hunter 2007). T.C.G. acknowledges an extensive use of TOPCAT (Taylor 2005) and Jupyter notebooks (Kluyver et al. 2016).

References

Appendix A: Selection on renormalised unit weight error (ruwe)

Lindegren et al. (2021b) introduced the renormalised unit weight error (ruwe) as a diagnostic to assess the reliability of the astrometric solution for individual sources. A maximum value of 1.4 was later proposed as a simple filter to reject ill-behaved sources and inconsistent astrometric fits caused by unresolved binaries. Rejecting sources with ruwe> 1.4 has become a common first step when selecting data from the Gaia archive (e.g. Fabricius et al. 2021; Maíz Apellániz et al. 2021; Gaia Collaboration 2023; Cruz Reyes & Anderson 2023; Medan & Lépine 2023; Quintana et al. 2023; Qin et al. 2023; Alfonso & García-Varela 2023; Lee et al. 2023; Ratzenböck et al. 2023; Gontcharov et al. 2023; Liao et al. 2023), but the selection biases introduced by this filtering have not been investigated to our knowledge.

The quantity ruwe was constructed so that its typical value was 1 across the whole sky and at all magnitudes. In regions of the sky that received many Gaia visits (usually more than 120), sources in the magnitude range G∼11 to 12.5 have significantly higher ruwe values (Fig. A.1), and removing sources with values above 1.4 creates an artificial dearth of sources, especially in the range G = 11.5 to 12 (Fig. A.2). We note that Golovin et al. (2023) reported a very similar effect over the same magnitude range (shown in their Fig. A.3) in their catalogue of nearby stars (ϖ > 40 mas).

Fig. A.1. Effect of magnitude and number of matched transits on ruwe. Top: Number of matched transits for sources with G < 13 over a large patch of sky, showing the imprint of the Gaia scanning law. Bottom: ruwe value for sources in this region in three different magnitude ranges. In areas with more than ∼120 Gaia visits, removing sources with ruwe > 1.4 effectively discards almost all sources with 11 < G < 12.

Fig. A.2. Effect of magnitude and number of matched transits on ruwe. Top: Median ruwe value in bins of 0.3 mag for sources with fewer and with more than 125 matched transits. The shaded area corresponds to the 25th to 75th percentile range. Bottom: G-magnitude distribution of the sources with ruwe < 1.4, and with fewer or more than 125 matched transits.

The reason for this quirk happening within a very narrow magnitude range is currently not understood, but understanding the origin of this effect is not necessary in order to account for it. We note that it coincides with the magnitude limit (G = 11.5) between the window classes 0A and 0B of the Gaia video-processing units, and it might be related to calibration issues across this boundary. The Gaia technical note GAIA-C3-TN-LU-LL-124-01 describes how ruwe is computed from the unit weight error uwe, which is itself obtained from non-linear relations to various indicators of astrometric quality. The calculation of ruwe accounts for the fact that the residuals behave differently in different regimes of colour and magnitude, but it does not account for the number of scans, which might introduce this glitch for sources with a unfortunate combination of G magnitude and transit numbers.

Appendix B: Profile fit using only the most complete bins

In Sect. 2.3.3 we commented on the metallicity dependence of the RC-GSP-Spec selection function and proposed a simple approximation to account for this dependence. The G-magnitude distributions shown in the bottom panel of Fig. 10 suggest that the drop in the selection function from ∼80% to 0% completeness is steeper at a given metallicity than when averaged over the entire sample. In order to sidestep this uncertainty on the completeness profile of individual metallicity ranges, we attempted to only perform the fitting of the density profile over the volume of space where each sample is considered to be over 80% complete. The results are shown in Fig. B.1.

Fig. B.1. Density profile for RC-GSP-Spec stars built upon any (top) and high completeness levels (80%; bottom). The top row is identical to the middle row of Fig. 13. The high-completeness requirement greatly reduces the coverage in galactocentric radius R, but the profiles remain statistically similar.

This requirement greatly reduced the coverage in galactocentric radius R, especially for the low-metallicity sample. The density profiles are shown in the bottom left panel of Fig. B.1. They essentially match the trends shown in Fig. 11, where we see a flat (at low metallicity) or monotonic increase towards the Galactic centre (at high metallicity) density in the solar neighbourhood.

The reduced coverage in R leads to very large uncertainties on the flaring scale length R_flare⁻¹. We remark that because the distance range covered by this restricted data is still several times larger than the vertical scale height, we are still able to accurately recover the values of h_Z, ⊙∼450 pc for the [M/H] < -0.3 dex population, ∼300 pc for -0.3 < [M/H] < -0.15, and ∼250 pc for the more metal-rich stars.

Appendix C: Vertical density distribution of the metal-rich population

In the analytic model presented in Sect. 3.2, the vertical density distribution ζ(Z|R) follows an exponential profile with scale length varying with galactocentric distance R,

$$\begin{array}{c}\hfill h_Z^{-1}(R)=h_{Z,\odot }^{-1}exp(R_{\mathrm{flare}}^{-1}[R-R_0]).\end{array}$$(C.1)

With this definition, a negative value of $R_{\mathrm{flare}}^{-1}$ causes $h_Z^{-1}$ to decrease with R, causing the scale height h_Z to increase. This effect is illustrated in Fig. C.1 for three chosen values of $R_{\mathrm{flare}}^{-1}$. The density profiles we obtain from Gaia data (Fig. 13) show a strong degree of flaring, with the most metal-rich populations reaching, $R_{\mathrm{flare}}^{-1}\sim$0.2 kpc⁻¹. This strong flaring is illustrated (qualitatively) in Fig. C.2, where the right panel shows that the vertical density drop is steeper in the inner than in the outer disc.

Fig. C.1. Vertical density distribution ζ for a flared exponential disc with hZ, ⊙ = 150 pc as a function of R and Z for three different flaring scale lengths. For the purpose of this plot, the densities in all three panels are normalised so that ζ(Z = 0|R⊙) = 1.

Fig. C.2. Left: Number of observed RC stars with [M/H]> 0.3 dex computed in (X, Y, Z) bins of 200 pc, corrected for incompleteness (where the completeness is above 70% only), and displayed as a function of R and Z. The Sun is at R = 8.15 kpc. Right: Log-density as a function of distance from the Galactic mid-plane |Z| computed in six ranges of R (colour-coded as shown above the left panel). The flattening of the density profile with increasing R is interpreted by our model as a sign of flaring. We reiterate that the observed counts should only be corrected for completeness for data exploration, and the model fitting procedure of Sect. 3.2 is instead performed by applying the selection function to the model.

Using APOGEE DR12 data, Bovy et al. (2016) found that the α-poor disc only exhibits moderate flaring with $R_{\mathrm{flare}}^{-1}$=-0.1 at all metallicities, but they pointed out that the footprint of the APOGEE survey is dominated by low-latitude fields that provide little leverage on the vertical density profile. The dense coverage offered by Gaia in the solar neighbourhood provides much stronger constraints.

Appendix D: Best-fit parameter table

All Tables

All Figures

	Fig. 1. Spatial distribution of the SDSS/APOGEE red clump stars and Gaia RC-GSP-Spec stars (see Sect. 2.2.1). Left: location projected on the Galactic plane. Middle: Galactic altitude Z vs. galactocentric distance R. Right: distribution of galactocentric distance R.
In the text

Fig. 2. Binning in intrinsic colour-magnitude space used by APOGEE (left) and APOGEE-2 (right) to select their Main Red Star Sample targets from the 2MASS catalogue. The exact magnitude limits between the three cohorts (short, medium, and long) depend on the field and are provided by the APOGEE collaboration in auxiliary files. The selection is computed as the ratio of successfully observed targets to available sources in each of the three (in APOGEE) or six (in APOGEE-2) regions.

In the text

	Fig. 3. Effective selection function for the APOGEE DR17 red clump stars shown in the Galactic (X, Y) plane at b = 0. Each pencil beam corresponds to a different APOGEE pointing. The difference in bright magnitude limit between APOGEE-1 (H = 7 mag) and APOGEE-2 (H = 10 mag) is visible in the fourth quadrant (lower left corner).
In the text

	Fig. 4. Definition of the Gaia red clump star sample (see Sect. 2.2.1). The shaded area between the dashed lines represents our selection from the Gaia DR3 photometry described in Sect. 2.2.1, based upon the APOGEE DR17 red clump candidates indicated by the black dots. The grey background indicates the density of all sources in the “gradient sample” of Gaia Collaboration (2023) for reference.
In the text

Fig. 5. Abundance-scaling relations of Gaia DR3 GSP-Spec and APOGEE DR17 ASCAP abundances for a sample of 2450 red clump stars in common. Left: offset between Gaia DR3 mh_gspspec metallicities and APOGEE DR17 ASCAP [M/H]. Right: offset between [α/Fe] (from GSP-Spec) and [Ca/Fe] (from APOGEE). In both panels, the dashed line indicates the correction we applied to calibrate Gaia on the APOGEE scale, given in Eqs. (4) and (5).

In the text

Fig. 6. Red clump stars in chemical abundance space. Left: distribution of the APOGEE RC stars in the [α/Fe]-metallicity plane. The α-rich stars are highlighted manually. Middle: same stars in the [Ca/Fe]-metallicity plane. The [Ca/Fe] abundances do not have the same discriminating power as [α/Fe], which accounts for O, Mg, Si, S, and Ca for APOGEE. Right: distribution of the Gaia RC-GSP-Spec stars in the [α/Fe]-metallicity plane. The contours enclose 99%, 95%, 90%, 80%, and 50% of the sample. The Gaiaα-abundances are dominated by Ca lines and do not allow us to reliably distinguish the α-rich from the α-poor disc stars.

In the text

	Fig. 7. Completeness of our RC-GSP-Spec sample with respect to all Gaia DR3 red clump stars computed in bins of 0.5 mag (black line). The grey lines, computed independently in the 12 HEALPixels at level 0, show the same quantity.
In the text

	Fig. 8. Selection function of the RC-GSP-Spec sample, constructed from count ratios in spatial (healpix level 3) and magnitude bins. The maps are in Galactic coordinates, centred on (ℓ,b) = (0, 0), with ℓ increasing to the left.
In the text

Fig. 9. Effect of distance and extinction on the sample completeness. Top left: cumulative extinction at Z = 0, computed from the Lallement et al. (2022) dust map. Top right: predicted apparent G magnitude for a red clump star, assuming an intrinsic MG = 0.5 mag. Bottom left: selection function SGaia chem(X, Y, Z) plotted for b = 0. Bottom right: XY distribution of red clump stars with |Z|< 15 pc in our RC-GSP-Spec sample (Sect. 2.2.2). The Galactic centre is to the right in every panel.

In the text

Fig. 10. Effect of metallicity on the selection function. Top left: distribution of metal-poor RC-GSP-Spec stars projected the Galactic (X, Y) plane for stars with |Z|< 200 pc alone. The Sun is at (X, Y) = (0, 0), and the Galactic centre is to the right. Top right: same for the most metal-rich stars in the sample. Metal-rich stars contained in the sample reach larger distances from the Sun. Bottom: G magnitude distribution in the entire sample (black) and in four ranges of metallicity. The offset ΔG used in Eq. (9) is indicated for the most metal-poor metallicity group.

In the text

Fig. 11. Completeness-corrected number of α-poor RC-GSP-Spec stars in bins of Galactic (X, Y) bins of 200 × 200 × 200 pc near the Galactic mid-plane (−200 < Z < 200 pc) as a function of galactocentric distance in parsecs. We selected α-poor stars as [α/Fe]< 0.1 dex. A dot represents the corrected density for regions in which the expected completeness is greater than 20%. The top and bottom panels correspond to the data displayed in the top left and top right panels of Fig. 10, respectively. We stress that the observed counts should only be divided by completeness for qualitative visualisation and data exploration purposes. Section 3.3 implements a forward-modelling approach to apply the selection function to a density model.

In the text

	Fig. 12. Effective selection function (shown at b = 0°) for the combination of our Gaia and APOGEE samples obtained from the effective selection functions displayed in Figs. 3 and 9 and Eq. (14). The completeness reaches over 80% in the solar neighbourhood dominated by the Gaia sample.
In the text

Fig. 13. Density profiles of red clump stars from APOGEE, Gaia RC-GSP-Spec, and their combined data sets (from top to bottom). The left panels show the radial density profiles Σ(R) obtained in six ranges of metallicity in the α-poor disc ([α/Fe]< 0.1 dex). The profiles are shifted by an arbitrary offset to be visually separated. The panels on the right represent the flaring scale length and vertical scale height predicted by our model for each metallicity bin. The shaded areas (left) and error bars (right) represent the 95% confidence interval. The Gaia RC-GSP-Spec star sample covers a smaller range in R, especially at low metallicity. We list the best-fit parameters in Table D.1.

In the text

Fig. A.1. Effect of magnitude and number of matched transits on ruwe. Top: Number of matched transits for sources with G < 13 over a large patch of sky, showing the imprint of the Gaia scanning law. Bottom: ruwe value for sources in this region in three different magnitude ranges. In areas with more than ∼120 Gaia visits, removing sources with ruwe > 1.4 effectively discards almost all sources with 11 < G < 12.

In the text

	Fig. A.2. Effect of magnitude and number of matched transits on ruwe. Top: Median ruwe value in bins of 0.3 mag for sources with fewer and with more than 125 matched transits. The shaded area corresponds to the 25th to 75th percentile range. Bottom: G-magnitude distribution of the sources with ruwe < 1.4, and with fewer or more than 125 matched transits.
In the text

	Fig. B.1. Density profile for RC-GSP-Spec stars built upon any (top) and high completeness levels (80%; bottom). The top row is identical to the middle row of Fig. 13. The high-completeness requirement greatly reduces the coverage in galactocentric radius R, but the profiles remain statistically similar.
In the text

	Fig. C.1. Vertical density distribution ζ for a flared exponential disc with hZ, ⊙ = 150 pc as a function of R and Z for three different flaring scale lengths. For the purpose of this plot, the densities in all three panels are normalised so that ζ(Z = 0|R⊙) = 1.
In the text

Fig. C.2. Left: Number of observed RC stars with [M/H]> 0.3 dex computed in (X, Y, Z) bins of 200 pc, corrected for incompleteness (where the completeness is above 70% only), and displayed as a function of R and Z. The Sun is at R = 8.15 kpc. Right: Log-density as a function of distance from the Galactic mid-plane |Z| computed in six ranges of R (colour-coded as shown above the left panel). The flattening of the density profile with increasing R is interpreted by our model as a sign of flaring. We reiterate that the observed counts should only be corrected for completeness for data exploration, and the model fitting procedure of Sect. 3.2 is instead performed by applying the selection function to the model.

In the text