The Nearby Evolved Stars Survey II: Constructing a volume-limited sample and first results from the James Clerk Maxwell Telescope

Evolved stars are typically analysed using two tracers: dust and molecular lines. NESS aims to build a coherent picture of mass-loss rates, by comparing these two mass tracers.

Dust-based mass-loss rates scale with the fraction of stellar light reprocessed from the optical to the mid-infrared (and derive the optical depth). This traces warm dust ($\sim$50–1500 K) close to the star. The scaling factor depends on the distance and luminosity of the star, and requires assumptions on the expansion velocity (e.g., from CO lines or OH masers) of the wind and its radial profile, how strong the dynamical coupling between dust and gas is, the wind’s geometry and clumpiness, the dust-to-gas ratio, and various geometrical, chemical and mineralogical properties of the dust grains (e.g. van Loon, 2006; McDonald et al., 2011c; Jones et al., 2014), as these factors all influence how stellar radiation is reprocessed by the outflowing dust.

By contrast, gas mass-loss rates from radio or (sub)millimetre wavelength observations can directly yield the wind’s velocity and density structure. This simplifies the scaling factors, requiring only the unknown but modellable excitation conditions of the observed molecule and its abundance in the wind. Of the abundant molecules, CO is expected to be the most stable, with a consistent abundance of $\sim$10^-4 with respect to molecular hydrogen across a range of mass-loss rates and chemistry (e.g. Schöier et al., 2002; Ramstedt et al., 2008; De Beck et al., 2010; Danilovich et al., 2015), though it will be lower in metal-poor environments (e.g. Leroy et al., 2011).

The gaseous envelope ejected over the last few millennia is cold (kinetic temperatures are a few $\times$ 10 K). With energies of $E_{u}/k=5-50$ K, (Pickett et al., 1998), low-$J$ pure-rotational transitions of CO are sensitive to the molecular content of the entire cold envelope, particularly the outer regions (Kemper et al., 2003). Many of these transitions fall in the (sub-)mm atmospheric windows, facilitating ground-based observation. Timescales of $\sim$10 000 years can be probed (Mamon et al., 1988; Groenewegen, 2017; Saberi et al., 2019), and time variability in mass-loss rates can sometimes be uncovered (e.g. Olofsson et al., 1990; Decin et al., 2006; Guélin et al., 2018; Maercker et al., 2010; Dharmawardena et al., 2019).

To generate dust-based estimates comparable to $\dot{M}$ from CO lines, we must therefore probe the bulk of the dust mass, which is likely to be at similarly cold temperatures to the gas seen in low-excitation CO lines. The sub-mm continuum traces emission from the coldest dust, at similar temperatures as the low-$J$ CO lines.

CO-based mass-loss rates have successfully been exploited by a wide range of studies. However, these are all either limited in scope or exhibit biases that make comparative statistics difficult.

Many existing studies focus on the brightest AGB stars with the highest mass-loss rates. “Extreme” mass-losing stars often dominate the dust budget of galaxies (e.g. Riebel et al., 2012; Boyer et al., 2012; Srinivasan et al., 2016), and are thus observed to explore these rapid and catastrophic evolutionary phases, and to constrain the dust budget of galaxies. However, strong mass loss ($\dot{M}\gtrsim 10^{-6}$ M_⊙ yr^-1) only occurs during the final stages of the AGB evolution of higher-mass AGB stars, hence observing them misses the mass loss that occurs in lower-mass AGB stars, or earlier on in AGB evolution, frustrating efforts to quantify these final stages in evolutionary models and biasing our understanding of chemical enrichment. In particular, very few AGB stars have been observed in the early phases of their dusty wind production, or before strong dust-production starts (Groenewegen, 2014; Kervella et al., 2016; McDonald et al., 2016; McDonald et al., 2018), despite these being both numerically the most numerous AGB stars in most galaxies (e.g. Boyer et al., 2015a; Boyer et al., 2015b) and the terminal phase of the AGB for low-mass stars (e.g. McDonald & Zijlstra, 2015a).

Other samples offer volume completeness, but only of specific types of targets, such as only carbon stars (Olofsson et al., 1993) or S-type stars (where C/O $\sim$ 1; e.g., Sahai & Liechti, 1995; Groenewegen & de Jong, 1998; Ramstedt et al., 2009). These samples can be compared (e.g. Ramstedt et al., 2009), but the different relative sizes, completenesses and volumes of these samples make it difficult to draw robust comparative statistics.

Most of these studies are relatively modest, comprising typically of tens of stars (e.g. Neri et al., 1998; Kemper et al., 2003; Teyssier et al., 2006; McDonald et al., 2018; Dharmawardena et al., 2018), and studies with larger samples of up to $\sim$100 stars (e.g. Kerschbaum & Olofsson, 1999; Olofsson et al., 2002; Loup et al., 1993; Kastner et al., 1993; Zuckerman & Dyck, 1986a; Zuckerman et al., 1986; Zuckerman & Dyck, 1986b, 1989) typically have tens of detections. Nyman et al. (1992) detected some 160 out of $\sim$500 IRAS-identified evolved stars in CO $J$=1–0, but did not target a complete sample or full sky coverage, and had relatively low detection rates due to a lack of sensitivity; while Nyman et al. (1992) had $T_{\rm sys}\sim 400-1000\,$K, the JCMT routinely observes with $T_{\rm sys}\lesssim 100\,$K enabling a factor of 3 or more improvement in sensitivity. This made it difficult for previous studies to accurately extract trends across the complex parameter space of AGB evolution, and to be robust against statistical outliers. This motivates the need for a volume-complete survey of nearby AGB stars, of sufficiently large scale to extract relationships that can improve comparative and evolutionary studies.

The optimal observational setup depends strongly on the spatial scale to be recovered. While the CO $J$=1–0 line has historically been the easiest to observe (e.g. Nyman et al., 1992), the large beam sizes and resultant contamination by cold interstellar gas and dust mean that the ¹²CO $J$=2–1 and 3–2 transitions (230.538 and 345.796 GHz, respectively) and sometimes the related ¹³CO transitions (220.399 and 330.588 GHz, respectively) are generally chosen to observe AGB stars (e.g. Fig. 1). Higher spatial resolution can be obtained with interferometers, but pressure on these instruments prohibits very large programmes at high sensitivity, and emission can be filtered out at large scales due to incomplete filling of the $uv$ plane. Large single-dish telescopes provide high sensitivity without filtering out emission. In particular, the James Clerk Maxwell Telescope (JCMT, 15 metre diameter, with resolutions of 21^′′/ 14^′′ at 230/345 GHz, respectively) and Atacama Pathfinder EXperiment (APEX, 12 metre diameter, 25^′′/ 17^′′), equipped with multi-beam detector arrays, provide high mapping speeds capable of observing hundreds of stars, efficiently revealing historic mass loss through extended emission, while the observing methods employed ensure that all the flux is captured.

By making this selection, the NESS project is performing a wide survey of a large number of AGB stars, $\sim$850, albeit at the cost of recovering limited spatial information about the envelope of each one. This places NESS as the widest project, with strong synergies to other, ongoing large observing programmes. At medium spatial resolutions (5$\aas@@fstack{\prime\prime}$5/3–4^′′) with only slight expense of $uv$ coverage (maximum recoverable scale, MRS $\sim$ 25^′′/18^′′) and sample size ($\sim$180 stars), is the Morita (Atacama Compact Array) survey DEtermining Accurate mass-loss rates for THermally pulsing AGB STARs (DEATHSTAR; Ramstedt et al. 2020³³3 http://www.astro.uu.se/deathstar.). Finally, at the highest resolutions (0.05–0.025^′′/0.03–0.016^′′) of the main Atacama Large Millimeter/submillimeter Array (ALMA), is the project “ALMA Tracing the Origins of Molecules forming dust In Oxygen-rich M-type stars” (ATOMIUM; Decin et al. 2020⁴⁴4https://fys.kuleuven.be/ster/research-projects/aerosol/atomium.), which includes 14 nearby AGB stars and red supergiants.

We require a large, volume-complete sample of Galactic evolved stars covering the largest possible range of $\dot{M}$, homogeneously observed in both CO lines and continuum, aiming to minimise biases. Defining the sample volume depends critically upon our ability to determine distances. Unfortunately, while distance estimates are rapidly improving, distances to Galactic evolved stars remain poor. We must make assumptions, and revisit the volume-completeness of the NESS survey in future, ultimately allowing us to complete our survey aims.

Where observed, maser parallax provides the most exquisite precision, but only evolved stars with the highest $\dot{M}$ exhibit the masing transitions required. We have included maser distances as our preferred choice where available (e.g. Orosz et al., 2017).

Convective motions and dust obscuration can move the astrometric centres of AGB stars in the optical, introducing significant astrometric noise to optical parallax measurements (e.g. Chiavassa et al., 2018). Similarly, AGB stars are poorly captured by the prior model used in Bailer-Jones et al. (2018). Consequently, data from Gaia DR2 (Gaia Collaboration et al., 2018) were not used (Gaia EDR3 is addressed below), as they rely solely on short-timescale Gaia data (see McDonald et al., 2018 for further discussion). Instead, we adopt the parallaxes from the Tycho–Gaia Astrometric Solution (TGAS) of Gaia Data Release 1 (DR1) (Gaia Collaboration et al., 2016) as our next preference, followed by the Hipparcos parallaxes (van Leeuwen, 2007), provided their fractional uncertainty is small ($\sigma_{\varpi}/\varpi<0.25$).

Other distance-determination methods typically only apply to a small subset of AGB stars; e.g., period–luminosity relationships can only be used if the pulsation mode and mean magnitude are known (e.g. Uttenthaler, 2013; Huang et al., 2018; Yuan et al., 2018; Goldman et al., 2019), as overlap between different sequences and scatter from single-epoch photometry can confuse measurements.

Instead, we take a statistical approach to missing distances, using the luminosity distribution of evolved-star candidates in the Large Magellanic Cloud (LMC), taken from the Spitzer Surveying the Agents of Galactic Evolution (SAGE) database (Meixner et al., 2006), accounting for the geometry and thickness of the LMC. The luminosity of each star is determined by integrating its SED from the optical to the far-infrared and the luminosities of all sources are combined to give a probability distribution, from which we extract the median and the central 68 per cent. This gives a “typical” luminosity and uncertainty of $6200^{+2800}_{-3900}$ L_⊙ for the LMC population.

Assuming that AGB stars in the solar neighbourhood have sufficiently similar luminosities to stars in the LMC, we can apply this same luminosity distribution to derive a probabilistic distribution of distances to individual Galactic sources. A galaxy’s luminosity distribution is mainly set by its star-formation history, as the final AGB luminosity depends on initial mass. Secondary effects include variation in mass loss, stellar evolution and the initial-mass function with metallicity. To test this assumption, we compare the AGB luminosity functions of the LMC and SMC to each other, and to the solar neighbourhood. No significant difference is seen (e.g., Fig. 15 in Srinivasan et al. 2016). The spread in median luminosity, 17 per cent (Boyer et al., 2015a; McDonald et al., 2017), indicates a systematic galaxy-to-galaxy error of 8.5 per cent in this distance-based method. In comparison, the intrinsic width of the luminosity distribution corresponds to a distance uncertainty of $\approx$25 per cent. Consequently, we expect any galaxy-to-galaxy differences to produce much smaller systematic uncertainties than the random uncertainties inherent in the measurement itself (see also Section 2.3.7). We therefore use this LMC luminosity distribution and photometry from infrared all-sky surveys to determine distances to Galactic sources where parallaxes from optical or maser astrometry are unavailable or imprecise.

Such a luminosity-based, volume-limited survey is subject to forms of Malmquist bias (Malmquist, 1925), whereby many sources from larger distances scatter into our distribution, while only a few sources from smaller distances scatter out. The $\dot{M}$ cutoffs that define our tiers (below) further re-enforce this bias, ensuring distant objects scattering into the sample are assigned anomalously high mass-loss rates, hence are placed into more extreme categories of mass loss. Furthermore, parallax errors generate a Lutz–Kelker bias (Lutz & Kelker, 1973), preferentially scattering sources into our sample for the same reasons. Hence, while the distances to our sources undergo frequent revision, we expect our samples to be relatively robust against obtaining larger datasets in future, though ultimately some sources may no longer meet our distance criteria. A small Malmquist bias in the LMC sample cancels out some of these effects, but AGB stars in the SAGE data are bright enough that this is largely negligible. In the following sections, we attempt to reduce these biases in our data.

The Infrared Astronomical Satellite (IRAS) point-source catalogue (Beichman et al., 1988) is our fiducial reference for source selection. IRAS photometry were matched against the 2MASS catalogue (Cutri et al., 2003) to find the nearest neighbour within 30$\arcsec$, providing homogeneous near- and mid-infrared photometry for all sources. The photometry were then integrated numerically to estimate the bolometric fluxes of the sources; we then compute the distribution of distances as that required to scale the bolometric flux to the LMC luminosity distribution. The 50th centile luminosity (6200 L_⊙) is chosen as representative, so that we have a point-distance estimate for each source following the constraints described below in Sects. 2.3.4 & 2.3.5. In order to restrict their sample to stars with mass-loss rates $>2\times 10^{-6}$ M_⊙ yr^-1, Jura & Kleinmann 1989 select sources with a minimum IRAS 60 $\mu$m flux of 10 Jy at 1 kpc. In their work as in ours, the mass-loss rates are computed from the dust-production rates assuming a gas-to-dust ratio of about 200. Their limit therefore corresponds to a DPR of $1\times 10^{-8}$ M_⊙ yr^-1. Extending this criterion to 2 kpc, we select sources brighter than 2.5 Jy, thus reducing the full point-source catalogue to $\sim$60,000 sources. This removes most extra-galactic sources and less-evolved stars, and is set sufficiently high that it avoids strong effects from Malmquist bias. We consider the effects of this choice on completeness in Section 2.3.8.

To this sample, we add evolved stars within 300 pc from McDonald et al. (2012) and McDonald et al. (2017), which would otherwise have been excluded by the IRAS 60 $\mu$m flux cut mentioned above. Stars were added if they have $T_{\rm eff}<5500$ K and $L>700$ L_⊙, and if their Hipparcos or TGAS parallaxes have fractional uncertainties of $\sigma_{\varpi}/\varpi<0.25$ (560 sources meet these criteria). This cutoff matches the expected accuracy of our luminosity-based distances and avoids strong effects from Lutz–Kelker bias.

Since $\sim$60,000 sources are beyond our capability to observe, we design a tiered system to observe representative samples of them, selected from a plane of DPR versus distance, using first-order estimates of the DPR outlined in Sect. 2.3.3, and cuts outlined in Sect. 2.3.4. This results in a set of 2277 potential targets, from which sources with simbad⁵⁵5http://simbad.u-strasbg.fr/ classifications that are inconsistent with dusty evolved stars are removed⁶⁶6To measure mass return, we are primarily interested in sources with non-negligible ongoing or past mass loss, i.e., for low- and intermediate-mass stars, those on the early-AGB up to and including planetary nebulae (PNe), and for high-mass stars red supergiants (RSGs), yellow hypergiants (YHGs), B[e] supergiants (B[e]SGs), luminous blue variables (LBVs) and Wolf-Rayet stars (WRs).. Of the 852 remaining sources in our final sample (defined below), only nine have maser distances available, and seven and 193 respectively have TGAS and Hipparcos measurements with uncertainties $<25$ per cent. All remaining sources use luminosity distances.

We fit the matched photometry in the 2MASS $J$, $H$, and $K_{s}$ bands and the IRAS 12, 25 and 60 $\mu$m bands with models from the Grid of Red supergiant and AGB ModelS (GRAMS; Sargent et al. 2011; Srinivasan et al. 2011) and extract dust-production rates ($\dot{M}_{\rm dust}$) for the entire sample. $\dot{M}_{\rm dust}$ is thus used as a criterion to define a tiered survey, limited in distance and $\dot{M}_{\rm dust}$ as described in Sect. 2.3.4. Our source lists are additionally divided into two groups, a large set to be observed in spatially unresolved modes (the “staring” sample) and a smaller group to be mapped in detail (the “mapping” sample).

GRAMS is tailored to the LMC, and the assumed dust properties reflect its metal-paucity (i.e. the relative abundance of different dust species may be marginally different in Galactic sources; e.g., Srinivasan et al., 2010). However, for the broadband photometry used here, this should provide sufficient precision to perform sample selection. Under certain assumptions of velocity scaling, dust opacity and dust-to-gas ratio, $\dot{M}_{\rm dust}$ should provide a close proxy to total mass-loss rate, $\dot{M}$: an assumption that NESS will test (Section 4). While GRAMS also returns a best-fit chemical classification (O- or C-rich), the seven bands of photometry used in this paper are insufficient to accurately constrain the dust chemistry. Given that M-stars are vastly more numerous in the Milky Way, we therefore assume that each source is oxygen-rich unless there is spectroscopic confirmation of C-rich dust chemistry from either the IRAS LRS or ISO SWS spectra, with the ISO SWS classification from Kraemer et al. (2002) taking precedence over IRAS LRS classifications from Kwok et al. (1997). This assumption applies only to the determination of initial DPRs and is necessary because the choice of dust chemistry can make a significant difference to the derived DPR. We demarcate these unclassified sources separately below.

Using the distances and dust-derived $\dot{M}$ derived above, we define five sub-samples of sources to be observed (Fig. 2). These samples are defined in terms of increasing distance and $\dot{M}_{\rm dust}$.

• •

  Tier 0 (or “very low” DPR sources) is a special addition of ten sources with $L\geq 1600$ L_⊙, $d<250$ pc and $\delta\geq-30\deg$, without limit on $\dot{M}_{\rm dust}$, drawn from McDonald et al. (2012) and McDonald et al. (2017). This tier explores mass loss from bright red giant branch (RGB) and AGB stars not producing dust. Nine more sources at $\delta\leq-30\deg$ meet these criteria, but have not yet been scheduled for observation by NESS. One (SX Pav) has been published already (McDonald et al., 2018).

• •

  Tier 1 ( “low”; 105 sources) includes all sources with $\dot{M}_{\rm dust}<10^{-10}$ M_⊙ yr^-1 at $d<300$ pc, and samples the AGB stars with the lowest $\dot{M}_{\rm dust}$. Sources were only included if there was a 3$\sigma$ dust excess in the GRAMS models, or if the source has an infrared excess in McDonald et al. (2012) or McDonald et al. (2017).

• •

  Tier 2 ( “intermediate”; 222 sources) contains all sources with $10^{-10}\leq\dot{M}_{\rm dust}<3\times 10^{-9}$ M_⊙ yr^-1 and $d<600$ pc, excluding the Galactic plane ($\left|b\right|<1.5$) for $d>400\,pc$;

• •

  Tier 3 ( “high”; 324 sources) consists of all sources with $3\times 10^{-9}\leq\dot{M}_{\rm dust}<10^{-7}$ M_⊙ yr^-1 and $d<1200$ pc, excluding the Galactic plane for $d>800$ pc; and

• •

  Tier 4 ( “extreme”; 182 sources) comprises all sources with $\dot{M}_{\rm dust}\geq 10^{-7}$ M_⊙ yr^-1 and $d<3000$ pc, excluding the Galactic plane for $d>2000\,pc$.

These tier limits result in each of Tiers 1–4 containing enough objects ($>$100) to adequately establish the typical range of properties of AGB winds within that tier, while the limit of $L\geq 1600$ L_⊙ approximates the typical luminosity at which dust-production begins for the lowest-mass stars (e.g. McDonald et al., 2011c; McDonald et al., 2011b; McDonald et al., 2014). Sources in Tier 4 are poorly classified in the literature, so are more likely to be contaminants. The NESS survey aims to improve on these objects’ classifications. Since the original selection, updated distances for some sources have shifted their locations in Fig. 2.

Of the entire sample of 852 sources, Table 1 gives the first ten entries, sorted by right ascension, to demonstrate the information included and its format. For each source, we report the distance used, the source of this distance, the GRAMS best-fit $\dot{M}$, and classifications based on IRAS LRS or ISO SWS spectra when available. The S-star classifications are from Sloan & Price (1998), Yang et al. (2007), and Hony et al. (2009), with the Hony et al. (2009) classifications taking precedence. The full interactive table, available at https://evolvedstars.space, represents the initial state of the NESS catalogue. Future data releases will attach the stellar parameters assumed and NESS data products to each source. Figure 3 illustrates the sky distribution.

We find that the CO studies discussed above in Sect. 2.1.2 detect only 204 of the 852 NESS targets (less than 25%) in either or both of the CO $J$=2–1 and $J$=3–2 transitions. Figure 4 shows that the NESS survey will extend these observations to the sources with least dust content, which are almost completely missing in previous studies. Furthermore, NESS will provide a four-to-five-fold improvement in the statistics at DPR $\geq 5\times 10^{-9}$ M_⊙ yr^-1 and at distances greater than 600 pc.

In addition, 46 of the nearer sources in the staring sample with significant dust production were selected for mapping on arcminute scales (black circles in Fig. 2). These comprise stars outside the Galactic Plane ($|b|\geq 1.5$ deg) with $\dot{M}_{\rm dust}>5\times 10^{-10}$ M_⊙ yr^-1 and $d<340$ pc, which are bright and whose CO envelopes have a large predicted angular size ($\gtrsim 20\arcsec$, based on Mamon et al. 1988). Although the GRAMS best-fit prediction for U Ant was lower than the cutoff, it is known to have an extended envelope in the FIR (e.g. Kerschbaum et al., 2010), such that GRAMS underestimates the DPR. It has been published separately as Dharmawardena et al. (2019). Seven sources satisfying these criteria were not selected into the mapping sample; they will be folded in as part of future proposals to complete this set.

Many stars in our sample have existing JCMT or APEX observations of the (2–1) and (3–2) transitions of ¹²CO and ¹³CO, or existing JCMT continuum data from SCUBA-2, and we use these archival data rather than re-observing these sources. Data from JCMT/SCUBA and APEX (both LABOCA and SABOCA) at the same wavelengths are not included because the fields of view are much smaller.

Gaia Early Data Release 3 (EDR3) occurred during the refereeing process of this paper. Its 34-month timespan (comparable to that of individual Hipparcos stars) significantly reduces astrometric noise compared to DR2, improving parallaxes for nearby evolved stars. While this improves the distances of many of our sources, we here retain the distances used to define the membership of our catalogue tiers for self-consistency, and defer decisions on individual distances to a dedicated catalogue paper (McDonald et al., in prep.).

Meanwhile, we compare the Gaia EDR3 distances to those used to select our sources. To identify Gaia EDR3 cross-matches to our IRAS sources, we selected Gaia sources within 3^′′ of the simbad co-ordinates, provided they had red Gaia colours ($B_{P}-R_{P}>1$ mag; 771 sources cross-matched), and compute distances for these by naïvely inverting the parallax (see comments below).

As Gaia EDR3 is largely independent of our source data (only a few sources were selected using TGAS parallaxes), this comparison tests the accuracy of our distance-selection methods. Figure 5 shows a good correlation overall. The scatter comes from the combined errors in the Gaia EDR3 and adopted distance data. We assume they are uncorrelated, thus

A long tail of sources towards large Gaia distances, caused by parallaxes with large fractional errors, skews the sample from a normal distribution. Some of these stars could be RSGs. However, this long tail continues into negative parallaxes for 31 sources (3 per cent), four of which are statistically significant at the 4–6$\sigma$ level. A further 47 sources have a Gaia cross-match with no determined parallax. Many sources in this long tail have astrometric noise exceeding the parallax itself, showing photo-centric motion is still a significant issue for AGB stars in Gaia EDR3.

A few sources also scatter to lower Gaia distances than expected, mostly in Tier 2. Some may be RGB stars, which are not contained in the prior luminosity function from the LMC. However, their inclusion requires detectable infrared excess, which is generally not observed in RGB stars (e.g. McDonald et al., 2011a). As the distribution is markedly non-Gaussian at both under- and over-estimated distances, we investigate the central 68th centile of the distribution of the three $\sigma$ terms in Equation 1, rather than their standard deviation. This avoids detailed treatment of stars with poor, zero and negative parallaxes, as these are rare or non-existent within the central 68th centile.

Where the Hipparcos parallax is the adopted distance measure, $\sigma$[Gaia / Hipparcos] = 0.16. These are typically close stars (median Hipparcos distance = 209 pc), with low astrometric noise, and the quadrature addition of fractional parallax errors (Equation 1) is only slightly more than the value of 0.12 expected from the quoted astrometric uncertainties of the Hipparcos and Gaia surveys.

Beyond $\sim$400 pc, imprecise historical parallaxes become replaced by luminosity distances. In tier 2 (including four sources with luminosity distances from tier 1), $\sigma$[Gaia / adopted] = 0.36, on average, for a median adopted distance of 477 pc. Similar numbers for tiers 3 and 4 are 0.45 (849 pc) and 0.84 (2121 pc). If we first assume $\sigma_{\rm adopted}/d_{\rm adopted}=0.25$, then we can approximate the true uncertainty in Gaia EDR3 for tiers 2, 3 and 4 as $\sigma_{\rm Gaia}/d_{\rm Gaia}$ = 0.26, 0.38 and 0.80, respectively.

While $\sigma_{\rm adopted,lum}$ should be largely invariant with true distance, $\sigma_{\rm Gaia}$ should increase linearly with true distance. We can enforce this by setting $\sigma_{\rm adopted,lum}=0.316$. However, the extremity of our sources (thus their astrometric noise) does increase with distance, while their brightness decreases as they become more optically obscured, and they become more crowded as they concentrate closer to the Galactic plane (Section 4.1). Consequently, $\sigma_{\rm Gaia}$ will increase with distance, $\sigma_{\rm adopted,lum}=0.316$ can be treated as an upper limit, and we estimate that our adopted value ($\sigma_{\rm adopted,lum}=0.25$) is approximately correct.

We can also assume that the ratio of Gaia EDR3 distances to luminosity distances should have a median of unity within each tier⁷⁷7We have not performed the parallax bias correction of Lindegren et al. (2020). Of the 771 Gaia EDR3 cross-matches, only 183 have solution types or colours within the range where these corrections are valid. For these remaining sources, the correction requires an accurate and invariant colour and magnitude, hence a detailed treatment of the epoch photometry, to be precise. The corrections for these 183 sources based on their given colours and magnitudes are fairly small ($<$52 $\mu$as and $<$9.6 per cent in 95 per cent of cases) and less than the parallax error in all but five cases. Consequently, we ignore the correction.. Where we have adopted maser distances, the median $d_{\rm adopted,maser}/d_{\rm Gaia}=0.80$ with a standard error of 0.17. Where we have adopted parallaxes, $d_{\rm adopted,plx}/d_{\rm Gaia}=0.97\pm 0.03$. Finally, for our luminosity-based distances in tiers 2, 3 and 4, respectively, $d_{\rm adopted,lum}/d_{\rm Gaia}$ = 0.98 $\pm$ 0.07, 0.85 $\pm$ 0.07, 1.02 $\pm$ 0.26. While this generates a 2.0$\sigma$ outlier for tier 3, the sample overall is in statistical agreement with $d_{\rm adopted,lum}/d_{\rm Gaia}$ = 1. Hence we consider there to be no discernable systematic offset between our adopted distances and the true distance, despite the varying and competing effects of Malmquist and Lutz–Kelker biases in the Galactic and LMC samples, and the differing properties of the two galaxies.

Since a substantial part of the infrared emission of high-mass-loss-rate AGB stars is from the dust, the completeness of IRAS to AGB stars at a given distance is a function of dust-production rate. The 2.5 Jy cut (Section 2.3.2) should retain all AGB stars brighter than the tip of the RGB out to $\sim$500 pc, hence we can consider our sample of AGB stars complete to these distances for all tiers. It should further retain the brightest objects and those with any measurable mass loss out to much greater distances, and we have checked that the distance distribution of each tier beyond $\sim$500 pc approximates the $d^{2}$ distribution expected for the Galactic disc, allowing for stochastic variance and the missing sample from $\left|b\right|<1.5$. Since AGB stars are intrinsically bright in the infrared, we do not expect significant source confusion either at small distances, or at larger distances when $\left|b\right|<1.5$.

If the 60 $\mu$m flux were the only constraint, our limiting DPR would be proportional to the square of the distance. Tiers 4 and 3 would then be more than 95% complete, but Tier 2 would be incomplete at large distances (missing some sources with lower DPRs). Tier 1 would remain unaffected thanks to the inclusion of nearby sources from McDonald et al. (2012); McDonald et al. (2017).

However, the 60 $\mu$m flux is not linearly dependent on the DPR because, below a DPR of roughly $0.3-3\times 10^{-9}$ M_⊙ yr^-1, the photosphere becomes the dominant contributor to the far-IR emission. Since our sample is complete to naked (zero-DPR) stars at the tip of the RGB out to at least 500 pc, under the worst-case scenario the outer 100 pc of Tier 2 would become incomplete. However, since the minimum DPR in Tier 2 is non-zero, the combined emission from dust and photosphere are expected to render our sample in Tier 2 complete as well. It is worth noting that changes in gas-to-dust ratio impact these calculations: our sample becomes more complete as the ratio increases.

Deriving DPR and $\dot{M}$ from fitting optical–mid-infrared SEDs (Sect. 2.3.3) best traces material ejected from the star in the last few decades, as this hot and warm dust emits in the mid-infrared. However, it is relatively insensitive to variations in mass-loss rate or to reservoirs of material ejected in the past (Sect. 2.1.1).

Resolving the envelopes at longer wavelengths isolates emission from cold dust emitted millennia ago (e.g. van der Veen et al., 1995; Dehaes et al., 2007; Ladjal et al., 2010; Dharmawardena et al., 2018). Cox et al. (2012) notably found extended far-infrared (FIR) emission in 49 out of 78 evolved stars¹⁴¹⁴14Some of the 29 “non-detections” may have extended emission, as the authors were looking for specific shapes of the extended emission., mostly drawn from the Mass-loss of Evolved StarS (MESS) Herschel Guaranteed Time Key Programme (Groenewegen et al., 2011), confirming cold dust is commmonplace.

However, sensitivity to the coldest dust requires (sub-)mm observations. Ladjal et al. (2010) found extended emission in four out of nine sources at 870 $\mu$m using the APEX Large Bolometer Camera (LABOCA); and Dharmawardena et al. (2018) resolved sub-mm emission at one or more wavelengths in 14 out of 15 evolved stars with JCMT/SCUBA-2, with extended emission accounting for an average 40 per cent of the sub-mm flux. Furthermore, by deriving time-averaged mass-loss rates from the resolved emission, Dharmawardena et al. (2018) show that the dust mass in the resolved observations is typically a factor of a few larger than predicted by the GRAMS models, which assume a constant $\dot{M}$ out to a fixed radius. This has potential to significantly increase the contribution of AGB stars to the interstellar dust reservoir, and/or their interaction with it, and allows us to probe longer-timescale variations in mass loss.

The deep SCUBA-2 mapping observations of NESS will more than double the existing sample of evolved stars with deep, wide-field imaging in the sub-mm, providing details on the range of contributions from extended emission. Meanwhile, the shallower continuum observations of the wider sample will provide a statistical overview of cold dust across a range of classes of evolved stars.

The far-infrared spectral slope, $\alpha$ (where $F_{\nu}\propto\nu^{\alpha}$) provides the nature of the FIR emission source. In particular, it can identify variations in the dust-opacity slope, $\beta$ (where $\kappa_{\nu}\propto\nu^{\beta}$; $\alpha\approx-\left(2+\beta\right)$ with some temperature dependence for very cold dust), which can depend on the composition and size of circumstellar dust grains (e.g. Hoogzaad et al., 2002; Dehaes et al., 2007). Many authors have shown that the $\beta$ depends strongly on grain size (e.g. Kruegel & Siebenmorgen, 1994; Rodmann et al., 2006; Ricci et al., 2011; Testi et al., 2014) and that, assuming an MRN-like power-law size distribution, this dependence is strongest for maximum grain sizes $\lambda/2\pi<a_{\rm max}<3\lambda$ (Draine, 2006). Our mapping observations are hence most sensitive to $a\sim 50~{}\micron-3\,mm$, significantly larger than predicted by dust-formation models (e.g. Höfner, 2008) but consistent with the largest pre-solar grains (e.g. Amari, 2014). Similarly, grains of different composition or structure can have very different spectral indices (e.g. Mennella et al., 1998; Jager et al., 1998), with laboratory and theoretical values for typical dust analogues covering the range $\beta=1-2.5$. Dharmawardena et al. (2018) found variations in $\beta$ in the outer envelope of IRC+10216, suggesting evolution of the dust as it is processed by the interstellar radiation field and integrated into the ISM.

Other sources contributing to sub-mm emission include diffuse backgrounds, and the star’s optical and radio photospheres. These can be disentangled by a combination of mapping observations and the spectral slope, $\alpha$.

Diffuse background emission is normally filtered, but inhomogeneities may add or subtract flux from the PSF core. Interstellar dust is expected to have $\alpha\approx-3.7$ (Planck Collaboration et al., 2014), but bright patches will not normally be spatially coincident with the AGB star in mapping observations, and diffuse emission this strong and inhomogeneous is likely to affect only a few stars (notably those in the Galactic plane) and be visible on large scales.

The star’s optical photosphere should have $\alpha=-2$, will be spatially compact, and is already included in the GRAMS model fit. Uncertainties in the GRAMS luminosity and temperature may lead to uncertainties in the extrapolated sub-mm flux, but will be small.

The star’s radio photosphere (unresolved emission arising from free electrons interacting with H i and H₂ in the stellar chromosphere) is opaque at millimetre wavelengths, with a wavelength-dependent radius (e.g. Matthews et al., 2018). These canonically have $\alpha\approx-1.86$ (Reid & Menten, 1997; Vlemmings et al., 2019): slightly shallower but insufficiently different to $\alpha=-2$ to distinguish from cold dust in the staring sample. The sub-mm photospheric radius is similar to the optical photosphere near the peak of the SED, but typically cooler ($\sim$2000 K), so the sub-mm flux does not greatly exceed that of the optical photosphere’s Rayleigh–Jeans tail (e.g. Kervella et al., 2016; Vlemmings et al., 2019).

Observations of IRC+10216 by Menten et al. (2006) also suggest the sub-mm contribution from both optical and radio photospheres is small: the radio flux can be extrapolated to the sub-mm following

for frequency $\nu$, distance from the source $D$, and stellar luminosity $L_{\rm IRC+10216}$. Given our detection limit of $\sim$ 11 mJy, an X-band (8.4 GHz) flux of $\sim 0.6$ mJy (Menten et al. 2006, Table 2), and assuming that IRC+10216 is roughly twice as luminous as the median of the LMC population, only sources within $\sim$200 pc would have detectable radio photospheric emission in our observations.

Stellar chromospheres can also contribute significant radio flux with a much shallower spectral dependence that can dominate stellar continuum emission at longer (radio) wavelengths (e.g. O’Gorman et al., 2020, for Antares). The radio photosphere canonically extends to $\sim$2–3 stellar radii (R_∗) at 22 or 43 GHz, and typically has a lower brightness temperature than the star’s effective photospheric temperature (e.g. Lim et al., 1998; O’Gorman et al., 2015, for Betelgeuse). In the case of Betelgeuse, the chromosphere only contributes 10–20 per cent more flux than the expected Rayleigh–Jeans tail at 850 µm (O’Gorman et al., 2017). Consequently, we expect little additional flux at these wavelengths from stellar chromospheres.

Other physical explanations for sub-mm excess include a low $\alpha$ due to variations in grain size and composition described above. However, mm-sized grains are not generally predicted to form in the outflow (the largest predicted and observed being a few $\mu$m in diameter; Höfner 2008; Norris et al. 2012; Scicluna et al. 2015), while few dust analogue materials have bulk spectral indices $\beta>-1$ (e.g. Mennella et al., 1998; Demyk et al., 2017a, b; Mutschke & Mohr, 2019). Brunner et al. (2018) investigated whether changes in grain structure, composition or size could explain similar excess emission in LABOCA observations of R Scl, but found all of these options unsatisfactory, ultimately suggesting that either additional dust components (such as polycyclic aromatic hydrocarbons) were required, or that the emissivities of the dust may exhibit a significant change in the sub-mm (see also Gordon et al., 2014, in the Magellanic Clouds). While data reduction may also play a role, Fig. 7 suggests a good agreement between the NESS and interpolated SPIRE fluxes at 450 $\mu$m: the only deviations being stars where we may underestimate the flux by a factor of up to two, whereas anomalous cold dust will produce an excess of flux at both 450 and 850 $\mu$m.

Consequently, we expect a marked excess in sub-mm flux in NESS should only be attributable to an anomalous cold dust reservoir, especially where the mapping sample resolves it from the point-source star.

We here present some early results of NESS continuum observations, based on the data reduction and fluxes presented in Sect. 3.4. At 850 µm our detection rate is 73 per cent, but varies significantly across our subsamples; the rate approaches 90 per cent for sources in Tier 2 (shown in green in Fig. 6), while the detection rate in Tier 4 (in red) is just above 50 per cent. The origin of this decrease is not yet clear, although increased confusion from diffuse interstellar dust at larger distances seems like an obvious contributor. Follow-up observations with a compact interferometer would reveal whether this is the case by filtering out the contamination. However, this requires care not to remove emission from historic dust mass loss.

The lower sensitivity at 450 µm  leads to a reduced detection rate of $\sim$30 per cent. A similar pattern exists throughout the subsamples, decreasing from $\sim 40$ per cent for Tier 2, to 24 per cent for Tier 4.

Fig. 11 compares our measured fluxes to the model predictions of GRAMS (Sect. 2.3.3). A 10 per cent calibration uncertainty has been added in quadrature to the statistical uncertainty. The GRAMS models clearly and systematically underestimate the sub-mm continuum flux, typically by a factor of 3–10, the factor being larger at 850 µm than 450 µm. This is not unexpected, given the limitations of the GRAMS grid (i.e., it is only intended to predict optical–mid-IR photometry and is tailored to the Magellanic Clouds), the uncertainties of the modelling (i.e., the models are fitted to relatively few bands in the near-/mid-infrared) and the uncertainties of the data reduction. However, similar to Dharmawardena et al. (2018), Dharmawardena et al. (2019) and Maercker et al. (2018), we find systematically higher sub-mm fluxes than the models predict, arguing against data-reduction uncertainty as the origin of the discrepancy. This might indicate that large reservoirs of cold dust are present, or that the dust properties are not well represented by those used in GRAMS. This figure also includes sources where extended dust emission was detected by Dharmawardena et al. (2018). These also all show systematic excess emission compared to the GRAMS models, hence physically linking a sub-mm excess flux (compared to GRAMS) with spatially extended cold dust emission.

Fig. 12 plots the 850 $\mu$m excess as a function of the DPR derived from GRAMS, roughly equivalent to plotting a mid-infrared colour against a mid-IR – sub-mm colour (e.g. [12]$-$[25] vs. [25]$-$[850]). A trend of increasing excess with DPR is visible, particularly across Tiers 2 and 3. To determine its statistical evidence, we use Bayesian model selection to calculate the posterior odds (Kass & Raftery, 1995), using dynesty to compute the evidence with Nested sampling (Skilling, 2004; Speagle, 2020), following Dharmawardena et al. (2020). This approach shows that a power-law relationship is preferred over no relationship by a significant margin ($\sim 10^{10\,000}$), while a power law with a constant flux above a break point is indistinguishable from the unbroken power law. While this framework is susceptible to some uncertainty, particularly from dependence on the choice of prior, it is difficult to envisage this overcoming a multi-googol preference.

The easiest interpretation of this trend relates to the amount of cold dust in the envelopes: stars with higher present-day DPRs will tend to have been producing dust for longer, and are therefore likely to have filled a greater fraction of the SCUBA-2 beam with dust. Due to specific circumstances of source distances the cold dust component does not have to be extended (spatially resolved), but by the dictates of physics the extended component has to be cold (assuming thermal equilibrium). As the GRAMS models use a fixed-density profile, dust composition and outer radius, they do not account for physical changes that may impact the sub-mm emission; detailed models which treat the sub-mm fluxes and radial profiles appropriately may better capture the behaviour. Stochastic changes in DPR or differences in dust composition may be responsible for the scatter.

We derive spectral indices from the above SCUBA-2 fluxes (Fig. 13 and Tab. 3), ensuring that the results are robust and upper limits are handled consistently, by adopting the method of Scicluna et al. (2020), to which we refer the reader for the details of the fitting and parameters for the Markov Chain Monte Carlo (MCMC). As expected, sources with a well-constrained spectral index (blue points) separate from lower limits (orange); the latter corresponding to sources detected at 850 $\mu$m but not 450 $\mu$m, which can have any positive spectral index. In general, the upper limits at 450 $\mu$m do not place strong constraints on $\alpha$, though the limits remain consistent with the rest of the sample. Sources with well-constrained $\alpha$ cluster around $\alpha\approx-2$, which manifests as the tight correlation in Figure 11, and consistent with a cold-dust reservoir (Sec. 4.4.2).

The two outliers (one at positive $\alpha$ and one at negative $\alpha$) are also the two outliers in Fig. 6. The source at $\alpha=0.3\pm 1.0$ is IC 418, a planetary nebula with significant free-free emission. The spectrum continues to rise into the radio, 1.72 Jy at 6 cm (Griffith et al., 1994). The source at $\alpha=-7.4\pm 1.1$ is AK Hya, a nearby AGB star with a relatively low mass-loss rate and high space motion (cf. McDonald et al., 2018). To eliminate the discrepancy, the true 450 $\mu$m flux would have to be lowered by factor of 20, which seems improbable even for a 4-$\sigma$ detection, as the expected frequency of such outliers is $\sim 0.00007$. If real, the origin of such an anomalous spectral index is an enigma; a small number of nearby debris discs have similarly steep spectra, attributed to anomalous dust compositions or size distributions (e.g. Ertel et al., 2012; Marshall et al., 2016).