© The Authors 2025

1 Introduction

Mass loss by massive stars in their red supergiant (RSG) phase is poorly understood, yet it is important to characterise this mechanism if we are to understand the late stages of stellar evolution and the behaviour and impact of their progeny, namely, core-collapse supernovae (SNe). Mass-loss rates retrieved from spectral-energy distribution (SED) modelling of RSGs in clusters cover a range 10⁻⁷-10⁻⁵ M yr⁻¹ (e.g. Beasor et al. 2020). In the case of some nearby extreme Galactic RSGs, mass-loss rates can reach up to several 10⁻⁴ M yr⁻¹ (e.g. Andrews et al. 2022; De Beck et al. 2010; Shenoy et al. 2016). Observations of VY CMa, the most-studied extremely dusty Galactic RSG, further suggest a complex mass-loss history, with a combination of steady and episodic mass loss over time and outflows with different velocity fields (e.g. Ziurys et al. 2007; Quintana-Lacaci et al. 2023).

Whereas the driving mechanism behind mass loss from lower-mass evolved stars in the asymptotic giant branch (AGB) phase is broadly understood as a combination of stellar pulsations, convection, shocks high up in the stellar atmosphere, and radiation pressure on dust particles formed close to the star (e.g. Höfner & Olofsson 2018), the underlying mechanism(s) for RSG mass loss is still unknown. Suggested ingredients include, among others, turbulent convective motions and radiation driving on molecules (Josselin & Plez 2007), giant convection plumes lifting material beyond the stellar gravitational potential (e.g. López Ariste et al. 2023), and magnetic fields (e.g. Tessore et al. 2017; Vlemmings et al. 2017).

RSG mass loss forms the circumstellar environment (CSE) that Type II SN expand into. Studies have reported a requirement for dense CSE to be present close to the progenitor star to explain the SN light curve properties (e.g. Khazov et al. 2016; Morozova et al. 2018). The dense CSE forms the basis for the build-up of a cold, dense shell where dust can efficiently form, an effect strengthened by clumpiness in the progenitor CSE (e.g. Inserra et al. 2011). Furthermore, it has been suggested that geometrical asymmetries found in the time-resolved observations of SN2023ixf, a Type II SN in M101, could be linked to strong asymmetries in the pre-explosion CSE (e.g. Singh et al. 2024; Shrestha et al. 2025). The progenitor to SN2023ixf has been described as a star with properties very close to those of VY CMa (Jencson et al. 2023), strongly motivating further highangular resolution observations of Galactic RSGs with known high mass-loss rates to help understand the CSE structure of progenitor stars to Type II SN.

NML Cyg is an extreme Galactic RSG with similar stellar properties to VY CMa. Molecular maser and thermal line emission profiles and dust emission observed towards NML Cyg are suggestive of a similarly complex circumstellar appearance and mass-loss behaviour as found for VY CMa (e.g. Etoka & Diamond 2004; Monnier et al. 2004; Schuster et al. 2006, 2009; Andrews et al. 2022; Singh et al. 2021, 2022). Characterising the CSE of this member of the class of rare, extreme Galactic RSGs will help us set constraints on the mass-loss mechanisms active in RSGs and set valuable constraints on the progenitors of Type II SNe.

In this paper, we describe observations of the Galactic RSG NML Cyg at λ = 1.3 mm (Sect. 2) and report on the properties of its circumstellar environment. We present the continuum emission and the derived dust properties in Sect. 3 and selected molecular-line emission and the derived outflow properties in Sect. 4. We discuss our findings in light of observations of other well-studied RSGs and constraints on SN progenitors and pre-SN mass loss in Sect. 5. We present our conclusions in Sect. 6.

2 Observations

Observations were gathered on 25 February and 4 and 9 March 2021 using the Band 3 receiver on the NOEMA Interferometer. The observations were taken in the A configuration with 11 antennas with baselines between 19.3-751.7 m. The observations gathered on 25 February were discarded due to the unstable atmospheric conditions. On 4 March, some data were flagged due to tracking issues in one antenna, leading to 2.6 h of on-source integration. On 9 March, 1.1 h of on-source integration was achieved. Calibration was carried out with the triple stellar system MWC349 as the flux calibrator and the quasar 2010+723 as the phase calibrator. The typical absolute flux calibration uncertainty was determined to be 10%.

The observations were obtained with a local oscillator frequency of 226.013 GHz resulting in measurements of the lower sideband (LSB) at 214.4-222.6 GHz and the upper sideband (USB) at 229.6-237.8 GHz. Both LSB and USB were observed at a spectral resolution of 2 MHz, corresponding to a velocity resolution of ∼3 km s⁻¹.

The data were reduced with the use of GILDAS (Bardeau & Pety 2020), with reduction and calibration applied via GILDAS/CLIC. The imaging and deconvolution were carried out with GILDAS/MAPPING, and additional analysis was performed with CASA Version 6.1.1.15 (McMullin et al. 2007). After the initial reduction and calibration, line identification was carried out across the spectra extracted for apertures of 8″, 1″, and 0.4″ around the central pixel, using the molecular databases provided by CDMS (Müller et al. 2001, 2005) and JPL (Pickett et al. 1998). Spectra in LSB and USB extracted towards the star and for a 1″ aperture are shown in Fig. A.1.

When imaging, HOGBOM cleaning was applied with robust weighting. Continuum images and spectral line cubes were generated with a pixel size of 0″.08. A noise limit was set to 20% of the expected sensitivity of the primary beam. Clean components were also selected with the use of the SUPPORT procedure within GILDAS/MAPPING to optimise the sensitivity of the final continuum images and spectral line cubes.

Final continuum images for each sideband were constructed by masking for all lines and artefacts across the total bandwidth observed. Self-calibration solutions were applied to improve the final fidelity of the results, increasing the signal-to-noise ratio (S/N) of the continuum images by more than a factor 3. Our final continuum observations (Fig. 1) achieved sensitivities of 0.4-0.5 mJy beam⁻¹. The synthesised beam size, noise level, and maximum recoverable scale (MRS) in the two observing bands are given in Table 1.

Spectral line cubes were extracted for each spectral line of interest after continuum subtraction was applied to each sideband. To provide additional information on large-scale extended structure for the CO, ¹³CO, SO₂, and SO lines, short-spacing observations from the IRAM 30m telescope using an on-the-fly map were combined with the NOEMA observations using the GILDAS/UVSHORT command within GILDAS/MAPPING. Selfcalibration solutions were transferred from the continuum emission of the relevant sideband in order to improve the fidelity of the results, doubling the S/N of the weakest line detections. The typical rms of the spectra extracted from the data with an aperture of 0″.4 is 2 mJy at 3 km s⁻¹ resolution and it is fairly homogeneous across the entire spectral range.

Fig. 1 Continuum emission from NML Cyg. (Top) Emission in LSB (colour map) and USB (white contours at 3, 5, 10, 30, 50, 100, 200 σ), with two continuum components labelled as A and B and the stellar position as a black star. The synthesised beam is shown in the bottom left corner (LSB: yellow, USB: white). Residuals in USB after subtraction of the IMFIT results are shown in green contours (at -10, -5, -3, 3, 5, 10 σ) with dashed contours for negative signal. (Bottom) Estimating the stellar position. Peak flux positions of selected molecular transitions with errors according to the 50 mas astrometric accuracy (filled circles) at systemic velocity overlaid on USB continuum contours (grey) from the top panel. Also indicated: continuum peak (purple X), stellar position extrapolated from Zhang et al. (2012b) (black star), peak position of HST 0.55 μm emission (red star), and LSB and USB beams (cyan and black ellipses).

3 Continuum emission

We detected two major continuum components, which we refer to as A and B in the remainder of this paper, peaking at ≈475σ and ≈145σ, connected by a band of low-brightness emission at up to 10-20σ (Fig. 1), consistent between the two sidebands. Both continuum components are unresolved in the current observations and component A contains the stellar contribution. We found no evidence of any other continuum emission across the 1′.25 × 1′.25 field of view at the current sensitivity. We discuss the possible configuration of the system and estimate basic properties for components A and B below.

3.1 System configuration

3.1.1 Stellar position

Using the previously determined stellar position and proper motion of the star reported by Zhang et al. (2012b), we calculated the stellar position at the epoch of the NOEMA observations to be α = 20^h46^m25.538^s ± 0.005^s, δ = 40°06′59″.320 ± 0″.005. This matches the emission peaks of a selection of high-S/N and rather compact molecular emission lines from SiO, ²⁹SiO, SiS, ²⁹SiS, and H²O well within the synthesised beam size of the observations (e.g., Figs. 1b and A.9) and considering the astrometric accuracy of the NOEMA observations of ≈50 mas. We expect these emission regions to be centred on the star as, for example, predicted by Menten & Melnick (1989) for the H²O maser at 232 GHz and confirmed by Quintana-Lacaci et al. (2023) for VY CMa. We conclude that the expected stellar position indeed likely coincides with the actual stellar position.

We find a small offset between the peak position of the HST 550 μm emission from Schuster et al. (2006, corrected for the proper motion), and the expected stellar position. However, we note that the astrometric accuracy reported for the HST observations is of the order 0″.3 (The DrizzlePac Handbook, Hoffmann et al. 2021), which does not exclude a match between this peak and the stellar position we derived. Furthermore, since those observations likely show dust-scattered stellar light, the emission peak does not necessarily coincide with the stellar position, as also discussed for the case of VY CMa by Zhang et al. (2012a).

3.1.2 Component A

We used the IMFIT procedure within CASA to retrieve the beam-deconvolved properties of the two continuum components approximated as 2D Gaussians (see Table 2). The main component, component A, is highly elongated (θ_min/θ_maj ≈ 0.35) along a roughly NE-SW axis (position angle PA≈35°) and encompasses the stellar position, but it peaks at a position offset from it.

Assuming the stellar parameters from Zhang et al. (2012b) and assuming a spectral index of 2, the star has an expected diameter of roughly 40 mas at the observed frequencies, that is, 10-15% of the synthesised beam size, and an expected stellar flux of only ≈7% of the total flux of component A (S_⋆,LSB = 21.1 mJy, S_⋆,USB = 24.2 mJy). Considering this, we do not subtract the stellar contribution to recover potential substructure in component A.

Using the stellar position described above, we find that the peak of component A is offset by ≈275 mas from the star (≈440 AU at 1.6 kpc), ata PA≈ −135°, that is, to the SW. This offset is not resolvable by the angular resolution of the observations, but is significantly larger than both the uncertainty on the stellar position (Zhang et al. 2012b) and the astrometric accuracy of the NOEMA observations (≈50 mas). We therefore suggest the presence of a bright dust clump physically offset from the stellar position in a SW direction.

Considering the similarity in size and orientation between our component A and the eastern component measured in the mid-infrared (MIR) by Schuster et al. (2009), we propose that these measurements might be tracing the same material. This is supported by the close match between our estimated dust properties (Sect. 3.2) and the results of Schuster et al. (2009). Given the considerations about the stellar position discussed above, our observations of continuum and thermal spectral line emission at 230 GHz thus imply that the assumption by Schuster et al. (2009), namely, of the star being located at their detected MIR peak brightness, is likely incorrect. This was also suggested by Zhang et al. (2012b) based on continuum and SiO and H₂O maser measurements at 43 GHz and 22 GHz. We discuss the match between our observations and the maser observations in Sect. 4.

We can use the molecular emission observed along the line of sight towards component A to constrain this component's location with respect to the plane of the sky. Pencil-beam spectra of selected transitions of SO₂, SO, H₂S, and PN, extracted towards the peak of component A reveal strong emission lines with a lowered emission and even absorption at blue-shifted velocities roughly in the range [−40, −16] km s⁻¹ (Fig. 2), and very clearly visible in the channel maps of CO, SiO, SO, and SO₂ at velocities ≤−22 km s⁻¹ in Figs. A.3, A.5, A.6, and A.7. We note here that not all of the detected SO₂ lines show absorption in this pencil beam towards the peak of A and we discuss this further in Sect. 4.2.1. The absorption features can be explained by the presence ofa warm body of dust behind colder gas. Although the velocity ranges covered by the absorption components are not identical across all molecular transitions, we propose that component A is located in front of the plane of the sky. The presence of emission at other (less blue-shifted) velocities along this line of sight suggests that the dust could be, at least partially, optically thin at the observed wavelengths. However, since the deconvolved size of component A is smaller than the synthesised beam, we would also detect emission originating to the sides of the component in the same beam, even if the dust were optically thick. We note that the position-velocity diagrams for a cut along PA = 65^° (discussed in Sect. 4) show a clear indentation for the SiO v = 0 and H₂S contours along the negative velocity axis at an offset of -0″.3. This corresponds to the location of component A and suggests that the dust could indeed be optically thick at these wavelengths, effectively blocking the emission coming from a very small area on the sky. Since we do not see a similar behaviour for the SiO v = 1 emission, which is likely more compact than the SiO v = 0 emission, we suggest that component A is located beyond the SiO v = 1 emission volume, relative to the star. Observations at a higher angular resolution and at significantly different wavelengths are needed to obtain a spectral index for component A to set stronger constraints on the dust properties.

Fig. 2 Pencil-beam spectra towards NML Cyg, extracted towards the stellar position (left) and towards the peak positions of components A (middle) and B (right). The spectra have been normalised to their peak and a vertical offset has been added. The transitions are labelled in each panel. Note: the spectrum of PN 5 - 4 towards component B has been scaled for the sake of visualisation. Shaded regions indicate velocity ranges of interest with respect to absorption and emission features, with the relevant boundaries indicated in the top of each panel.

3.1.3 Component B

The secondary continuum component, component B, is roughly 40-45% as bright as and somewhat larger than component A. The IMFIT results suggest a slight elongation along a roughly NNE-SSW direction (PA≈13°) and an offset with respect to the stellar position of ≈985 mas (≈1570 AU) at PA≈−65° (i.e., to the NW of the star).

For comparison, Schuster et al. (2009) also found an elongated structure in the 8.8 μm MIR observations, with emission from warm dust (200-500 K) reaching up to 300-500 mas to the NW of the star, placing this warm dust between the star and our component B, where our data reveal a band of continuum emission connecting components A and B.

We found emission from several molecular species roughly coincident with the projected location of the peak emission of component B (CO, SO, and H₂S; see Sect. 4), with a clear prevalence of red-shifted emission, mostly in the range [17, 23] kms⁻¹. The morphology of the emission from H₂S and H${}_2^{33}$S traces the location of component A, the ‘bridge’ from A to B, and, finally, the location of component B when moving from blue-shifted to red-shifted velocities. A similar co-location of H₂S and continuum emission has been reported in observations of VY CMa (Quintana-Lacaci et al. 2023). Hence, we propose a location for component B and possibly also the dust excess measured by Schuster et al. (2009), behind the plane of the sky. As opposed to our results for component A, we found no absorption in the spectra extracted towards component B.

3.1.4 Geometry

Multiple mass ejections and/or dust shells have been invoked to explain a number of previous observations of NML Cyg (e.g. Monnier et al. 1997; Blöcker et al. 2001). The presence of multiple continuum components was also suggested by Schuster et al. (2006, 2009) based on the irregular morphology found in the NIR/MIR observations and is now solidified by our mm-wavelength observations. We note that component B reaches far beyond the 0″.1 dust destruction radius suggested by Schuster et al. (2006, 2009).

After subtraction of the two 2D Gaussian components, the IMFIT procedure reveals residual emission (Fig. 1a) close to the peak positions of components A and B, at the NE tip of component B, and in the region connecting components A and B. Several negative residual components are also visible. Clearly, the approximation of the observations as a combination of two 2D Gaussian components is insufficient. Therefore, we consider it likely that higher-angular resolution observations will reveal substructure within the here presented continuum emission, similar to what was found for VY CMa (O'Gorman et al. 2015; Kaminski 2019; Asaki et al. 2020).

Components A and B are located along a NW-SE axis and are, at the same time, both elongated along a roughly perpendicular direction (NE-SW). We note that OH, H₂O, and SiO maser observations have revealed a similar NW-SE geometrical axis, leading to the suggestions of a possible bipolar outflow from NML Cyg, with a receding NW lobe and an approaching SE lobe, and a rotating SiO maser shell (Richards et al. 1996; Etoka & Diamond 2004; Boboltz & Marvel 2000). This orientation is in line with our suggestion that component B is located behind the plane of the sky. We refer to Sect. 4 for a further discussion on the envelope morphology based on our NOEMA observations and the high-angular resolution maser results from the literature.

3.2 Dust temperature and mass

As argued above, dust is the dominant contributor to the observed continuum emission. To estimate the dust mass responsible for the observed flux densities, S_V, we assume the dust is optically thin and that the Rayleigh-Jeans approximation is valid to obtain the dust mass, M_d, from $$S_{\nu }=\frac{3M_{\mathrm{d}}{Q}_{\nu }{T}_{\mathrm{d}}k{\nu }^2}{2a_{\mathrm{g}}{\rho }_{\mathrm{g}}{c}^2{d}^2},$$(1)

where Q_V is the grain emissivity, T_d is the dust temperature, k is the Boltzmann constant, a_g and ρ_g are the radius and mass density of the dust grains, respectively, c is the speed of light, and d is the distance to the star (Knapp et al. 1993). In the case of component A, S_ν is the measured flux density of the component corrected for the expected stellar flux, S_⋆,v (see Sect. 3.1.2).

We assume that the grain emissivity behaves as Q_ν ∝ V^β and S_V ∝ ν^α with α = 2 + β. For the dust around NML Cyg, we expect β-values in line with those found for other evolved stars. Previous models of the dust around VY CMa have used typical values of β = 0.7-0.9 (Knapp et al. 1993; O'Gorman et al. 2015; Kaminski 2019). Our continuum detections and the non-detection of dust at 43 GHz by Zhang et al. (2012b) are indeed in line with these literature values. We proceed to derive dust properties under the assumption that β = 0.9.

Optically thin dust at a radius r_d, heated by the incident stellar radiation field yields a dust temperature $$T_{\mathrm{d}}={\left(\frac{L_{\star }{T}_{\star }^{\beta }}{16\pi \sigma r_{\mathrm{d}}^2}\right)}^{1/(4+\beta )}.$$(2)

Assuming a stellar luminosity of $L_{\star }\approx 2.7\times {10}^5{L}_{\odot }$, stellar temperature T_⋆ = 2500 K (Zhang et al. 2012b), and distances, r_d, corresponding to the offsets Δ listed in Table 2, we derived T_d ≈ 445 K for component A and T_d ≈ 265 K for component B. Assuming typical dust grain properties for M-type evolved stars (a_g = 0.2 μm, ρ_g = 3.5 g cm⁻³; Knapp et al. 1993), we derived dust masses M_d ≈ 11 × 10⁻⁴ M_⊙ for component A and M_d = 8-9 × 10⁻⁴ M_⊙ for component B.

The assumed distances r_d are lower limits, meaning that the calculated dust temperatures, T_d, are upper limits and the calculated masses, M_d, are lower limits (under the assumptions on optical depth and dust emissivity). An angle to the plane of the sky of 45^° or 80^° would lead to larger spatial offsets from the star and a consequent ≈15% or ≈50% decrease in dust temperature, respectively. The derived masses, in turn, would increase by ≈10% and a factor of ≈2, respectively. Unfortunately, we cannot, currently, constrain these inclination angles; we refrained from deriving an inclination angle based on an assumed expansion velocity considering the outflow complexity (see Sect. 4) and that we do not have information on the clump geometry.

4 The molecular outflow

We detected 72 spectral features of which we identified 67 to pertain to 26 different isotopologues. Here, KCl and its isotopologue K³⁷Cl have been identified for the first time around an M-type (carbon-poor) supergiant. An overview of the spectra is given in Fig. A.1.

In this paper, we use the first spatially resolved observations of thermal emission from CO, SiO, SO₂, SO, and H₂S to trace the extent, morphology, and kinematics of NML Cyg's outflow. We investigate their relation to the NOEMA continuum measurements presented in Sect. 3.2 and literature results, in particular, from maser emission. An in-depth presentation of the additional molecular features in this data set will be made in upcoming publications.

4.1 Extent

We detect emission from ¹²CO and ¹³CO and for the first time also from C¹⁸O towards NML Cyg. Figures 3, 4, A.3, and A.4 show channel maps and position-velocity diagrams of the ¹²CO and ¹³CO(2 − 1) emission that we base the rest of the discussion on.

The CSE shows significant asymmetry, with the largest CO extent measured to the SE, the smallest to the NW. We detect ¹²CO(2 - 1) emission at >10σ out to 8″.6 from the stellar position in the SW direction, 6″.8 in the NE direction, 6″.1 in the NW direction, and 11″.4 in the SE direction (Fig. 4). This largest extent in the SE direction is also traced by a 3σ component in the ¹³CO emission. We note that the PV-diagrams in Fig. 4 show ¹²CO emission >3σ beyond 10″ along all four of these directions (SW, SE, NE, NW), but we do not consider these here, as there are artefacts in the images at similar levels (closer to the edge of the FoV) and because we cannot easily disentangle the CSE emission from the ISM contribution in the channels at 5-8 kms⁻¹ (Fig. A.2), where the largest CO extent is seen in Fig. 4.

The detected CO angular extents of 11″.4 and 6″.1 translate to linear, radial (projected) distances from the star of 2.7 × 10¹⁷ cm and 1.5 × 10¹⁷ cm, respectively (at 1.6 kpc; Zhang et al. 2012b). We expect these scales to be somewhat smaller than the CO photodissociation radius, by up to a factor of a few, as a consequence of the excitation and detectability of CO(2 - 1), as also described by Ramstedt et al. (2020) for the densest AGB outflows. Assuming the CO abundance typically assumed for an M-type AGB star (CO/H₂ = 2 × 10⁻⁴) and extrapolating the results of Saberi et al. (2019, see their Table B.1) to a massloss rate of 5 × 10⁻⁴ M_⊙yr⁻¹ (which would be an upper limit for NML Cyg; Gordon et al. 2018) leads to a half-abundance radius r_1/2 ≈ 10¹⁸ cm. This would be at the high end of the range we would find likely to be the actual photodissociation radius for NML Cyg's outflow. Considering the proximity to Cyg OB2, NML Cyg plausibly resides in a more intense radiation field, which would cause a smaller molecular outflow, but likely by less than 15% because of the high mass-loss rate (Saberi et al. 2019). Furthermore, the smaller extent to the NW could possibly be due to irradiation from Cyg OB2, although the observed SE-NW axis (PA ≈ -35°) is not quite aligned with the direction to the cluster (PA ≈ -63°). We conclude that the CO extent around NML Cyg fits the expectations rather well for an assumed mass-loss rate and outflow velocity and that the molecular CSE extends significantly beyond the limits suggested by Schuster et al. (2006). We also note that according to Dharmawardena et al. (2018), the transition from the stellar to the interstellar radiation field dominating grain heating occurs at ≈0.19 pc ≈ 24″ in the case of NML Cyg, that is, well beyond where we detect CO(2 - 1) emission.

We note that asymmetries, with a predominant SE-NW orientation, are also seen in the red-shifted SO and SO₂ emission. These show a clear extension to the SE, most clearly visible in the 14 km s⁻¹ channels (Figs. A.6 and A.7). An extension along a similar SE-NW axis - but on different spatial scales - is traced by SiO, OH, and H₂O maser emission (PA ≈ 140° - 150°; Boboltz & Marvel 2000; Richards et al. 1996; Etoka & Diamond 2004) and dust emission (Blöcker et al. 2001). The observed asymmetry could be a consequence of extrinsic factors, for example, a less intense UV irradiation of the molecular gas in the SE compared to the NW direction, connected to the location of NML Cyg relative to the Cyg OB2 cluster, and/or it could be a consequence of factors intrinsic to the NML Cyg system, such as a highly irregular mass loss.

Although we refrain from discussing in detail the chemical networks in the CSE of NML Cyg in this paper, we wish to highlight some findings from a first examination of the extent of the molecules other than CO presented in this paper. We find that the H₂S emitting region is strongly correlated with the continuum emission presented in Sect. 3 and only extends somewhat beyond that, out to a distance of ≈1″.4 from the star, which is significantly smaller than the 1000R_* ≈ 9″ suggested by Singh et al. (2022). Their modelling efforts, based on single-dish observations, led to the suggestion that CO, HCN, and H₂S have the most extended molecular abundance distributions. Furthermore, while the abundance profiles modelled by Singh et al. (2021, 2022) suggest an SiO envelope that is similar to or smaller than those observed for SO and SO₂, we find that the emitting region of SiO (v = 0) reaches angular distances from the star of 2″.8-4″.5 and is slightly larger and more circularly symmetric than the overall SO and SO₂ emitting regions, which exhibit more substructure along the SE-NW axis and extend out to ∼3″.This direction-dependent difference in extent is best seen when comparing Figs. 4a and 4b.

Fig. 3 Blob of CO emission to the SW of the star. The panels show of the 12CO J = 2 - 1 emission (colour) and 13CO J = 2 - 1 emission (white contours at [3, 5, 10]σ) across the velocity range [17; 26] kms−1. The extent of the colour scale has been limited for the sake of visibility. The stellar position is indicated as a black cross (+), the NOEMA beam is plotted in the bottom left corner as a white filled ellipse.

Fig. 4 Position-velocity diagrams. Cuts are centred on the assumed stellar position (marked as X) and at position angles of 65° to cut through the SW clump (top) and perpendicular to that, at 155° (bottom). Contour levels are drawn at [3,5,10,15, 20,40,80,100] × σ with σ = 2.0 mJy/beam; for 12CO, we only show the 3σ contour (red). Note: the CO emission at 5kms−1 and 8 km s−1 is contaminated by ISM contributions.

4.2 Morphology

The CO emission traces a complex outflow morphology, as shown in the channel maps and position-velocity diagrams in Figs. A.2–A.4, and 4. We find significant substructure, visible as arcs and blobs, also in the emission of SiO, SO, SO₂, and H₂S (Figs. A.5–A.9).

4.2.1 Molecular emission co-spatial with the continuum emission

We find a lack of emission and even find absorption features in CO, SiO, SO, SO₂, H₂S, and PN at blue-shifted velocities ≲16 km s⁻¹ in the line of sight towards the peak of component A and rather see a ring of emission surrounding this spatial region (see spectra in Fig. 2 and channel maps in Appendix A). Considering that we can see a similar effect for a variety of molecules, but not all detected lines, and that it is limited to the blue-shifted part of the spectrum, we suggest that this is the result of absorption in the (colder) gas against a background of hot dust (see also Sect. 3.1.2), rather than optically thick dust obscuring gas located behind it. Another possibility is that the affected lines are optically thick in these velocity ranges, obscuring the underlying dust emission, and that continuum subtraction (performed across all channels) causes an erroneous absorption feature, as suggested by Humphreys et al. (2024) in the case of some molecular lines observed towards VY CMa. However, since we do not attempt to constrain detailed properties of the dust or gas based on these features and since component A is smaller than our synthesised beam and possibly made up of multiple smaller components, we did not explore this option further.

We observe that the detected SO₂ emission shows a diverse range of behaviour: some transitions show the absorption feature towards the peak of component A mentioned above, whereas others show distinct emission in this same line of sight and velocity range (Fig. 2). This suggests that the excitation of SO₂ could be very strongly coupled to the radiation field in line with the findings of Danilovich et al. (2016), emphasising the importance of considering the possible complexity of excitation when interpreting molecular emission maps of only selected transitions. Cernicharo et al. (2011) predicted and confirmed absorption and maser behaviour for a number of long-wavelength rotational transitions of SO₂ towards cold (<40 K) molecular clouds.

In contrast to what we find for component A, the spectra towards the peak of component B do not exhibit strong absorption features. We find that the SiO emission shows a depression in a red-shifted velocity range around 20km s⁻¹ where several other molecular transitions, including ¹³CO, H₂S, SO, and SO₂ conversely show rather distinct emission components. This could be suggestive of a local depletion of SiO.

We find that the ¹²CO, ¹³CO, and H₂S emission trace a clear ‘bridge' between components A and B across the velocity range [-7, 23] kms⁻¹, coincident with a ‘bridge' of significant continuum emission (Figs. 1, A.3, A.4, A.8). The velocities at which we see this emission structure suggest that the connection from the stellar position to B extends behind the plane of the sky. At the same time, the SiO, SO₂, and SO emission do not trace this ‘bridge'. The SiO emission is mostly limited to the NE and SE quadrants of the maps, whereas the SO and SO₂ emission regions do reach significantly beyond the location of component B in the NW quadrant across the velocity range of [-14, 13] kms⁻¹. The SO emission maps, however, show a distinct gap at the location of the ‘bridge' at redshifted velocities (most clearly seen for the 5₆ - 4₅ transition at [8, 17] kms⁻¹; Fig. A.6).

The morphology of the detected H₂S emission is closely related to that of the continuum, with emission components clearly spatially correlated with components A and B across velocity space and extending only little beyond these (Fig. A.8). This is in line with the finding for VY CMa that the presence of H₂S appears strongly connected to the presence of warm dust (Quintana-Lacaci et al. 2023). Higher-angular resolution observations are needed to resolve the structure of both the gaseous and dusty component to establish how close to the star H₂S actually is formed and how this relates to the presence of dust. In the case of AGB stars, H₂S is thought to be formed at high densities in the inner layers of the outflows of oxygen-rich stars with high mass-loss rates (Danilovich et al. 2017). According to the models by Van de Sande & Millar (2022), the presence of H₂S and a shell-like distribution of SO, possibly quite similar to what we detected for NML Cyg, rule out a very strong UV field from a companion star.

4.2.2 SW blob

We find an emission blob in the ¹²CO and ¹³CO emission to the SW of the star extending from roughly 5″ to 6″.7 from the star, at PA = 245^° with respect to the stellar position, covering an opening angle of roughly 30°. It is most apparent at [17, 26] kms⁻¹, but likely already present at less red-shifted velocities (Figs. 3 and A.2). We note that the 29 km s⁻¹ channel bears no sign of this clump at all beyond two 3σ spots of only ≈0″.3 across.

The SW blob is not measurably traced by any of the other emission lines in our data set and we find no spatially coincident continuum emission at our current sensitivity. We note that the PA = 245^° is similar to that of component A (PA = 225°), but we recall that the blue-shifted absorption of spectral line emission led us to derive a location in front of the plane of the sky for component A, whereas the SW blob in CO is found at significantly red-shifted velocities, suggesting a location behind the plane of the sky.

Because we detected the blob in only one transition of ¹²CO and ¹³CO (J = 2 - 1) and the emission is highly complex, we cannot reliably estimate the temperature or (CO) mass in the SW blob. We are not aware of any previously reported observations tracing this blob.

4.2.3 Arcs and shells

We used position-velocity (PV) diagrams to highlight some of the substructures and (a)symmetry in the outflow. Figure 4 shows PV-diagrams for cuts through the stellar position along PA = 65^°, cutting also through the CO clump to the SW of the star, as well as perpendicular to that, at PA = 155^°, that is, roughly along the SE-NW axis identified in earlier maser and dust observations (Blöcker et al. 2001; Boboltz & Marvel 2000; Etoka & Diamond 2004; Richards et al. 1996). The cut along PA = 65^° and through the stellar position results in the best alignment in the PV-diagram of the CO emission with a central ellipse representing a spherical shell of ejected mass out to a radius of 1″ and at a velocity of ≈35 km s⁻¹. Cuts along other position angles (through the stellar position) lead to worse aligned CO emission in the PV-diagrams, mostly noticeable as spatial offsets at the extreme negative velocities. This is suggestive of some deviation from spherical symmetry in the CSE already within about 1″ from the star, but our current observations do not provide the necessary angular resolution to characterise this structure in detail.

Figure 4a clearly shows the SW blob, extending (to the top left) at red-shifted velocities and negative spatial offsets. In addition to this, we see multiple arcs in the PV-diagrams of the CO emission, suggestive of (partial-) shell-like substructure in the CSE and reminiscent of structures seen in the outflows of several AGB stars (e.g. Guélin et al. 2018; Decin et al. 2020; Randall et al. 2020). From a manual exploration of the PV-diagrams, we find that the contrast between the arc and inter-arc emission ranges from a factor of a few to almost ten. However, the complexities in the PV-diagrams are too many to get accurate contrast values. We propose the presence of multiple, likely incomplete shells of material having been ejected with possibly different velocities during episodes of enhanced mass loss. The tell-tale arcs are seen in the PV-diagram of CO for the cut at PA = 65^°, where they appear reaching out to offsets of roughly -4″.5, -2″.5, -1″.2, 1″.0, 2″.5, 4″.0, and 7″.5. In Fig. 5, we show what spherical shells of a given radius and moving at a given expansion velocity would look like in the PV-diagram. The largest of the arcs is seen only for positive offsets, that is, corresponding to a NE offset from the stellar position.

The CO emission is too complex and the angular resolution too low to determine with certainty whether the shells are concentric. Similarly, it is difficult to constrain the velocities corresponding to the respective arcs with high precision, but a by-eye fitting in the PV-diagrams indicates that the largest shells are noticeably slower than the smaller ones, with ≈22 km s⁻¹ at 7″.5, ≈25 km s⁻¹ at 4″, ≈30 kms⁻¹ at 2″.5, and ≈35 km s⁻¹ for 1″.5. This significant decline in expansion velocity with increasing shell radius can be seen as a ‘tightening' of the arcs in velocity space (Fig. 5). This gradient in velocity could be a consequence of increasing wind velocities or deceleration of the material over time. The former could be linked to, for example, an increase in the stellar luminosity, the onset of a steadily increasing superwind, a change in the dust properties, or a change in the underlying ejection mechanism(s). On the other hand, it is not clear what would cause a deceleration of the CSE.

An additional complexity arises when we consider the PV-diagrams constructed along different position angles. In the PV-diagram for the cut at PA = 155^° (Figs. 4c and 4d), we found an arc extending out to roughly 10″, but moving at the same velocity as the 7″.5 shell at PA = 65^°. A cut at an intermediate position angle leads to a PV-diagram with an arc at an intermediate shell radius, but again with a similar velocity. This change in the size of the shell with position angle suggests that we are seeing ellipsoidal shells rather than spherical shells, with the major axis aligned roughly at PA = 155^°. This is broadly in line with the previously identified dominant orientation of the infrared and maser emission (Blöcker et al. 2001; Boboltz & Marvel 2000; Etoka & Diamond 2004; Richards et al. 1996). A similar behaviour is also seen for the smaller shells. The uninterrupted structure in the PV-diagrams suggests that it is likely a continuous process that created these structures, but since the shell material seems to be moving at the same velocities (within the uncertainties of the by-eye fitting to the PV-diagram features), the ejection could not have taken place at the same time. It is therefore possible that we are seeing the imprint of a three-dimensional spiral, formed by mass being ejected while NML Cyg moves on a significant orbit or by a binary companion affecting the outflow with an orbital plane along PA = 155°. Additional data at higher angular resolution are necessary to confirm this and obtain more information on a potential orbital configuration.

A first estimate leads us to suggest a temporal spacing between the shells of the order of several hundred years, however, tightening by roughly a factor of two between consecutive detected arcs with estimated average ages of ∼2800 yr, ∼1400 yr, ∼650 yr, and ∼220 yr.

Fig. 5 Same as Fig. 4a, but omitting SO2 contours and adding spherical shells centred on the star with radii 1″.0, 2″.5, 4″.0, 7″.5, and radially outward directed velocities of 35kms−1, 30kms−1, 25 km s−1, and 22 km s−1 , respectively (dashed cyan lines).

4.3 Kinematics

The NOEMA observations and, in particular, the spectra and PV-diagrams as presented above, reveal significantly higher velocities in the wind of NML Cyg than the ≈30 km s⁻¹ previously reported (Andrews et al. 2022; Singh et al. 2021). We detected CO emission across the velocity range [-45; 54] kms⁻¹, with possibly some low-level emission even extending all the way to 65kms⁻¹. The SiO v = 0, J = 5 - 4 emission is seen across the range [-62; 55] kms⁻¹, whereas the SiO v = 1, J = 5 - 4 emission is limited to the range [-37; 46] km s⁻¹. We find SiS v = 0, J = 12 - 11 emission across [-53, 46] km s⁻¹ and SiS v = 1, J = 12 - 11 emission across [-32, 26] kms⁻¹. In all of these cases, the highest-velocity emission originates from within one synthesised beam offsets from the stellar position. The 30 km s⁻¹ ‘threshold', however, is crossed at offsets of ≈ 1″ from the star for SiS v = 0, at ≈ 2″ for SiO v = 0, and for CO even at ≈ 3″.5. This would mean that the large-scale outflow is faster than previously assumed and that there might be a very fast inner-wind component with speeds up to ≈ 60kms⁻¹. The latter would be in line with our earlier suggestion ofa wind-speed gradient across the arcs traced in the PV-diagrams (Sect. 4.2.3) with increasing wind speed for decreasing arc radius and, hence, some intrinsic change in the wind parameters over time, possibly connected to a change in stellar parameters. However, we note that neither the H₂O nor the SiO maser observations were reported to show components at velocities beyond 30 km s⁻¹ (Richards et al. 1996; Zhang et al. 2012b), although both originate from the inner regions of the outflow, within one synthesised beam in the NOEMA observations.

The high-velocity blue-shifted (≤-22kms⁻¹) SiO v = 0 emission is absent towards the SW portion of component A. Also, SiO v = 1 seems to be absent in the directions of the SW part of component A and all of component B, but it is generally not detected at equally large blue-shifted velocities (≤-40kms⁻¹) either. We conclude (see also the discussion in Sect. 4.2.1) that part of the high-velocity SiO is in front of and colder than the dust in component A. We deem the possibility that an optically thick component A obscures a compact highvelocity SiO cloud unlikely because (1) we detect emission in this direction at other velocities (and in other spectral lines) and (2) we would expect emission coming from small projected offsets from the peak of A to still enter the beam and be smeared out across the peak position and the SW portion of component A. However, this puts the high-velocity component at a significant spatial offset from the star, as the projected offset of the peak of component A is Δ ≈ 0″.3 (Table 2), which corresponds to a lower limit to the distance to the star of -4800 AU. Additional observations are needed to spatially resolve this inner region and better constrain both the dust and gas properties.

5 Discussion

5.1 Comparison to maser emission

To get additional insight in the detailed structure of the CSE of NML Cyg, we can compare the high-resolution maser emission reported in the literature to the continuum and line emission in our NOEMA observations. With our assumed stellar position (see Sect. 3.1.1), the 22 GHz H₂O maser component to the NW of the star reported by Richards et al. (1996) is located at the NW edge of component A (peak-PA ≈ -45°; projected offset ≈290 mas; v = 5.5 km s⁻¹). All other 22 GHz H₂O maser components reported by Richards et al. (1996) and Zhang et al. (2012b) are located so close to the stellar position that they coincide with the northern portion of component A and the region where we detect thermal H₂O and SiO emission.

For lack of absolute positioning calibration, Etoka & Diamond (2004) assumed the stellar position to correspond to the brightest maser spot in the blue-shifted OH 1612 MHz maser emission. Considering the overall alignment of the OH 1612 MHz maser emission with the SO₂ emission (both red-shifted and blue-shifted; Fig. 6), this assumption appears to be valid. This positioning places one bright component of the blue-shifted OH 1612 MHz emission mostly concentrated at the NE tip of component A, at the stellar position, and one to the SE of the star. The blue-shifted OH 1612 MHz emission has no particularly bright or compact component at the projected position of component B. The blue-shifted OH 1612 MHz emission matches the blue-shifted SO₂emission well in both extent and shape, with the exception of the northernmost SO₂ emission reaching beyond the OH emission. Similarly, the overall structure of the red-shifted OH 1612 MHz maser emission matches that of the red-shifted SO₂ emission well, with strikingly similar hook-like structures to the NW of component B. However, we find that the SO₂emission extends in an arc-like shape beyond the OH maser emission to the SE of the star. We further also see an apparent avoidance in the red-shifted OH 1612 MHz emission of the peak of component B, which we previously discussed and suggested it is strongly correlated with red-shifted molecular emission.

The blue-shifted OH 1665 MHz emission from Etoka & Diamond (2004) is distributed in compact components corresponding to the stellar position, three positions to the E and NE, and three weak, but significant components to the W, of which two overlap with the on-sky position of component B. The location of all but the stellar components is coincident with a ring-like structure seen in the SO₂emission (Fig. 6f).

Considering the strikingly similar offsets between our continuum components and the two strongest components in the blue-shifted OH 1612 MHz maser maps, we tested the viability of matching those pairs of emission, effectively putting the very brightest maser spot close to the peak of component B, rather than at the stellar position. Based on a worse agreement between the maser maps and several of the substructures identifiable in our CO and SO₂maps, we are fairly certain that the original assumption is correct. This, in turn, means that we find some symmetry around the star between 1) continuum component B to the NW of the star (peak-PA ≈ -65°; projected offset -985 mas), likely behind the plane of the sky (velocities of around 20kms⁻¹) and 2) the OH maser component SE of the star (peak-PA-127°; projected offset -940 mas; Etoka & Diamond 2004) and in front of the plane of the sky (velocities of around -20kms⁻¹). Our observations hence provide further support for the finding of a receding NW lobe and an approaching SE lobe as traced by the H₂O masers on a smaller scale and the suggestion of a bipolar-like outflow embedded in a larger CSE (e.g. Richards et al. 1996).

Fig. 6 Comparison of OH maser emission and our NOEMA observations. The colour maps show moment-0 maps of CO J = 2 - 1 (top) and SO2 JKa,Kc = 111,11 - 100,10 (bottom) line emission (colour maps) constructed for the velocity ranges specified in each panel's lower right corner, in accordance with the ranges covered by the maser emission. Contours show the NOEMA continuum emission (white) and the OH maser emission from Etoka & Diamond (2004) in cyan. These are, from left to right: red-shifted OH 1612 MHz, blue-shifted OH 1612 MHz, and blue-shifted OH 1665 MHz. The beams of the line emission measured with NOEMA are given in the bottom-left corner of each panel as white ellipses.

5.2 Mass-loss history

The interpretation of the continuum components A and B as dust clumps implies directional, enhanced mass loss that has been sustained over significant timescales, combined with very efficient dust formation. Assuming a typical terminal velocity of 30 km s⁻¹ and keeping in mind that we only measure a projected distance, such outflows could have occurred over timescales of 60 and 250 years for components A and B, respectively. However, as indicated in the previous section, the observations indicate a higher outflow velocity close to the star, making such estimates highly uncertain. At the same time, we note that the MIR emission presented by Schuster et al. (2009) partially coincides with our component A and the ‘bridge' between our components A and B, hinting at a potentially more continuous process.

The timescales connected to convective cells, Alfvén waves, and hot spots, where lifetimes of variation fall between 150 days and 3-6 years (Hartmann & Avrett 1984; López Ariste et al. 2018; Kervella et al. 2018; López Ariste et al. 2022), are far shorter than would be required to explain the formation of the observed components A and B. A strong magnetic field, as observed for VY CMa and NML Cyg, could possibly support the formation of such large clumps (O'Gorman et al. 2015; Vlemmings et al. 2002, 2017). Of course, higher-angular resolution observations might reveal these components to be made up of multiple smaller dust clumps and properties and formation timescales will need to be revised.

In contrast to these localised (limited opening angle) massloss enhancements, the likely presence of arcs in the outflow (see Sect. 4.2.3) is suggestive of repeated large-scale density or massloss enhancements, that is, covering a large part of the stellar surface. As mentioned, these might be a consequence of binary interaction, but based on the currently available data set, we are not able to constrain the properties of such a potential binary system.

5.3 Considering the origins of arcs and shells

In Sect. 4.2.3 we estimated a declining temporal spacing between arcs, ranging from 2800 yr to only 220 yr. Although very long periods seem unlikely for a binary system containing an RSG (e.g. Patrick et al. 2022), the potential tightening of a binary orbit, ultimately resulting in a binary engulfment and a merger/common-envelope event is considered in theoretical models. Recently, Landri & Pejcha (2024) explored the effect of a binary companion (2 M_⊙) on an eccentric orbit around an RSG (20 M_⊙), which grazes the RSG and finally plunges into it after a number of orbits. Although VY CMa was considered a reference object for this study because of the high level of complexity in its CSE, observations of this object have not revealed similar structures to those produced in the simulations thus far. The large-scale asymmetries of the CSE and the arcs and shells revealed by our observations of NML Cyg, on the other hand, do exhibit reasonable similarity to the simulated density structures. However, we did not trace the same velocity structure as that simulated by Landri & Pejcha (2024, their Fig. 8). This could be a consequence of the achieved angular resolution, preventing us from detecting more tightly wound arcs, or an indicator that this simulation is not an adequate representation of the NML Cyg system. The separation between the detected arcs or shells of the order of a few hundred years is indeed considerably different from the 20-40 year orbit of the simulated system and such a binary separation might induce different density and velocity structures. However, the tightening of the density enhancements we tentatively detect in our NOEMA observations are in line with a continuously shrinking orbit. We might hence be witnessing the inspiral of a binary companion into the RSG. The angular resolution of the observations (∼0″.4) translates to a physical size scale of ∼600 AU, or a timescale of 50-100 yr, for a 35-60 km s⁻¹ outflow velocity. This means that there could be additional arcs in the CSE following the same tightening pattern as described above, but which we cannot spatially resolve. Higher-angular resolution observations will be critical to test the plausibility of this ongoing binary inspiral scenario.

As an alternative to the impact of a binary companion, we could consider the mass-loss enhancements as intrinsic to the star, similarly to detached shells in AGB stars (Olofsson et al. 1996; Maercker et al. 2024). These are brief mass-loss enhancements that arise as a consequence of thermal pulses. In these AGB stars, only one detached shell is ever seen at a time in CO, because of the long inter-pulse periods (of the order of 10^4-5 yr). In models for super-AGB stars, on the other hand, the inter-pulse periods go down to afew 10² yr (Siess 2010; Doherty et al. 2014), with decreasing inter-pulse periods for an increasing initial stellar mass. However, the 25 M_⊙ initial mass estimate for NML Cyg (Zhang et al. 2012b) is at odds with the super-AGB suggestion, given that the current framework of stellar models puts an upper limit to the initial mass of super-AGB stars to roughly 10M_⊙ (Doherty et al. 2015). Furthermore, the age of NML Cyg for an assumed 25 M_⊙ initial mass (8 Myr; Zhang et al. 2012b) is already significantly higher than the estimated age of the Cyg OB2 cluster (2-3 Myr). This problem is only aggravated if the star had a much lower initial mass. In addition, the most massive super-AGB models by Doherty et al. (2015, for 9.8 M_⊙ initial mass) reached a maximum luminosity of 1.11 × 10⁵ L_θ, whereas Zhang et al. (2012b) estimated L_⋆ = 2.7 × 10⁵ L_θ for NML Cyg. In conclusion, we deem it unlikely that NML Cyg is a superAGB star rather than an RSG if the distance and luminosity estimates by Zhang et al. (2012b) are correct. Thus, we reject the characterisation of the arcs as detached-shell analogues.

5.4 Comparison to other RSGs - VY CMa in particular

Previous studies of RSGs have revealed complex outflows around VY CMa, μ Cep, Betelgeuse, and Antares (Ohnaka 2014; Kervella et al. 2011, 2018; Kaminski 2019; Montargès et al. 2019; Humphreys & Jones 2022; Murakawa et al. 2003; Tabernero et al. 2021). The similarities between the stellar properties ofVY CMa and NML Cyg (Zhang et al. 2012a,b) call for a direct comparison between their circumstellar environments, in particular.

The amount of dust present around NML Cyg appears to be similar to that around VY CMa. Similarly to the case of VY CMa, the continuum emission towards NML Cyg is dominated by circumstellar contributions that reveal substructures asymmetrically offset from the star. O'Gorman et al. (2015) reported 4 × 10⁻⁵ M_⊙ of hot (∼1000 K) dust close to the star and 2.5 × 10⁻⁴ M_⊙ of warm (<450 K) dust at a distance of at least 400 AU from the star, based on observations at ∼325 GHz and ∼659 GHz. Adding measurements at ∼178 GHz, Vlemmings et al. (2017) reported a dust mass of at least 1.2 × 10⁻³ M_⊙ with an upper limit to the dust temperature of ∼100 K. Humphreys et al. (2024) reported ALMA detections at ∼250 GHz of ≈ 2.5 × 10⁻⁴ M_⊙ of dust at temperatures in the range 300-500 K, which are not detectable in the optical or infrared. The dust-mass estimates we report for NML Cyg suggest at least the same amount of dust close to the star, but possibly considerably larger masses in case of optically thick dust (see Table 2).

The extent of NML Cyg's CO envelope also exceeds that of VY CMa. Recent observations with ALMA of VY CMa (Singh et al. 2023) have revealed an arc of emission of CO(2 - 1) and HCN(3 - 2) located at roughly 9″ offset NE from the star, considerably further out than the roughly 5″ traced by the CO(3 - 2) emission presented by Kaminski et al. (2013), which might have been affected by the filtering out of some large-scale emission. As reported in Sect. 4.1, we find a maximum CO extent for NML Cyg of ≈ 3 × 10¹⁷ cm in the SE direction. Taking into account the distance differences between VY CMa and NML Cyg, we traced CO emission out to distances up to a factor 1.7 larger than for VY CMa. Similarly, the NML Cyg CO emission-region size significantly exceeds the location of the furthest clumps of CO emission reported by Montargès et al. (2019, ≈4 × 10¹⁶ cm) for μ Cep.

Recently, high-angular resolution ALMA observations of VY CMa revealed collimated Hubble-type outflows accelerated to ≈ 100kms⁻¹, carving through inhomogeneous large-scale, slower (30-40kms⁻¹) ejecta (Quintana-Lacaci et al. 2023). The current observations of NML Cyg have not revealed equally collimated streams at such high velocities, but they do indicate that there are components in the circumstellar environment moving at higher velocities than the roughly 30 km s⁻¹ previously reported for this star (Singh et al. 2022; Andrews et al. 2022). Higher-angular resolution follow-up is needed to further localise and characterise these high-velocity components. We note that in Sect. 4.3, we report a velocity gradient across the detected arcs around NML Cyg, with the velocity increasing towards the star. We are not aware of similar results around other RSGs.

5.5 Impact on SN explosion

A dense circumstellar medium surrounding the collapsing RSG is required to explain the observations of many SNe II (Morozova et al. 2018; Yaron et al. 2017). Furthermore, asymmetries in the pre-SN circumstellar medium have a significant effect on the early detectability of the SNe via the probability of X-ray emission escaping the system and on the SN light curve. For example, Singh et al. (2024) and Shrestha et al. (2025) reported the need for asymmetry in the dense circumstellar matter and clumpiness in the large-scale circumstellar environment of SN 2013ixf to explain the evolution of the polarisation of the SN emission in the first ≈80 days after core collapse. Jencson et al. (2023) described the progenitor to SN 2013ixf as an RSG of 10⁵ L_⊙, with a pulsation period of ≈1120 days and a pulsation amplitude of 0.6mag at 3.6 μm, no bright pre-SN outbursts, $T_{\mathrm{e}\mathrm{f}\mathrm{f}}={3500}_{-1400}^{+800}$ K, M_init = 17 ± 4M_⊙, and a pre-SN mass-loss rate of 3 × 10⁻⁵ - 3 × 10⁻⁴ M_⊙ yr⁻¹ for an assumed wind velocity of 10 km s⁻¹ between 3 and 19yr before explosion. These stellar and wind properties are very close to those of VY CMa and NML Cyg (Zhang et al. 2012a,b).

The model presented by Singh et al. (2024) includes a confined dense circumstellar component with an r⁻² density profile, created by a short-lived mass loss at a rate of 10⁻² M_⊙ yr⁻¹ in the inner 5 × 10¹⁴ cm. However, there is no observational evidence for pre-SN episodic bursts (Singh et al. 2024, and references therein). At larger radii, out to 10¹⁶ cm, the material can be modelled with an r⁻³ density profile, ejected from the star at an average 10⁻⁴ M_⊙ yr⁻¹ . Even though we observe somewhat higher expansion velocities in the material surrounding NML Cyg, with velocities possibly reaching up to ≈60kms⁻¹, the average gas mass-loss rate of a few 10⁻⁴ M yr⁻¹ (Andrews et al. 2022) is in line with this pre-SN model. In addition, we find strong asymmetries in the CSE of NML Cyg, in the form of rather massive, localised dust clumps, which could eventually give rise to a polarisation behaviour that is similar to what has been observed for SN 2013ixf.

We note that Singh et al. (2024) mentioned the possible scenario of a binary component impacting the circumstellar morphology and the pre-SN mass stripping from the RSG. We suggest that we could be seeing early stages of this in the case of NML Cyg, considering that the arcs reported and discussed in Sects. 4.2.3 and 5.3 are possibly a consequence of a binary companion interacting with the star and/or the stellar outflow.

6 Conclusion

We present NOEMA observations of the red supergiant NML Cyg around 230 GHz (λ = 1.3 mm) at an angular resolution of ≈0″.4. Our observations reveal the presence of two bright continuum components that we derived to be the emission from warm dust close to, but significantly offset from the star. We estimated dust masses of roughly 10⁻³ M_⊙ for both components and temperatures of the order of 440 K and 265 K. We find that the cooler component is likely to be located in the red-shifted outflow and, hence, behind the plane of the sky that cuts through the star - and vice versa for the warmer component.

We detected ≈70 spectral lines across the observed spectrum (214.2-222.3 GHz and 229.7-237.8 GHz) and present maps for some transitions of selected molecules, such as CO, SO₂, and H₂S. The molecular lines trace multiple components in the outflow, including blobs, arcs, and shells and reveal higher-velocity components than previously identified for NML Cyg. A comparison to previously published OH maser maps confirms our assumption with respect to the stellar position and highlights a strong SE-NW orientation in the system. We propose that some of the observed features could be linked to an interaction with a hitherto unknown binary companion located in an orbital plane oriented along that same SE-NW direction.

It is likely that higher-angular resolution observations will reveal a more complex substructure of both the dust and gas in the circumstellar environment of this source and that masses, temperatures, and timescales will need to be revised. Observations at other wavelengths will help constrain the dust properties.

Acknowledgements

This work is based on observations carried out under project number w20bi001 with the IRAM NOEMA Interferometer as well as complementary single-dish observations carried out with the 30m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). We thank Ka Tat Wong for his work as the local contact astronomer and his help and advice in the reduction and calibration of the observations. EDB and HA acknowledge financial support by Chalmers Gender Initiative for Excellence (Genie). G.Q. would like to thank funding support from Spanish Ministerio de Ciencia, Innovación, y Universidades through grant PID2023-147545NB-I00.

Appendix A Supporting material

Across the two covered frequency ranges 214.2 - 222.3 GHz (LSB) and 229.7 - 237.8 GHz (USB), we detect 67 spectral features from 15 species (26 species when counting isotopologues separately), of which about half are new detections compared to the observations presented by Singh et al. (2021) and Andrews et al. (2022). Figure A.1 shows the spectra across the full frequency coverage towards the star and for a 1″ aperture centred on the star. Five features in these spectra could not be unambiguously identified with a carrier molecule and are marked as unidentified (U) in Fig. A.1.

Figures A.2–A.8 show channel maps of selected emission lines of ¹²CO, ¹³CO, SiO, SO, SO₂, and H₂S and Fig. A.9 shows moment-0 maps of selected emission lines of SiO, SiS, and H₂O that support the discussions in the paper.

Fig. A.1 NOEMA spectra of NML Cyg: continuum-subtracted LSB (top panels) and USB (bottom panels) spectra extracted for the central beam (red) and for a 1″ aperture centred on the star (black). Molecular species responsible for the detected spectral lines are indicated in blue at the respective transitions' rest frequencies, unidentified features are marked in red as “U”.

Fig. A.2 Channel maps showing emission from CO J = 2 - 1 across the full field of view. The flux scale is cut off for the sake of visibility. Note: the channels at 5 km s−1 and 8 km s−1 are contaminated by ISM contributions spread across the entire field of view. The stellar position is indicated with a black +.

Fig. A.3 Channel maps showing emission from 12CO J = 2 - 1 on a smaller scale than in Fig. A.2. The flux scale is cut off for the sake of visibility. Contours overlaid show the continuum emission (white) at [5, 10, 20, 30, 40, 60, 80, 100]σ. The stellar position is indicated with a black +.

Fig. A.4 Channel maps showing emission from 13CO J = 2 - 1. The flux scale is cut off for the sake of visibility. Contours overlaid show the continuum emission (white) at [3, 5, 10, 30, 50, 100]σ. The stellar position is indicated with a black +.

Fig. A.5 Channel maps showing emission from SiO v = 0, J = 5 - 4. The flux scale is cut off for the sake of visibility. Contours overlaid show the continuum emission (white) at [3, 5, 10, 30, 50, 100]σ. The stellar position is indicated with a black +.

Fig. A.6 Channel maps showing emission from SO J = 56 - 45. The flux scale is cut off for the sake of visibility. Contours overlaid show the continuum emission (white) at [3, 5, 10, 30, 50, 100]σ. The stellar position is indicated with a black +.

Fig. A.7 Channel maps showing emission from SO2 J = 111,11 - 100,10. The flux scale is cut off for the sake of visibility. Contours overlaid show the continuum emission (white) at [3, 5, 10, 30, 50, 100]σ. The stellar position is indicated with a black +.

Fig. A.8 Channel maps showing emission from H2S J = 22,0 - 11,1. The flux scale is cut off for the sake of visibility. Contours overlaid show the continuum emission (white) at [3, 5, 10, 30, 50, 100]σ.

Fig. A.9 Moment-0 maps for selected emission lines. From left to right: overall moment-0 map, blue-shifted, central, and red-shifted emission, covering the velocity ranges indicated in the bottom right corner of each panel. The stellar position is overlaid as a white star, and continuum emission plotted as white contours at [3, 5, 10, 30, 50, 100]σ.

References

All Tables

All Figures

Fig. 1 Continuum emission from NML Cyg. (Top) Emission in LSB (colour map) and USB (white contours at 3, 5, 10, 30, 50, 100, 200 σ), with two continuum components labelled as A and B and the stellar position as a black star. The synthesised beam is shown in the bottom left corner (LSB: yellow, USB: white). Residuals in USB after subtraction of the IMFIT results are shown in green contours (at -10, -5, -3, 3, 5, 10 σ) with dashed contours for negative signal. (Bottom) Estimating the stellar position. Peak flux positions of selected molecular transitions with errors according to the 50 mas astrometric accuracy (filled circles) at systemic velocity overlaid on USB continuum contours (grey) from the top panel. Also indicated: continuum peak (purple X), stellar position extrapolated from Zhang et al. (2012b) (black star), peak position of HST 0.55 μm emission (red star), and LSB and USB beams (cyan and black ellipses).

In the text

Fig. 2 Pencil-beam spectra towards NML Cyg, extracted towards the stellar position (left) and towards the peak positions of components A (middle) and B (right). The spectra have been normalised to their peak and a vertical offset has been added. The transitions are labelled in each panel. Note: the spectrum of PN 5 - 4 towards component B has been scaled for the sake of visualisation. Shaded regions indicate velocity ranges of interest with respect to absorption and emission features, with the relevant boundaries indicated in the top of each panel.

In the text

Fig. 3 Blob of CO emission to the SW of the star. The panels show of the 12CO J = 2 - 1 emission (colour) and 13CO J = 2 - 1 emission (white contours at [3, 5, 10]σ) across the velocity range [17; 26] kms−1. The extent of the colour scale has been limited for the sake of visibility. The stellar position is indicated as a black cross (+), the NOEMA beam is plotted in the bottom left corner as a white filled ellipse.

In the text

Fig. 4 Position-velocity diagrams. Cuts are centred on the assumed stellar position (marked as X) and at position angles of 65° to cut through the SW clump (top) and perpendicular to that, at 155° (bottom). Contour levels are drawn at [3,5,10,15, 20,40,80,100] × σ with σ = 2.0 mJy/beam; for 12CO, we only show the 3σ contour (red). Note: the CO emission at 5kms−1 and 8 km s−1 is contaminated by ISM contributions.

In the text

	Fig. 5 Same as Fig. 4a, but omitting SO2 contours and adding spherical shells centred on the star with radii 1″.0, 2″.5, 4″.0, 7″.5, and radially outward directed velocities of 35kms−1, 30kms−1, 25 km s−1, and 22 km s−1 , respectively (dashed cyan lines).
In the text

Fig. 6 Comparison of OH maser emission and our NOEMA observations. The colour maps show moment-0 maps of CO J = 2 - 1 (top) and SO2 JKa,Kc = 111,11 - 100,10 (bottom) line emission (colour maps) constructed for the velocity ranges specified in each panel's lower right corner, in accordance with the ranges covered by the maser emission. Contours show the NOEMA continuum emission (white) and the OH maser emission from Etoka & Diamond (2004) in cyan. These are, from left to right: red-shifted OH 1612 MHz, blue-shifted OH 1612 MHz, and blue-shifted OH 1665 MHz. The beams of the line emission measured with NOEMA are given in the bottom-left corner of each panel as white ellipses.

In the text

	Fig. A.1 NOEMA spectra of NML Cyg: continuum-subtracted LSB (top panels) and USB (bottom panels) spectra extracted for the central beam (red) and for a 1″ aperture centred on the star (black). Molecular species responsible for the detected spectral lines are indicated in blue at the respective transitions' rest frequencies, unidentified features are marked in red as “U”.
In the text

	Fig. A.2 Channel maps showing emission from CO J = 2 - 1 across the full field of view. The flux scale is cut off for the sake of visibility. Note: the channels at 5 km s−1 and 8 km s−1 are contaminated by ISM contributions spread across the entire field of view. The stellar position is indicated with a black +.
In the text

	Fig. A.3 Channel maps showing emission from 12CO J = 2 - 1 on a smaller scale than in Fig. A.2. The flux scale is cut off for the sake of visibility. Contours overlaid show the continuum emission (white) at [5, 10, 20, 30, 40, 60, 80, 100]σ. The stellar position is indicated with a black +.
In the text

	Fig. A.4 Channel maps showing emission from 13CO J = 2 - 1. The flux scale is cut off for the sake of visibility. Contours overlaid show the continuum emission (white) at [3, 5, 10, 30, 50, 100]σ. The stellar position is indicated with a black +.
In the text

	Fig. A.5 Channel maps showing emission from SiO v = 0, J = 5 - 4. The flux scale is cut off for the sake of visibility. Contours overlaid show the continuum emission (white) at [3, 5, 10, 30, 50, 100]σ. The stellar position is indicated with a black +.
In the text

	Fig. A.6 Channel maps showing emission from SO J = 56 - 45. The flux scale is cut off for the sake of visibility. Contours overlaid show the continuum emission (white) at [3, 5, 10, 30, 50, 100]σ. The stellar position is indicated with a black +.
In the text

	Fig. A.7 Channel maps showing emission from SO2 J = 111,11 - 100,10. The flux scale is cut off for the sake of visibility. Contours overlaid show the continuum emission (white) at [3, 5, 10, 30, 50, 100]σ. The stellar position is indicated with a black +.
In the text

	Fig. A.8 Channel maps showing emission from H2S J = 22,0 - 11,1. The flux scale is cut off for the sake of visibility. Contours overlaid show the continuum emission (white) at [3, 5, 10, 30, 50, 100]σ.
In the text

	Fig. A.9 Moment-0 maps for selected emission lines. From left to right: overall moment-0 map, blue-shifted, central, and red-shifted emission, covering the velocity ranges indicated in the bottom right corner of each panel. The stellar position is overlaid as a white star, and continuum emission plotted as white contours at [3, 5, 10, 30, 50, 100]σ.
In the text