© ESO 2018

1. Introduction

Despite their rarity, massive (> 20 M_⊙) stars are major agents of galactic evolution via the deposition of chemically enriched material, mechanical energy and ionising radiation, while dominating integrated galactic spectra in the UV and IR windows (the latter via re-radiation by hot dust). Despite this central role in galactic astrophysics, their short lives are poorly understood in comparison to stars such as the Sun. Unlike lower-mass stars, heavy mass-loss has long been recognised as a critical factor, along with rotation and the presence (or otherwise) of magnetic fields, in governing their evolutionary pathway (Ekström et al. 2012) and the nature of their demise, core-collapse or pair-production supernova (SN), Gamma-ray burst or prompt, quiet collapse to black hole (BH). As might be anticipated from this, the nature of the stellar corpse, either neutron star or BH, is intimately related to the pre-demise mass-loss history, which is of particular interest given the detection of gravitational waves from coalescing black holes (Abbott et al. 2016).

For single stars it has long been thought that line-driven radiative winds serve as the mechanism by which massive stars lose mass. Unfortunately, recent studies have reported observationally derived mass-loss rates of OB stars which are discordant by up to a factor of 10 (e.g. Puls et al. 2006; Massa et al. 2003; Fullerton et al. 2006; Sundqvist et al. 2011), due to uncertainties related to the degree of structure (“clumping”) present in O and B star winds (e.g. Prinja et al. 2010; Prinja & Massa 2013; Sundqvist et al. 2011; Šurlan et al. 2012). In order to emphasise the far-ranging consequences of this ambiguity, we highlight that it is no longer clear that radiatively-driven winds are able to drive sufficient mass-loss to transition from H-rich main sequence (MS) to H-depleted post-MS Wolf Rayet (WR) star. A popular supposition is that the consequent mass-loss deficit is made up by the short-lived transitional phase between these evolutionary extremes. Observations of ejection nebulae associated with both hot (luminous blue variable; LBV) and cool (yellow hypergiant and red supergiant; YHG and RSG) transitional stars imply phases of instability in which extreme, impulsive mass-loss may occur (Ṁ ≥ 10⁻⁴M_⊙ yr⁻¹). However it is not clear which stars encounter such instabilities, nor whether the duration of this phase is sufficient to strip away the H-rich mantle to permit the formation of WRs. Observations of apparently “quiescent” LBVs show that they support dense winds with mass-loss rates comparable to WRs (Clark et al. 2014b) outside of eruption. However direct mass-loss rate determinations for very luminous cool transitional stars are few and far between (de Jager et al. 1988); a critical weakness of current stellar evolutionary codes (e.g. Ekström et al. 2012).

An alternative evolutionary channel has been suggested by the recent recognition that ∼70% of massive stars are found within binary systems (de Mink et al. 2014; Sana et al. 2012, 2013). Interaction between both components may lead to extreme mass-loss from the primary (e.g. Petrovic et al. 2005), allowing for the formation of WRs and favouring the production of NSs over BHs (e.g. Ritchie et al. 2010), while the mass gaining undergoes rejuvenation. In extreme cases binary merger may also lead to the production of a massive blue straggler.

In both mass transfer and merger scenarios for binary evolution, substantial modification of the stellar initial mass function may be anticipated, while binary merger may provide a viable formation route for very massive stars (> 100 M_⊙; Schneider et al. 2014, 2015). Moreover, even outside of evolutionary phases dominated by thermal- and nuclear-timescale mass-transfer, massive binaries provide valuable insights of the properties of stellar winds via observations of the resultant wind collision zones. Indeed both the high energy (thermal) X-ray and low energy (non-thermal synchrotron) mm/radio emission that results from shocks within colliding wind binaries (CWBs) provide additional binary identifiers.

Extending this paradigm, when found within stellar aggregates, the collective action of stellar winds and supernovae (SNe) yield powerful cluster winds (e.g. Stevens & Hartwell 2003) that can either disperse or compress their natal giant molecular clouds, respectively inhibiting or initiating subsequent generations of star-formation. Moreover, SNe and their interaction with cluster winds (and the shocks within CWBs) have also been implicated in the production of Galactic cosmic rays and attendant very high energy γ-ray emission (e.g. H.E.S.S. Collaboration 2018; Abramowski et al. 2012; Bykov et al. 2015).

Given the above considerations, young massive stellar clusters (YMCs) form ideal laboratories for the study and resolution of these issues due to their co-eval stellar populations of a single metallicity. Consequently YMCs have recently received increased attention as (near-IR) Galactic surveys have yielded numerous further examples. Discovered over half a century ago, Westerlund 1 (Wd1; Westerlund 1961) subsequently escaped detailed investigation due to significant interstellar extinction. The serendipitous detection of multiple radio sources associated with cluster members (Clark et al. 1998; Dougherty et al. 2010, henceforth Do10) sparked renewed observational efforts which revealed Wd1 to host a uniquely rich and diverse population of massive stars (Clark & Negueruela 2002; Clark et al. 2005). Specifically, Wd1 appears to be co-eval and at an age (∼5 Myr; Negueruela et al. 2010; Kudryavtseva et al. 2012) where cool supergiants and hypergiants may co-exist with WRs; the uniquely rich population of both (Clark et al. 2005; Crowther et al. 2006a), as well as > 100 OB supergiants (Negueruela et al. 2010; Clark et al., in prep.) makes Wd1 a powerful laboratory for the study of massive stellar evolution.

As such Wd1 has received attention across the electromagnetic spectrum from radio (Kothes & Dougherty 2007, Do10), IR and optical (e.g. Brandner et al. 2008, and refs. above) through to X-ray (Muno et al. 2006a,b; Clark et al. 2008, henceforth Cl08) and higher energies (GeV and TeV; Ohm et al. 2013; Abramowski et al. 2012, respectively). Multi-epoch spectroscopic radial velocity (RV) surveys have identified a rich population of binaries (Ritchie et al. 2009a, and in prep.), with tailored modelling of individual systems clearly revealing the influence of binarity on massive stellar evolution (e.g. Ritchie et al. 2010; Clark et al. 2011, 2014b).

In order to better understand the nature of the massive single and binary stellar populations of Wd1 we undertook Atacama Large Millimetre Array (ALMA) Band 3 (100 GHz) continuum observations of Wd1 in 2015. The millimetre waveband is a uniquely powerful diagnostic of mass-loss and, in conjunction with observations at other wavelengths to constrain the continuum spectral energy distribution (SED), may determine the radial run of wind clumping via thermal Bremsstrahlung emission (e.g. Blomme et al. 2003). Moreover, CWBs may be identified by millimetre-radio observations due to either the presence of a non-thermal continuum component from synchrotron emission originating in the wind collision zone, or excess mm emission if such shocks are optically thick (Pittard 2010, 2011).

Fenech et al. (2017) presented results for the brightest continuum source within Wd1, the supergiant B[e] star and interacting binary Wd1-9 (Clark et al. 2013). In this paper we present a census of the remaining 3-mm sources and discuss their association with stellar counterparts where appropriate. To enable ease of comparison to the radio study of Wd1 by Do10, we choose to mirror the structure of that work in this paper, with departures where necessary due to novel results or synergies between radio and mm datasets. Finally, significant work has been undertaken to establish a distance to Wd1, with most estimates falling between 4 and 5 kpc (Clark et al. 2005; Crowther et al. 2006a; Kothes & Dougherty 2007; Brandner et al. 2008; Negueruela et al. 2010; Kudryavtseva et al. 2012; Clark et al., in prep.), for the purpose of this paper we adopt a distance, d = 5 kpc to Wd1.

2. Data acquisition, reduction and analysis

2.1. Observations and imaging

ALMA was used to observe Wd1 on the 30th June and the 1st July 2015 (Project code: 2013.1.00897.S) covering the central approximately 3.5 sq. arcmin area with 27 pointings. The observations were made at a central frequency of ∼100 GHz with a total usable bandwidth of 7.5 GHz over four spectral windows centred on 92.5, 97.5, 102.5 & 104.5 GHz respectively. Each spectral window contains 128 channels with a channel frequency width of 15.625 MHz (only 120 channels are usable). The array consisted of 42 antennas with baselines ranging from 40 to 1500 m and a total on-source integration time per pointing of 242.3 s (∼4 min). The data were calibrated using the standard ALMA pipeline procedures in Common Astronomy Software Applications (CASA; pipeline version 4.3.1: r34044) and included the application of apriori calibration information as well as flagging of erroneous data. Observations of J1617-5848 were used to perform the phase and bandpass calibration and observations of Titan and Pallas were used to amplitude calibrate the data with assumed flux densities of 228.96 and 82.01 mJy respectively (at 91.495 GHz).

As will be discussed in detail in Sects. 3–6 a number of the detected stars in these observations are resolved. In order to ensure that this is not the result of potential phase errors in the data associated with the pipeline water vapour correction performed, the data were re-calibrated applying slightly less flagging than in the original pipeline calibration and also making use of the REMCLOUDs package. This is used to correct for the effect of any water in the form of fog or clouds present in certain weather conditions. As each section of data (30th June and 1st July) were calibrated independently, we found the addition of REMCLOUD marginally improved the phase calibration for the 30th June but made no discernible improvement for data from the 1st July. We therefore proceeded using re-calibrated 30th June data with REMCLOUD and re-calibrated 1st July data without REMCLOUD.

Following initial calibration (both with initial and re-calibrated data), several iterations of phase self-calibration were applied to the data. At each stage images of the full mosaic field were produced. A selection of the brightest compact sources across the field were cleaned to produce the desired model. A selection of pointings, mainly centred around Wd1-9 were used with the model by GAINCAL to calculate the calibration solutions. APPLYCAL was then used to apply the corrections to all pointings. A final round of amplitude and phase calibration was also performed incorporating information from both Wd1-9 and Wd1-26 in the model.

The data were initially imaged using the mosaicing and multi-frequency synthesis functions in the CASA CLEAN software. Cleaning was done only on the bright compact component of Wd1-9 (which has a flux density of 153 mJy) in order to subtract this component from the data and limit high dynamic range problems in subsequent imaging. The final images were produced using the mosaicing and multi-scale capabilities within CLEAN including data from all pointings to produce a single wide-field map of Wd1 (see Figs. 2 and 3). This utilised Cotton-Schwab cleaning and natural weighting with scale sizes set at 1, 14, 12, 37 & 119. The final image has a fitted beam of 750 × 570 mas. The final full-field image was primary beam corrected following cleaning to account for the change in sensitivity across the primary beam (PB). Example images of individual sources will be shown in Sects. 4–6 and have been taken directly from the PB-corrected wide-field image. A complete catalogue of images of each source is also included in Appendix C. For subsequent data analysis, the images were transferred to the AIPS (Astronomical Image Processing Software) environment.

To aid cross-correlation and identification, the absolute astrometry of the final images were checked utilising the 8.6 GHz radio image and the FORS R-band image (which had previously been aligned to the radio image) as presented in Do10. Following the procedure from Do10, we determined any offset between the ALMA and FORS image by first using the AIPS task HGEOM to make the geometry of the two images consistent and then performing a Gaussian fit (using the AIPS task JMFIT) to the compact component of Wd1-9 in both images. A comparison of the fitted peak positions revealed a small offset between the images. As we performed self-calibration on the ALMA data, which causes absolute positional information to be lost, we chose to shift the ALMA image in order to bring the peak positions of Wd1-9 (as determined by the Gaussian fit) into alignment. The aligned FORS and ALMA images are presented in Fig. 2, the ALMA image showing identified stellar sources is shown in Fig. 3 and the ALMA and radio 8.6 GHz image is shown in Fig. 4. Both the original image and the primary-beam corrected image were shifted in this way and used for performing all subsequent analysis.

2.2. Identification of 3 mm stars

In order to identify and catalogue the sources present in Wd1, the SEAC source extraction software was used (see Peck 2014; Morford et al. 2017, for further details). This utilises a floodfill algorithm to search the map and identify discrete “islands” of emission by locating pixels above a given (seed) threshold and subsequently adding neighbouring pixels to the “island” down to a lower (flood) threshold. Thresholds of 5 σ and 3 σ were set for the initial use of SEAC on the ALMA wide-field image, where σ is a regionalised noise level calculated by determining the rms level in individual cells of a 6 × 6 grid across the whole image. This was necessary to account for the change in noise level across the image especially towards the strong emission from Wd1-26 and Wd1-9. SEAC was used on both the original and primary-beam corrected images. Sources determined by SEAC with a seed threshold of 5 σ in the original image (i.e. not primary beam corrected) are considered as detections.

This resulted in the detection of 98 sources in the Wd1 field. These detections were cross-correlated with previously published catalogues (including Clark et al. 2005; Crowther et al. 2006a; Negueruela et al. 2010; Do10 and Cl08 as well as the VPHAS point sources catalogue for the region) in order to identify the observed sources. Following this cross-correlation process, the sources identified at 3 mm fall broadly into three categories; those that are previously known and identified at multiple wavelengths including optical and radio, those that are detected in the radio observations from Do10 though have no other counterparts and those that are detected only in these ALMA observations. Table 1 lists a summary of the number of sources detected at 3 mm within these categories along with those for relevant stellar spectral types where possible.

There are also a number of optically identified sources that do not appear in these ALMA images. In particular there are several sources that have the same or similar spectral classifications to those that have been detected. In order to attempt to locate any weak emission from these sources, further SEAC runs were performed with seed thresholds of 4.5 and 4 σ. This identified a further three sources albeit with lower detection thresholds, which have also been included in the final source list. The total 3 mm catalogue list therefore contains 101 sources. Each source has been given an ALMA 3 mm catalogue number, designated as FCP18, in order of right ascension. All information for these sources is contained in Tables 2 and 3 for the stellar sources and Table A.1 for the ALMA-only sources.

3. The millimetre emission in Wd1

The millimetre emission observed in Wd1 shows both distinct isolated sources as well as more extended regions of diffuse emission which broadly coincide with similar regions seen in the radio images (see Fig. 4).

3.1. 3-mm star characteristics

Of the total 101 objects detected in these observations, 50 have been identified as stellar sources that have previously observed optical and/or radio counterparts. Their identification, flux density as measured from the primary-beam corrected image alongside the radio-mm spectral index are presented in Table 2. In addition to SEAC, Gaussian fitting was performed on each of the sources (see Sect. 3.1.1). The flux densities presented are either those taken from the Gaussian fitting (for unresolved or Gaussian-like sources) or the integrated flux density of the source “island” returned by SEAC. The errors on the flux density are those from calculating the integrated flux density (either returned from the Gaussian-fitting or SEAC) and the expected absolute amplitude calibration error for ALMA which is taken to be 5% (see ALMA Cycle 2 Technical Handbook) combined in quadrature. For sources with a clearly defined compact component with extended structure, flux densities and where possible spectral indices have been provided for both the compact and extended components separately. Radio observations have shown extended nebulae associated with several sources including Wd1-4 and Wd1-20. This extended emission appears in some cases as multiple components within these higher resolution ALMA images presented here. In particular, four of the FCP18 catalogue sources are part of the diffuse extended emission associated with the nebula of Wd1-20 identified by comparison with the radio emission, as well as three sources which are associated with the nebula of Wd1-4. These components have been listed individually in the final catalogue (Table A.1) for completeness. However they have been treated as part of Wd1-4 and Wd1-20 for calculating both flux densities and spatial extent.

3.1.1. Spatial resolution

A large number of the stellar sources appear to be resolved and in order to determine their spatial extent two approaches were used. Primarily a Gaussian fit was performed using the AIPS task JMFIT. Estimates of the peak flux density and position from SEAC were used and only pixels >3 times the rms local to the source were included in the fit. Three dimensional plots of the source structure and the resulting model fit were used to visually assess the fit (an example of which can be seen in Fig. 1). The final convolved and deconvolved source sizes are presented in Table 3 which lists the spatial information for each source, including the positional offset from the observed ALMA peak to the catalogued source position taken from the literature (see the caption of Table 3 for details). Where the Gaussian fit was deemed unreasonable i.e. when the source structure was distinctly non-Gaussian, a largest angular size (LAS) was measured and is listed in Table 3 (the LAS sizes are not deconvolved). The IRING AIPS task was also used to perform integrated annular profiles of each source centred on the position of the peak brightness. The errors included for source sizes measured using Gaussian fitting are those directly returned by the fitting procedure. For the convolved sizes, these represent the errors calculated in the fit. For the deconvolved sizes the errors represent the potential minimum to maximum range the deconvolved size could have (based on the convolved size errors), as a result of subtracting the image restoring beam, and therefore likely overestimate the true error.

A total of 27 of the stellar sources have been determined to be resolved i.e. have sizes larger than the image restoring beam and have convolved and deconvolved sizes listed in Table 3. A further 17 appear to be partially resolved (i.e. those that have only one dimension deconvolved listed in Table 3), while the remaining 7 sources are unresolved. The sources which are partially or fully resolved separate into two main categories: those that appear point-like i.e. are generally representable by a Gaussian and those that have extended nebular emission. The former include stars of every spectral type found within the cluster. Given that they are in general relatively faint detections, the errors associated with the WRs and OB super-/hypergiants are systematically larger in relation to source size than those associated with the cool super-/hypergiants. Despite this, the majority of WRs and OB super-/hypergiants appear extended, with objects such as Wd1-243 (LBV; Sect. 6.2) and WR A and L (Sect. 4) robustly so and hence serving as exemplars. The latter mainly appear to be extended nebulae associated with the cool hypergiants e.g. Wd1-26.

One final feature of note is that a number of compact sources for which it is possible to fit a Gaussian appear to show excess emission in the shoulders of their profiles (see Fig. 7 and Sect. 4.2.3). This is clearly present in the profiles of the compact sources associated with the YHGs Wd1-4, -12a and 16a but also fainter sources such as those associated with e.g. WR D, J, K and Q, albeit with less confidence.

We defer discussion of these sources for subsequent sections, where their extent can be compared to predictions derived from their varied wind properties as well as the properties of ejection nebulae associated with other examples.

Comparison of the mm and radio flux densities of those sources associated with the extended emission reveal an apparent flux density deficit in our ALMA observations if, as suggested by the radio observations, the emission mechanism is optically-thin free-free from ionised circumstellar ejecta. This strongly suggests that some of the mm emission has been resolved out i.e. occupies larger spatial scales than those effectively sampled by ALMA during these observations. As a result we refrain from determining ionised ejecta masses for these sources, favouring the estimates presented in Do10.

3.1.2. Spectral indices

Where possible the spectral index between the 8.6 GHz (3.6 cm) radio and 3 mm ALMA flux densities are presented in Table 2. The majority of stellar sources appear to have spectral indices consistent with partially optically-thick free-free emission. However a number appear to show flatter or inverted spectral indices. The divergence of the spectral index from the canonical value of α = 0.6 expected for a partially optically-thick wind (Wright & Barlow 1975) can occur from a number of reasons, such as changes in the run of ionisation or wind geometry with radius, or the presence of both optically-thick and thin components (such as optically-thick clumps in a structured wind; Ignace & Churchwell 2004). In long-period (P _orb ≳ 1 yr) colliding wind binaries, an additional non-thermal component caused by wind interaction outside the respective radio/mm photospheres can lead to a reduction or flattening of the spectral index (e.g. Chapman et al. 1999, Do10). Conversely, in more compact binaries thermal emission from material in the wind collision region (WCR) may come to dominate the spectrum at mm or radio wavelengths (Stevens 1995; Pittard et al. 2006; Pittard 2010; Montes et al. 2015). For instance, optically-thin thermal emission (α ∼ −0.1) from an adiabatic WCR may dominate emission at long (cm) wavelengths, leading to a flatter than expected spectrum (Pittard et al. 2006). Alternatively, for shorter period systems, cooling of material in the WCR becomes increasingly important (Stevens et al. 1992) leading to optically-thick thermal emission (Pittard 2010) which could dominate the spectrum at mm wavelengths, resulting in a steeper mm-radio continuum spectral index.

For all eventualities it is important to recall that the time-dependent line of sight through the complex geometry of the WCR and stellar winds imposed by orbital motion may also be reflected in variability in the spectral index. This is of particular concern for very compact (contact) binaries in which interactions occur before appreciable wind acceleration, where computational limitations preclude accurate modelling of the resultant wind and the WCR geometries and interactions, and hence quantitative predictions for the emergent spectrum.

3.1.3. Mass-loss determination

There is considerable advantage to using free-free mm/radio fluxes for determining mass-loss for massive stars in that, unlike Hα and UV, the emission due to electron-ion interactions in their ionised winds arises at large radii, where the terminal velocity will have been reached. Therefore the interpretation of the mm/radio flux densities is more straightforward and is not strongly dependent on details of the velocity law, ionization conditions, inner velocity field, or the photospheric profile. Though the greater geometric region and density squared dependence of the free-free flux makes these continuum observations sensitive to clumping in the wind, there is evidence that clumping decreases in the outer wind regions (e.g. Runacres & Owocki 2002; Puls et al. 2006).

The mass-loss rate is related to the observed free-free emitted radiation as$$S_{\nu }=2.32\times {10}^4{\left(\frac{\dot{M}\sqrt{f_{\text{cl}}}}{\mu v_{\infty }}\right)}^{4/3}\frac{1}{D^2}{\left(\gamma g_{\text{ff}}\nu \overline{Z^2}\right)}^{2/3},$$(1)

where, S_​ν is our observed radio flux in mJy measured at frequency ν in Hz; Ṁ is in M_⊙ yr⁻¹; the terminal velocity v_∞ is in km s⁻¹; D is the distance in kpc (see e.g. Wright & Barlow 1975). The quantities μ, Z, and γ are the mean molecular weight per ion, ratio of electron to ion density, and mean number of electrons per ion. The Gaunt factor, g _ff, can be approximated by$$g_{\text{ff}}\approx 9.77\left(1+0.13 \log \left(T_{\text{ e}}^{3/2}/\nu \sqrt{(\overline{Z^2})}\right)\right),$$(2)(e.g. Leitherer & Robert 1991).

Since the free-free emission process depends on the density-squared, it is affected by wind clumping. Equation (1) includes a simple account of this, such that all clumps are assumed to have the same clumping factor given by f _cl = ⟨ρ²⟩/⟨ρ⟩², where the angle brackets indicate an average over the volume in which continuum emission is formed. A given flux can therefore be interpreted as a certain mass-loss rate for a smooth wind (f_cl = 1), or correspondingly as a lower mass-loss rate in a clumped wind (f_cl > 1). Note, with an assumption of no inter-clump material, f_cl here is related to the reciprocal of the volume filling factor. Regarding the remaining terms in the above relations we adopt T_e = 0.5T_eff (e.g. Drew 1989) and that for OB stars hydrogen is fully ionised, He⁺ dominates over He²⁺, and a helium abundance of n_He/n_H = 0.1 (i.e. μ ∼1.4, Z = 1, γ = 1). We note that for hot stars the dominant H ionization stage is controlled essentially by the wind density through mass-loss rate and clumping. However in the case of strong clumping in the mm formation region, recombination would be enhanced and thus favour He⁺. For the chemistry of the WR star winds we have adopted here the following generalisation for μ based on Leitherer et al. (1997); μ = 4.0 for WN6 and earlier types; μ = 2.0 for later-type than WN6; μ = 4.7 for WC8 and WC9. Once more, He⁺ is assumed to be the most prevalent in the WR mm-emitting region and we have assumed Z = 1 and γ = 1. Following Leitherer et al. (1995) we assume the B hypergiant (Ia⁺) candidates in Wd 1 to have a chemical composition akin to LBVs and adopt μ = 1.6, γ = 0.8 and Z = 0.9. For the fundamental stellar parameters (T_eff, terminal velocity v_∞) we have adopted the (observational) compilations of Crowther et al. (2006b) and Searle et al. (2008) for B0-B5 supergiants; Crowther (2007) and Sander et al. (2012) for WN stars; and Sander et al. (2012) for WC stars.

Our mass-loss estimates assume the mm fluxes are not affected by binarity, either via synchrotron emission and/or enhanced thermal emission from a (strongly radiative) WCR (e.g. Pittard 2010). After excluding two sources with spectral indices strongly indicative of non-thermal emission we present the resultant mass-loss rates for a total of 27 stars in Table 4 and further discuss these in the following sections.

3.2. Mm sources lacking counterparts

The remaining 51 detected sources which have no catalogued optical or radio counterpart are listed in Table A.1. Their positions in Wd1 are shown in Fig. B.2 and images of each source are available in Fig. C.1. Seven of these sources have been determined to be knots of emission associated with the extended nebulae of Wd1-4 and Wd1-20. Of the 44 > 5 σ detected sources that are unique to these data, at least five appear to be obviously associated with extended emission, in most cases these are clearly low surface-brightness areas with brighter points of emission accounting for the source detections. A further 31 sources are more identifiable as potential sources appearing as well-defined components within diffuse mm emission. However, it is not possible from these observations alone to determine if they are for example stellar sources or again merely brighter components of the more diffuse emission. Five of the ALMA-only sources with no previous counterpart do appear to be bright discrete sources in isolation and not associated with any diffuse emission.

Large surveys of the background galaxy population at mm wavelengths are now making use of high sensitivity observations to probe the faint source population. Whilst some of these surveys are performed at 3-mm (100 GHz), they tend to be less sensitive (e.g. Mocanu et al. 2013, use observations at 95, 150 and 220 GHz though only probe down to ∼4 mJy). The vast majority utilise observations at 1.1–1.3 mm and the most sensitive of these using ALMA Band 6 or 7 (e.g. Ono et al. 2014; Carniani et al. 2015; Oteo et al. 2016; Hatsukade et al. 2016). Using the cumulative source counts from these surveys, it is however possible to estimate the number of background sources expected within the field observed for Westerlund 1.

Massardi et al. (2016) study the spectral evolution of sources from mm to radio wavelengths and show that the majority of the population display a down-turning spectral index with indices ranging from 0 to ∼ −1.8 for the 100–200 GHz frequency range (3–1.5 mm) with an average value around −0.5. The noise level in the central regions of the Westerlund 1 field is ∼28 μJy beam⁻¹ resulting in a 3-σ value of ∼85 μJy beam⁻¹. Assuming a −0.5 spectral index is representative, this would equate to ∼56–51 μJy beam⁻¹ at band 6–7. Umehata et al. (2017) used ALMA and the Aztec camera to measure the cumulative source counts and find S _{> 0.4 mJy} = ${9800}_{-2200}^{+5100}$ deg⁻². The observations of Wd-1 cover an area of approximately 7.3 arcmin², providing an estimate of approximately 20 background sources within our field. Likewise the theoretical models from Hayward et al. (2013) give S _{ALMA6 > 0.4 mJy} = 6897 deg⁻², estimating of ∼14 sources in the field. Both of these however, lack some of the faint source population (< 0.4 mJy) which should be detectable in our observations. More recently, Oteo et al. (2016) use ALMA Band 6/7 calibration survey observations to estimate source counts down to a flux density of 0.20 mJy, equating to ${11}_{-9}^{+10}$ sources in our field-of-view above this threshold. They also highlight the increased level of uncertainty in calculating source counts due to potential spurious source contamination, when employing lower detection thresholds often used to probe the fainter source population.

Interestingly, Shimizu et al. (2012) use hydrodynamical simulations to predict the background submillimetre galaxy population and present results for a number of ALMA bands including 3, 6 and 7. Their results would predict < 1 background sources in our field directly at 3-mm. However assuming the previous correction for spectral index and calculating for Band 6 would instead predict ∼20 background sources for this study. As discussed in Carniani et al. (2015), without knowing the intrinsic spectral energy distribution of the sampled sources (and their redshift distribution) it is difficult to accurately perform large extrapolations in flux density. Whilst Mocanu et al. (2013) have suggested that the 3-mm models under-predict the expected source counts, it is likely that there is also large uncertainty in predicting the background source counts in this manner.

3.3. Extended continuum emission

There are several areas of low surface brightness emission in the central region of Wd1. These are spatially coincident with the extended regions identified in the radio images presented in Do10 as A1–A8. The higher resolution of these ALMA observations is “resolving-out” some of the more diffuse emission on the larger spatial scales. This results in the significantly more patched or clumped background emission and the breaking down into more discrete components within these regions seen in the ALMA images. This is most clearly evident in Fig. 4 where the 3 mm emission (in colour-scale) appears to be associated with knots of emission in the more extended structures seen in the radio (blue contours). Five of the extended regions are readily associated with individual stars namely the RSGs Wd1-20 and Wd1-26, the BSG D09-R1, the sgB[e] star Wd1-9 and WR B, with these stars showing as more discrete sources in the ALMA observations. The remaining three extended regions are not associated with any currently known stellar counterpart, however, several of the ALMA-only sources identified in these observations are found in these regions. In some cases, such as around Wd1-20, the ALMA-only sources can be identified as low brightness diffuse emission associated with the larger nebula. In others such as near WR B, the other ALMA-only sources appear bright, compact and very similar to the other known stellar sources.

4. The stellar sources I. The Wolf Rayets

A key finding of our survey is that almost half of the 50 3-mm continuum sources associated with stellar counterparts within Wd1 are Wolf-Rayet (WR) stars (Table 2); we detect all such stars within our field-of-view. This is an important result since it raises the possibility of determining their mass-loss rates. Given the key role of WR-phase mass-loss in stellar evolution, e.g. changes in the mass-loss rates by a factor of 2–3 are sufficient to distinguish between black hole and neutron star formation post-SN (Nugis & Lamers 2000; Wellstein and Langer 1999), undertaking such a study for a sample of known distance and for which the properties of the progenitor population may be determined is particularly attractive.

Understandably, given the above, a number of previous studies have investigated the mm- and radio-continuum properties of WR stars. Consequently, in order to place our results in context, we first briefly review prior findings before discussing our results.

4.1. Previous mm and radio surveys

Abbott et al. (1986, and references therein) provide the first comprehensive survey of the radio properties of WR stars, presenting a distance (≤3 kpc) and declination (δ > −47^∘) limited 4.9-GHz radio survey of 42 of the 43 stars so selected. They report 28 detections, with spectral types ranging from WN5-8, WC4-9 and a single WN/C star. Complementing these data, radio observations of southern hemisphere WRs were obtained by (Leitherer et al. 1995 1997); Chapman et al. (1999) who report 16 detections (from 36 targets), with spectral sub-types ranging from WN6h-8h and WC5-8¹. Cappa et al. (2004) returned to study WRs north of δ > −46^∘, detecting 20 of 34 stars and yielding 16 detections of stars missing from the above studies at 8.46 GHz. Spectral sub-types of detections range from WN5-9h, WC6-9 and a single example of a WN/WCE transitional star. Regarding the nature of sources, Abbott et al. (1986) reported non-thermal emission for ∼21% of their detections, Cappa et al. (2004) 30% and Leitherer et al. (1995,1997), Chapman et al. (1999) at least 40% of their sample. Trivially, while the percentage of non-thermal emitters at radio wavelengths appears significant, the differences between these values likely reflect the different methods imposed on the authors by the nature of their respective datasets, making a final fraction difficult to ascertain.

After exclusion of apparent non-thermal sources, mass-loss rates may be inferred for the remainder of the stars. Willis (1991) re-interpreted the observations of Abbott et al. (1986) using refined wind terminal velocity and ionisation/abundance values, finding a mean mass-loss rate Ṁ ~ 5.3 ± 2.3 × 10⁻⁵ M_⊙ yr⁻¹ for the 24 stars considered. Leitherer et al. (1997) reported a mean Ṁ ~ 4 × 10⁻⁵ M_⊙ yr⁻¹ for all spectral sub-types with the exception of WC9 stars, for which it appears lower by a factor of > 2, while Cappa et al. (2004) find a mean Ṁ ~ 4±3 × 10⁻⁵ M_⊙ yr⁻¹ for the WN stars and Ṁ ~ 2±1 × 10⁻⁵ M_⊙ yr⁻¹ for the WC8-9 stars.

Mass-loss rates for individual stars from these samples, broken down by spectral sub-type and assuming f = 1 (i.e. no wind-clumping), are presented in Fig. 6 alongside those from these observations (see Sect. 4.2.2). Note that a number of confirmed binary systems² demonstrate apparently thermal radio emission and hence have been included in the plot. In such cases contamination from emission from the WCR (Sect. 3.1.2) cannot be excluded, which would lead to an artificially elevated mass-loss rate.

Finally, although not a systematic survey, 42-GHz observations centred on Sgr A^* by Yusef-Zadeh & Morris (1991) identified a large number of stars within the Galactic Centre cluster. Given the age of the cluster (6–7 Myr; Martins et al. 2007) it is expected that these stars will be less massive than those in Wd1, while no cut for non-thermal sources is possible with extant data. However, given the presence of a large number of the hitherto unrepresented WN9h sub-type we include these values in Fig. 6. In doing so we immediately highlight that the mass-loss rates of the WN9h stars appear systematically lower than those of other spectral sub-types, which we discuss further in Sect. 4.2.2.

To the best of our knowledge, the only systematic mm continuum survey of WRs was undertaken by Montes et al. (2015), who report 1.2-mm fluxes for all 17 sources observed. Targets appear to have been selected on the basis of previous radio detections; hence this compilation will automatically be biased towards brighter continuum sources. Combining these detections with radio data revealed that all sources showed a positive (thermal) mm-radio spectral index, even for stars such as WR79a and 105, which Cappa et al. (2004) suggest are possible non-thermal sources. Montes et al. (2015) refrain from determining mass-loss rates from these data and given the substantially larger sample size, we choose to compare mass-loss rated derived for our Wd1 WR cohort to those from radio observations.

4.2. WRs within Wd1

Crowther et al. (2006a) provide a census of the WR population of Wd1, finding spectral types spanning WN5 to WN10-11h and WC8-9. We detected all of the 21 WRs within the observational field – ranging in sub-type from the WN5 stars WR J and R through to the WN9-11 stars WR S and Wd1-13 as well as WC stars of sub-types WC8 (WR K) and WC9 (e.g. WR M). Only WR N, T, W and X were external to the field. Of the detections, 8 WN (WR A, B, D, G, L, O, U and Wd1-13) and 5 WC stars (WR C, E, F, H and M) display current binary signatures (Table 2)³, with the WN10h/BHG star WR S postulated to have evolved via binary channel prior to disruption via a SN (Clark et al. 2014b). We caution that the six WN detections for which no evidence for binarity currently exists (WR I, J, P, Q, R and V) have not been the subject of an RV survey due to a lack of appropriate emission lines that can function as RV diagnostics. Finally the WC8 star WR K shows no evidence for RV shifts (Ritchie et al. 2009a), although this could simply reflect a face-on inclination or a wide and/or highly eccentric orbit, with the latter also potentially explaining the lack of secondary binary diagnostics, which are known to be transient phenomena in some such systems, in this and other cluster WRs.

4.2.1. mm-radio spectral indices

Employing the 3.6-cm radio observations of Do10 we calculate the two-point mm-radio spectral index, α, for each star. Where no flux density is reported for a source in D010, we find a limiting spectral index based on a 3.6 cm flux density limit of 170 μ Jy beam⁻¹ (see Do10 for details). We first examine those sources with both mm and radio detections since their spectra are most constrained and they provide exemplars for the remaining objects.

WR A (WN7b+OB?, P _orb = 7.62d binary) has a flat radio spectrum (α _radio ∼ 0.0) that is best interpreted as a composite of thermal (wind) and non-thermal (WCR) components (Do10), although this is a slightly surprising result given that the short orbital period should place the WCR well within the radio photosphere. The mm-radio spectral index, α_mm ∼ 0.85 is unambiguously thermal, suggesting that the stellar wind dominates at mm wavelengths, although this value is slightly steeper than expected for a canonical stellar wind. Following the discussion in Sect. 3.1.2 this could be the result of a highly structured wind or optically-thin emission from the WCR. The latter would be expected for a short period binary system but would appear difficult to reconcile with the apparent non-thermal emission component at radio wavelengths. The canonical α = 0.6 spectral index assumes the wind to be at terminal velocity. The spectral index can however be non-linear if the wind is accelerating, with a steeper spectral index at higher frequencies. For WR A (and the other WR stars), the characteristic radius where the free-free optical depth is equal to 1 is ∼100 R_⋆ and is therefore significantly beyond the wind acceleration zone for a canonical ( β = 1) velocity law.

WR V (WN8o) likewise has a flat radio spectrum (α_radio ∼ 0.0) suggestive of a composite thermal+non-thermal origin although, unlike WR A, no corroborative evidence for binarity exists. As with WR A α_mm ∼ 0.46 suggests the increased dominance of the stellar wind at mm-wavelengths.

WR B (WN7b+OB?, P_orb = 3.52d binary) appears comparable to WR A. It too has a flat radio spectrum (α_radio ∼ 0.04); however when considering the total flux of the mm source it has a moderately negative spectral index (α_mm ∼ −0.10). Examination of both radio and mm data shows that it is significantly extended at both wavelengths (Fig. 5 and Table 3). Following the reasoning in Do10, we consequently refrain from interpreting the emission beyond noting that the mm continuum flux is comparable to other cluster WRs, suggesting that it is likely dominated by emission from the stellar wind. Finally we note the twin “lobes” of emission equidistant from the star and apparently aligned upon a single axis (Fig. 5). While this is suggestive of the jet lobes that are sometimes associated with high mass X-ray binaries such as SS433, the double eclipse in the light-curve clearly indicates a normal stellar companion (Bonanos 2007).

WR L (WN9h+OB?, P_orb = 54d binary) has a radio spectrum (α_radio > 0.5) consistent with the canonical value for a partially optically-thin stellar wind, although current limits would also accommodate an additional contribution from a non-thermal component. The mm-radio spectral index (α_mm ∼ 0.90) is directly comparable to that of WR A and so similar conclusions apply here.

WR F (WC9d+OB?, P_orb = 5.05d binary) has weaker constraints on the radio spectrum (α_radio > 0.0) but a mm-radio spectral index, α_mm ∼ 0.64 entirely consistent with emission from a stellar wind.

WR S (WN10-11h/BHG) is thought to be a single star (Clark et al. 2014b) and the radio spectrum (α_radio > 0.0) is consistent with such an hypothesis however, the mm-radio spectral index (α _mm ∼ 0.22) is unexpectedly flat for such a scenario. We discuss this below.

The final 15 stars only have lower limits for α_mm due to a lack of radio detections (Table 2). None of these demonstrate spectral indices that imply purely non-thermal emission. Examining the remaining stars by spectral sub-type and all five of the WN5 and WN6 stars have mm-radio spectra consistent with thermal emission from a stellar wind, but which could also accommodate contributions from optically-thin thermal or non-thermal emission; in this regard we note that binarity is suggested for WR O, R and U on the basis of their X-ray properties. Of the WN7 stars, WR D, G and P, all appear to have emission dominated by the stellar wind, despite binarity being suggested for WR D and G. As with WR L, lower limits to the spectral index of WR I (WN8) are higher than expected for a canonical stellar wind, while the WN11h star Wd1-13 replicates the properties of, and conclusions drawn for, the WN5 and WN6 stars.

Despite it being likely that all four WC9 stars (WR C, E, H and M) are binaries, we find no compelling evidence for a non-thermal emission component in their spectra, with a similar conclusion drawn for the apparently single WC8 star. Unfortunately, the lack of period determinations for the WC9 cohort precludes us commenting on the likelihood of an additional optically-thick continuum contribution from the WCRs (cf. Pittard 2010, Sect. 3.1.2).

In summary, we find the mm-radio spectral indices of all the WRs to be consistent with partially optically-thin thermal emission, with no compelling a priori reason to suspect a significant contribution from a source other than the stellar wind for any star, with the exception of WR V and S. The non-thermally emitting fraction derived from mm observations is therefore clearly lower than derived from extant radio surveys (Sect. 4.1), which may be expected given the greater flux expected from the stellar wind at mm wavelengths due to the ν^+0.6 dependence.

4.2.2. Mass-loss rates

As a consequence of the above discussion we utilise the mm-fluxes in order to infer mass-loss rates for all detections following the methods outlined in Sect. 3.1.3. We present the results in Table 4 and Figs. 6 and 12, where we also compare our values to previous radio surveys (Sect. 4.1). While we defer individually tailored non-LTE model-atmosphere analyses for a future work, extant modelling for both WR F and S (Clark et al. 2011, 2014b) shows an encouraging consistency in mass-loss rates to within a factor of ∼2. Utilising solely the errors in flux density we calculate example ranges for the mass-loss rates for four representative sources. WR O, E, D and U have ranges of 2.90–3.56 × 10⁻⁵, 5.09–5.67 × 10⁻⁵, 1.89–2.18 × 10⁻⁵ and 2.11 − 2.97 × 10⁻⁵ M_⊙ yr⁻¹ respectively, implying errors of the order ±0.2 dex in Fig. 6. The mm mass-loss rates are compared directly with those calculated from radio measurements as published in the literature (see Sect. 4.1). Utilising directly the radio flux densities and re-calculating the radio mass-loss rates assuming the same chemistry, v_∞ and T_eff, as for the mm mass-loss rates, results in a change in the range of ∼0.05–0.3 × 10⁻⁵ M_⊙ yr⁻¹ for each object, with the WC mass-loss rates experiencing the smallest shift.

Our survey has a number of advantages over those performed previously. Specifically, our observations detect a much higher percentage of sources (100% versus 44–65%; Sect. 4.1). As a consequence we can be confident that the range of mass-loss reported accurately reflects the underlying distribution, whereas values from previous surveys will be biased due to the non-detection of fainter objects supporting lower mass-loss rates. Moreover, all of our objects are co-located, minimising uncertainties in the distance estimate that adversely affect surveys of isolated field objects. Additionally, accurate determination of cluster age and hence progenitor masses for the WRs greatly improves the utility of our mass-loss rates when employed to test theoretical stellar evolutionary predictions and also gives confidence that scatter in the values observed is intrinsic to the stars themselves.

It is immediately obvious that, with the exceptions of the WN10-11h stars WR S and Wd1-13, mass-loss rates for cluster WRs span a limited range of log(Ṁ) ~ −4.1 → −4.8 or ∼1.6–7.4 ×10⁻⁵ M_⊙ yr⁻¹, a scatter which, following the proceeding discussion, is at least in part likely intrinsic. As is apparent from Fig. 5, mean mass-loss rates for the whole ensemble (Ṁ ~ 3.6 ± 0.4 × 10⁻⁵ M_⊙yr⁻¹; excluding WR S and Wd1-13) and both WN (Ṁ ~ 3.4 ± 0.5 × 10⁻⁵ M_⊙yr⁻¹) and WC cohorts (Ṁ ~ 4.2 ± 0.7 × 10⁻⁵ M_⊙yr⁻¹) are entirely consistent with previous radio determinations (Sect. 4.1) as well as values derived from spectroscopic model-atmosphere analysis (Hamann et al. 2006; Sander et al. 2012). Strikingly, we find no evidence of the systematically lower mass-loss rate for the WCL stars that had been previously supposed (Sect. 4.1).

The WNVLH/BHG hybrids WR S and Wd1-13 are (past) members of close binaries for which we hypothesise that binary interaction has resulted in the stripping of their outer layers (Ritchie et al. 2010; Clark et al. 2014b). As such it is not immediately apparent that they are directly comparable to the wider population, appearing intermediate between objects such as the WN9h star WR L and non-cluster early-B hypergiants such as ζ¹ Sco. WR L itself has a mass-loss rate in excess of those reported for the lower-mass WN9h stars within the Galactic Centre cluster, but which is directly comparable to those of other WRs within Wd1.

Finally in Sect. 3.1.2 we highlighted the potential blending on WR B with extended emission and the unexpectedly flat mm-radio spectra of WR S and V. Despite this, in each case we find that our mass-loss rate estimate is comparable to those other members of the same spectral sub-type, supporting our implicit assumption that the mm-continuum emission in each case remains dominated by the stellar wind.

4.2.3. Spatially resolved emission

As reported in Sect. 3.1.1 Gaussian fitting to the sources implies that excluding WR C, J and U the emission associated with the remaining 18 WRs appears spatially extended (Table 3; Fig. 5)^4,5. As expected the magnitude of the associated errors broadly correlates with source luminosity; as a consequence we cannot comment on any possible relation between source extent and luminosity, WR sub-type or binary status with confidence. Despite this we emphasise that the errors associated with the more luminous sources are sufficiently small that the conclusion that a subset of WRs are indeed resolved appears robust, with e.g. WR A, E, F and L resolved at a significance of ∼5.3σ, ∼3.6σ, ∼4σ and ∼3.3σ respectively.

As described in Sect. 3.1.1 an additional emission component, visible as a “shoulder” in the wing of the Gaussian profile, appears present in a number of sources (e.g. WRs J, K and M; see Fig. 7) in a similar manner to the cool super-/hypergiants, although at a lower significance given the fainter nature of the WRs. As a consequence we refrain from quantitative modelling of this feature, although we speculate that these deviations imply a composite structure for these sources, with a bright core and surrounding low surface-brightness halo. Pre-empting Sect. 5.2.2, possible corroboration of this hypothesis is provided by the structure of the central component of the RSG Wd1-20, which is sufficiently resolved to directly observe a similar potential core + halo structure in the image and its corresponding annular profile (also shown in Fig. 7), which displays a similar (albeit more pronounced) deviation from a Gaussian morphology.

Utilising the mm-fluxes and adopted wind and stellar parameters of individual stars (Sect. 3.1.3) we may adopt the formalism of Wright & Barlow (1975) to infer the radius of the mm-photosphere of each object. These range from 0.1 (WR K) to 3.8 (WR L) milliarcsec; even for WR L this is a factor of ∼20× smaller than the minor axis (Table 3). Comparison of such an analytic estimate to one derived from non-LTE model-atmosphere analysis for WR S (Clark et al. 2014b) shows coincidence to within a factor of ∼2. As a consequence we may reject the hypothesis that we are resolving the stellar wind in any of these sources. Likewise, while (field) WR stars are commonly associated with wind-blown bubbles, such structures are orders of magnitude larger than we observe here. We defer further discussion and interpretation of these phenomena until Sect. 8.2.

5. The stellar sources II. The cool super- and hyper-giants

We next turn to the YHGs and RSGs within Wd1. Such stars are of particular interest since, despite the brevity of the phase, it is thought that they shed mass at sufficient rates to profoundly affect evolutionary pathways. However current models (e.g. Ekström et al. 2012) rely on sparse empirical constraints for mass-loss rates in this phase (de Jager et al. 1988), which often have to be inferred via secondary diagnostics (e.g. via dust emission which then requires the adoption of an uncertain dust:gas ratio). Wd1 has the potential to be a critical test-bed for understanding mass-loss from stars for three reasons. Firstly it is sufficiently massive that it contains the richest population of RSGs and YHGs known in the Galaxy; moreover we find it at an age (∼5 Myr) at which such stars first occur and, as a consequence, these are also likely to be the most extreme objects (M_initial ∼ 35–40 M_⊙, log(L/L_⊙)∼5.8) permitted by nature at solar metalicities (Ekström et al. 2012).

Radio observations by Do10 reveal the presence of ionised circumstellar gas around a number of cluster members, with recent, deeper, observations detecting and resolving such ejecta around nine of the ten YHGs and RSGs (missing only Wd-8; Andrews et al. 2018). Given their photospheric temperatures are insufficient to ionise this material this leaves the presence of a hot companion or the diffuse cluster radiation field as possible physical agents (cf. discussion in Do10). More importantly this opens the possibility of directly determining time averaged mass-loss rates for such systems. Being more sensitive to extended structure, the radio observations of Do10 provide more accurate constraints on the masses of the circumstellar nebulae while, with their enhanced spatial resolution, the ALMA observations provide complementary information on the geometry; both prerequisites if mass-loss estimates are to be determined. Unfortunately, we cannot easily combine both datasets to obtain spectral index information, but Do10 conclude that the spatially extended emission associated with each star is consistent with an optically-thin thermal origin.

5.1. The yellow hypergiants

We detected and resolved four of the five cluster YHGs within the field of view of our ALMA observations – Wd1-4, 12a, 16a and 32 (Table 2 and Figs. B.1 and B.3), with Wd1-8 remaining undetected at both mm and radio wavelengths.

The nebulae associated with Wd1-4 (F2 Ia⁺), Wd1-12a (A5 Ia⁺) and W16a (A2 Ia⁺) are clearly asymmetric at 3-mm, comprising a dominant (and resolved) point-like source embedded at the apex of more extended, trailing nebulosity, with an arc- or arrow-head like morphology (Figs. 8, and Table 3). These sources are in turn located at the head of elongated radio nebula of greater extent, with an overall cometary morphology comprising a resolved, compact nucleus coincident with the YHG and an extended tail. While the nebulosity associated with Wd1-32 is too compact to reveal such a configuration, the radio continuum observations of Wd1-265 by Andrews et al. (2018) show that it, too, is associated with a spectacular cometary nebula. While such behaviour has never before been observed for YHGs it is reminiscent of the nebulae associated with the RSGs NML Cyg and GC IRS7 (with the latter being more elongated), which are thought to be shaped via an interaction with the nearby Cyg OB2 association and Galactic Centre cluster respectively (Schuster et al. 2006; Serabyn et al. 1991; Yusef-Zadeh & Morris 1991). Given that the cometary tails of the Wd1 YHGs are also all orientated away from the cluster core we suggest that a similar physical process is in operation here; we return to this below.

What is the origin of the circumstellar material? In interpreting both radio and ALMA data it is important to recognise that individual cluster YHGs may be in either a pre- or post-RSG phase and hence may have very different physical properties (e.g. L/M ratio, Ṁ, etc.). Nevertheless the comparable angular sizes and fluxes (Tables 2 and 3) of the extended emission components and, where available, the radio-determined nebular masses⁶ point to a common physical origin. Ejection nebulae have been associated with the post-RSG YHGs IRC +10 420 (Tiffany et al. 2010; Shenoy et al. 2016) and IRAS 17163-3907 (Lagadec et al. 2011; Hutsemékers et al. 2013; Wallström et al. 2017) but are physically more extended and contain orders-of-magnitude more mass than those considered here. Likewise, while the YHGs within Wd1 exhibit pulsations (Clark et al. 2010), none have (yet) demonstrated the long-term secular evolution or giant eruptions that characterise post-RSG examples such as ρ Cas, IRC +10 420 and HR8752 (Lobel et al. 2003; Nieuwenhuijzen et al. 2012; Oudmaijer et al. 1996; Oudmaijer 1998). Finally, none show the rich emission spectra of both IRC +10 420 and IRAS 18357-0604, which appear to be caused by extreme post-RSG mass-loss rates (Oudmaijer 1998; Clark et al. 2014a); instead they share a continuous morphological sequence with the mid-late BHGs within Wd1 (Clark et al. 2005) which we argue are in a pre-RSG phase (Sect. 5.2).

Given the above, if the majority of Wd1 YHGs are indeed in a pre-RSG phase and therefore are not yet prone to the large scale instabilities that characterise post-RSG stars, then the circumstellar material appears most likely to result from the accumulation of a quiescent stellar wind. Adopting a wind velocity of ∼200 km s⁻¹, we infer a wind-mediated mass-loss rate of Ṁ ~ 10⁻⁵ M_⊙ yr⁻¹ for the thermal component of Wd1-4a (when scaled to 5 kpc; Do10). For a freely expanding wind this would correspond to a dynamical age of ∼10² yr for the nebulae and hence an accumulation of ∼10⁻³ M_⊙ over this period. Given the nebular morphologies of Wd1-4, 12a 16a and 265 clearly suggest wind confinement in the hemisphere nearest the cluster core this is better interpreted as a lower limit to nebular mass; indeed estimates from Do10 and Andrews et al. (2018) are in excess of this value.

5.2. The red supergiants

The ALMA observations resolved all four RSGs within Wd1, including Wd1-75 for the first time (Figs. 9–11 and Tables 2, 3). To within observational uncertainties, all extended nebular components have radio spectral indices consistent with optically-thin thermal emission (Do10). Intriguingly, their nebular masses are larger than inferred for the YHGs (Do10; Andrews et al. 2018), suggesting either a longer duration for their mass loss or that a differing physical process may be responsible for their formation. Determination of current mass-loss rates via mid-IR fluxes is precluded by the saturation of the stars in all available flux-calibrated datasets⁷. Nevertheless the detection of SiO and H₂O maser emission associated with both Wd1-26 and −237 argues for very high mass-loss rates which are thought to occur during the latter stages of the RSG phase, when stars have the most extreme L/M ratio (Davies et al. 2008; Fok et al. 2012).

Observations of (predominantly lower-luminosity) field RSGs show that, when present, their circumstellar nebulae show a diverse set of properties. Few determinations of nebular masses are present in the literature, but the nebulae associated with Wd1-20, Wd1-26 and Wd1-237 appear broadly comparable to that surrounding the extreme RSG VY CMa, a star for which transient mass-loss rates in excess of ∼10⁻³ M_⊙ yr⁻¹ have been inferred (Shenoy et al. 2016; Smith et al. 2001, 2009a). Nebular geometries range from quasi-spherical through to elongated (VX Sgr and S Per respectively; Schuster et al. 2006), highly aspherical and clumpy (VY CMa; Smith et al. 2001) and cometary (GC IRS7; Serabyn et al. 1991; Yusef-Zadeh & Morris 1991). A number show evidence for interaction with their immediate environment (e.g. GC IRS7 and Betelgeuse; Noriega-Crespo et al. 1997), while others provide evidence for time-variable mass-loss rates (e.g. VY CMa and μ Cep; Shenoy et al. 2016). Given the physical information encoded in the nebular geometries we chose to group and discuss the cluster RSGs on this basis.

5.2.1. Wd1-75 and Wd1-237

At mm-wavelengths Wd1-75 presents as a compact, elongated source (Table 3 and Figs. B.1 and B.3). In terms of spatial extent it is similar to the approximately spherical nebula detected around μ Cep in the mid-IR by de Wit et al. (2008) ⁸, although additional structure is observed at both larger and smaller radii at different wavelengths around the latter star (Schuster et al. 2006; Shenoy et al. 2016; Cox et al. 2012). The new radio continuum observations of Andrews et al. (2018) show that this emission is embedded within a more extended nebula; however given the compact nature of the nebulosity at both wavelengths no conclusions regarding morphology may be drawn.

Located in the southern extremities of the cluster, at mm-wavelengths Wd1-237 appears similar to Wd1-75, albeit with a slightly greater extent (Table 3). However radio observations suggest a much more extended circumstellar nebula enveloping this structure (∼11.2 × 8.5 arcsec; Do10). The observations of (Andrews et al. 2018, also Fig. 9 here) clearly resolve this emission, revealing a compact central nebula co-incident with the ALMA source and star itself, which is offset from the centre of a larger, quasi-spherical nebula which appears brighter on the hemisphere facing the cluster core. This morphology is strikingly reminiscent of the mid-IR nebula associated with μ Cep (Shenoy et al. 2016). As with the YHGs this configuration immediately suggests interaction with the cluster proper. Such nested configurations have been associated with other RSGs (cf. preceding discussion in Sect. 5.2) and have historically been interpreted as arising from variations in the mass-loss rate of the star. While this is an obvious hypothesis for Wd1-237 it is not the only explanation; an issue we return to below.

Given the apparent sculpting of the nebula by the cluster, a robust determination of the dynamical age will require detailed hydrodynamical simulations. Nevertheless simply taking the displacement of the central core from the outer arc of emission facing the cluster core (∼0.12 pc for a distance of 5 kpc) and an outflow velocity of 10 km s⁻¹ suggests a minimum age of ∼10⁴ yr, which in turn would imply a time-averaged mass loss rate of ∼10⁻⁵ M_⊙yr⁻¹ via the nebular mass (∼0.07 M_⊙; Do10, Andrews et al. 2018). We regard these as at best order of magnitude estimates given uncertainties in nebular age, outflow velocity and the possibility of an additional neutral component to the nebula which would be invisible to current observations.

5.2.2. Wd1-20 and Wd1-26

Radio observations of the two remaining cluster RSGs have already revealed them to be associated with pronounced cometary nebula (Do10). Such linearly extended structures have previously only been associated with the Galactic Centre star GC IRS7, although bow shocks have been detected around a further three RSGs; Betelgeuse (Noriega-Crespo et al. 1997), μ Cep (e.g. Cox et al. 2012) and IRC-10414 (Gvaramadze et al. 2014). Clearly however, they bear close resemblance to the nebulae of the cluster YHGs, while that associated with Wd1-237 appears to be a less collimated analogue.

We first turn to Wd1-26. Analysis of optical spectroscopy reveals the gaseous element of the nebula to be composed of N-enriched material (Mackey et al. 2015), while mid-IR imaging reveals a co-spatial dusty component (Clark et al. 1998). The Hα image presented by Wright et al. (2014); over-plotted on ALMA data in Fig. B.1) appears to show Wd1-26 at the centre of a clumpy ring nebula. A second, triangular nebula is seen to the north-east, with both structures potentially connected by a thin, faint bridge of emission. Our ALMA observations clearly replicate this latter feature, implying a physical connection between both components. While the southern nebular component is obviously clumpy, the higher spatial resolution afforded by ALMA (and the lack of contamination by background stars that afflicts the Hα image) suggests that Wd1-26 does not sit at the centre of a ring nebula. Rather it is located at the apex of one of two bright elongated blobs of emission that run parallel to one another. These are also aligned with the major axis of the nebula, which in turn is directed towards the cluster core. Leading these two sub-components are a further two emission hot-spots, with all four components in turn embedded within more extended lower surface-brightness emission. Comparison to the 2-cm radio image of GC IRS7 (Yusef-Zadeh & Morris 1991) reveals a striking similarity, with both nebulae demonstrating “comet-like” morphologies of comparable sizes (∼0.38 pc for Wd1-26 versus ∼1 pc for GC IRS 7), with significant substructure visible in both (by analogy) “coma” and “tail” components.

As with Wd1-26, radio observations of Wd1-20 reveal the nebular tail to be orientated away from the core of the cluster and, at shorter wavelengths, it too appears to show similar substructure. Likewise both mass and spatial extent (Table 3; Do10 and Andrews et al. 2018) are comparable. While the nebular tail is not detected as a single coherent structure in the ALMA data, isolated hot-spots coincident with features in the radio data are seen; however these shortcomings are compensated for by the detailed view afforded of the immediate environment of the RSG itself (Fig. 10). Specifically a compact aspherical nebula (∼0.02 × ∼0.01 pc; Table 3) is coincident with the star. To the north of this is a further arc of emission ∼0.14 pc in extent with a projected separation between the apex of this structure and the RSG of ∼0.036 pc. We note that this feature appears absent in the nebula surrounding Wd1-26, it is not clear whether this corresponds to a real physical difference or instead is due to differing orientations/lines-of-sight through otherwise identical nebulae.

An obvious interpretation for the arc of emission is that it is a bow-shock, with the nebular morphologies of both Wd1-20 and Wd1-26 being shaped by interaction with the wider cluster environment⁹. Potential physical agents for this include (i) relative motion of the RSGs through the intra-cluster medium/wind (detected at both radio and X-ray wavelengths; Do10, Muno et al. 2006b), (ii) photoionisation by the hot massive stellar cohort and/or (iii) interaction with the stellar winds of individual cluster members and/or a recent supernova(e).

Of these, motion of the stars relative to the cluster as a whole may be disfavoured by the observation that the cometary nebulae associated with both RSGs and YHGs are all orientated away from the core region which, under such a scenario, would require all four objects to be falling towards this region. Likewise one might anticipate that velocities significantly higher than the cluster virial velocity would be required for bow shocks to form via interaction with a static intra-cluster medium; essentially all such objects would be “run-towards” rather than the more normal “runaways”. As a consequence we favour an interaction between the star and a dynamic stellar/cluster wind(s) as the most likely causal agent for the cometary nebulae.

Mackey et al. (2014, 2015) investigated the shaping of RSG winds by external agents providing, respectively, detailed simulations of the bow-shock around Betelgeuse and the nebula of Wd1-26. The latter study suggests ram confinement of the stellar wind by the nearby sgB[e] star Wd1-9 as a potential physical origin for the cometary morphology. More sophisticated simulations (e.g. incorporating more detailed models of the cluster environment) would be invaluable in order to (i) further test this hypothesis, (ii) constrain the mass-loss rates Wd1-26 (and Wd1-20) and (iii) determine whether the “clumpy” sub-structure is a result of aspherical mass-loss (cf. VY CMa Smith et al. 2001) or post-ejection hydrodynamical instabilities (cf. Cox et al. 2012).

Nevertheless, one important insight arising from these simulations is the prediction of a massive static shell interior to the bow-shock. This forms via the stalling of the neutral RSG wind as a result of external photoionisation and as a consequence of this the accumulation of mass lost by the star. In this regard the observation of such compact nebulae surrounding both Wd1-20 and −237 is of particular interest.

6. The stellar sources III. OB stars – supergiants, hypergiants, LBVs and sgB[e] stars

The last cohort we turn to comprise the hot, early-type stars within Wd1; a somewhat more heterogeneous grouping than those previously considered. We detect a total of 21 mm-continuum sources coincident with cluster members (Tables 2–3 and Figs. B.1 and B.3). Subject to the classification of two objects (D09-R1 & R2)¹⁰ detections range from O9 Iab (Wd1-25) to B9 Ia⁺ (Wd1-42a) plus the sgB[e] star Wd1-9 and the LBV Wd1-243. Despite their presence in large numbers, no stars of lower luminosity class were detected (Clark et al., in prep.). In this regard we note that based on radio data (Do10), the O9 Ib star Wd1-15 should have been detected but was not, implying that it is variable – a potential signature of non-thermal emission.

Three further stars – D09-R1 and R2 (OB SGs) and Wd1-17 (O9 Iab) – have radio properties suggestive of non-thermal emission (Do10). Our mm-data are consistent with such a conclusion for Wd1-17 but imply flat or moderately positive mm-radio spectra for D09-R1 and R2 (Table 2). However, the radio and mm point sources associated with each star are further coincident with extended continuum emission; consequently determination of absolute fluxes and hence the physical nature of the emission is uncertain. While we currently favour non-thermal emission for D09-R2 and Wd1-17 – and hence identification, along with Wd1-15, as CWBs – no corroborative evidence for such a classification is available at optical or X-ray wavelengths (Table 2; Cl08, Ritchie et al., in prep.).

Of the remaining objects the emission associated with the LBV Wd1-243 and the central component of Wd1-9 is entirely consistent with arising in a partially optically thick stellar winds. Lower limits to the radio-mm spectral indices for Wd1-21, 8b, 11, 52, 56a and 61a do not constrain the emission mechanism, while those for Wd1-7, 23, 25, 28, 30, 33, 42a, 43a, 46a and 71 are consistent with partially optically thick thermal emission, though the presence of an additional optically thin and/or non-thermal component may not be excluded. For the purposes of this paper we proceed under the assumption that these are indeed thermal wind sources for which mass-loss rates may be determined (Table 4 and Fig. 12).

Finally, we find a number of sources, including both super- and hypergiants, to be partially resolved (Table 3), but with an extent greater than expected on the basis of mass-loss rates derived below; an exemplar being the LBV Wd1-243 (0.15 ± 0.01 arcsec); we discuss this phenomenon in Sect. 8.2.

6.1. The OB super- and hypergiants

Post-MS late-O and B stars within Wd1 follow a smooth progression in spectral morphologies leading to a corresponding evolutionary passage from O9 III through O9-9.5 Iab to B0-4 Ia (Negueruela et al. 2010; Clark et al. 2015, and in prep.). The status of the B5-9 Ia⁺ stars Wd1-7, 33 and 42a is less clear-cut since such classifications have also been applied to post-RSG objects such as HD168607 and HD168625. The simulations of Groh et al. (2014) and Martins & Palacios (2017), suggest massive (> 40 M _⊙) late O/early B supergiants evolve directly into cooler hypergiants (consistent with the empirical results of Clark et al. 2012). Alongside their close spectral similarity to the cluster B supergiants this does however argue for these objects to be a direct extension of this O9III – B0-4 evolutionary sequence, and consequently for all these stars to be in a pre-RSG phase.

Accurate classifications are unavailable for D09-R1 and -R2 and, while the spectral variability and X-ray properties of Wd1-30a reveal that it is a massive O+O star interacting binary, determination of physical parameters is non trivial (Clark et al., in prep.). Following the discussion above both Wd1-15 and Wd1-17 appear likely non-thermal colliding wind sources. No observational constraints on the binary properties of either system are available and consequently it is not possible to disentangle contributions from the components of either system (stars and wind collision zone). Excluding these stars, as well as the LBV Wd1-243 and the sgB[e] star Wd1-9, which we discuss later, leaves a total of 15 cluster members with modern classifications for which mass-loss rate determinations may be made.

Of this subset, three stars, Wd1-23a (B2 Ia + B Ia?), 43a (B0 Ia + ?) and 52 (B1.5 Ia + ?) (Negueruela et al. 2010), show compelling evidence for binarity¹¹. Conversely, a number of known and suspected massive OB binaries are not detected within the cluster (e.g. Wd1-27 and 53a; Bonanos 2007; Clark et al. 1998) suggesting that their properties do not uniformly favour discovery at mm-wavelengths (cf. Sect. 3.1.2). Indeed, the mass-loss rate inferred for Wd1-52 is entirely comparable to other (apparently single) B1.5Ia stars and those of Wd1-23a and 43a to the wider supergiant cohort and ∼0.6 dex lower than the apparent CWB Wd1-17 (O9.5Ia) if we were to make the questionable assumption that its 3mm flux derives solely from wind emission.

Mass-loss rate estimates for this subset are given in Table 4, with upper-limits to the continuum-flux and mass-loss rates for non-detections in Table A.2; these are presented graphically in Fig. 12¹². We supplement these data with the inclusion of the WNVLh/BHG hybrids Wd1-5 and −13¹³.

The main finding of our study is the preference for detecting super-/hyper-giants of spectral type B1.5 and lower (11 from 13 including the B1-1.5 Ia star Wd1-56a) versus earlier spectral types (4 from > 60 stars), with no cluster O giants detected. Mass-loss rates vary from ∼1.5–4.0 × 10⁻⁶ M_⊙ yr⁻¹ for the B1.5-9 super-/hypergiants and ∼4.0–8.3 ×10⁻⁵ M_⊙ yr⁻¹ for the earlier supergiants. However, given the upper limits inferred from the large population of undetected supergiants of early spectral type we are currently unable to conclude that their mass loss rates are systematically greater than their descendents of later spectral types. Indeed, inspection of Fig. 12 suggests that the ≤B1.5 supergiants detected represent the extreme tail of a wider distribution of mass-loss rates for such stars; taken as a whole both detections and non-detections suggest a significant intrinsic scatter in the mass-loss rates of stars of a given spectral type (e.g. ∼0.4 dex and ∼0.6 dex for B0 and B1 supergiants respectively).

How do these results compare to other studies? Despite an absence of mm-continuum observations of OB stars in the literature, a number of investigations have been undertaken at radio wavelengths (Benaglia et al. 2007; Bieging et al. 1989; Lamers & Leitherer 1993; Leitherer et al. 1995; Puls et al. 2006; Scuderi et al. 1998). After exclusion of non-thermal sources (Benaglia et al. 2007) provide mass-loss rates for 19 O8-B5 supergiants which we reproduce in Fig. 12, supplemented with the late-type BHG HD160529 (Leitherer et al. 1995); our ALMA observations compare very favourably to these in terms of both numbers of detections and distributions of spectral types. Upon consideration of both datasets it is apparent that the mass-loss rates derived by both studies are broadly comparable. However non-LTE model-atmosphere analysis of spectroscopic data for field B0-9 super-/hypergiants (Crowther et al. 2006b; Searle et al. 2008) suggest that the majority of stars so analysed support lower mass-loss rates than those derived from these radio and mm continuum observations. Fortuitously seven stars have both radio continuum and spectroscopic mass-loss rate determinations¹⁴ which reveal that the radio detected stars comprise the subset of objects with the highest spectroscopic mass-loss estimates. This implies that the radio detections reported in Benaglia et al. (2007) are biased towards the high radio luminosity subset of a larger population of predominantly fainter objects with lower mass-loss rates, as we infer for our ALMA sample.

How, then, are these observations to be interpreted? A synthesis of the spectroscopic studies suggests that supergiants of later spectral type support lower mass-loss rates (cf. Fig. 9 of Clark et al. 2014b), but this likely reflects the systematically lower luminosities (and hence masses) of the (field) stars considered; this bias would not be expected to be present in the co-eval population of Wd1. Indeed a more complicated relationship between mass-loss rate and terminal wind velocity and spectral type/stellar temperature might be expected as a result of the recombination of Fe IV to Fe III – the so called “bi-stability jump” – around 21kK/spectral type B1 (Pauldrach & Puls 1990; Vink et al. 1999, 2000, 2001; Benaglia et al. 2007). Lamers et al. (1995) show that the ratio of terminal-wind to escape-velocity v _∞/v_escape) drops from ∼2.6 to ∼1.9 around spectral type B1. Given the ${\upsilon }_{\infty }^{-4/3}$ dependence of mm-continuum flux (Eq. (1)); in the absence of any other changes to wind properties one would therefore anticipate that it would be easier to detect stars of later spectral type, as appears to be borne out by our observations. In practice one might also anticipate changes in mass-loss rate and ionisation structure of the wind with spectral type; however quantitative analysis of such affects must await tailored non-LTE analysis of individual stars.

6.2. The LBV Wd1-243

The LBV Wd1-243 is a strong, resolved source in our ALMA data. Spectroscopic observations show that it remained in a long-lived cool-phase between 2002 and 2009 (approximate spectral type of A3 Ia⁺; Ritchie et al. 2009b), with unpublished data suggesting this state persisted until 2015, despite the continued presence of low level photospheric pulsations. The mm-radio spectral index suggests a thermal spectrum with mass-loss rates broadly comparable between both radio (log(Ṁ/M_⊙ yr⁻¹) ∼ − 5.0; data obtained between 1998–2002) and mm determinations (log(Ṁ/M_⊙ yr⁻¹) ∼ − 4.6; data obtained in 2015). Intriguingly however, both values are an order-of-magnitude greater than that determined from non-LTE model-atmosphere analysis (log(Ṁ/M_⊙ yr⁻¹) ∼ − 6.1; Ritchie et al. 2009b), with the caveat that the latter analysis was unable to simultaneously replicate both the prominent H I and He I emission features and photospheric lines evident in the spectra.

Even outside giant eruptions, LBV mass-loss rates are clearly variable (Clark et al. 2009; Groh et al. 2009), potentially providing an explanation the moderate discrepancies between the radio and mm-continuum values. In any event both values are fully consistent with the range of values quoted for LBVs in the literature (Clark et al. 2014b, and refs. therein) and in particular the radio derived value for the post-RSG LBV HD160529 (nominal classification of B8 Ia⁺ and log(Ṁ/M_⊙ yr⁻¹) ∼ − 4.87 Leitherer et al. 1995).

6.3. The sgB[e] star Wd1-9

While ALMA observations of Wd1-9 have been discussed in detail in Fenech et al. (2017) the improved sensitivity of the mosaic images presented here reveals that the two filaments of emission extending to the south-east and south-west of the source, which were identified with low S/N in Fenech et al. (2017) are much more clearly seen (Fig. 13). This strengthens the hypothesis that Wd1-9 shares a comparable nebular morphology to the B[e] star MWC349 (White & Becker 1985; Gvaramadze & Menten 2012). However, in the case of Wd1-9 we don’t see filamentary structure at mm-wavelengths to the north of the nebular. Intriguingly, initial results of an analysis of Very Large Telescope Interferometer/MIDI data suggest the presence of compact, elongated mid-IR continuum source aligned along the N-S axis of the radio and mm-nebula (Hummel et al., in prep.), with a size comparable to that of the circumbinary torus inferred from the IR spectral energy distribution of Wd1-9 (Clark et al. 2013). If interpreted as the circumbinary torus this would provide powerful corroboration of the hypothesis that the mm-continuum emission extended in the NS direction resides in the orbital plane of the binary, with the fainter filamentary structure resulting from the disc wind.

7. The magnetar CXO J164710.2-455216

Only a handful of out-bursting magnetars have been found to exhibit (transient, pulsed) radio emission (Halpern et al. 2005; Rea et al. 2012; Szary et al. 2015). Unlike normal pulsars they appear to have shallow or flat spectral indices, making detection at mm wavelengths viable, with Torne et al. (2015) demonstrating that the magnetar SGR 1734-2900 is visible up to 225 GHz despite a distance of 8.3 kpc. Although the origin of radio emission from magnetars is still unclear, both Rea et al. (2012) and Szary et al. (2015) predict that given the physical properties of CXO J164710.2-455216 (viz. a ratio of quiescent X-ray luminosity to spin down power greater than unity and its residual temperature) we should not expect to detect it at mm- or radio-wavelengths. Moreover CXO J164710.2-455216 has yet to exhibit the intense flares that were associated with the transient radio emission of SGR1900+14 and SGR1806-20 (Halpern et al. 2005, and Refs. therein). Unsurprisingly, therefore, our ALMA observations yield only an upper limit to the continuum flux of 65 μJy beam⁻¹.

8. Discussion

8.1. Post-MS mass-loss rates and stellar evolution

Following the discussion in Sect. 1 a key omission in our description of the physics of massive stars is an accurate determination of mass-loss rates as a function of initial mass, metallicity and age. These ALMA observations of Wd1 provide multiple detections of every post-MS phase and offer the potential of uniformly determined, unbiased mass-loss rate determinations for stars of M_init ∼ 30–50 M_⊙.

In the absence of tailored modelling of individual stars, detailed comparison of mass-loss rates to predictions from stellar-evolution models (e.g. Brott et al. 2011; Ekström et al. 2012; Martins & Palacios 2017) are premature, although we may make limited qualitative comparisons. Each of these codes utilises the theoretical mass-loss recipe of Vink et al. (2001) for OB stars (although Martins & Palacios 2017 scales this by a factor of a third to account for recent observational findings) and the empirical relation of Nugis & Lamers (2000) for field WRs, which together imply a significant increase in mass-loss rates as stars transition between these phases. Such behaviour is indeed observed for Wd1, with a simple comparison of Figs. 6 and 12 indicating that WRs span log(Ṁ/M_⊙ (yr⁻¹) ∼ − 4.1 → −4.8 and B super-/hypergiants log(Ṁ/M_⊙(yr⁻¹) ∼ − 5.1 → −5.8 (with the two WNVLh/BHG hybrids intermediate between both groups). However we caution that the large number of OB supergiants non-detections currently biases the sample to more extreme stars and so one might expect an extension of this range to lower mass-loss rates; deeper observations to detect the remaining cohort would be of great value. This is especially true since the range of spectral types exhibited by early super-/hypergiants within Wd1 (O9-B9) spans the predicted location (∼B1/21 kK) of the bi-stability jump due to the transition from Fe IV to Fe III as the dominant ion in the stellar wind (Sect. 6.1).

Mass-loss rates for the cool super-/hypergiants utilised by codes are again empirical in origin and for RSGs span multiple orders of magnitude (de Jager et al. 1988; Nieuwenhuijzen & de Jager 1990; Sylvester et al. 1998; van Loon et al. 1999). Do10 report log(Ṁ/M_⊙ (yr⁻¹) ~ −5 for the partially optically-thick thermal component associated with the YHG Wd1-4a and infer a similar time-averaged mass-loss rate from the nebular properties of the RSG Wd1-237. Such values are in excess of those exhibited by the OB supergiants within Wd1 and are comparable to the quiescent mass-loss rates of the handful of other similarly extreme RSGs (log(Ṁ/M_⊙ (yr⁻¹) ∼ − 5 (quiescent) → − 3 (outburst); Blöcker et al. 1999, 2001; Smith et al. 2001). However they are not as high as those assumed in stellar evolution codes which, for the Geneva models for a M_init ∼ 40 M_⊙ star, may exceed log(Ṁ/M_⊙ (yr⁻¹) ∼ − 4. However we caveat this with the observation that the mass-loss histories of such stars within Wd1 are clearly complex, as evidenced by their nebular morphologies, and further detailed modelling will be required to provide more accurate mass-loss rates. Moreover, observations to investigate the possibility of a reservoir of neutral material within the circumstellar environment, to which our observations would be insensitive, would be particularly valuable since the presence of such gas would render current estimates of time-averaged mass-loss rates from nebular properties as lower limits.

Finally we turn to the apparent binary Wd1-9. Current evolutionary codes do not incorporate binary interaction, but simulations of potentially comparative systems (Petrovic et al. 2005)¹⁵ suggest time averaged mass-loss rates of log(Ṁ/M_⊙ (yr⁻¹) ≥ −4 over the 4–6 ×10⁴ yr period of fast case-A mass-transfer. The current mass-loss rate inferred for Wd1-9 (log(Ṁ/M_⊙ (yr⁻¹) ~ − 4.2Fenech et al. 2017) is indeed broadly consistent with such predictions. Moreover it is directly comparable to the highest mass-loss rates determined for the cluster WRs, greater than that of the “quiescent” LBV Wd1-243 and provisional estimates for the cool super-/hypergiants and most striking of all, a factor of ≳20× larger than found for the most “extreme” B super/hypergiants within the cluster. Placing Wd1-9 into such an evolutionary context therefore provides powerful confirmation of the hypothesis that binary interaction plays a central role in the lifecycle of massive stars.

Similar comparative studies of the evolution of mass-loss rates through distinct post-MS phases are possible for both the Arches and Galactic Centre clusters. Martins et al. (2007, 2008) present the results of a comprehensive model-atmosphere analysis of near-IR spectro-photometric data for members of both clusters, from which we find qualitative agreement with our results in the sense that mass-loss rates are also found to increase with evolutionary state. Specifically the mass-loss rates for the mid-O supergiants within the Arches cluster are approximately an order-of-magnitude lower than those of the more massive and evolved WN7-9h stars, with the mid-O hypergiants intermediate between both extremes (Martins et al. 2008). Likewise the mass-loss rate obtained from analysis of a mean spectrum derived from observations of ten late-O/early-B supergiants within the Galactic centre cluster is over a magnitude lower than the WNLh and WC groups (Martins et al. 2007).

Unfortunately, we may not quantitatively compare mass-loss rates between Wd1 and the two galactic centre clusters since results from the two different methods employed appear prone to a systematic discrepancy. Specifically, Martins et al. (2008) demonstrate that the clumping-corrected ($\dot{M}/f_{\text{cl}}^{0.5}$) spectro-photometric mass-loss estimates for the Arches cluster members are significantly larger than the radio determinations of Lang et al. (2005) for the subset of stars for which such a comparison is possible. We arrive at a similar conclusion after comparison of the results of Yusef-Zadeh & Morris (1991) and Martins et al. (2007) for the Galactic Centre cluster¹⁶.

Martins et al. (2008) suggests that the discrepancy may in part arise due to differences in the clumping fractions in the regions of the inner and outer wind responsible for, respectively, the near-IR spectral diagnostics and mm-/radio-continuum. Such an evolution in the run of clumping with radius has also been suggested by Puls et al. (2006) and Runacres & Owocki (2002); we return to this issue in Sect. 8.2.

8.2. Resolved emission associated with the OB super/hyper-giant and WR cohorts

As detailed in Sects. 4.2.3 and 6 a significant number of the hot stars within Wd1 (both WRs, and OB super-/hyper-giants) appear to be moderately extended in our continuum observations. Stellar wind sources have been spatially resolved at mm and radio wavelengths in the past; examples being the BHG Cyg OB2#12 (Morford et al. 2017) and the LBV P Cygni (Skinner et al. 1997), in both cases facilitated by the combination of relative proximity of the star, low terminal wind velocity and high mass-loss rate. However as discussed in Sect. 4.2.3, analytic and numerical estimates for the sizes of the mm-/radio photospheres of the WRs indicate they should not be resolved with the current observations, with an identical conclusion following for the OB super-/hypergiants.

High resolution radio observations of the CWB WR140 has resolved the WCR between both stellar components (Dougherty et al. 2005) and for WR146 both binary components and the WCR (Dougherty et al. 2000). However despite the distance to both being lower than to Wd1 and their orbital separation much greater than e.g. WR A and B, the resolved emission has a smaller angular scale than observed for the majority of stars within Wd1; we therefore conclude that it unlikely that we have resolved the (putative) WCRs of these objects.

Both binary driven mass-loss (cf. Wd1-9; Sect. 6.3) and impulsive ejection events associated with the LBV phase (e.g. Duncan & White 2002) can give rise to resolved ejection nebulae. However not all of the resolved stars within Wd1 currently show evidence for binarity while, to the best of our knowledge, LBV-like instabilities have not been associated with early OB supergiants nor, with only a couple of exceptions, with WN stars earlier than WN9 or WC stars. Moreover, given the similarity in size of the resolved sources (Table 3) one would require multiple stars to synchronise the timing and duration of any such interaction or instability in order to replicate this observational finding.

Similarly, we note that a number of resolved sources demonstrate mm-/radio spectral indices inconsistent with optically-thin thermal emission, as might be expected for ejection nebulae. Additionally, mass-loss rates derived for the sources are also fully consistent with expectations for WR and OB super-/hypergiants respectively (Sects. 4.2.2 and 6.1). Both observations suggest sources dominated by emission from stellar winds; it seems highly improbable that either binary or LBV ejection mechanisms could conspire to replicate the fluxes expected from the stellar winds of such stars as a function of their spectral type.

Given the results of profile fitting (Sects. 3.1.1 and 4.2.3), it appears possible that these sources are composite, with an unresolved core that dominates the emergent flux due to the stellar wind and an extended, low surface brightness halo surrounding this. If confirmed via higher-resolution observations, what would be the physical origin of such a configuration? After eliminating the above possibilities, and following our findings for the cool super-/hyper-giant cohort, an obvious explanation would be the confinement of the stellar winds leading to the formation of a compact wind-blown bubble.

Since, to the best of our knowledge, such a phenomenon has not been observed for isolated massive stars it is attractive to attribute this to their membership of a YMC and we highlight that Lang et al. (2005) report spatially resolved radio-continuum emission of comparable extent associated with a number of stars in the Quintuplet cluster. Prospective physical agents for this would include pressure confinement by the intra-cluster medium or wind/wind interaction with the stellar cohort and/or one or more recent SNe. Both X-ray (Muno et al. 2006b), mm and radio continuum observations (Do10; Sect. 8.3) argue for the presence of an intra-cluster medium or wind associated with Wd1. Mackey et al. (2015) conclude that pressure confinement is unlikely to explain the properties of the nebula surrounding Wd1-26; given the greater wind momenta of, for example, the WR cohort such a conclusion would appear to apply to such stars as well. This would appear to favour interaction between stellar and/or a cluster wind as a physical mechanism. However, given the highly structured nature of the intra-cluster medium/wind within Wd1 (as revealed by mm and radio continuum emission; Sect. 3) and the asymmetrical distribution of massive stars within Wd1, an assessment of both possibilities would required detailed hydrodynamical simulation which are clearly beyond the scope of this work.

8.3. Stellar and cluster wind interaction

It has long been understood that feedback from clusters via a stellar-wind- and SN-driven cluster-wind contributes to chemical evolution of the interstellar medium (Krumholz et al. 2014). However the properties and the mechanisms by which such winds are initiated and maintained are currently uncertain, the later due to uncertainties regarding the efficiency by which stellar kinetic energy is converted to thermal energy, which in turn drives the outflow (e.g. Stevens & Hartwell 2003). Previous observations of Wd1 have shown that the stars are embedded in a highly complex cluster medium, with both hot (Muno et al. 2006b) and cool components present (Do10). The resolution of our observations confirm that the cluster medium is highly structured, with multiple “clumps” of characteristic scale present in the core regions (Sect. 3.3), while also allowing us to study the interaction and entrainment of the winds of individual stars in unprecedented detail (Sects. 4–6).

In particular the pronounced cometary nebulae of the cool super-/hypergiant cohort and the compact “bubbles” associated with both OB super-/hypergiant and Wolf-Rayet stars imply that the circumstellar environments of stars embedded within young massive clusters may be profoundly different from their isolated brethren. One important consequence of this interaction may be in modifications to the phenomenology of the resultant SNe.

There is a wealth of observational evidence for a subset (∼10%) of SNe where the explosive ejecta interacts with pre-existing circumstellar ejecta/material (Smartt 2009). Although comparatively infrequent, such type IIn events can be highly luminous (e.g. Smith et al. 2007), potentially allowing their detection in the early Universe. Consequently the nature of the progenitors has been subject to considerable interest. The apparent requirement for the ejection of large amounts of circumstellar material in the past has led to the suggestion that LBVs may be the direct progenitors of such SNe (Smith et al. 2007; Gal-Yam & Leonard 2009). Multiwavelength studies by Fox et al. (2013) show a diversity of behaviour amongst type IIn SNe, suggesting a range of progenitor systems, while Smith et al. (2009a) suggested that extreme RSGs such as VY CMa may also be viable progenitor systems (see also Yoon & Cantiello 2010). Moreover spectroscopic signatures revealing the presence of circumstellar material has also been observed in a wide range of other SNe sub-types¹⁷.

However eruptions supporting greatly enhanced mass-loss rates just prior to SNe are not the only route for the production of dense circumstellar shells. The possibility of the retention of significant circumstellar material close to irradiated red supergiants in the absence of impulsive mass-ejection events has already been highlighted by Mackey et al. (2014, 2015); our observations corroborate this hypothesis, while potentially extending it to yellow hypergiants. Equally exciting is the identification of compact emission halos surrounding the WRs which, given the age of Wd1, are expected to be the immediate progenitors of core-collapse SNe. Flash spectroscopy of the type II event SN2013fs (Yaron et al. 2017) and the type IIb SN2013cu (Gal-Yam et al. 2014) suggest circumstellar material confined to within < 10¹⁵ cm of the progenitor; directly comparable to the sizes of the extended emission associated with the cluster WRs. While SN2013fs has been attributed to the explosion of a RSG, the wind signatures identified for SN2013cu suggest instead a WN(h) progenitor (Gal-Yam et al. 2014). Further observations to quantify the amount of material contained within these structures to compare to the estimate obtained for SN2013cu would be of obvious interest.

Finally we note that multiple lines of evidence suggest that the pre-SN circumstellar environment may be aspherical, with both clumpy (Smith et al. 2009b) or disc-like configurations proposed (Levesque et al. 2014) for a subset of SNe, of which SN1987A is the most obvious exemplar. Clearly the circumstellar nebula surrounding e.g. Wd1-26 satisfies the former condition but we also highlight the apparent toroidal geometry of the supergiant B[e] star Wd1-9, which in this case is thought to arise from binary interaction. Nevertheless, taken as a whole the properties of the circumstellar environments of the massive stars within Wd1 suggest that they, and potential SNe arising from them, may not be fully understood in isolation.

9. Concluding remarks

In this work and that of Fenech et al. (2017) we have presented the results of a Band 3 (3 mm) continuum survey of the Galactic YMC Wd1. This target was chosen due to a combination of proximity, co-evality, age and mass; the combination of the latter three parameters resulting in a rich population of post-MS stars characterised by the co-existence of both cool super-/hypergiants and Wolf Rayets. An unprecedented total of 50 stellar sources were detected, comprising examples of every phase of massive post-MS evolution present within the cluster. The full exploitation of this exceptional dataset in order to constrain the physics of (radiatively-driven) stellar winds – e.g. the presence of “clumping” – requires quantitative analysis (e.g. non-LTE model-atmosphere analysis and/or hydrodynamical simulations of nebular emission) which is beyond the scope of this work. Nevertheless we may immediately draw two broad conclusions from the provisional analysis included here.

Firstly we were able to utilise canonical stellar- and wind-properties and the formulation of Wright & Barlow (1975) to infer mass-loss rates for both the hot (WR, LBV, sgB[e] and OB Ia/Ia⁺) and cool (YHG and RSG) cohorts, with additional estimates for the latter objects being drawn from their bulk nebular properties (Do10). As a consequence we were able to follow the evolution of the mass-loss rate through all post-MS phases to the brink of SN; an essential observational constraint for quantifying the physics of massive stars. As expected we found that the mass-loss rates of the (brighter) OB super-/hypergiants are approximately an order of magnitude lower than those of the WRs, with both appearing broadly consistent with determinations from field populations. Where available, estimates for the sole LBV and the YHGs/RSGs within the cluster were comparable to the WR population, but not as extreme as measured for other such examples (e.g. Ṁ > 10⁻⁴ M_⊙ yr⁻¹ for such stars during outburst). We highlight that the distribution of spectral types of the OB supergiants is particularly fortuitous since it spans the predicted location of the bi-stability jump induced by the recombination of Fe IV to Fe III; Wd1 therefore offers a powerful test of predictions for wind properties either side of this discontinuity.

In general we find no evidence that massive binaries are systematically brighter than single stars due to an additional contribution from the wind-collision zone. However the mass-loss rate inferred for the sgB[e] star Wd1-9 is exceptionally high (Fenech et al. 2017), being directly comparable to the most extreme WRs present (WR F (WC9) and WR A (WN7b)). Multiple lines of evidence point to Wd1-9 being a massive interacting binary currently exhibiting rapid (case A) mass transfer/loss. Our observations suggest this has led to the formation of a massive circumstellar torus which drives a bipolar wind. Since Wd1-9 is the only star within Wd1 to exhibit this phenomenon it suggests that the duration of this violent phase is comparatively short, although if the current mass-loss rate were maintained for only the ∼5–10 × 10⁴ yr predicted by theory (Petrovic et al. 2005) its influence on the subsequent evolution of both components would be profound.

Secondly, we highlight the degree to which the environment in which the stars are located affects their circumstellar properties. This is most evident in the manner in which the ejecta associated with both the RSGs and YHGs appears to be sculpted during interaction with the wider cluster, although the mechanism(s) involved – cluster wind and/or radiation field – remain uncertain. While this finding was anticipated the resolution of the sources associated with a number of the OB supergiants and WRs was not. Intriguingly, while such a phenomenon has not been associated with isolated stars, Lang et al. (2005) report a similar extension to the radio sources associated with a subset of the O supergiants and WRs within the Quintuplet. Following Sect. 8.2 the least objectionable hypothesis is that these represent compact wind blown bubbles confined via one of the preceding physical effects.

If correct this would have important ramifications both for the initiation and driving of the cluster wind via the combined action of individual stellar winds and SNe and also for the nature of the SNe themselves. We are long accustomed to the idea that the binary nature of a SN progenitor will affect the nature of the explosive endpoint via stripping of the outer H-rich stellar layers and/or the subsequent shaping of the pre-explosion circumstellar material – with Wd1-9 providing a case-in-point. However if the circumstellar envelope of a SN progenitor may be shaped by its immediate environment, potentially leading to the formation of a dense envelope in the absence of an episode(s) of enhanced impulsive mass-loss, then one must accept that the diversity of SNe morphologies does not solely result from the properties of the progenitor and instead originates in part from the wider context of explosion. Given the expectation that the majority of massive stars are born in a clustered environment, this is potentially a far reaching conclusion.

To conclude, ALMA appears a uniquely powerful tool for the quantitative analysis of mass-loss rates for massive stars, as well as their interaction with their wider environment. As such it has the potential to revolutionise our understanding of such objects, providing a window into short-lived, violent phases of their evolution and opening a window onto quantitatively new phenomena.

---

Acknowledgments

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00897.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. D. Fenech wishes to acknowledge funding from a STFC consolidated grant (ST/M001334/1). IN is partially supported by the Spanish Government Ministerio de Economía y Competitivad (MINECO/FEDER) under grant AYA2015-68012-C2-2-P. The authors would also like to thank Dr. Roger Wesson (UCL) for providing the VPHAS data for comparison.

References

Appendix A: Tables

Appendix B: Figures

Appendix C: Supplementary images

All Tables

All Figures