Main

At a distance of only 8.2 kpc from the Sun (ref. 9), the Galactic Centre provides a unique laboratory for studying the complex physical processes that occur within a galactic outflow. The ‘Fermi bubbles’2,10, two giant lobes extending up to ~10 kpc from the Galactic plane, are thought to outline the current boundaries of the Milky Way’s nuclear wind. Several hundred neutral gas clouds have been found recently within this volume through observations of the atomic hydrogen (H i) line at a wavelength of λ = 21 cm (refs. 7,8). Figure 1 shows a column density map of H i clouds in the nuclear wind8 detected with the Green Bank Telescope (GBT). Although the bulk of the cloud population lies within the boundaries of the Fermi bubbles (green dashed line11), it has not been established whether this outflowing H i gas arises from the same event that generated the Fermi bubbles. These clouds were identified through their anomalous line-of-sight velocities, which are incompatible with Galactic rotation and can instead be described using a biconical wind model in which clouds accelerate from the Galactic Centre, reaching a maximum velocity of 330 km s−1 after about 2.5 kpc (refs. 8,12). To assess whether outflowing H i structures carry molecular gas, we targeted two objects (hereafter, MW-C1 and MW-C2), highlighted by red boxes in Fig. 1, in the 12CO(2 → 1) emission line at 230.538 GHz with the 12-m Atacama Pathfinder Experiment (APEX) telescope. These two clouds have relatively high H i column densities (>1019 cm−2) and show an elongated head-to-tail morphology along the direction pointing away from the Galactic Centre. We mapped both clouds in 12CO(2 → 1) emission over a 15′ × 15′ field centred on the peak of the H i emission, at a spatial resolution of 28″ (full-width at half-maximum, FWHM), corresponding to ~1 pc at the distance of the Galactic Centre, and a spectral resolution of 0.25 km s−1. These data revealed molecular gas outflowing from the centre of our Galaxy.

Fig. 1: Atomic hydrogen gas outflowing from the Galactic Centre.
figure 1

The colour scale shows the column density of anomalous H i clouds in the Milky Way’s nuclear wind, detected with GBT8. The green dashed line is the boundary of a volume-filled model for the Fermi bubbles11. The two H i clouds observed in the 12CO(2 → 1) line with APEX are marked by red boxes.

Full size image

Figure 2 shows H i column density maps (Fig. 2a, c) from GBT observations and integrated brightness temperature maps (Fig. 2b, d) from the 12CO(2 → 1) line obtained with APEX for MW-C1 and MW-C2. Higher-resolution H i data from the Australia Telescope Compact Array (ATCA) for MW-C2 are also overlaid as contours on the CO map. CO velocity fields and three representative spectra across each field are presented in Fig. 3. CO emission is detected in both H i clouds, with substantial morphological and kinematical differences between them. MW-C1 shows five distinct compact clumps of molecular gas concentrated towards the part of the H i cloud that faces the Galactic Centre (arrows in Fig. 2). At least three clumps have a velocity gradient along the direction pointing towards the tail of the H i cloud. All the CO emission in MW-C1 lies in the local-standard-of-rest (LSR) velocity range VLSR ≈ 160–170 km s−1. Typical FWHM line widths are ~2–3 km s−1 (see spectra in Fig. 3). By contrast, in MW-C2 most of the CO emission is distributed along a filament-like structure, with some fainter and more diffuse clumps in the region away from the Galactic Centre. CO emission is spread over a larger velocity range than in MW-C1, spanning 30 km s−1 over VLSR ≈ 250–280 km s−1, and the velocity field does not show any clear ordered motion. 12CO(2 → 1) line profiles in MW-C2 are much broader than in MW-C1, with an FWHM ranging from ~5 km s−1 to 12 km s−1.

Fig. 2: Atomic hydrogen and molecular gas in two clouds in the Milky Way’s nuclear wind.
figure 2

a, c, H i column density maps from GBT data8 at an angular resolution of 570″ for MW-C1 (a) and MW-C2 (c). Black arrows point towards the Galactic Centre (GC). Red boxes highlight the 15′ × 15′ fields observed with APEX. Contour levels are at (0.2, 0.5, 1, 2, 4) × 1019 cm−2. b, d, 12CO(2 → 1) integrated brightness temperature maps from APEX data at 28″ resolution for MW-C1 (b) and MW-C2 (d). H i contours at (4, 8, 16, 24) × 1020 cm−2 from ATCA data at 137″ resolution are overlaid on the MW-C2 map in d. The circles at the top left of each panel show the angular resolution of the telescopes. RA, right ascension; dec., declination; WCO, integrated brightness temperature.

Full size image
Fig. 3: Molecular gas kinematics in MW-C1 and MW-C2.
figure 3

a, b, Velocity fields derived from a Gaussian fit to the 12CO(2 → 1) data for MW-C1 (a) and MW-C2 (b). c, d, 12CO(2 → 1) spectra for MW-C1 (c) and MW-C2 (d) at the positions labelled in a and b. We note the differences in velocity spread and line shape between MW-C1 and MW-C2. Tb is the line brightness temperature.

Full size image

The observed features indicate that cold gas in MW-C2 is interacting and mixing with the surrounding medium more efficiently than in MW-C1, resulting in a more turbulent molecular gas. An interpretation of the differences in the morphokinematics of the molecular gas in the two clouds is that we are witnessing two evolutionary stages of a cold cloud being disrupted by interaction with a hot flow. Our idealized biconical wind model12 with a maximum wind velocity of 330 km s−1 places MW-C1 at a distance of 0.8 kpc and MW-C2 at a distance of 1.8 kpc from the Galactic Centre, implying that MW-C2 may have been within the nuclear outflow twice as long (7 Myr, versus 3 Myr for MW-C1). Our model also predicts that MW-C2 is moving faster than MW-C1 (~300 km s−1 versus ~240 km s−1). In the classical picture, in which cold gas is entrained in the hot wind, MW-C1 may therefore represent an early stage of the interaction with the surrounding medium, at which molecular gas is still relatively intact and undisturbed near the initial dense core; whereas molecular gas in MW-C2 could have been stripped off from its core, resulting in a disordered morphology/velocity field and broader linewidths. However, the observed characteristics of the two clouds may also be explained in terms of different local conditions of the hot outflow. A larger and more complete sample of molecular gas detections in outflowing clouds is needed to provide a more robust picture.

The two clouds analysed in this work have atomic gas masses of Mat ≈ 220M☉ (MW-C1) and Mat ≈ 800M☉ (MW-C2), as derived from the GBT H i data (M☉, mass of the Sun). All mass measurements from observations are scaled by a factor of 1.36 to account for the presence of helium. It is not straightforward to estimate the mass of molecular matter, because the gas may have considerable opacity in the 12CO(2 → 1) line and the appropriate CO-to-H2 conversion factor XCO in the Milky Way’s wind is unknown. We used the observed CO integrated brightness temperatures, cloud radii and line widths to constrain the acceptable XCO values by means of chemical and thermal modelling of a cloud undergoing dissociation by photons and cosmic rays. We found that XCO for the 12CO(2 → 1) transition in our clouds lies in the range (~2–40) × 1020 cm−2 (K km s−1)−1. The lowest value, XCO = 2 × 1020 cm−2 (K km s−1)−1, is consistent with the Galactic conversion factor13, and was used to derive lower limits to the molecular gas mass Mmol. We obtained Mmol ≥ 380M☉ for MW-C1 and Mmol ≥ 375M☉ for MW-C2, implying molecular-to-total gas mass fractions of fmol = Mmol/(Mmol + Mat) ≥ 0.64 and fmol ≥ 0.32, respectively. We emphasize that these values are lower bounds and the molecular gas mass may be higher by a factor of ten. As a consequence, the total mass of molecular gas in the nuclear wind of the Milky Way is large. Under the conservative assumption of an average fmol ≈ 0.3–0.5 for all outflowing H i clouds in the GBT sample, and using an atomic outflow rate8 of \({\dot{M}}_{{\rm{at}}}\approx 0.1{M}_{\odot }\,{{\rm{yr}}}^{-1}\), we estimated an outflow rate of \({\dot{M}}_{{\rm{a}}{\rm{t}}}\ge (0.05\mbox{--}0.1){M}_{\odot }\,{{\rm{y}}{\rm{r}}}^{-1}\) in molecular gas. This value is of the same order of magnitude as the star formation rate (SFR) of the Central Molecular Zone14 (CMZ), implying a molecular gas loading factor \(\eta ={\dot{M}}_{{\rm{mol}}}/{\rm{SFR}}\) at least of the order of unity at a distance of ~1 kpc from the Galactic plane, similar to that estimated in nearby starburst galaxies15. This cold outflow affects the gas cycle in the inner Galaxy and may constitute an important mechanism that regulates the star formation activity in the CMZ.

From a theoretical point of view, such a large amount of high-velocity molecular gas is puzzling16. It is believed that cool gas in a disk can be lifted and accelerated by both drag force from a hot outflow17 and by radiation pressure18. This requires a source of strong thermal feedback and/or radiation feedback. The Milky Way does not currently have an active galactic nucleus (AGN), nor is the SFR of the inner Galaxy comparable to that of starburst galaxies with known molecular winds (for example, NGC253)15. Current simulations of AGN-driven winds have focused on very powerful AGNs19,20 and there have been no investigations studying whether a relatively small black hole like Sagittarius A* could expel large amounts of cold gas, even if it had undergone a period of activity in the recent past. On the other hand, the current SFR of the CMZ is not large enough to explain the estimated outflow rate of cold gas21, and no observational evidence so far suggests a sizable change in the SFR of the CMZ in the last few million years22. A scenario in which the star formation in the CMZ is episodic on a longer cycle23,24 (10–50 Myr) and is currently near a minimum might help to partly reconcile the observed and predicted cool gas mass loading rates, although our wind model suggests that the lifetimes of cold clouds are shorter than 10 Myr. Cosmic rays are also believed to contribute to the pressure on cold gas25, but their role is only just starting to be understood and needs observational constraints. Moreover, in either an AGN- or a starburst-driven wind, the extent to which cold gas survives under acceleration is a matter of debate17,26, and several different mechanisms have been investigated to extend the lifetime of cool gas in a hot wind (for example, magnetic fields27 and thermal conduction28). An alternative scenario has been recently proposed in which high-velocity cool neutral gas (temperature T < 104 K) forms directly within the outflow as a consequence of mixing between slow-moving cool clouds and the fast-moving hot wind29,30. This mechanism overcomes the problem of accelerating dense material without disrupting it, and may explain the high velocities observed in cool outflows. However, current simulations cannot trace the gas down to the molecular phase.

In conclusion, this detection of ouflowing cold molecular gas in the Milky Way is a challenge for current theories of galactic winds in regular star-forming galaxies, because none of the above processes seems able to easily explain the presence of fast molecular gas in the Milky Way’s wind. Targeted observations of molecular gas tracers in the Milky Way’s nuclear wind are expected to contribute considerably to our understanding of these fascinating phenomena.

Methods

Observations and data reduction

Observations of the 12CO(2 → 1) emission line at 230.538 GHz were made with the 12-m APEX antenna31 using the PI230 heterodyne receiver (ESO project 0104.B-0106A; principal investigator E.M.D.T.). The spectrometer32 covers a bandwidth of 8 GHz at a spectral resolution of 61 kHz, corresponding to a velocity resolution of about 0.08 km s−1 at 230 GHz. The beam size at this frequency is 27.8″ (FWHM), the main-beam efficiency is 0.72 and the jansky-to-kelvin conversion factor is 40 ± 3. We observed our targets in on-the-fly position-switching mode, integrating for 1 s every 9″. Both fields were 15′ × 15′ wide, centred at (RA, dec.)J2000 = (17 h 56 min 34.0 s, −32° 29′ 14″) for MW-C1 and at (RA, dec.)J2000 = (17 h 18 min 22.2 s, −27° 56′ 28″) for MW-C2. The observed regions are shown in red boxes in Figs. 1, 2. The total integration time was approximately 25 h for each field. Throughout the observing session (September to November 2019), the precipitable water vapour varied between 0.6 mm and 3 mm.

We reduced the data using the Continuum and Line Analysis Single-dish Software (CLASS) from the GILDAS package33. A first-order baseline was subtracted from the calibrated spectra by interpolating the channels outside the velocity windows in which we expected to see the emission based on the H i observations. The spectra were then smoothed in velocity and mapped onto a grid with a pixel size of 9″ and a channel width of 0.25 km s−1. The root-mean-square noise (σrms) in the final data cubes was 65 mK and 55 mK for MW-C1 and MW-C2, respectively, in a 0.25 km s−1 channel.

Atomic gas and molecular gas mass

The H i GBT data8 and the 12CO(2 → 1) APEX data were analysed to estimate the atomic and molecular gas masses, respectively. First, the three-dimensional source finder DUCHAMP34 was applied to the data cubes to identify regions of sizable emission. During this process, we set a primary threshold to identify emission peaks at 5σrms and reconstructed sources by adding pixels down to a secondary threshold of 2.5σrms.

The column density at a given position (x, y) on the sky can be written as:

$${N}_{{\rm{H}}}(x,y)=C\int {T}_{{\rm{b}}}(x,y,v){\rm{d}}v,$$
(1)

where the integral considers pixels in only one detection, Tb is the line brightness temperature, dv is the channel width (1 km s−1 for GBT and 0.25 km s−1 for APEX) and C is a constant. For the H i line, under the assumption that the gas is optically thin, the constant is35 C = 1.82 × 1018 cm−2 (K km s−1)−1. For CO lines, this constant is also known as the CO-to-H2 conversion factor XCO (ref. 13). Because the conversion factor in the nuclear wind cannot be constrained with existing data, we used the value estimated in molecular clouds in the disk of the Milky Way36, XCO = 2 × 1020 cm−2 (K km s−1)−1. We checked this XCO value against the predictions of radiative-transfer models, described in the next section, and found that the Galactic value is probably a lower limit for clouds in the nuclear wind.

The total mass of gas can be calculated as:

$$M=1.36m{D}^{2}\int {N}_{{\rm{H}}}(x,y){\rm{d}}x{\rm{d}}y,$$
(2)

where the factor 1.36 takes into account helium, D ≈ 8.2 kpc is the adopted distance to the clouds, m is the mass of atomic/molecular hydrogen for atomic/molecular gas, dx and dy are the pixel sizes in radians (105″ for GBT, 9″ for APEX). The observed properties and estimated masses are summarized in the Extended Data Table 1.

Radiative-transfer models

We used the chemistry and radiative-transfer code DESPOTIC37 to constrain the CO-to-H2 conversion factor of the clouds. DESPOTIC computes the chemical and thermal state of an optically thick cloud given its volume density and column density. The turbulent velocity dispersion of the gas was assumed to be 1–5 km s−1 (see Fig. 3) in our modelling. The chemical equilibrium calculation uses solar abundances for dust and all elements in the H−C−O chemical network38, whereas the thermal equilibrium calculation includes heating by cosmic rays, the grain photoelectric effect, cooling by the H i, C i, C ii, O i and CO lines, and collisional energy exchange between dust and gas. Level populations were calculated using an escape probability method, with escape probabilities estimated using the spherical geometry option of DESPOTIC.

We investigated different combinations of the interstellar radiation field χ and the cosmic-ray ionization rate ζ through a set of DESPOTIC models using log(χ/G0) = [−1, 0, 1, 2], where G0 is the solar radiation field39 and log[ζ (s−1)] = [−16, −15, −14]. The interstellar radiation field was varied between subsolar (χ ≈ 0.1G0) and highly supersolar (χ ≈ 100G0) values, representative of a highly star-forming environment like the CMZ. The cosmic-ray ionization rate ranges from the value measured in the solar neighbourhood40 (ζ ≈ 10−16 s−1) to the estimated upper limit for the CMZ41 (ζ ≈ 10−14 s−1). We stress that our CO clouds lie at about 1 kpc from the Galactic plane and that both the interstellar radiation field and the cosmic-ray ionization rate are expected to drop with distance from the disk. Therefore, although the estimated values of χ and ζ in the CMZ are orders of magnitude higher than in the solar neighbourhood, models with intermediate interstellar radiation fields and cosmic-ray ionization rates should be more representative of the conditions high in the Milky Way’s wind.

For each model, DESPOTIC returned the 12CO(2 → 1) integrated brightness temperatures (WCO) as a function of the number density (nH2) and column density (NH2) of molecular hydrogen. We only considered solutions consistent with the observed integrated brightness temperature (1–5 K km s−1; see Fig. 2) and observed cloud radius of R = 0.75nH2/NH2 = 1–5 pc, and we calculated the expected CO-to-H2 conversion factor XCO = NH2/WCO for the 12CO(2 → 1) transition. We found that there are no acceptable solutions for a strong interstellar radiation field (log(χ/G0) ≥ 1), which indicates that molecular clouds with the observed properties cannot exist in the presence of a CMZ-like radiation field. Instead, models with solar and subsolar radiation fields returned solutions compatible with the observational constraints for any cosmic-ray ionization field. An interstellar radiation field weaker than the one produced in the CMZ is therefore more representative of the environment at 1 kpc above the Galactic Centre. The predicted XCO varies by an order of magnitude, ranging between ~2 × 1020 cm−2 (K km s−1)−1 and ~4 × 1021 cm−2 (K km s−1)−1, depending on the combination of radiation field and cosmic-ray ionization rate. The value of XCO = 2 × 1020 cm−2 (K km s−1)−1 that is commonly assumed in the Milky Way disk13 and used in this study is consistent with the smallest values returned by our radiative-transfer models, obtained with a weak, subsolar radiation field and a solar-like cosmic-ray ionization rate of ζ = 10−16 s−1. As a consequence, the molecular gas masses calculated in this work probably represent lower limits to the real cold gas mass in our CO clouds.

Wind kinematic model

To estimate the position, velocity and lifetime of MW-C1 and MW-C2, we used a biconical wind model8,12 calibrated on the full population of H i clouds. This model is based on the assumption that clouds were launched from a small region close to the centre of the Galaxy and are moving with a purely radial velocity Vw(r), where r is the distance from the Galactic Centre. For simplicity, we considered models of the form:

$${V}_{{\rm{w}}}(r)=\left\{\begin{array}{cc}{V}_{{\rm{i}}}+({V}_{max}-{V}_{{\rm{i}}})\frac{r}{{r}_{{\rm{s}}}} & {\rm{f}}{\rm{o}}{\rm{r}}\,r < {r}_{{\rm{s}}}\\ {V}_{max} & {\rm{f}}{\rm{o}}{\rm{r}}\,r\ge {r}_{{\rm{s}}}\end{array}\right.,$$
(3)

where Vi is the initial velocity at r = 0 and rs is the scale distance at which the maximum velocity Vmax is reached. Equation (3) describes a kinematic model in which clouds are subjected to a constant acceleration up to rs and maintain a constant velocity at distances r ≥ rs. Although equation (3) is purely empirical and chosen to reproduce the H i data, recent hydrodynamical simulations of starburst-driven winds have found qualitatively similar trends for the velocity of the cool gas with distance30. The LSR velocity VLSR of a cloud travelling in the wind and seen at Galactic coordinates (l, b) can be written as:

$${V}_{{\rm{L}}{\rm{S}}{\rm{R}}}(l,b,r)={V}_{{\rm{w}}}(r)[\sin \phi \sin b-\cos \phi \cos b\cos (l+\theta )]-{V}_{0}\sin l\sin b,$$

where the polar angle φ and the azimuthal angle θ can be easily written as a function of (l, b, r) (ref. 8) and V0 = 240 km s−1 is the rotation velocity of the LSR around the Galactic Centre42. In our model, clouds are restricted inside a bicone with half-opening angle φmax. We constrained the four free parameters of this model, Vi, Vmax, rs and φmax, by matching the LSR velocity distributions predicted by our model with that observed from the H i cloud population8,12. Our fiducial model is a biconical wind with opening angle φmax = 70°, where clouds accelerate from an initial velocity of Vi = 200 km s−1 to a maximum velocity of Vmax = 330 km s−1 at rs = 2.5 kpc. According to this wind model, MW-C1 and MW-C2 have travelled a distance of 0.8 kpc and 1.8 kpc from the Galactic Centre in about 3 Myr and 7 Myr, and their current outflow velocity is about 240 km s−1 and 300 km s−1, respectively.