Molecular gas in MGE 042.0787+00.5084

A slowly expanding torus associated with the candidate LBV MGE 042.0787+00.5084

Cristobal Bordiu, J. Ricardo Rizzo, Alessia Ritacco
Centro de Astrobiología (INTA-CSIC), Ctra. M-108, km. 4, 28850 Torrejón de Ardoz, Madrid, Spain
Institut de RadioAstronomie Millimétrique (IRAM), Granada, Spain
Accepted XXX. Received YYY; in original form ZZZ

The luminous blue variable (LBV) phase is a poorly understood stage in the evolution of high mass stars, characterized for its brevity and instability. The surroundings of LBV stars are excellent test beds to study their dense stellar winds and eruptive mass-loss events. Aiming to improve our knowledge of the LBV phase, we observed the and lines of CO and CO in a field of around the recently identified candidate LBV MGE 042.0787+00.5084, using the IRAM 30-m radio telescope. We report the first detection of molecular emission associated with this source, tracing a structure with an evident circumstellar distribution. Morphology and kinematics of the gas can be explained by an expanding torus, a structure that may have originated from stellar ejecta or the action of stellar winds onto the parent molecular cloud. We derive the physical properties of the gas by means of LTE and non-LTE line modelling, obtaining densities of H in the order of cm and kinetic temperatures below 100 K. In addition, we build a kinematic model to reproduce the structure and velocity fields of the gas, which is in good agreement with the observations. We estimate a total molecular gas mass of M and a dynamical age of years, leading to an average mass-loss rate of 0.8–1.2 M yr.

stars: mass-loss – stars: massive – stars: evolution – ISM: molecules – ISM: kinematics and dynamics – ISM: clouds
pubyear: 2018pagerange: A slowly expanding torus associated with the candidate LBV MGE 042.0787+00.5084B

1 Introduction

High mass stars strongly disturb the composition and structure of the interstellar medium by delivering large amounts of thermal and mechanical energy into their surroundings.

Before exploding as core-collapse supernovae, high mass stars evolve through a series of short-lived and violent phases (Langer1994). Among these transitional stages, the so-called luminous blue variable phase (hereafter LBV) is possibly the most puzzling due to its intrinsic instability. LBV stars are luminous and hot supergiants characterized by dense winds that sustain high mass loss rates –sometimes reaching  M yr(VanBoekel2003)–. This steady mass loss is occasionally accompanied by violent outbursts that eject the outer stellar layers into the ISM, like the eruptions of Car in the 19th century (Davidson1989). Consequently, most LBV stars are surrounded by expanding structures of dust and gas. These circumstellar structures typically have very different shapes, from simple spherical shells to more complex, axisymmetric arrangements, including bipolar nebulae and rings (Marston2006). With characteristic times of years (Humphreys1994), the galactic population of LBV stars is extremely scarce, with less than twenty confirmed members (Clark2005) and a few candidates which still require confirmation (Stringfellow2012). Circumstellar envelopes are also indirect probes of the existence of these sources, because they may contain information about past mass-loss events. (Mizuno2010; Gvaramadze2012).

LBV stars have been widely studied at optical and infrared wavelengths, with most of the research focusing on the determination of stellar parameters (temperatures, luminosities, metallicities and variabilities) and, to a lesser extent, the properties of the associated nebulae (neutral and ionized gas, and warm dust, Umana2011; Ingallinera2014). However, most of the underlying physical processes during the LBV phase are far from being totally understood. Little is known about the mechanisms that drive its variability and violent outbursts, and its actual role in the late evolution of high mass stars is unclear, with recent observations suggesting that LBV stars may be (under certain conditions) direct progenitors of type II supernovae (Groh2013). For these reasons, the current picture of the LBV phase is rather incomplete. As many LBVs are surrounded by complex and multi-phase envelopes, targeting other components associated with these stars may provide valuable information to tackle some of these questions. In this regard, molecular gas is specially worth of attention. This component has been historically overlooked, as molecules were thought to be rapidly destroyed by the strong UV radiation in the outskirts of LBV stars. Despite of this, a series of pioneering studies have resulted in the detection of CO and NH around a few LBVs. These detections prove that significant amounts of molecular gas can form and survive for long time-scales around sources of this type.

The best example is G79.29+0.46, a LBV nebula where Rizzo2008 reported nested shells of warm CO, strikingly attached to the infrared nebula and tracing successive mass-loss events. In a further study, Rizzo2014 found several NH warm clumps clearly associated with the star, proving that the chemistry in the environs of LBV stars is more complex than initially assumed. These results underlined the potential of dynamical studies of molecular gas associated with LBV stars to achieve a complete vision of this elusive phase.

A particularly interesting but still unexplored case is that of MGE 042.0787+00.5084, a recently identified candidate LBV star. This source exhibits a nearly spherical circumstellar envelope with intense infrared emission, that was initially discovered in the Spitzer MIPSGAL 24 m survey (Mizuno2010). The central star was later catalogued as a candidate LBV by means of near-infrared spectroscopy, which revealed strong hydrogen emission lines (including Pfund and Brackett series) and several metal lines (\ionMgii, \ionNaii, \ionFeii) usually found in other LBV spectra (Flagey2014). On the other hand, recent 6-cm VLA observations (Ingallinera2016) unveiled a clumpy radio nebula strongly coupled to the infrared shell, with hints of a differentiated spectral index, as in other well known LBV sources such as HR Car (Buemi2017).

In this paper we report the first detection of molecular material around MGE 042.0787+00.5084, using IRAM 30-m observations of CO and CO and lines. In Sect. 2 the observing setup and procedures are detailed. In Sect. 3 we present the results of the observations. In Sect. 4 we estimate the physical parameters of the gas and provide an interpretation for the observed structure based on the derived mass-loss rate, which is compatible with the LBV nature of the object. We summarize our findings in Sect. 5.

2 Observations

We observed MGE 042.0787+00.5084 with the 30-m radio telescope of the Institut de RadioAstronomie Millimétrique in Granada (Spain). Observations took place the night of 2017 July 24, under good summer conditions (opacity , corresponding to 4.2 mm of precipitable water vapour at 225 GHz), with a total observing time of 3 hours.

Telescope was configured to operate in dual-receiver mode to map the distribution of CO and CO at 1.3 and 3 mm simultaneously, while covering other molecules of interest such as CO and CO. EMIR receivers E090 and E230 were used with two spectral setups: a first setup addressing CO , CO and CO , and a complementary setup specifically targeting CO . The selected backend was the FTS spectrometer, providing an instantaneous bandwidth of 4 GHz per polarization and a velocity resolution of 0.26 and 0.52  km s at 1 and 3 mm respectively.

On-the-fly (OTF) mapping in zig-zag mode was used to make maps around the source, using a position-switching strategy with a common reference for all subscans. Two orthogonal mapping directions were used for each setup, with 21 parallel subscans per direction. Subscan length was set to 90″ with a mapping speed of 3″ s, using a Nyquist-compliant spacing of 43 between subscans. The reference was located 10 away of the source to minimize background/foreground contamination from the galactic plane. Map resolution corresponded to the half-power beam width of the telescope, which is and at 1 and 3 mm respectively. In addition, a deep integration towards the star’s position was done to search for other, less abundant molecules.

We used Saturn and the quasar 1749+096 for pointing correction, and G34.3+0.2 for line calibration. Pointing was regularly checked, while calibrations were run every fifteen minutes and between consecutive maps. The pointing accuracy was always better than 5″.

Line Freq V HPBW rms
[GHz] [ km s] [] [mK]
CO 115.27 0.51 21.3 25
CO 110.20 0.53 22.3 10
CO 230.54 0.26 10.7 30
CO 220.40 0.27 11.2 20
Table 1: Observational parameters

A summary of the observational parameters for the target lines is provided in Table 1. Individual spectra were baseline-subtracted and combined to produce data cubes for each rotational transition. The data reduction and analysis was carried out with gildas, a software package especially tailored to deal with native 30-m data (Guilloteau2000). Hereafter in this work we adopt the following conventions: (1) data is presented in main beam temperature ( scale unless specified otherwise. Conversion from antenna temperature () is trivial via the relation , with


where and denote the antenna main beam and forward efficiencies respectively, expressed as a function of the frequency for the IRAM 30-m telescope (SanchezContreras2015); (2) velocities are expressed with respect to the local standard of rest (); (3) positions are given as offsets relative to the J2000 equatorial coordinates of the source, which are (RA, Dec.) = (19°06′2457, +08°22′019); and (4) position angles are measured from north to east.

3 Results

After a thorough search in the whole data cubes we found significant emission of CO and CO at different velocities in the range from 10 to 80  km s. Figure 1 displays the spectra of the four transitions spatially averaged over the whole field. CO was detected only towards the southeast of the source, and CO was not detected at any position.

A detailed analysis of the spatial distribution of the emission reveals that, in the four CO and CO lines, most of the velocity components are present over a significant part of the observed field, without a clear pattern with respect to the star; these components are presumably not related to the source, being most likely foreground/background contamination arising from molecular clouds in the line of sight. In contrast, the velocity component in the range (+13,+19)  km s, marked with vertical lines in the spectra, depicts certain symmetry around the central star, in the form of a rather elongated ring-like structure. This component is represented as channel maps in Figs. 2 to 5, in steps of 0.5  km s. Several arcs more or less centred at the star position are noted; these arcs present different sizes at different velocities, moving from southwest to northeast as velocity increases. At 15–15.5  km s  the emission is closed in elliptical features, especially noted in the CO line.

Figure 1: CO and CO spatially averaged spectra towards MGE042.0787+00.5084. CO spectra is scaled by a factor of 6 for easier comparison. Hanning smoothing has been applied to reduce the noise. The short vertical line on top of each spectra marks the component in the range (+13,+19) km s.

Figure 6 shows the corresponding velocity-integrated intensity maps, as contours superimposed on the 24 m Spitzer MIPSGAL image of the infrared nebula (taken from Mizuno2010). These maps reveal a rather clumpy ring-like structure centred at the star position, displaying a remarkable elongation in the approximate southeast-northwest direction. By fitting an ellipse to the emission peaks as shown in the figure, we derive an angular size of with a position angle of 150.

Figure 2: CO emission towards MGE 042.0787+00.5084 in the velocity range 12.75–19.35  km s. Absolute is shown in the top left corner of each panel. Colour scale represents in K. Contours start at 1 K with a spacing of 1 K. Crosses at (0,0) indicate the position of the star.
Figure 3: Same as Fig. 2 for CO , with a channel binning of 2 (resolution 0.5  km s) for easier comparison. Contours start at 1 K with a spacing of 1 K.
Figure 4: Same as Fig. 2 but for CO . Contours start at 100 mK with a spacing of 100 mK.
Figure 5: Same as Fig. 2 but for CO , with a channel binning of 2 (resolution 0.5  km s) for easier comparison. Contours start at 200 mK with a spacing of 200 mK.

The structure is mostly unresolved so we are not able to determine its actual thickness. Never the less, the 30-m angular resolution is sufficient to observe a certain degree of morphological differentiation, allowing us to identify a number of regions that might deserve further attention:

  1. The main ring-like feature, which completely surrounds the infrared nebula and accounts for most of the emission.

  2. The inner cavity, a region in which molecular material is likely depleted. The extent of this cavity strikingly matches the infared nebula.

  3. The partially isolated emission blob towards the northwest, hereafter the northwestern clump, projected from the star. This region is particularly prominent in the images.

  4. The southeast region, from which arises the most intense and widespread emission.

Figure 6: Velocity-integrated maps of CO and CO towards MGE042.0787+00.5084 in the range (+13.4,+17.4)  km s, shown as contours, superimposed on the Spitzer MIPSGAL 24 m image of the infrared nebula (Ingallinera2016). Contour starting values and steps are 6.5 and 0.6; 6.5 and 1.2; 1.35 and 0.15; and 1.5 and 0.25 K  km s  for panels a, b, c and d respectively. The grey circle in the bottom right corner of each panel represents the beam size. Dashed lines in panel b indicate the cuts for the position-velocity diagrams in Fig. 9, and the ellipse used for the fitting is shown for reference.

This southeast region, where the overall symmetry of the structure is lost, may be seriously contaminated by emission from a molecular cloud at a close . Figure 7 compares the spectra of CO and CO averaged over three different parts of the map: the clump, the stellar position and the southeast region. CO is only detected in the latter, as a weak and narrow component in the velocity range (+25.6,+27.8) km s. The non-detection of CO in other positions suggests the existence of a chemical differentiation, indicating that this widespread emission may be –at least partially– disaggregated from the main ring-like structure.

Figure 7: Overabundance of CO. Left panel: Contour maps representing the emission of CO integrated in the velocity range (+25.6,+27.8) km s (solid) and CO in the velocity range (+13.4,+17.4) km s (dashed), superimposed on the Spitzer MIPSGAL 24 m image of the infrared nebula in gray scale. Contours starting values and steps are 0.4 and 0.1 and 1.35 and 0.15 respectively, in units of K km s. The labelled squares represent regions of where the spectra in the corresponding right panels have been spatially averaged. Right panels: averaged spectra of CO (black) and CO (filled) over the region with the corresponding label. The former has been scaled by 2 for easier comparison. The velocity components represented by the contour maps are marked with vertical lines. Intensity scale in units of K.
Figure 8: Hints of substructure in the CO gas: (a) the two velocity components as contours superimposed to the 24 m image of the infrared nebula in gray scale. In dashed contours, the component at = 15.4  km s, with levels starting at 6.5 in steps of 1.2 K km s; in solid contours, the component at = 17.8  km s, with levels starting at 2.5 with steps of 0.6 K km s. (b) spectra at position (20″, 10″). (c) spectra at star’s position (0″, 0″). (d) spectra at position (-20″, -10″). A gaussian fitting for each line is plotted.

As shown in Fig. 8, a closer inspection of spectra in individual positions reveals that this structure is the blending of two different components at = 15.4 and 17.8  km s, depicted in dashed and solid contours respectively. The component at 15.4  km s is relatively ubiquitous and completely surrounds the source in all directions, but the component at 17.8  km s is highly localized and only visible in the CO lines. It is prominent towards the East, with a peak temperature comparable to the other component, but almost non-existent in other positions (e.g. southwest). Similarly, its motion with respect to the star is unclear. This makes extremely difficult to confirm whether this component is related to the source in any way: it could be part of an interfering cloud or a mere stratification of the main component. Regardless of its nature, the contribution of this component is negligible, as the overall emission is dominated by the component at 15.4  km s. For this reason, we limit our analysis to the latter, adopting its central velocity as the systemic velocity of the structure.

Figure 9: Position velocity diagrams of the CO line in the slices defined in Fig. 6. The direction of the cut is indicated in the top left corner of each panel, and the vertical dashed line represents the position of the central star.

For a better understanding of the velocity structure of the CO and CO emission, we studied the kinematics of the gas by means of position-velocity (PV) diagrams along the strips indicated with dashed lines in Fig. 6, which were selected according to the axes of the ellipse previously used for the fitting. The first slice goes from SE to NW (P.A. 150°), following the direction of maximum elongation, and the second is orthogonal, from NE to SW (P.A. 60°). As seen in Fig. 9, we find emission highly concentrated in two isolated blobs either side of the star, each of them spanning for 2  km s. The velocity difference between the blobs is small regardless of the direction of the cut, in the range 0.5–1.5  km s from peak to peak.

Figure 10: Intensity-weighted velocity maps of CO and CO towards MGE042.0787+00.5084 shown as contours. Gray scale represents the intensity of each transition at the original resolution, in units of K  km s. Dashed blue contours indicate and solid red contours . Contour spacing is 0.15 and 0.25 K  km sfor and lines respectively. Border of the maps have been masked out to minimise noise.

4 Discussion

4.1 Morphology and kinematics

Despite the contribution from the southern cloud, most of the emission arises from a circumstellar ring-like feature that encloses the infrared shell. The PV diagrams in Fig. 9 showed that material in the ring moves at different velocities either side of the exciting star, but a deeper analysis is required to disentangle the motion of the structure. We resorted to study the first-order moment of the emission, i.e. the intensity-weighted velocity distribution, presented in Fig. 10, which approximates the velocity field of the gas. The four CO and CO lines share a common pattern, a smooth velocity gradient in the NE-SW direction that roughly splits the field across the semimajor axis of the structure. Therefore, the northeast and southwest regions show opposite motions, represented by two sets of contours in the maps: gas in the northeast moves away (redshifted emission) while gas in the southwest moves towards the observer (blueshifted emission).

The spatial distribution of the emission can be explained by a torus seen at a certain angle, but the velocity gradient observed is compatible with a structure undergoing expansion or contraction. However, near-infrared spectroscopic observations by Flagey2014 proved that MGE 042.0787+00.5084 is a massive evolved star, allowing us to confidently rule out the contraction scenario –i.e. material collapsing/accreting onto the star–, and therefore conclude that the gas is distributed in an expanding toroidal structure.

If the ring is expanding isotropically in a uniform medium, the observed elongation may be interpreted as a mere projection effect, which allows us to infer the inclination angle by fitting ellipses to the emission maxima as described in Sect. 3, and calculating the minor to major axis ratios. With this method we derive an inclination of . However, the overall morphology of the emission is heavily altered by the southern widespread component, and a local density gradient in the direction of the major axis may also contribute to the elongation, so we could be overestimating the inclination of the structure.

In this regard, we note similar velocity shifts in the two PV diagrams from Fig. 9, showing just a slightly larger difference in the NE-SW cut, which is likely more affected by projection according to the viewing geometry. This could be taken as an indicator of a low inclination. Therefore we cannot merely rely on the ellipse fitting to constrain the inclination angle of the structure, and in Sect. 4.3.4 we considered a larger range of possible inclinations.

We can use the NE-SW cut to measure the projected expansion velocity as half the velocity difference between the emission blobs. We estimate a projected velocity of  km s, which relates to the expansion velocity as


where is the inclination angle and the increment in velocity. Using the previously estimated inclination angle, we obtain 1.2–1.4  km s. This is an extremely low expansion velocity, an order of magnitude lower than those measured in other LBV sources (e.g. the main CO shell in G79.29+0.46 expands at 14  km s, Rizzo2008), but it is important to bear in mind that this result is just a coarse lower limit, strongly affected by the inclination angle. Lower inclinations would be compatible with higher velocities.

4.2 Physical parameters of the gas

Being MGE 042.0787+00.5084 a high-mass evolved star, mass-loss processes –both heavy steady winds and/or eruptive mass ejections– could have originated the structure seen in CO and CO. The simultaneous observation of two rotational lines from two isotopologues allows us to derive some physical parameters of the gas.

To do so, we performed LTE (local thermodynamic equilibrium) and non-LTE analyses in two separate regions, to account for the spatial differentiation observed in the structure: the torus as a whole, and the northwestern clump. The two regions are depicted in Fig. 11. We excluded most of the the southeastern region to minimize the contribution from the widespread emission.

Figure 11: Regions in which physical parameters were derived independently: the northwestern clump (blue) and the main torus (magenta). For reference, CO emission in the velocity range (+13.4, +17.4)  km s is shown in greyscale.

4.2.1 LTE analysis

As a first approach we determined the global physical parameters of the gas in the LTE scenario, using CO, the least intense and presumably least abundant observed species. By assuming optically thin emission, the column density in the lower level is given by


where is the frequency of the line, is the spontaneous decay rate, is the degeneracy of each rotational level and the integral term depends on the excitation temperature through the relation , which is true in the optically thin regime (Wilson2009). In LTE, levels are populated following a Boltzmann distribution described by a single . In order to relate to the total population of all energy levels in the molecule we use a rotational partition function such that,


where is the energy of the upper level and is the excitation temperature. The latter can be estimated from the line ratio by solving the relation


with the peak temperature of the transition. The derived excitation temperatures –between 9 and 12 K– and column densities – all in the range of cm– are summarized in Table 2. The column densities of each molecule do not vary significantly between the torus and the clump.

(CO) N(CO) (CO) N(CO)
[10 cm] [10 cm]
Clump 11.2 1.0 (0.2) 11.6 6.1 (0.6)
Torus 9.8 1.0 (0.2) 12.1 6.1 (1.1)
Table 2: LTE column density estimates.

4.2.2 Non-LTE analysis

To study a more realistic scenario we ran simulations with the non-LTE excitation and radiative transfer code RADEX (Tak2007). RADEX uses the escape probability formulation with the Large Velocity Gradient (LVG) approximation (Sobolev1960) to constrain the excitation conditions of the gas according to the observed line intensities. We involved the four lines in the fitting, convolving the maps to the resolution to obtain true line ratios.

The calculations were made taking the CO and CO collisional coefficients from LAMDA database (Schoier2005), with a background temperature of 2.73 K. First, we ran RADEX for a range of kinetic temperatures () from 10 to 250 K in steps of 10 K, solving the radiative transfer in a two-parameter space defined by –the H volume density– and , with the aim of finding the best fit for the observed CO intensities and ratios. Later we cross-checked the resulting values with the CO line intensities and ratios in order to determine (CO) and restrict the range of possible .

[ cm] [10 cm] [10 cm]
10 7.2 (1.4) 1.2 (0.1) 0.05 0.11 8.3 (1.2) 0.35 0.7
30 1.5 (0.2) 1.0 (0.2) 0.02 0.09 7.1 (0.9) 0.12 0.51
50 1.0 (0.1) 1.0 (0.2) 0.01 0.09 6.8 (1.0) 0.09 0.49
70 0.8 (0.1) 1.0 (0.2) 0.01 0.08 6.8 (1.0) 0.08 0.49
[ cm] [10 cm] [10 cm]
10 4.1 (1.2) 1.3 (0.2) 0.06 0.14 7.6 (1.7) 0.39 0.71
30 1.0 (0.2) 1.2 (0.2) 0.04 0.12 6.2 (1.8) 0.19 0.56
50 0.7 (0.2) 1.1 (0.2) 0.03 0.12 6.2 (1.8) 0.18 0.56
70 0.5 (0.1) 1.1 (0.2) 0.03 0.12 6.2 (1.5) 0.17 0.55
Table 3: RADEX results.

Simulation outcomes are compiled in Table 3 and some RADEX runs are presented in Fig. 12. We find CO column densities in the range 1.0–1.3 cm, which only differ by 10–20% from those of the LTE analysis. There is little variation between the column densities of the clump and the torus, with CO opacities supporting the optically thin assumption. Similarly, the CO and CO column densities are relatively insensitive to temperature, which is consistent with optically thin emission as well. Considering the uncertainties involved, (CO) is only 5–7 times higher than (CO), which indicates a surprisingly low isotopic ratio. The implications of this finding are discussed in Sect. 4.3.2.

Figure 12: CO RADEX fittings computed for K (left), K (centre) and K (right), in the clump (thick blue line) and the torus (thin magenta line); see Fig. 11. Solid curves represent the observed CO intensity, and dashed curves the CO ()/() ratio. Intersections between two sets of contours indicate the most likely values of and .

The resulting gas densities are quite low, of a few 10 cm, and we observe significant differences between the torus and the clump, with the latter being 1.5–2 times denser. These low densities are consistent with the non-detection of other high-density tracers, and somehow limit the range of possible kinetic temperatures, with solutions up to 70 K for the torus and 100 K for the clump. Above these values, falls well below the critical density of the transition.

In the case of CO we find higher values, suggesting that the observed intensities may be slightly affected by opacity, specially at low temperatures. In addition, we observe a significant dilution in the clump after convolution, which contrasts with the prominence of this region in the original maps. This may be related to a beam-filling issue, leading to a potential underestimation of the true densities. However, this effect is difficult to calibrate as the whole structure is unresolved. Interferometric observations would reveal the actual angular extent of the structure, allowing for a more accurate estimation of its physical parameters.

4.3 Origin of the molecular gas

The morphology and kinematics of the structure raise the question of the origin of the gas. While tori are a quite unusual arrangement for circumstellar material around evolved high mass stars, a few examples are found in the literature, such as the dusty torus around the LBV candidate HD168625 (OHara2003). There are two possible explanations for the formation of such a structure around a LBV star, namely: (1) the observed molecular gas is stellar ejecta from a process of non-isotropic mass loss. In this scenario the observed CO and CO formed when the stellar wind got colder and coagulated; and (2) the molecular material is the relic of the molecular cloud where the star was born, accumulated and compressed by the steady action of a non-spherical wind from a previous stage.

The total mass of molecular gas present in the CO structure would shed some light on its origin, since very different orders of magnitude are expected in each scenario. This mass can be derived from the column density of H via the relation


where is the distance to the object, the mass of the hydrogen atom, is a correction factor to account for the hydrogen relative abundance, and the last term is the integral of the measured H column density over the solid angle subtended by the source. Most of these factors are known or can be estimated with confidence, but the distance and the true H column density are subject to a great uncertainty.

4.3.1 About the distance to MGE 042.0787+00.5084

There are no available estimates in the literature for the distance to MGE 042.0787+00.5084. By taking the systemic velocity of the structure of 15.4 km s, we calculated two kinematic distances, namely 1.2 and 11.4 kpc. Details on how we determined these distances are provided in Annex A. We can expect moderate deviations from galactic rotation that would alter these estimates, but even departures of 20 km s would increase the near distance just by a factor of 2, which is not critical for the subsequent interpretation.

Although the two distances are equally possible, we found a number of qualitative arguments that support the near distance: MGE 042.0787+00.5084 is by far the strongest object in the sample studied by Flagey2014, being remarkably bright in bands J, H and K, and the symmetry and homogeneity of the nebula suggest a short age, which is more compatible with closer distances.

In order to narrow the range of possible distances, we explored the recently released data from ESA’s Gaia mission (Gaia DR2, GaiaCol2016, GaiaCol2018), finding a match within 2″ of the source position. This object (SourceID 4307460432455822976) has a parallax of mas according to Gaia. With such a large uncertainty, distances in the range 4–20 kpc are possible within the error margin. Still, it is possible to deal with inaccurate parallax measurements by means of the Bayesian inference method described by BailerJones2018, who applied this statistical approach on the whole DR2 catalog to determine geometric distances according to the following probability density function


where and are the parallax and its uncertainty respectively, is the parallax zero point and is a prior that depends on the galactic coordinates of the source.

This method has proven to be adequate for narrowing Gaia distances for other LBVs, reaching results in good agreement with literature values even when large uncertainties are involved (Smith2018). The inferred geometric distance for MGE 042.0787+00.5084 is thus kpc, but we must note that this result is based on a very inaccurate parallax and hence must be regarded with caution.

In the following we adopt a conservative distance of 2.5 kpc, as a compromise between the two approaches and taking into account possible departures from galactic rotation. None the less, we admit that a larger distance is possible (e.g. 10 kpc). The implications of such a long distance on the global parameters of the structure are discussed in Sect. 4.4.

Figure 13: Isotopic ratio maps shown in colour scale: (a) the CO/CO ratio and (b) the CO/CO ratio. The colour scale is the same for the two maps, with contours at 1, 3, 5 and 7. In each map the two lines have been convolved to the same beam size, represented as a gray circle in the bottom right corner. Noisy pixels have been masked out and spatial smoothing has been applied to avoid artefacts.

4.3.2 Isotopic ratios

Another essential factor when determining masses is (CO), the CO to H relative abundance, which allows to translate (CO) into (H). This parameter strongly depends on the C abundance. We observe [CO/CO] line intensity ratios as low as 5–7 in large regions of the structure, as shown in the line ratio maps of the and transitions (Fig. 13). According to the LVG results, opacity by itself is not enough to explain an order of magnitude decrease in the [CO/CO] ratio with respect to the typical ISM values of 70–90 (Wilson1992). This leaves just two possible explanations for such a low isotopic ratio: a CO depletion or a CO enhancement.

In the first scenario, chemical fractionation could account for the destruction of CO. Similarly, a CO enhancement could be a direct consequence of stellar evolution. It is known that the CNO cycle favours the production of C in intermediate-mass stars (Berdyugina1994). The stellar mass range in which this overproduction takes place is not completely constrained, but some studies have measured low [CO/CO] ratios in a number of massive stars (e.g. Ori, Lambert1984). If this process occurs in MGE 042.0787+00.5084 as well and the low [C/C] ratio is the result of the central star nucleosynthesis, it would be widespread in the surroundings, therefore indicating that a significant amount of molecular material is of stellar origin.

However, the non-detection of CO in the structure supports the CO overproduction hypothesis. We detected this isotopologue in the southeastern cloud, where abundances are close to the canonical values (e.g. [CO/CO] 5.5, see Fig. 7), but not in the structure, despite having CO and CO intensities and rms comparable to those in the cloud. This may imply an excess of CO. Never the less, other processes could be playing a role: high [CO/CO] ratios have been proposed as signposts of high FUV fields acting upon the outermost layers of molecular clouds (e.g. photon dominated regions) and causing selective photodissociation of some isotopes (Shimajiri2014). In this case, this effect would explain the apparent absence of CO, but should only be noticeable locally, in the densest parts of structure, particularly the northwestern clump.

In the map (panel b) we note a closed minimum toward the source. The spatial extent of this minimum is comparable to the beam, and the levels of CO and CO in this region are above 5 and 3 respectively. Therefore this minimum is probably a meaningful feature.

Unfortunately, the angular resolution of the 30-m radio telescope does not allow to study this feature in detail and discriminate between the two proposed scenarios. Interferometric observations would be valuable to assess if the low isotopic ratio observed is the result of stellar evolution or a localized effect of the strong UV field of the star. Taking into account the strong uncertainty that affects the (CO) factor, we use (CO) to estimate the H column density, adopting a canonical relative abundance of (CO) = 10.

4.3.3 Kinematic model

Using the global physical parameters as constraints, we developed a simple model to test expanding torus interpretation. The model aimed at reproducing the observed morphology and kinematics following a phenomenological approach. We used the pseudo-Monte Carlo non-LTE radiative transfer code LIME v1.8 (Brinch2010). LIME simulates photon propagation through unstructured 3D Delaunay grids by applying the SimpleX algorithm (Ritzerveld2006) and assuming ballistic transport. For a given molecule and transition LIME calculates level populations and generates synthetic spectral cubes, based on the LAMDA collisional rates. We processed these cubes with a custom Python pipeline to simulate the observation with the IRAM 30-m telescope, obtaining maps directly comparable to the original data. We used the CO line because it provides the highest angular resolution and signal-to-noise ratio.

We modelled the structure as an expanding torus characterized by a set of geometrical and physical parameters. A complete list of the parameters is provided in Table 4 and the model is described in detail in Annex B. Geometry is controlled by distance, inclination and position angle. We defined as a scaling radius coincident with the inner edge of the torus to parametrize density and temperature. Density is described by a toroidal function where the parameter , the half-opening angle, determines the oblateness of the structure and the dimensionless parameters and dominate the steepness of the density gradient in the radial and angular directions. Since the torus is unresolved in the 30-m data we do not know its actual thickness, so we set some of these parameters manually for convenience. For the temperature profile we used a segmented function, which remains constant at a reference temperature in the cavity and behaves as a radial power law for . Regarding the velocity field, we considered a purely radial expansion.

4.3.4 Model results

We adopted a distance of 2.5 kpc with a relative abundance (CO) = as explained above. The free parameters of the model are the reference density –related to –, the reference temperature –related to the kinetic temperature of the gas–, the expansion velocity and the inclination angle . To account for the uncertainty in the determination of the inclination angle and the possibility of a nearly face-on structure discussed in Sect. 4.1, we built three types of models, A, B and C, corresponding to inclinations of 15, 30 and 45° respectively. We used reference temperatures in the range 10–70 K and densities of H in the range 0.8–4 cm, in agreement with the LVG results (see Table 3). We also let the expansion velocity vary from 1 to 5  km s to cover possible scenarios where a low inclination somehow compensates for a higher velocity.

We explored the parameter space iteratively looking for the best fit in the velocity range (+13.4, +17.4) km s, based on a channel-wise criteria described in Annex B. The best-fitting set of parameters for each kind of model is presented in Table 5 and a comparison of the main observables (line intensity maps, velocity fields and position-velocity diagrams) is shown in Fig. 14. The three models show a preference towards warm, low density solutions, as they reach the best fit for cm and = 70 K. It is important to note that these are reference values, so the average and in the structure may be slightly lower. In any case, the best-fitting density is not particularly high, which could explain the nondetection of other less abundant molecules in the structure. Other denser and colder models fail to match the observed intensities, even by an order of magnitude.

Model A, which corresponds to an inclination of 15°, provides the best match, with a quite convincing reproduction of the morphology and kinematics of the gas. The expansion velocity of 3  km s is 2 times higher that the prior estimates from the PV diagrams in Fig. 9, but still much lower than typical values. The integrated intensity map is slightly weaker than the data in most of the structure, but this difference is an expected side effect of the simpleness of our torus model and the uniformity of the density distribution. Never the less, model A is able to reproduce with high fidelity the observed velocity gradient everywhere but in the cavity, which in our model is highly depleted. Models B and C, corresponding to higher inclinations, fail to reproduce the observed morphology as the resulting structure is exceedingly elongated. Similarly, while the velocity gradient approaches that of model A, the velocity extent of the emission blobs in the PV diagrams is smaller than expected, of just 0.75–1  km s, which means that the material in the line of sight is confined in a narrower velocity range.

4.4 Mass and mass-loss

In view of these results, we conclude that the structure of molecular gas seen around MGE 042.0787+00.5084 is a low-density, nearly face-on torus that expands slowly. Adopting a conservative of 50 K, we derive a total mass for the structure of M ( M in the torus and 0.1 M in the clump). This mass is probably rather low to be the remnant of a molecular cloud piled up by the stellar wind, but it is compatible with masses of molecular gas found around other LBV stars (Rizzo2008).

Therefore, assuming all the gas formed from stellar ejecta, a crude estimate of the mass loss rate of MGE 042.0787+00.5084 can be provided simply by calculating the dynamical age of the circumstellar structure. The torus deprojected radius is 15″ on average, which at 2.5 kpc translates into 0.18 pc. For 3  km s, the dynamical age of the structure is


The LBV phase typically lasts for a few years, so this dynamical age is compatible with a structure formed during this period. However, the age derived by this method must be regarded as a coarse upper limit, as the expansion velocity of the gas could have been significantly higher in the past, implying a younger structure.

Never the less this time-scale is useful to estimate the average mass-loss rate of the star during the LBV phase, yielding a value of 0.8–1.2  M yr. This result is compatible with evolutionary models (Langer1994), laying in the lower end of the LBV spectrum, where stars exhibit high temperatures and moderate mass-loss rates (Lamers1989).

Still, the mass-loss rate estimate provided here must be interpreted carefully. It is a time average over a period of years, so the star could have experienced brief periods of severe mass-loss in the past and we cannot use it to infer its current state. In addition, we only considered molecular gas in the calculation, which most likely is not the only contributor to the total mass. Other significant components (e.g. neutral and ionized gas) would increase the mass-loss rate and were not taken into account.

Finally, it is worth to mention that the derived time-scale for the molecular gas structure also constrains the age of the inner nebula. While the torus is the result of steady stellar winds, the infrared nebula may trace a much more recent eruptive event. The possibility of simultaneously observing two successive mass-loss events traced by different phenomena makes MGE 042.0787+00.5084 a very interesting target for further studies.

Figure 14: Comparison between data and best model for each inclination angle (see Table 5). Top row: CO integrated line intensity maps in the velocity range (+13.4, +17.4)  km s. Colour scale in units of K  km s, with contours starting at 5 with steps of . Middle row: Intensity-weighted velocity maps in the same velocity range. Bottom row: Position-velocity diagrams in the velocity range (+10, +20)  km s for a slice with a P.A. of 60°(see Fig. 6).

In our analysis we have adopted a conservative estimate of 2.5 kpc based on a series of reasonable assumptions from the scarce data available. We find that larger distances, in the range of 10 kpc or more, are still possible –according to Gaia parallaxes and admitting feasible departures from galactic rotation– but somehow less likely, as they lead to scenarios difficult to reconcile with the observations.

Size and age of the structure increase with while mass scales with . Therefore, at 10 kpc the structure would be years old for  km s, which places its origin into the main sequence stage. The total mass of molecular gas in the torus would then be M. Groh2014 analyzed the evolution of a non-rotating 60 M star in terms of mass-loss, finding that the mass lost during the 3 Myr prior to the LBV phase is less than 10 M in total. Considering that molecular material represents just a fraction of the mass ejected, MGE 042.0787+00.5084 would need to rival the luminosity of Car and should be an extremely efficient molecule producer in order to form such a massive structure. In addition, with the contribution of the inner nebula, the circumstellar material around this star would probably exceed the mass of any other known LBV nebula. On the other hand, the density of the gas is proportional to ; as a consequence, at 10 kpc the gas would be diffuse, with an average of cm. Under these conditions, LTE modelling would require extremely high kinetic temperatures far from the star (above 400 K at 0.75 pc) to reproduce the observed CO line intensities.

These are only indirect arguments that do not allow us to rule out far distances but stress the need for an accurate distance determination, which would permit to trace back the evolution and mass-loss record of this source.

Parameter Meaning Value Units
d distance 2.5 kpc
inner radius 0.18 pc
p, q density steepness 2
temperature steepness 1.5
half-opening angle 15 deg
i inclination 15–45 deg
temperature at 10–80 K
density at 0.8–4 cm
expansion velocity 1–5  km s
Table 4: LIME model parametrization. Free parameters are expressed as a range.
Parameter A (i = 15°) B (i = 30°) C (i = 45°)
[cm] 800 800 800
[K] 70 70 70
[ km s] 3 1 1
Table 5: Best-fitting parameters for each LIME model.

5 Conclusions

We report the first detection of CO gas associated with the candidate LBV MGE 042.0787+00.5084 forming a circumstellar structure that encloses the infrared nebula. This detection was obtained by means of 1.3 and 3 mm observations in a region of in the field of the source carried out with the IRAM 30-m telescope. After a detailed study of the morphology and kinematics of the structure, we estimate the global physical parameters of the gas by LTE and non-LTE modelling. Finally we provide an interpretation for the formation of the molecular gas consistent with two scenarios in the LBV phase. The key findings in this paper are:

  1. We found significant emission of CO and CO and in the velocity range (+13.4, +17.4) km s, tracing a structure with a clear circumstellar distribution. We interpreted this structure as an expanding torus according to the observed morphology and kinematics.

  2. We studied the physical conditions in two separate regions of the structure, namely the northwestern clump and the main torus. We constrained density, column density and kinetic temperature under LTE and non-LTE assumptions, finding that the clump is denser than the main torus by a factor of 1.5–2 regardless of the temperature, with H volume densities of a few cm. On the other hand, only kinetic temperatures in the range 10–70 K for the torus and 10–100 K for the clump are consistent with the observed line intensities and ratios.

  3. The measured [CO/CO] ratio was rather low, in the range from 5 to 7 in most of the structure. This result was particularly intriguing as no CO was detected, except in the southern cloud. We proposed a CO overproduction as an explanation for such an abnormal isotopic ratio, discussing the implications of this enhancement as the result of stellar evolution or selective photodissociation.

  4. The formation of a toroidal structure around a LBV star can be explained by two scenarios: (1) molecular ejecta from the star via non-spherical mass-loss, or (2) compression of the parent molecular cloud due to the effects of the strong and steady stellar wind. The determination of the total amount of molecular material could pinpoint the most likely scenario, but the distance to the source is poorly constrained. By adopting = 2.5 kpc, we derived a mass of molecular gas of 0.6 M, consistent with an average mass-loss rate of 0.8–1.2  M yrduring the LBV phase.

  5. We developed a simple kinematic model with the aim of reproducing the morphology and velocity fields of the structure. We studied the influence of density, temperature, inclination and expansion velocity. The model with the lowest inclination (=15°) and a rather low expansion velocity (3 km s) provided the best fit, being able to accurately reproduce the observed velocity gradient in most of the structure, as proven by first-order moments and PV diagrams. This supports the interpretation of the structure as an expanding torus of CO with a nearly face-on viewing geometry.

The characterization of the molecular material associated with MGE 042.0787+00.5084 constitutes an excellent example of how studies of this kind provide useful information on the mass-loss processes during the LBV phase. This source has shown a great potential to learn about the interplay between evolved high mass stars and the ISM, deserving follow-up studies to better constrain its physical parameters and evolutionary record. This can be achieved by means of interferometric observations and near-infrared high-resolution spectroscopy. Most importantly, this successful detection proves that the existence of molecular material in the vicinity of evolved massive stars is not an exceptional finding but likely a common feature, therefore laying the foundations for future studies involving similar sources and less abundant molecules.


This work is based on observations carried out under project number 043-17 with the IRAM 30-m telescope. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). We acknowledge the IRAM staff from Pico Veleta for the help provided during the observations. J.R.R. acknowledges the support from project ESP2015-65597-C4-1-R (MINECO/FEDER).


Appendix A Computation of kinematic distances

Assuming MGE 042.0787+00.5084 is in circular rotation around the Galactic Centre, its angular velocity is given by


where and denote the galactic longitude and latitude of the source and and are the galactocentric distance and angular velocity of the Sun respectively. Considering the fit for the rotation curve of the Milky Way provided by Fich1989, the galactocentric radius of an object with angular velocity is


This galactocentric radius relates to the distance through the expression


which is a regular second degree equation and thus admits two possible solutions. Taking the canonical values for the local galactic constants kpc and  km s, we obtain a galactocentric radius of R = 7.63 kpc for MGE 042.0787+00.5084, which corresponds with kpc and kpc.

Appendix B LIME model

The model is an ideal, axisymmetric torus defined in a cylindrical coordinate system . Geometry of the model is defined by three parameters: distance, position angle and inclination angle of the structure with respect to the line of sight.

We used ortho- and para- H as the main collision partners for photon propagation, with an ortho-para ratio of 3, and grains covered by thin ice mantles for dust opacity (Ossenkopf1994). A CO relative abundance of was used. No magnetic field was included in the model.

The model assumes gas and dust are thermally coupled, following the same torus-like density distribution, such that


where is a scaling radius equal to the inner radius of the torus, is a reference density value, and are dimensionless parameters that determine the radial and angular density gradients respectively, and is the half-opening angle of the torus. To model the inner cavity we added a depletion factor that weights the reference density when . Accordingly, points with (i.e. outside the torus) are similarly depleted.

We modelled the temperature profile with a simple radial power law controlled by the dimensionless parameter ,


for , remaining constant at reference temperature within the cavity, which is consistent with the absence of material to absorb radiation.

Regarding the torus kinematics, we adopted a velocity field consistent with a structure expanding at a constant rate, as discussed in Sect. 4.1. Consequently the velocity vector is completely radial in the model space, with .

To measure the goodness of the fitting we used a channel-wise criteria in the velocity range of interest (+13.4, +17.4  km s), defined as


where and are the main-beam temperature per pixel for model and data respectively, is the data rms scaling factor and denotes the velocity channel. By doing so we account for the morphology and velocity extent of the emission simultaneously.

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description