Feedback in W49A diagnosed with radio recombination lines and models

TL;DR

This study uses radio recombination lines to analyze the ionized gas dynamics in W49A, revealing a feedback-driven shell structure that influences star formation, modeled with WARPFIELD and indicating re-collapsing feedback processes.

Contribution

First application of RRL imaging to W49A combined with feedback modeling using WARPFIELD, demonstrating a feedback-driven, re-collapsing shell structure influencing star formation.

Findings

Ionized gas forms a shell around W49A's central OB cluster.

The shell structure is consistent with feedback-driven re-collapsing models.

Feedback processes may have triggered subsequent star formation.

Abstract

We present images of radio recombination lines (RRLs) at wavelengths around 17 cm from the star-forming region W49A to determine the kinematics of ionized gas in the THOR survey (The HI/OH/Recombination line survey of the inner Milky Way) at an angular resolution of 16.8"x13.8". The distribution of ionized gas appears to be affected by feedback processes from the star clusters in W49A. The velocity structure of the RRLs shows a complex behavior with respect to the molecular gas. We find a shell-like distribution of ionized gas as traced by RRL emission surrounding the central cluster of OB stars in W49A. We describe the evolution of the shell with the recent feedback model code WARPFIELD that includes the important physical processes and has previously been applied to the 30 Doradus region in the Large Magellanic Cloud. The cloud structure and dynamics of W49A are in agreement with a feedback-driven shell that is re-collapsing. The shell may have triggered star formation in other parts of W49A. We suggest that W49A is a potential candidate for star formation regulated by feedback-driven and re-collapsing shells.

Figures

Tables

Table 1. Table 1: Cloud properties of W49A.

Quantity	Assumed values	Comments	Reference
Cloud radius	$R\approx 6\mathrm{pc}$	Approximate distance from Cluster 1 to W49SW and W49S.	Peng et al. (2010)
Molecular gas mass (W49A; r$\lesssim 6$ pc)	$M_{\mathrm{gas}}\approx 2_{-1}^{+2}\times {10}^5{M}_{\odot }$	Assumed uncertainties of a factor of 2; used for modeling.	Galván-Madrid et al. (2013), Urquhart et al. (2018)
Molecular gas mass (W49 molecular cloud; r$\lesssim 60$ pc)	$M_{\mathrm{gas}}\approx 1\times {10}^6{M}_{\odot }$		Galván-Madrid et al. (2013)
Stellar cluster mass (all infrared sub-clusters)	$M_{\ast }\approx 5-7\times {10}^4{M}_{\odot }$		Homeier & Alves (2005)
Stellar cluster mass (Cluster 1)	$M_{\ast }\approx 1\times {10}^4{M}_{\odot }$	Used for modeling, assuming that mass is possibly up to a factor of two higher.	Homeier & Alves (2005)
Derived molecular gas density (W49A; $r\lesssim 6\mathrm{pc}$)	$n\approx 4_{-2}^{+4}\times {10}^3{\mathrm{cm}}^{-3}$	Assuming homogeneity and spherical symmetry.

Table 2. Table 2: Datasets used in this work.

Data	Frequency	Angular res.	Spectral res.	Noise	Survey	Reference
	[GHz]		[km $s^{-1}$]
Stacked RRLs (H151$\alpha$-H156$\alpha$; H158$\alpha$)	1.780a	$16\stackrel{\prime \prime }{.}8\times 13\stackrel{\prime \prime }{.}8$	5	$\sigma \left(F\right)\sim 1.5$ mJy beam-1	THORb	This work
13CO(3-2)	330.587	$15″$	0.5	$\sigma \left(T_A\right)\sim 0.6K$	CHIMPSc	Rigby et al. (2016)
C18O(3-2)	329.331	$15″$	0.5	$\sigma \left(T_A\right)\sim 0.8K$	CHIMPSc	Rigby et al. (2016)
13CO(1-0)	110.201	$46″$	0.2	$\sigma \left(T_A\right)\sim 0.1K$	GRSd	Jackson et al. (2006)

Equations

$$ϵ=M_{\ast }/M_{\mathrm{cloud}},$$

Full text

1 Introduction

Radiative and mechanical feedback from UV radiation and winds from OB stars and supernova explosions alter the potential of molecular gas to fragment and form future generations of stars, but their precise role in star formation is yet to be clarified. Recent investigations of 30 Doradus in the Large Magellanic Cloud (LMC) indicate that feedback from the older stellar population in its main star cluster NGC 2070 has been unable to destroy its parent molecular cloud (Rahner et al., 2018a). Rather, feedback-driven shells of ionized and molecular gas from these older cluster members may have re-collapsed due to the gravitational attraction of the star cluster and the self-gravity of the shell, forming a second generation of stars. Only with this second generation can the combined stellar feedback become strong enough to disperse the parent molecular cloud.

While star formation induced by re-collapsing shells may not change the total number of stars formed in a giant molecular cloud (GMC), such a process would naturally imply that star formation in at least some GMCs may occur episodically. As pointed out by Rahner et al. (2017), re-collapse of a feedback-driven shell can occur only in the most massive and densest GMCs. It is therefore an open question whether this also takes place in GMCs in the Milky Way.

In our search for a 30 Doradus analog in the Milky Way, we focus this study on W49A. It is one of the most massive and luminous young star-forming regions in the Galaxy, with one of the highest luminosity-to-molecular-mass ratios, indicating a high star formation efficiency (it is responsible for 12% of all luminosity in the Galaxy; Urquhart et al., 2018). While regions like W49A lie within the statistical distribution of luminosity-to-molecular-mass ratios of the Galaxy-wide sample of star-forming regions, the nature of the physical processes that cause this extreme star formation is still not fully understood. This makes it an ideal region to search for re-collapse of feedback-driven shells.

W49 was first discovered as a radio source by Westerhout (1958) and is located at a distance of 11.1 kpc from the Sun (Zhang et al., 2013) in the Perseus spiral arm (e.g., Moore et al., 2012), at a Galactocentric radius similar to that of the solar system. It is associated with a GMC of mass $\sim 10^{6}\,{\rm M_{\odot}}$ (Galván-Madrid et al., 2013), which contains the massive H ii region complex W49A (Mezger et al., 1967) with a molecular gas mass of $\sim 2\times 10^{5}\,{\rm M_{\odot}}$ (Urquhart et al., 2018; Galván-Madrid et al., 2013). An overview of W49A in Galactic coordinates is shown in Fig. 1, with 8 $\mu$m emission (GLIMPSE; Benjamin et al., 2003; Churchwell et al., 2009) as a tracer of ongoing star formation (e.g., Stock et al. 2014), and 870 $\mu$m emission (from ATLASGAL; Schuller et al., 2009) indicating the presence of cold and dense gas – the sites of future star formation. The main physical parameters are summarized in Table 1.

The W49A cloud harbors multiple ultracompact (UC) H ii regions (e.g., Dreher et al., 1984; De Pree et al., 1997), with a central condensation in a ring-like structure (Welch et al., 1987). The region has the highest concentration of UC H ii regions in the Galactic disk (18 compact and UC H ii regions; Urquhart et al., 2013). Four star clusters are detected by infrared imaging (Alves & Homeier, 2003), amounting to a total stellar mass of $5-7\times 10^{4}\,{\rm M_{\odot}}$ (Homeier & Alves, 2005). The so-called ‘Cluster 1’ (Homeier & Alves, 2005) in the central part of W49A (r$\lesssim 2.5\,{\rm pc}$) contains about 50 O stars with $M_{*}\geq 20\,{\rm M_{\odot}}$, and a total of $\sim$270 O stars are found in and around W49A (Homeier & Alves, 2005). Several of these candidate cluster members have masses between $20$ and $250\,{\rm M_{\odot}}$, as determined from near-infrared spectroscopy and photometry (Wu et al., 2016, Fig. 8).

The molecular gas of W49A shows a complex velocity structure, with blue- and red-shifted components with respect to the systemic velocity ($\varv_{\rm LSR}=8.6$ km ${\rm s}^{-1}$; Quireza et al. 2006a): There are two main velocity components at 4 km ${\rm s}^{-1}$ and 12 km ${\rm s}^{-1}$(or at $-$4.6 km ${\rm s}^{-1}$  and 3.4 km ${\rm s}^{-1}$ with respect to the systemic velocity) observed in many gas tracers on scales ranging from the entire W49A region (Mufson & Liszt, 1977; Miyawaki et al., 1986; Simon et al., 2001; Miyawaki et al., 2009) to the innermost parts, towards W49N in CS (e.g., Serabyn et al., 1993). While CO transitions may be optically thick, the double-peaked structure appears to reflect the true dynamics of the molecular cloud, as it is persistent in optically thin, high-density tracers such as H13CO+(1-0) and CS(2-1) (Fig. 3 in Galván-Madrid et al. 2013; also Miyawaki et al. 2009). The velocity components were attributed to background and foreground clouds, respectively, by Serabyn et al. (1993), as indicated by ${\rm H_{2}CO}$ absorption (Goss & Tilanus, 1985). This spatial location would mean that the two clouds are moving towards each other, either as the collision of two clouds, or the inside-out collapse of one cloud (Serabyn et al., 1993; Welch et al., 1987). Alternatively, the velocity components could be due to expanding motions (Peng et al., 2010). Other studies have suggested that the complexity of the ${\rm H_{2}CO}$ absorption, as well as of CS and HCO+ emission in this region, is connected to a complex arrangement of H ii regions surrounded by dense and diffuse molecular gas, outflows, and infall motions (Dickel & Goss, 1990; Williams et al., 2004).

The number of ionizing photons from the stellar population in W49A (Alves & Homeier, 2003), with the largest contribution being from the most massive of these stars (Wu et al., 2016), is sufficient, or may exceed that required, to produce the observed ionized gas within W49A, as determined from low- (Kennicutt, 1984) and high-resolution continuum observations (e.g., De Pree et al., 1997). The fraction of ionized to molecular gas mass is low ($\sim 1\%$; Galván-Madrid et al. 2013), indicating that the molecular gas in W49A is not yet penetrated by this radiation. Estimates of the radiation pressure on dust particles indicate that radiative feedback alone from the star cluster is not yet strong enough to disperse W49A (Galván-Madrid et al., 2013; Reissl et al., 2018), while models with other feedback mechanisms (thermal pressure from the H ii region and shocked winds, protostellar jets, radiation pressure; but not wind-momentum) find feedback to be approximately strong enough to disrupt the cloud (Murray et al., 2010). It was indeed recently noted that comprehensive modeling that includes all sources of feedback is needed to reliably predict the impact of radiative feedback on cloud evolution (Rahner et al., 2017).

This poses the question of how, if at all, feedback from Cluster 1 affected star formation in W49A. The H ii regions and O/B associations within W49A have different ages, with some having already dispersed some of the molecular gas around them (the O/B stars discovered in Alves & Homeier 2003), and others still deeply embedded (e.g., the UC H ii regions in the Welch ring; e.g., Dreher et al. 1984; De Pree et al. 1997). Both single O stars and UC H ii regions appear to be spread across W49A in small groups. While many studies focussed on the formation of the UC H ii regions in the Welch-ring/W49N dust clumps (e.g., Welch et al., 1987; Serabyn et al., 1993; Williams et al., 2004), two main scenarios have been invoked to explain the star formation in W49A as a whole. One interpretation is causally unrelated, sub-clustered star formation across W49A (Alves & Homeier, 2003) along filamentary inflows of molecular gas from a larger reservoir of gas (Galván-Madrid et al., 2013). Alternatively, star formation and the dynamics of the region are causally connected to the feedback of the central star cluster. The central cluster of O stars drives expanding shells of molecular gas (Peng et al., 2010), which may have triggered star formation in the Welch ring (Alves & Homeier, 2003). Feedback may also be responsible for the formation of W49S and W49SW, if they are interpreted as ejecta from the central star cluster. While previous interpretations of sub-clustered star formation may be a suitable description of the star formation activity in W49A, the dynamics show imprints of feedback (Peng et al., 2010).

In this work, we investigate the observational signatures of stellar feedback using emission from radio recombination lines (RRLs). These lines have been used in numerous works to determine the physical conditions of ionized gas in H ii regions (e.g., Hjellming & Davies 1970, Shaver 1980; see also review by Roelfsema & Goss 1992), in particular also for W49A (e.g., van Gorkom et al., 1980; Balser et al., 1999). They have also been mapped for different star-forming regions in the Milky Way (e.g., Pankonin et al., 1979; Lang et al., 2001; Balser et al., 2001). In this work, we use high-$n$ RRLs (1.6-1.9 GHz) to trace the kinematics of ionized gas in W49A. The continuum emission at these frequencies may be optically thick, especially towards the location of UC H ii regions such as the Welch-Ring. As the continuum optical depth affects the line emission, this caveat is addressed and the results are compared to previous RRL studies at higher frequencies where available. We use the RRL data from the THOR survey (The H i/OH/Recombination line survey of the inner Milky Way; Beuther et al. 2016) as constraints for models of the past and future evolution of a feedback-driven shell in a cloud with average physical properties corresponding to those of W49A.

After presenting the observations in Sect. 2, the RRL emission is analyzed and compared to molecular gas emission in Sect. 3 to characterize the morphology and velocity structure of different gas phases in W49A. We identify shell-like RRL emission at the interface between hot ionized gas and neutral gas. With the radius and velocity of the shell, we aim to constrain its evolution with the 1D code, warpfield111https://bitbucket.org/drahner/warpfield/ (Rahner et al., 2017) in Sect. 4. While acknowledging that a 1D code does not capture all the complex dynamics of W49A, it includes a refined treatment of the interplay of the relevant types of stellar feedback and allows us to easily probe a large parameter space in density and star formation efficiency. This provides general conclusions on the evolution of this feedback-driven shell in W49A. The results are discussed in the context of the previous interpretations of star formation in W49A, and compared to the feedback-driven star formation history of 30 Doradus in Section 5. In Sect. 6, we summarize our main findings.

2 Observations and Methods

2.1 Radio recombination lines at 1.6–1.9 GHz

The RRLs allow us to study the kinematics and structure of the ionized gas in W49A. We choose RRLs between orders of n=151 and n=158 (1.6-1.9 GHz) from the THOR survey (The H i/OH/Recombination line survey of the inner Milky Way; Beuther et al. 2016) to extract spatially resolved kinematic information of ionized gas. As these lines may in principle be affected by high optical depth in the continuum, we discuss the optical depth in W49A in Sect. 3.1 and compare our data in Sect. 3.3.1 to kinematic information obtained from RRLs at higher frequencies in the literature.

The H151$\alpha$-H156$\alpha$ and H158$\alpha$ RRLs were observed towards W49A within the THOR survey. Observations of 5$-$6 minutes per pointing were conducted with the VLA C-configuration in L-band. Each line was observed with a bandwidth of 2 MHz and a spectral resolution of 15.63 kHz. This corresponds to a velocity resolution of 2.5 km ${\rm s}^{-1}$ and 2.8 km ${\rm s}^{-1}$ for the highest and lowest frequency transitions, respectively (the H151$\alpha$ line at 1.891 GHz and the H158$\alpha$ line at 1.651 GHz; see also Table 2 in Beuther et al. 2016). Using CASA222http://casa.nrao.edu; version 4.2.2, the absolute flux scale and the bandpass were calibrated on the quasar 3C386 and the complex gain with the quasar J1925+2106.

With a significant detection of all RRLs, the data were gridded to a spectral resolution of 5 km ${\rm s}^{-1}$ as opposed to 10 km ${\rm s}^{-1}$ in the first data release of the THOR survey. All observations were continuum-subtracted, imaged, and combined into mosaics of $3.75\degree\times 2.5\degree$, and deconvolved as described in Beuther et al. (2016). We only consider the seven lines in the frequency interval between 1.6 and 1.9 GHz to maintain higher spatial resolution (while higher-order transitions from lower frequencies were included by Beuther et al. 2016). The angular resolution of the observations of the individual transitions is between $12\aas@@fstack{\prime\prime}0\times 14\aas@@fstack{\prime\prime}6$ (H151$\alpha$) and $16\aas@@fstack{\prime\prime}8\times 13\aas@@fstack{\prime\prime}8$ (H158$\alpha$). To improve the signal-to-noise ratio, the different Hn$\alpha$ lines were stacked in velocity after regridding each individual line to a common angular resolution (using CASA task imsmooth). Our image of the stacked continuum-subtracted emission of the RRLs has an angular resolution of $16\aas@@fstack{\prime\prime}8\times 13\aas@@fstack{\prime\prime}8$. The RMS noise of the stacked image in a line-free channel is 1.5 mJy beam*-1*. The characteristics of the RRL data are summarized in Table 3.

This study focusses on the continuum-subtracted RRL emission associated with W49A ($\varv_{\rm LSR}=8.6$ km ${\rm s}^{-1}$; Quireza et al. 2006a). A likely unassociated component at 60 km ${\rm s}^{-1}$, which has been found in other RRL observations (e.g., Liu et al. 2013; Liu et al., in prep.), is not significantly detected in our dataset and not discussed further here (faint emission (¡3 $\sigma$) at 60 km ${\rm s}^{-1}$ is seen in some of the spectra in Fig. 3). The integrated emission between $-$40 and 60 km ${\rm s}^{-1}$ (moment 0) is presented in Fig. 2. Individual spectra are shown in Fig. 3, which were extracted towards three positions indicated by the labels “Pos. 1”, “Pos. 2” and “Pos. 3” in Fig. 2 and in the upper-left panel of Fig. 4. Images of all channels between $-$10 and 25 km ${\rm s}^{-1}$ are shown in Fig. 4.

2.2 Archival CO observations

To investigate the interplay between ionized and molecular gas, we use 13CO(3-2) and C18O(3-2) emission (CHIMPS444The 13CO/C18O (J = 3$\rightarrow$2) Heterodyne Inner Milky Way Plane Survey. survey; Rigby et al. 2016) to trace the morphology and kinematics of molecular gas in W49A. Although 13CO may be optically thick at the high column densities found towards W49A (Galván-Madrid et al., 2013), lines from high-density tracers like CS and H13CO+, as well as from the high-column-density tracer C18O, show similar shapes in the same study. The ratio of main-beam brightness temperatures of the two CO isotopologs remains close to constant at $\sim 0.11$ across the peaks and troughs of the spectra, roughly in agreement with a ratio of $0.14$ from Wilson & Rood (1994) in the solar neighborhood. This is expected as W49A has a similar Galactocentric radius to that of the Sun. Similar conclusions have been obtained by Miyawaki et al. (2009) based on investigations of 13CO (J=1-0) and C18O (J=1-0) at 17″ resolution. Hence, we assume 13CO emission to be optically thin and use it with the C18O(3-2) emission for the kinematic analysis of the molecular gas in W49A.

Table 3 contains a detailed description of the CO data used, which includes archival 13CO(1-0) data from the GRS555Galactic Ring Survey. survey (Jackson et al., 2006) at lower angular resolution. Spectra of all tracers are shown in Fig. 3 towards selected positions in W49A. Figure 4 shows emission contours of C18O(3-2) after smoothing to a spectral resolution of 5 km ${\rm s}^{-1}$.

2.3 Gaussian line-fitting of the RRL data

To extract the kinematic properties of the RRL emission, maps of velocity centroids and line widths are derived for each tracer. The RRLs typically have a full width at half maximum (FWHM) of 25–45 km ${\rm s}^{-1}$ in W49A. The RRL data are well described by a single Gaussian. We fit the stacked RRLs at each pixel with a Gaussian profile with the CASA task specfit to obtain maps of the peak velocity and FWHM (Fig. 5). Despite the large channel width of 5 km ${\rm s}^{-1}$, the nominal uncertainty on the center of the Gaussian is $\lesssim 2.5$ km ${\rm s}^{-1}$ in all cases, and $<1$ km ${\rm s}^{-1}$ towards the RRL emission peak.

3 Results

3.1 Spatial distribution of the RRL emission in comparison to other gas tracers

The emission in the stacked Hn$\alpha$ data between 1.6 and 1.9 GHz is spread over the entire W49A star-forming complex. The spectral indices of the continuum emission between 1 and 2 GHz ($\alpha$) fall between $\sim 0$ and $\sim 1$, indicating thermal, partially optically thick emission (Wang et al., 2018, figure 10). Towards W49N, the spectral index approaches values of $\alpha\sim 0.8$, which can indicate ionized outflows (e.g., De Pree et al., 1997, and references therein). It is similarly extended in W49A as the 8 $\mu$m continuum (Fig. 1) and the 870 $\mu$m continuum emission (Fig. 2).

The emission peaks of the different tracers are slightly different. The integrated 1.6–1.9-GHz RRLs peak offset from the 870 $\mu$m emission (Fig. 2). The offset between the peaks is $\sim$25″ or $\sim$1.5 pc.

The peak in 870 $\mu$m emission corresponds to W49N, which harbors the Welch-Ring of UC H ii regions (Welch et al., 1987), as indicated in Fig. 1 by the 5-GHz emission. The fact that the 1.6–1.9-GHz RRL emission (and the continuum emission at the same frequency) peaks slightly offset is likely due to optically thick continuum from the high-density ionized gas in the very compact UC H ii regions at the location of the Welch-ring that lowers the flux at our longer wavelengths. A similar effect could be causing the relatively weak RRL emission towards W49S and W49SW, which also host UC H ii regions (e.g., Dreher et al., 1984; De Pree et al., 1997).

The peak of the 1.6–1.9-GHz RRL emission occurs towards the more extended radio sources “L” and “M” at the edge of the Welch ring (e.g., Dreher et al., 1984, see Fig. 2 and inlay to Fig. 1). Since it is also offset from the dust emission peak, this H ii region may have already cleared parts of its cocoon of gas and dust.

3.2 Shell-like distribution of RRL and CO emission around Cluster 1

At the center of W49A, the emission of stacked RRLs appears to have a shell-like morphology. The channel maps (Fig. 4) show that towards the position of Cluster 1 (Pos. 2 in Fig. 3), the RRL emission is suppressed. North and south of Cluster 1 (with respect to Galactic coordinates), the RRL emission increases, in parts resembling a ring-like structure. This emission is potentially indicative of an ionized bubble around Cluster 1 with a radius of $\sim 2\pm 1$ pc, as indicated by the inner white ring in Fig. 4.

The radius of the bubble is in agreement with the shells found in 4.5 and 8.0 $\mu$m emission by Peng et al. (2010, 2-3 pc diameter); the RRL emission here agrees better with their shell 1, though this cannot be said with certainty due to the lower angular resolution of the RRL image. Extended emission at this position is seen in the broadband 1.6 GHz continuum image (Fig. 1) as well as in radio continuum at higher frequencies (e.g., De Pree et al., 1997). The spectral index between 1 and 2 GHz in this particular region of W49A lies between $\alpha=-0.1$ and 0.3, and increases to $\alpha=0.8$ towards the Welch-Ring (Wang et al., 2018, figure 10).

Confirmation of this structure is seen in the C18O(3-2) data (Fig. 4), but it is less pronounced in 13CO(3-2) and (1-0) emission (Fig. 9). In Fig. 4 at 5 km ${\rm s}^{-1}$, the C18O emission resembles an arc around Cluster 1. Channel maps at 10 km ${\rm s}^{-1}$ and 15 km ${\rm s}^{-1}$show elongated emission edges around Cluster 1, together with a lack of emission towards the center of the bubble (see also spectrum of Pos. 2 in Fig. 3). To a lesser degree, this can also be seen in 13CO(3-2) and (1-0) emission in Fig. 9. At 10 km ${\rm s}^{-1}$ and 15 km ${\rm s}^{-1}$ the emission is lower towards Cluster 1 than to the north or south of it.

We do not observe a full ring in RRL or CO emission. Parts of the RRL emission may be due to outflows from the surrounding UC H ii regions. RRL emission appears towards the north, south, and west of Cluster 1 only (in the Galactic coordinate frame). The CO observations show similar structure. Together with arcs in 8 $\mu$m emission reported east of Cluster 1 by Peng et al. (2010, see also Fig. 1), we attribute the RRL emission to a ring-like structure indicative of a bubble around Cluster 1.

3.3 Kinematics of RRLs and 13CO

3.3.1 Velocity distribution

The 13CO(3-2) emission shows two components at similar strength towards Cluster 1 (Pos. 2). Blue-shifted with respect to the systemic velocity of W49A ($\varv_{\rm LSR}=8.6$ km ${\rm s}^{-1}$; Quireza et al. 2006a) lies a broad component with a peak at $\sim$4 km ${\rm s}^{-1}$ (or multiple blended narrow components between 1$-$7 km ${\rm s}^{-1}$). The other component is a narrow, red-shifted emission peak at $\sim$13 km ${\rm s}^{-1}$. The RRL emission at this position is located between the CO components, although the peak occurs closer towards the red-shifted CO component, with a fitted peak of 13.6$\pm$2.7 km ${\rm s}^{-1}$, which has high uncertainties due to the weak RRL emission at Pos. 2.

Across W49A, the stacked RRL emission shows variations in peak velocity from the north to the south. The peak velocity distribution of the RRL data (Fig. 5) indicates a peak velocity at the latitude of W49N (Pos. 1 and higher) of 6.2$\pm$0.7 km ${\rm s}^{-1}$. At lower latitudes towards Pos. 3, velocities of 12.6$\pm$1.2 km ${\rm s}^{-1}$ are observed. While RRL emission at low frequencies of $1-2$ GHz may in principle be affected by optically thick continuum emission, we find that for large parts of the bubble, the spectral index is between $-$0.1 and 0.3, which is only slightly higher than what would be expected for optically thin emission. The velocities of the H92$\alpha$ RRL from De Pree et al. (1997) are comparable to the velocities of the higher-order RRLs shown here in the case of the sources “JJ” and Pos. 3 (14.3$\pm$0.3 km ${\rm s}^{-1}$ and 12.6$\pm$1.2 km ${\rm s}^{-1}$, respectively). The velocity is lower for sources “M” and “L” (3.8$\pm$0.4 km ${\rm s}^{-1}$ and 1.8$\pm$0.13 km ${\rm s}^{-1}$) as compared to the high-order RRL emission in Pos. 1 (6.2$\pm$0.7 km ${\rm s}^{-1}$), which, however, may also be influenced by source “J”, for which the H92$\alpha$ RRL has a velocity of 9.3$\pm$0.5 km ${\rm s}^{-1}$. The closeness of the RRL velocities from this work to those of RRL observations at higher frequencies indicates that potentially high optical depth in the continuum at 1-2 GHz is not the governing factor determining the appearance of the kinematics.

The 13CO(3-2) emission extends between 0 and 20 km ${\rm s}^{-1}$ in W49A. Towards Pos. 1, a red-shifted peak dominates the emission, while towards Pos. 3 it is a blue-shifted peak that dominates. In both positions, the 13CO peak velocity appears anti-correlated with the peak velocity observed in the RRL emission at the same positions. At the RRL peak velocities, the 13CO emission seems to be suppressed. The C18O(3-2) emission shows similar profiles and confirms the location of the 13CO peaks. This confirms the assumption that the 13CO emission profiles are not affected by optical depth effects. The lower-angular-resolution GRS 13CO(1-0) emission shows two smooth peaks towards all three positions, which are mainly the result of averaging over variations within the 46″ GRS beam.

To summarize, towards the position of Cluster 1, the RRL emission is located in between, or aligned with one of, the two 13CO emission components. The two 13CO components are separated by $\sim 9$ km ${\rm s}^{-1}$. If these are attributed to fore- and background parts of the shell bubble (Peng et al., 2010), the expansion or collapse velocity of this shell would be approximately between 5 and 10 km ${\rm s}^{-1}$, depending on whether both parts of the shell are moving, or one of them is at rest with respect to the star cluster. An alternative explanation for the kinematics of the region is, for example, ionized outflows from the surrounding UC H ii regions. To the north and south of the center of the shell bubble (Peng et al., 2010), both the RRL and CO emission show one main component, with the peaks at each position anti-correlated with respect to the rest velocity of W49A ($\varv_{\rm LSR}=8.6$ km ${\rm s}^{-1}$; Quireza et al. 2006a). This indicates that, at least in parts of the ring, the emission may be influenced by dynamics of individual places of star formation in the region. For the following modelling (Sect. 4) we concentrate on the large-scale dynamics. For the ionized bubble, we assume a spherical structure with a shell radius of 1$-$3 pc from the RRL emission data and a relative velocity of the shell components of 5$-$10 km ${\rm s}^{-1}$.

3.3.2 Distribution of the FWHM of the RRL emission

It is intrinsically difficult to see expansion or infall signatures in the RRL line profiles, due to their large line widths of up to 50 km ${\rm s}^{-1}$. The large line width itself, especially towards Cluster 1 (note the increase of line widths up to $45$ km ${\rm s}^{-1}$ in Fig. 5), however, may point to unresolved dynamical motions. Contributions to the line width are thermal broadening, pressure broadening, and nonthermal motions (dynamical broadening). The contribution of pressure broadening may be significant for high-order RRLs, while spatially unresolved dynamical motions can also contribute to the line width, especially in very active star-forming regions.

We evaluate the relative importance of pressure and dynamical broadening for two positions. Assuming an electron temperature for W49A of $T_{e}\sim 8300\,{\rm K}$ (Quireza et al., 2006b), the thermal line width is approximately 20 km ${\rm s}^{-1}$. An estimation of the pressure broadening requires observations of lower-order RRLs, as pressure broadening depends on the order of the energy level of the transition (n) to the seventh power (e.g., Keto et al., 1995). Observations of millimeter(mm)-H$\alpha$ lines (n=$39-42$) are available at an angular resolution of $\sim 20\arcsec$ towards the ATLASGAL dust emission peaks in W49A (Kim et al., 2017). For W49N (source AGAL043.166+00.011), a line width of 37.0$\pm$0.2 km ${\rm s}^{-1}$ is reported for the mm-RRLs, while the centimeter(cm)-RRLs presented here have a line width of 31.4$\pm$3.3 km ${\rm s}^{-1}$. Since the widths are larger or comparable at most, we conclude that dynamical broadening dominates the emission, which is expected due to the large number of UC H ii regions at this location. For this particular position, optical depth effects may also influence the kinematics of the cm-RRLs. Towards Pos. 3 (source AGAL043.178-00.011), the line width of the cm-RRLs is larger than the mm-RRLs (33.5$\pm$4.0 km ${\rm s}^{-1}$ and 29.7$\pm$0.5 km ${\rm s}^{-1}$, respectively). Using equations 4 and 5 from Nguyen-Luong et al. (2017) and a thermal line width of 20 km ${\rm s}^{-1}$, we derive line widths of dynamical and pressure broadening of 23 km ${\rm s}^{-1}$ and 6 km ${\rm s}^{-1}$, respectively. Again, this indicates that dynamical broadening dominates the RRL line width over pressure broadening. While this is likely to be true for parts of W49A which contain high star formation activity (as UC H ii regions indicate in many positions), pressure broadening is certainly present and may be dominant towards some positions. Especially for large line widths in Fig. 5 (e.g., Pos. 2), observations of RRLs with lower quantum numbers at the same angular resolution as the cm-RRLs are necessary to verify the role of pressure broadening (see, e.g., Nguyen-Luong et al., 2017).

An investigation of the spectra for double-peaked structures, which would be indicative for relative motions of the front and back of a shell, is inconclusive: Figure 4 shows slightly more emission in Pos. 2 at 5 km ${\rm s}^{-1}$ and 15 km ${\rm s}^{-1}$, compared to 10 km ${\rm s}^{-1}$. The variation, however, is only at a 2-$\sigma$ level and the spectrum is not fitted well by two Gaussians. Also, two distinct Gaussian components with such a small velocity separation are not expected to be distinguishable in the RRL data from a single Gaussian due to their large intrinsic, thermal line widths.

4 Stellar feedback models for W49A

How has star formation progressed in W49A? We now use the presented data to answer important questions about W49A: What influence has the formation of high-mass stars had on the subsequent evolution of the cluster? To explore possible scenarios of the evolution of W49A, we employ models of an expanding shell around a star cluster to compare with our observations.

In the following, we compare our data to models of an expanding shell of a star cluster inside a molecular cloud to constrain the evolution of the shell around Cluster 1, which was discussed in Section 3.2. Given its complexity in terms of currently ongoing star formation and density distribution, we do not attempt here to reproduce the RRL emission in W49A directly. Rather we aim to use the radius at which the RRL emission peaks to constrain the evolutionary history of the shell. As discussed below, we use literature values of the stellar mass in Cluster 1, the molecular cloud mass, and the average density of W49A (see Table 1) as initial conditions, and constrain the evolutionary stage of the shell with its observed radius and the age of the stellar cluster.

Rahner et al. (2017) have developed feedback models (warpfield666The modelling in this paper is based on warpfield 1 (Rahner et al., 2017). The recent update of the code (warpfield 2; Rahner et al. 2018b) was tested to yield consistent results for this analysis.) describing the expansion of a hydrostatic, spherical shell due to radiative feedback from young stars, including the energy and momentum input from stellar winds and supernovae, while also accounting for the effects of the gravity of the cluster and the self-gravity of the cloud. Models with the physical prescription for the hydrostatic shell, which is used in warpfield, have been successfully compared to observations in previous works (e.g., M17; Pellegrini et al., 2007). In the models by Rahner et al. (2017), the expansion of the shell is initially governed by adiabatic expansion due to wind-shocked gas from the stellar winds of O/B stars in the massive star cluster. Once cooling becomes comparable to the wind-energy input, a second phase starts, and momentum-driven expansion from radiation and winds dominates. A potential third phase of expansion is reached as the entire molecular cloud is swept up and the system freely expands into the (low-density) ambient medium.

These models are applied here due to their advantage in efficiently probing large parameter spaces in molecular cloud mass, density, and star formation efficiency. For the computation of multiple models to remain computationally feasible, the models assume a 1D geometry (i.e., spherical symmetry of the H ii region). While this is certainly a limitation, especially for determining detailed dynamics and morphologies of H ii regions, it allows us to trace the general evolution of feedback-driven shells. The models also incorporate all relevant physical aspects of feedback, in particular for very massive and luminous regions in which the initial mass function should be fully sampled, with masses of individual stars of up to 120 ${\rm M}_{\odot}$. In the sub-region of W49A relevant for this investigation, Cluster 1, the most massive star has a mass of $M=130\pm 30\,{\rm M}_{\odot}$ (Wu et al., 2016, Fig. 8), which is within the uncertainties of the highest mass of sampled stars in the models. Therefore, since W49A is one of the most massive and luminous regions in the Milky Way, it is an ideal target for such an investigation (Urquhart et al., 2018).

4.1 Initial conditions

The input parameters for the models are listed in Table 1. Galván-Madrid et al. (2013) find a molecular gas mass of $M_{\rm gas}\approx 2\times 10^{5}\,{\rm M_{\odot}}$ for W49A within a radius of 6 pc, which we adopt for the modeling in the following. The value is in agreement with the cloud mass derived from dust emission (Sievers et al., 1991; Urquhart et al., 2018). As a side note, W49A is embedded in a larger molecular cloud complex, which has a total mass (including W49A) of $1\times 10^{6}\,{\rm M_{\odot}}$ within 60 pc (Galván-Madrid et al., 2013). Hence, free expansion of the feedback-driven shell into the low-density ISM would only occur far beyond a radius of 6 pc. Since the shells in the models considered here do not expand beyond this point during the first expansion (see Sect. 4.2), we treat the molecular cloud mass as fixed at $M_{\rm gas}\approx 2\times 10^{5}\,{\rm M_{\odot}}$ in the models. Low-density channels, in which free expansion would occur earlier, could exist in the true initial density distribution. Such channels could influence the expansion of the shell in the energy-driven phase, i.e., in the initial stages of the expansion. The momentum-driven expansion phase of the shell is not significantly affected.

Assuming a homogenous medium and spherical geometry for $r<6$ pc, the average initial molecular density of W49A is estimated as $n\sim 4\times 10^{3}\,{\rm cm^{-3}}$. To address uncertainties of at least a factor of two in the molecular gas mass, multiple models are computed to explore a parameter space in density777Mass densities chosen as input to warpfield are: $\rho\approx(8-40)\times 10^{-21}\,{\rm g\,cm^{-3}}$, in steps of $4\times 10^{-21}\,{\rm g\,cm^{-3}}$. of $n=(2-10)\times 10^{3}\,{\rm cm^{-3}}$, in steps of $1\times 10^{3}\,{\rm cm^{-3}}$. The initial radial gas density profile of the models is assumed to be flat for simplicity. Implications of radially decreasing density profiles and inhomogeneities in general are discussed in Sects. 4.2, 5.4, and 5.5. It must be noted that in reality, the cloud density might evolve as a function of time due to inflow of gas onto the cloud (e.g., Fukui et al., 2009; Seifried et al., 2017; Ibáñez-Mejía et al., 2017). This is currently not taken into account in warpfield.

With literature constraints on the mass of the stellar clusters in W49A (Alves & Homeier, 2003) and the above assumption on the molecular cloud mass, we constrain the star formation efficiency (SFE). The mass of the central O/B-star cluster, Cluster 1 (Fig. 1), is $M_{*}=1\times 10^{4}\,{\rm M_{\odot}}$ (Homeier & Alves, 2005), but may potentially be higher due to the large visual extinction in the region. Therefore, we choose SFEs of $\epsilon=0.05$ and 0.09 to account for star cluster masses of $M_{*}\simeq(1-2)\times 10^{4}\,{\rm M_{\odot}}$.

For these parameter ranges in density and SFEs the cluster does not disperse its cloud in the first 2 Myr, assuming a cloud mass of $M_{\rm cl}=2\times 10^{5}\,{\rm M_{\odot}}$ (see Table 1; Galván-Madrid et al. 2013; Urquhart et al. 2018), nor does it drive feedback shells beyond the outer radii of W49A (see Sect. 4.2). The total cloud mass is only of dynamical importance once free expansion into the ambient ISM sets in. Since this stage is not reached in any of the models in the first 2 Myr, the important parameters are the density and SFE. We choose not to vary the initial molecular cloud mass of $M_{\rm cl}=2\times 10^{5}\,{\rm M_{\odot}}$, which is the mass in W49A within $r<6\,{\rm pc}$. The shell will encounter the ambient ISM only at $r>60\,{\rm pc}$, since W49A is embedded in a larger molecular cloud complex. Therefore, the assumption of a fixed initial cloud mass is reasonable, as the phase of free expansion of the shell is likely to start only at larger radii, and after accumulating more mass than assumed here.

4.2 Models of stellar feedback indicate collapse of an initial expanding shell

The time evolution of the shell radius in each model is shown in Fig. 6, with $t=0$ at the beginning of the expansion of the shell. For simplicity, we neglect the (ultra)-compact H ii-region phase of the high-mass stars, which is assumed to occur before the shell expansion and is considerably smaller than the stellar ages ($3\times 10^{5}$ yr; Mottram et al. 2011). The shell is expected to have a similar age to the O stars, that is, 1-2 Myr (Wu et al., 2016). This time range is highlighted in Fig. 6 in gray, together with the adopted current expansion radius of 1-3 pc (in comparison, the angular resolution of the RRL images is $\sim 0.8$ pc), in agreement with both the RRL emission (see Sect. 3.2) and the 8 $\mu$m emission (Peng et al., 2010). For reference, the approximate radius of the entire W49A region is shown as the gray-shaded region C in Fig. 6.

These conditions are met by several models: For a large range of densities, shells fall within the limits of the radius after an evolution of 1-2 Myr. In all models in Fig. 6, the observed shell size agrees with the simulated shell size at two distinct times – once during collapse and again during re-expansion. For three models ($n=2.5,4.0,7.5\times 10^{3}\,{\rm cm}^{-3}$), the radius and range in time at which the observational constraints are met are indicated schematically in Fig. 6 (by a colored bar which is labeled “now”). If the double-peaked velocity profile in the 13CO data is attributed to relative motions of the back and the front of the shell, it would be expanding or collapsing at a speed of 5$-$10 km ${\rm s}^{-1}$ (see Sect. 3.3). This agrees well with the modeled velocities of the re-collapsing shells. For each model, the formation period of the O/B-stars is highlighted as being between 1 and 2 Myr before $t=$ “Now” (Wu et al., 2016). Similarly, the presence of UC H ii regions in W49A indicates a more recent episode of star formation (a few $\times 10^{5}$ yr ago or even more recent; Kawamura & Masson 1998; Wood & Churchwell 1989; Mottram et al. 2011). The formation period of the UC H ii regions is marked as a time span of $3\times 10^{5}$ yr before the current time in Fig. 6.

All three models discussed here show the same qualitative evolution: A shell, triggered by the formation of a star cluster (left panel of Fig. 7), expands to a certain radius inside the molecular cloud (middle panel of Fig. 7), at which point the self-gravity of the shell and the gravity of the star cluster start to dominate the force balance. The shell would re-collapse to form a new star cluster. Using the observed shell radius and star cluster age as constraints, we infer that the shell is either at the end of re-collapse or at the beginning of a new expansion (right panel of Fig. 7).

All models predict that Cluster 1 is not powerful enough to disperse the molecular cloud. Rather the dispersal would only happen after repeated collapse, which would trigger more star formation. Additionally, feedback-triggered star formation in the outskirts of W49A (i.e., fragmentation of the shell or shell collision with outer dense clumps) could increase the total amount of stellar feedback which might eventually become sufficient enough to disperse the cloud.

As can be seen from Models 1-3 in Fig. 6, the shell expands farther into surroundings with lower densities, and will also take a longer time to re-collapse for models with lower initial cloud densities. For a given stellar cluster mass with the same strength of feedback, this trend is expected, since the shell has piled up more material at a given radius in a denser medium. At even lower densities, feedback from the first cluster would be sufficient to destroy the cloud. However, such densities would be lower than the range of cloud densities estimated from observations. Similarly, increasing the stellar cluster mass (within the masses probed here) allows the model to push to larger radii, while the timescales for re-collapse are not strongly affected. More massive stellar clusters provide stronger feedback forces, which accelerate the shell to higher velocities and therefore cause the larger expansion. This is shown in Model 4 in Fig. 6, for $M_{*}=2\times 10^{4}\,{\rm M_{\odot}}$ at a constant density of $n=4\times 10^{3}\,{\rm cm^{-3}}$.

The radius of maximum expansion varies among the different models between the currently observed radius (gray-shaded region B in Fig. 6; Model 2) and the projected radius of W49A (gray-shaded region C; Model 3). Radial, monotonically decreasing density profiles would imply denser material in the center of W49A. The acceleration of the shell would be lower within the maximum expansion radii of the models considered here, since the material which is swept up in the shell grows faster with radius. Also, feedback due to thermal pressure in phase 1 (see above) would be lower, as dense gas cools faster. Both effects would lead to smaller radii of maximum expansion.

For all densities probed here, the models confirm that shells driven by Cluster 1 may have at least affected the Welch ring and possibly triggered star formation there (in agreement with Alves & Homeier, 2003; Peng et al., 2010). They have in common that the feedback-driven shell underwent re-collapse. On the lower end of the density range justified by the observations, shells have pushed to the outskirts of W49A in the past, and may have affected star formation all across W49A.

5 Discussion

5.1 Dynamics of re-collapse in W49A

The stellar feedback models predict a shell in contraction, which is about to induce or may have already induced a new star-formation event. From an observational perspective, the double-peaked emission of molecular gas towards Pos. 2 (Cluster 1) is consistent with the presence of a moving shell. Both peaks are $\sim$4.5 km ${\rm s}^{-1}$ offset from the systemic velocity of $\varv_{\rm LSR}=8.6$ km ${\rm s}^{-1}$(Quireza et al., 2006a). This indicates expansion or contraction of $\sim$4.5 km ${\rm s}^{-1}$, if the star cluster, which drives the shell, is at the systemic velocity. It is not possible to determine from CO and RRL emission alone, whether the shell is expanding or collapsing. The RRL emission appears to be associated more with the high-velocity peak. However, the exact peak position is uncertain, due to the low signal-to-noise ratio of the RRL emission at this position. Depending on the rest velocity of the star cluster, it may also be possible that expansion or collapse are asymmetric and that one part of the shell is at rest.

Anti-correlation between molecular and ionized emission in velocity is seen towards Pos. 1 and Pos. 3, offset from Cluster 1. At the peak velocities of the ionized gas emission in both positions, the strength of molecular gas emission is weaker in comparison to emission at other velocities. If both gas phases at a given velocity are spatially connected, this emission may be weaker as the molecular gas is destroyed.

The emission structure at both positions can also be affected by local star formation and does not have to be dominated by the shell dynamics alone. Position 3 harbors the low-cm continuum emission source “JJ”, while Pos. 1 contains the sources “M” and “L” (see Fig. 2; all three sources are characterized with low-cm continuum emission in De Pree et al., 1997). Position 1 is also located at the edge of the Welch-Ring which harbors multiple UC H ii regions and is dynamically complex (see e.g., Serabyn et al., 1993; Williams et al., 2004). The velocities from the high-order RRLs presented here are similar to the velocities of the individual UC H ii regions (see Sect. 3.3.1), which are determined at significantly higher angular resolution ($<1$″) from the H92$\alpha$ RRLs (De Pree et al., 1997). The velocity signature therefore may not only be shaped by large-scale motions, such as the expanding bubble, but also by the small-scale dynamics, which are spatially unresolved in the RRL observations presented in this work.

The precise time of the formation of the second generation of stars postulated by the models is difficult to constrain since the shell is either close to the end of the re-collapse or at the beginning of the second expansion phase. Towards the center of Cluster 1, we do not see evidence of UC H ii regions, many young stars or enhanced dust opacities (compared to the rest of W49A), which would be indicative of new stars having been formed there recently. However, high foreground extinction towards W49A ($A_{V}>30$ mag; Alves & Homeier 2003) and high uncertainties in stellar ages make it impossible to rule out that a second generation of stars has already formed at this position.

Furthermore, as W49A is clearly not spherically symmetric, it is also possible that Cluster 1 did not form at the center of mass and the shell may collapse to a different center. In the idealized case of a 1D geometry the cloud is forced to re-collapse onto the same cluster. Depending on the stellar mass distribution, which is intrinsically difficult to measure for the region, the re-collapsing gas may fall more towards W49N, which is $\sim$3 pc offset from Cluster 1 in projection (Alves & Homeier, 2003) and therefore may be within the maximum radius of expansion for large parts of the range of densities in agreement with the observations. The youngest clusters dominate feedback. So long as the shell expands to a radius which encompasses both populations the models will be valid in the future. However, without clear evidence of the presence or absence of a new population of stars at the center of the shell, the molecular cloud is either shortly before, after, or even within a new episode of star formation.

Most likely, W49A is currently in the phase of the formation of this second star cluster. Indication of this is the large number of compact and UC H ii regions surrounding Cluster 1, including the Welch-Ring. These may have been triggered during the re-collapse of the shell. Since some of them are extended (i.e., source ”L”; De Pree et al. 1997), they may already be contributing to the feedback of the bubble. These sources may provide the feedback to halt the re-collapse of the shell in the future, or may have already started the second expansion period of such a shell.

5.2 Effect of feedback on the molecular gas in W49A

The warpfield models were used in Sect. 4.2 to characterize the expansion of the shell radius. This section discusses to what extent other parts of W49A were affected by feedback and where this feedback may have triggered star formation, for example via triggering of clumps of molecular gas into gravitational collapse (e.g., Elmegreen & Lada, 1977; Whitworth et al., 1994; Preibisch et al., 2002; Preibisch & Zinnecker, 2007).

For all models highlighted in Fig. 6, the expanding shell remains confined to the inside of W49A (¡6 pc). The radius of maximum expansion is between 2.5 and 5 pc. This implies that feedback of Cluster 1 has affected molecular gas in W49N. The Welch-ring harbors UC H ii regions, which have an approximate lifetime of $3\times 10^{5}$ yr (Mottram et al., 2011). Therefore, these need to have formed after the O/B stars in Cluster 1. Feedback may have triggered the star formation in the Welch ring, either at the point of maximum expansion (Model 2; $n\sim 7.5\times 10^{3}\,{\rm cm}^{-3}$), or during contraction at re-collapse (Models 1,3,4; $n\lesssim 7.5\times 10^{3}\,{\rm cm}^{-3}$). Also, as discussed in Sect. 5.1, these UC H ii regions may already be part of the second episode of star formation after re-collapse.

Towards the lower bound of the considered range of densities (Model 3; $n\approx 2.5\times 10^{3}\,{\rm cm}^{-3}$), feedback could have pushed to the outskirts of the cloud. During the expansion, it may have triggered star formation in molecular clumps in the swept-up molecular gas, which could be the reason for the formation of O stars and UC H ii regions outside of Cluster 1, as well as larger star forming sites such as W49S and W49SW (Fig. 1). Shell fragmentation could be a possible cause for these clumps to be decoupled from the bulk collapse of the shell. This may occur especially when the expansion of the shell has considerably slowed down (e.g., McCray & Kafatos, 1987, Eq. 14 therein). This is expected at the expansion maximum and could explain why we observe them today at projected distances from Cluster 1 of $\sim 5-6$ pc.

We note that for Models 1 and 2, the molecular clumps W49S or W49SW at the edge of W49A, would not have been affected by the star formation in Cluster 1; they were of independent origin. Their formation may be connected to global motions, like cloud-cloud collisions (Serabyn et al., 1993) or fragmentation of filamentary inflow (Galván-Madrid et al., 2013). This may fit in the picture of sub-clustered star formation, which does not require a causal connection between individual events of star formation in W49A (e.g., Alves & Homeier, 2003).

It must be stated that the models do not predict any ejecta at distances $d>12\,{\rm pc}$, as discussed by Peng et al. (2010). In the picture of re-collapsing shell models, these would need to be formed independently, at least if feedback were driven only by Cluster 1. In principle, strong inhomogeneities in the molecular cloud could have allowed feedback to channel out of the surrounding molecular cloud. However, as mentioned by Peng et al. (2010), high energies would be necessary to drive these ejecta (few $\times 10^{50}\,{\rm erg}$).

The previous discussion highlights that the expansion radius strongly depends on the density structure. For the models in Sect. 4, we assumed a homogeneous density profile. If the density instead increased towards the molecular cloud center, the evolution of the shell would have a smaller radius of maximum expansion (see Sect. 4.2). However, the models would still predict a re-collapsing shell similar to Models 1 and 2. Since most of the expansion occurs in the momentum-driven phase (see Sect. 4.1), an inhomogeneous density distribution may affect the expansion of the shell locally. The expansion of the shell would necessarily be asymmetric, that is, expanding further in directions of lower densities and less in directions of higher densities as compared to the expansion in the case of an average, homogeneous density distribution as delineated here. Modeling the exact shape and the precise impact on the cloud of the feedback-driven shell requires 3D simulations, as further discussed in Sect. 5.5, and requires us to make assumptions regarding the initial density distribution of the cloud. Keeping in mind the limitation of such an analysis, the models presented here give an estimation of the evolution of the feedback-driven shell using the simplest assumption on a density profile, namely that it is homogeneous.

To summarize, the warpfield models show that the feedback shell has affected different subsections of the molecular cloud, depending on the assumed density distribution of the molecular cloud. The compression of gas by the feedback-driven shell potentially triggers gravitational collapse of dense gas clumps and the formation of new stars. While widespread star formation is observed around Cluster 1, the effectiveness of triggering remains disputed in the literature (e.g., Dale & Bonnell, 2011; Krumholz et al., 2014), and is yet to be demonstrated to alter the star formation efficiency of molecular clouds in general and of W49A in particular. Nonetheless, we speculate that this can constitute an additional hypothesis for the high star formation activity at least of parts of W49A. Star formation in regions of the cloud that were not affected by feedback needs to have occurred in different ways. Therefore, these models do not necessarily rule out other hypotheses on the formation scenario of stars in W49A, such as sub-clustered star formation or cloud-cloud collisions, which have already been discussed in the literature.

5.3 Star formation rate and multiple generations of stars

Re-collapsing shells predict episodic events of star formation in a molecular cloud. While the previous section discussed potential triggering of star formation during the expansion of the shell into the molecular cloud of W49A, here we focus on the formation of a new star cluster when the shell re-collapses. At this time, the accumulated molecular gas of the shell-material is concentrated inside a small radius and forms a new star cluster (Rahner et al., 2017, 2018a).

At each re-collapse event, stars are formed at a certain instantaneous star formation efficiency (SFE),

[TABLE]

with $M_{*}$ the mass of the star cluster originating from a single star formation event, and $M_{\rm cloud}$ the cloud mass prior to the formation of a cluster (e.g., Murray, 2011). If star clusters continue to form with a fixed SFE for each star formation event (we adopt values of $\epsilon=0.05$ and 0.09), it would take multiple re-collapse events (each of which leads to a new star formation event) to disperse the cloud. Since the main focus of the model comparison was to follow the first expansion cycle of the feedback shell from Cluster 1, we assume here for simplicity that the SFE at future collapse events is the same as at the initial cluster formation. Changes in the SFE, however, are expected (as, e.g., in 30 Dor; Rahner et al. 2018a), since the SFE depends on the physical conditions of the re-collapsing gas, which are likely to be different from those during the first episode of star formation. The recently formed stars will eventually power a common feedback shell together with the already existing star cluster, which will affect the evolution of this shell and change the time-scale between re-collapse events.

The re-collapse events have implications for the age distribution of the future star cluster. Different generations of stars would be expected in W49A, induced by recurring collapse events and the subsequent formation of a new star cluster. Their age difference in the case of W49A will, however, be small compared to other star-forming regions, such as 30 Doradus (see Sect. 5.4). The models for Cluster 1 predict an age difference of only 1-2 Myr. The wide-spread star formation in W49A suggests that stars may also form in between (see Sect. 5.2). The age differences are therefore likely to be lower than 1 Myr, which is well within typical uncertainties of the ages derived for O stars, and therefore intrinsically difficult to differentiate observationally.

5.4 W49A in comparison to the massive star-forming region 30 Doradus

One of the most prominent prototypes for an extreme high-mass star-forming region can be found in 30 Doradus, which is located in the LMC. It has the observational advantage of being visible in the optical, which is not possible for W49A due to the obscuration of the dust by the Milky Way. There are obvious differences between both regions in molecular gas mass (depending on method and radius: $10^{6-7}\,{\rm M_{\odot}}$ for 30 Doradus; e.g., Dobashi et al. 2008; Faulkner 1967; Sokal et al. 2015; vs. $10^{5.3-6}$ for W49A; e.g., Galván-Madrid et al. 2013) and stellar cluster mass ($7\times 10^{4}-5\times 10^{5}$ ${\rm M}_{\odot}$ for 30 Doradus; Selman et al. 1999; Bosch et al. 2001, 2009; Cignoni et al. 2015; vs. $5-7\times 10^{4}$ ${\rm M}_{\odot}$ for W49A; Homeier & Alves 2005). Here, we focus the discussion on feedback-related star formation in both star-forming regions.

The star cluster NGC 2070 in 30 Doradus contains a sub-cluster, R136, which is younger ($\sim 1\,{\rm Myr}$) than the older stellar population of NGC 2070 ($\sim 5\,{\rm Myr}$; Brandl et al. 1996; Walborn & Blades 1997; Massey & Hunter 1998; Selman et al. 1999; Sabbi et al. 2012; Cignoni et al. 2015; Crowther et al. 2016). Rahner et al. (2018a) found that the formation of R136 can be naturally connected to the formation of NGC 2070 by the re-collapse of an expanding shell from the first star formation event, and the subsequent formation of R136. While the stars in W49A are of similar age to R136, there is no distinct older population. Rather, there seems to be only an even younger population of UC H ii regions present in W49A, with an age difference to the O-star population of 1-2 Myr (see Sect. 5.2). The model confirms this difference in star formation history: while for 30 Dor, the models find re-collapse events 3-4 Myr after the initial event of star formation, the re-collapse is predicted earlier for W49A, within 2 Myr after initial star formation.

The larger re-collapse interval in 30 Doradus is due to a combination of different cloud density, stellar cluster mass, and metallicity. The models require higher initial molecular gas densities in W49A to satisfy the observational constraints ($n\gtrsim 2.5\times 10^{3}\,{\rm cm}^{-3}$) than have been found for 30 Doradus ($n\gtrsim 5\times 10^{2}\,{\rm cm}^{-3}$). The mass of the first cluster of stars is a factor of a few higher in 30 Doradus and its metallicity is lower at $Z\approx 0.5\,{\rm Z_{\odot}}$ (Lebouteiller et al., 2008; Choudhury et al., 2016).

The lower metallicity in 30 Doradus changes the shell evolution in a nontrivial manner. While Lopez et al. (2011) suggest weaker winds in low-Z systems, Rahner et al. (2017) finds that the longer cooling time can lead to winds being more important, as the expansion of the shell stays in phase 1, the energy-driven phase, for a longer time. Pellegrini et al. (2011) finds that the shell structure is over-pressured near the cluster by X-ray bubbles over radiation pressure, consistent with the results by Rahner et al., and in opposition to Lopez et al. (2011). In turn, accelerations of the shell in phase 2, the momentum-driven phase, may decrease. Low metallicities imply weaker stellar wind luminosities (e.g., Lopez et al. 2011) and affect the coupling of radiation to the shell by changing the shell structure (Rahner et al., 2017): The confining pressure from the winds of the shell is weaker, which leads to lower densities in the shell (see Rahner et al., 2017, eq. 14). In turn, these lower densities result in a decrease in radiation pressure on a given part of the shell, as it absorbs less ionizing radiation than a high density shell, because the absorption of ionizing radiation by hydrogen depends on the recombination rate and is proportional to $n^{2}$. The effect of metallicity is therefore complicated and depends on the details of the investigated model whether there is faster or slower expansion. Further, Pellegrini et al. (2011) found significant fractions of the 30 Doradus region to be dominated by optically thin (fully ionized) gas leading to a decoupling between radiation and the gas, a result that can only be explained by fully accounting for feedback and ISM coupling and cooling, and not from simple scaling relations.

These models provide a framework with which to understand how feedback (both positive and negative) regulates star formation and the associated timescales. The fact that models with re-collapse of a feedback shell are compatible with observations in W49A indicates that re-collapse may not be a phenomenon that is unique to 30 Doradus. This would imply that feedback is not only responsible for cloud dispersal, but could also be the means by which star formation regulates the rate of stellar birth by introducing a timescale for subsequent cluster-formation events. The timescale would be much shorter (¡2 Myr) in W49A than in 30 Doradus, and consequently it would also be harder to observe (see Sect. 5.3). Other clusters besides 30 Doradus and potentially W49A may be affected by feedback-regulated star formation episodes. According to Rahner et al. (2018a), these need to have a SFE of $\lesssim 10\%$ and average densities of $n\geq 5\times 10^{2}\,{\rm cm^{-3}}$. Investigation of clusters that meet these criteria is necessary to decipher whether or not this is a common scenario in the formation of massive star clusters.

5.5 Future investigations

Detailed observations are needed of the morphology of the molecular gas and the stellar content as seen in infrared wavelengths in order to shed more light on the physical processes dominating the dynamics and evolution of W49A. While much work has been done on W49N, more studies focussing on the connection of the large-scale dynamics inside and outside of the region are necessary. With the Atacama Large Millimeter/sub-millimeter Array (ALMA888http://www.almaobservatory.org), kinematics of molecular and ionized gas can be mapped at significantly higher resolution than the observations of ionized and molecular gas applied here. Due to the high extinction towards W49A of $A_{V}>30$ mag (Alves & Homeier, 2003), soft, diffuse X-ray emission from wind bubbles (Townsley et al., 2003, 2014, 2018) are difficult to observe. Chandra observations of hard X-rays for W49A exist, and have been discussed for W49N, including both wind-driven bubbles and other possible formation scenarios (Tsujimoto et al., 2006). An interpretation of these in the light of wind-driven feedback is therefore nontrivial and remains to be explored for Cluster 1. The James Webb Space Telescope (JWST999https://www.jwst.nasa.gov) will yield high-sensitivity imaging of many star-forming regions in the Galaxy. This will provide better insight into the most embedded populations of O stars in regions like W49A.

The 1D models from Rahner et al. (2017) are well suited for exploring a large parameter space in molecular cloud mass, star formation efficiency, and density. While these models contain all the necessary physics to describe the expanding shell, modeling a shell in three spatial dimensions is necessary. Three-dimensional models can better account for escape routes for feedback in the energy-driven phase at the beginning of the expansion, as well as for instabilities in the shell, which may alter the effectiveness of feedback (Dale et al., 2014), even though the impact of the shell structure on the escape fraction of ionizing radiation is accounted for in the 1D models. More importantly, it is the density structure that has the greatest effect on the geometry and the timescales of the expansion. As seen in W49A, many star-forming regions appear not to be spherically symmetric, and expansion will not occur isotropically into a uniform medium. These models are naturally more computationally intensive. Such studies have been carried out (e.g., Howard et al., 2016; Peters et al., 2010), although these have not yet included all the relevant physics. However, given that with the 1D models we can explore the parameter spaces very well, the results from this work may provide the average initial cloud parameters for more detailed modeling.

6 Conclusions

The study presented here investigates hydrogen radio recombination line emission of the H151$\alpha$-H156$\alpha$ and H158$\alpha$ transitions at frequencies between 1.6 and 1.9 GHz with the VLA in C-configuration towards W49A from the THOR survey (Beuther et al., 2016), as well as archival CO observations from the CHIMPS and GRS surveys (Rigby et al., 2016; Jackson et al., 2006). We find shell-like RRL emission around the infrared star cluster at the center of W49A, Cluster 1. We find double peaked 13CO(3-2) emission towards Cluster 1, again indicative of emission from shell-like geometries. The RRL emission is rather broad over all velocities between the two CO components, but peaking towards the redshifted one. Towards positions at the edges of the shell-like bubble, the velocities of the ionized and molecular gas emission are anti-correlated. This anti-correlation may be a signature of interaction of the ionized shell with the molecular envelope.

One-dimensional models of expanding shells around massive star clusters (warpfield, Rahner et al., 2017) are used to investigate the evolution of the shell-like bubble in W49A. Given the observational initial conditions of stellar cluster mass, age of the O-star population, and the current radius of the shell, all models predict re-collapse of the shell after the first star formation event. Feedback of the first formed cluster is therefore not strong enough to disperse the cloud in the first instance.

The evolution of the shell however strongly depends on the assumed densities. In some of the models, the shell has expanded to the outskirts of W49A. With this, it could have affected molecular gas in the entire region of W49A. In this case, a causal connection between feedback from cluster 1 and star formation all across W49A is possible. However, for most of the models in the range of the assumed densities, the feedback shell did not expand to the outskirts of W49A. Hence, it is more likely that only limited parts of W49A were affected by feedback from the central stellar cluster, while stars in the outer parts of W49A formed independently. To what extent this feedback altered the physical conditions in the surroundings and nurtured (i.e., triggered) star formation, is left to future studies.

This modeling presents another alternative to sub-clustered star formation and cloud-cloud collision models, without ruling out either. In comparison to the star-forming region 30 Doradus, indications for re-collapse are less clear due to difficulties in identifying two distinct stellar populations. On the other hand, as W49A is younger, it is likely to continue to form stars, possibly in an episodic mode.

Acknowledgements.

We would like to thank the referee for the detailed, helpful, and insightful comments, which considerably improved the paper. M.R.R. is a fellow of the International Max Planck Research School for Astronomy and Cosmic Physics (IMPRS) at the University of Heidelberg. H.B., M.R.R., Y.W., J.S. and J.C.M. acknowledge support from the European Research Council under the Horizon 2020 Framework Program via the ERC Consolidator Grant CSF-648505. M.R.R., D.R., H.B., E.W.P, S.C.O.G. and R.S.K. acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) via Sonderforschungsbereich (SFB) 881 “The Milky Way System” (sub-projects B1, B2 and B8). S.C.O.G., E.W.P. and R.S.K. further acknowledges support from the DFG via Priority Program SPP 1573 “Physics of the Interstellar Medium” (grant numbers KL 1358/18.1, KL 1358/19.2, and GL 668/2-1) and from the European Research Council via the ERC Advanced Grant STARLIGHT (project number 339177). The research was carried out in part at the Jet Propulsion Laboratory, which is operated for NASA by the California Institute of Technology. R.J.S. acknowledges support from an STFC Ernest Rutherford fellowship. S.E.R. acknowledges support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement # 706390. F.B. acknowledges funding from the European Union’s Horizon 2020 research and innovation programme (grant agreement No 726384 - EMPIRE).

Appendix A Additional plots