# Insights into the collapse and expansion of molecular clouds in outflows from observable pressure gradients

Kalliopi M. Dasyra, Georgios F. Paraschos, Thomas G. Bisbas, Francoise

Combes, Juan Antonio Fernández-Ontiveros

► To cite this version:

Kalliopi M. Dasyra, Georgios F. Paraschos, Thomas G. Bisbas, Francoise Combes, Juan Antonio Fernández-Ontiveros. Insights into the collapse and expansion of molecular clouds in outflows from observable pressure gradients. Nature Astronomy, 2022, 6, pp.1077-1084. 10.1038/s41550-022-01725-9. insu-03851616

HAL Id: insu-03851616

https://insu.hal.science/insu-03851616v1

Submitted on 7 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Insights into the collapse and expansion of molecular clouds in outflows from observable pressure gradients

Kalliopi M. Dasyra ID 1,2✉, Georgios F. Paraschos ID 1,2,3, Thomas G. Bisbas2,4, Francoise Combes ID 5 and Juan Antonio Fernández-Ontiveros ID 6

The jets launched by actively accreting black holes can generate massive outflows in galaxies, which could suppress or enhance star formation by rarefying or compressing clouds of molecular gas. To study the stability of such jet-impacted clouds, we performed astrochemical, thermally balanced, radiative transfer modelling of the CO and HCO⁺ emission of the galaxy IC 5063. We found that jet-related mechanical heating and cosmic rays contribute to the molecular gas heating rate and could even individually sustain it. Clouds excited by these mechanisms have temperatures and densities reflecting an order-of-magnitude increase in their internal pressure. Variations of their external pressure, deduced from [S II] and [N II] ionized gas emission, further reveal that some clouds are undergoing rarefaction and others compression. Our work shows a new viewpoint on plausible links between galactic outflows and star formation conditions: that of observable pressure gradients. It also emphasizes the role of cosmic rays in contributing to these gradients.

In the past decade, numerous detections of outflows, that is, winds of molecular gas, revealed that a high fraction (often 1–20%) of the gas content of dense clouds in galaxies can be kinematically and dynamically perturbed<fnref n="1" />,². These detections raised questions about the suppression³⁻⁵ or enhancement⁶⁻⁸ of star formation in galactic regions with outflows due to the rarefaction or compression of clouds. An excellent laboratory to investigate for hydrodynamical changes of clouds in outflows is the galaxy IC 5063. The collision of radio-emitting jets with dense interstellar medium (ISM) clouds near the nucleus of IC 5063 results in the initiation of several discrete, spatially offset outflows along or near the jet trail. These outflows have been seen in [Fe II], H₂, CO and other lines⁹. Even the expansion of a jet-inflated cavity (or ‘cocoon’) has been proposed based on HⅠ (ref. ¹⁰) and CO (ref. ¹¹) data. A first study of the molecular gas excitation based on the CO(4–3)/CO(2–1) flux ratio revealed a notable temperature increase from ambient (not jet-impacted) to jet-impacted regions. Some molecular gas in the outflow has temperatures >100 K (ref. ⁴). It is thus plausible that, along with the change in the excitation state of jet-impacted clouds, a change in their dynamical state also occurred.

The stability, expansion or contraction of clouds is determined by the interplay of their self gravity with changes in their internal and external pressure. Despite the importance of pressure in this process, pressure maps for the molecular gas and its surrounding medium in regions with outflows are lacking from the literature. They have only been presented for the X-ray emitting gas delineating the surfaces of jet-inflated bubbles, such as the cavities in Perseus A¹² or the Fermi bubbles in the Milky Way¹³. Here we present pressure maps of the tenuous ionized medium surrounding molecular clouds in the nucleus of IC 5063 from Very Large Telescope (VLT) Multi Unit Spectroscopic Explorer (MUSE) and Spectrograph for INTEGRAL

Field Observations in the Near Infrared (SINFONI) data. We also present pressure maps of the molecular gas itself, computed from the deepest available CO and HCO⁺ Atacama Large Millimeter/ submillimeter Array (ALMA) cubes, which we reconstructed from all available data in the ALMA archive. Rotational lines of CO from J=1–0 up to J=4–3 and HCO⁺ from J=1–0 and J=4–3 were fitted for the gas density, temperature and pressure derivation in a spatially resolved, pixel-by-pixel manner (using the spectrally collapsed images of Fig. 1). The resulting pressure maps depend little on the mechanism of energy deposition into the dense ISM, whether it be mechanical heating or cosmic ray (CR) heating or both. This result owes to similar kinetic temperature, $T_{\text{kin}}$, times H₂ volume density, $n_{\text{H}_2}$, products for both mechanisms.

Mechanical heating is a widely examined mechanism that encapsulates a series of physical processes: radiation pressure drives winds of tenuous ionized gas, which then interact with molecular clouds via ram pressure, hydrodynamic and thermal instabilities that assist ram pressure by cloud evaporation, and cloud uplifting done by the work of adiabatically expanding plasma bubbles¹⁴⁻¹⁶. These mechanisms also apply to any wind initiated by a radio jet¹⁷. Another mechanism present in any wind, and in particular in a jet-driven wind, are CRs. CRs are less frequently discussed in the literature as a mechanism of gas heating, acceleration and wind launching¹⁸, even though they can be energetically important: in the solar neighbourhood, their energy density is comparable with the thermal energy density of the ISM¹⁹. CRs penetrate clouds, directly interacting with the dense gas at great cloud depths via ionizing collisions that can even change the gas chemistry²⁰⁻²². Hydrodynamical simulations indicate that CR-driven winds around supernovae are cooler and denser than thermally driven winds²³. The analysis that follows shows that both mechanical and CR heating are capable of

<footnote n="1">¹Section of Astrophysics, Astronomy, & Mechanics, Department of Physics, National and Kapodistrian University of Athens, Athens, Greece. ²Institute for Astronomy, Astrophysics, Space Applications, and Remote Sensing, National Observatory of Athens, Penteli, Greece. ³Max-Planck-Institut für Radioastronomie, Bonn, Germany. ⁴I. Physikalisches Institut, Universität zu Köln, Cologne, Germany. ⁵Observatoire de Paris, LERMA, Collège de France, CNRS, PSL University, Sorbonne University, Paris, France. ⁶Centro de Estudios de Física del Cosmos de Aragón, Teruel, Spain. ⁸e-mail: kdasyra@phys.uoa.gr</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | Observational maps used for the spatially resolved SLED fitting. The top row shows line fluxes (left to right): $F_{CO(1-0)}$, $F_{CO(2-1)}$, $F_{CO(3-2)}$ and $F_{CO(4-3)}$. Only pixels with $S/N > 3$ for CO (1-0), (2-1) and (3-2) are shown in this and in all following relevant figures. For CO(4-3), the dotted line traces the $S/N = 3$ contour. Upper limits were used for pixels outside this dotted contour. Solid contours mark the jet trail, as traced mainly by the [O III] emission in a Hubble Space Telescope Wide Field Planetary Camera 2 narrow-band image. Crosses mark the position of the radio core and lobes from the ALMA continuum. The bottom row shows velocity (left) and velocity dispersion (middle) averaged for all CO lines, and HCO$^+$ (4-3)/(1-0) ratio (right). In the bottom right panel, white contours trace the (1-0) emission at $S/N = 5, 10$ and grey contours trace the (4-3) emission at $S/N = 5, 20$. In all panels, the coordinates are relative to the galaxy's centre in the sky.</figcaption></figure>

reproducing the observed CO and HCO$^+$ spectral line energy distribution (SLED) in IC 5063.

To model SLEDs, we ran radiative transfer (RT) models that were carried out in a self-consistent manner for the available gas heating sources and that only considered dynamically meaningful cloud states with respect to their volume density, column density and kinetic temperature combination (see Methods and Extended Data Fig. 1 for more details). For this purpose, we used the three-dimensional photodissociation region (3D-PDR) code (ref. 24), which performs thermal balance calculations for the input heating sources and solves the astrochemical network equations before carrying out RT and outputting line fluxes. We used two heating sources of galactic origin in all our models: stellar ultra-violet (UV) photons and stellar-evolution-related CRs. CRs are included in 3D-PDR via the ISM ionization rate that they cause, $\zeta_{CR}$. We also used UV photons related to the active galactic nucleus (AGN) and two heating sources related to the jet: more CRs and mechanical heating. The additional CR flux is a highly attractive yet little-explored candidate for the dense gas acceleration and excitation in radio galaxies. The mechanical heating is associated with a widely implemented feedback mode in hydrodynamical simulations of cosmological scales$^{25}$ or galaxy scales$^{26,27}$. Indeed, a contribution of all these mechanisms to the ISM heating rate, $\Gamma$, was indicated by our best-fit model.

In our best-fit model, UV radiation is mainly important for ambient regions far from the radio lobes (for example, in the northern spiral arm of Fig. 2). In these regions, UV heating rates of a few $10^{-22}$ erg s$^{-1}$ cm$^{-3}$ are achieved thanks to unattenuated UV fluxes ($f_{UV}$) that are 5–10 times greater than those in the solar neighbourhood$^{28}$. At the northwest lobe, $\Gamma_{CR}$ increases to $6 \times 10^{-21}$ erg s$^{-1}$ cm$^{-3}$ with $\zeta_{CR}$ reaching $2 \times 10^{-14}$ s$^{-1}$ thanks to CRs from the jet. Additionally, the heating rate due to CRs ($\Gamma_{CR}$) is mainly attributed to CRs of galactic origin and is only at a level of $\sim 10^{-24}$ erg s$^{-1}$ cm$^{-3}$, despite a typical ionization rate similar to that of the solar neighbourhood ($\zeta_{CR} \approx 5 \times 10^{-17}$ s$^{-1}$). However, the role of CRs becomes important at the radio lobes, where the most intense jet-dense ISM interaction takes place (Fig. 2). At the northwest lobe, $\Gamma_{CR}$ increases to

$6 \times 10^{-21}$ erg s$^{-1}$ cm$^{-3}$, with $\zeta_{CR}$ reaching $2 \times 10^{-14}$ s$^{-1}$, thanks to CRs from the jet. Still, mechanical heating is the dominant mechanism throughout the entire radio jet trail, with the mechanical heating rate $\Gamma_m$ reaching values as high as $2 \times 10^{-20}$ erg s$^{-1}$ cm$^{-3}$. Within one beam from the northwest radio lobe, mechanical heating contributes ~3 times more than CRs to the heating budget. Taking into account all mechanisms, $\Gamma$ overall increases by a factor of 200 from ambient regions to the jet vicinity. The extra heating results in an increase of the jet-impacted gas temperature: $T_{kin}$ increases from 20 K in ambient regions to 130 K at the northwest lobe (Fig. 2), with a corresponding increase of the CO(1-0) excitation temperature, $T_{ex}$, from 10 K to 110 K. The heated molecular gas is also denser (with $n_{H_2} \approx 20,000$ cm$^{-3}$ instead of $3,000$ cm$^{-3}$), either because it was denser before the jet passage or because its deeper and denser layers got heated or exposed after the jet passage. The detection of HCO$^+$ (4-3) in the same region$^{29}$, which requires densities $\gtrsim 10^7$ cm$^{-3}$ for collisional excitation, suggests that the jet heats the clouds all the way down to their cores. The line-of-sight-integrated, intrinsic $H_2$ column density of clouds, $N_{H_2}^{LOS}$, there exceeds $10^{23}$ cm$^{-2}$. It is above the threshold for collapse and star formation$^{30}$ and two orders of magnitude greater than that of the ambient ISM. Even the densest gas phase that could form stars thus gets heated by the jet.

While our best-fit model points to a joint contribution of CR and mechanical heating to the energy budget, either mechanism can reproduce the CO and HCO$^+$ SLEDs of jet-impacted regions alone, as demonstrated by two more models (explained in Methods and shown in Extended Data Fig. 1): one with $\zeta_{CR}$ free and $\Gamma_m = 0$ and another with $\Gamma_m$ free and Galactic $\zeta_{CR} = 10^{-16}$ s$^{-1}$. Each mechanism individually leads to a similar increase in $\Gamma$ (Extended Data Figs. 2 and 3) for a reasonable increase of $\Gamma_m$ (in agreement with the outflow kinetic luminosity) or of $\zeta_{CR}$ (in agreement with values in the Galactic Centre or in the shock fronts near supernovae remnants$^{31}$). This result underlines the importance of incorporating CRs into outflow studies.

Differences between the models also exist. The CO emission of CR-heated gas originates from cooler and denser regions, with $T_{ex} \approx 75–80$ K and $n_{H_2} \approx 10,000–16,000$ cm$^{-3}$, than that of mechanically
<figure><img src="image_4.png" /><figcaption>**Fig. 2 |** Spatially resolved CO and HCO$^+$ SLED fitting results. Mechanical heating rate, CR heating rate, kinetic temperature and H$_2$ volume density maps are shown from left to right. Masks, contours and crosses are as in Fig. 1.</figcaption></figure>

heated gas, with $T_{\text{ex}} \approx 90-130$ K and $n_{\text{H}_2} \approx 7,000-10,000 \text{ cm}^{-3}$.
Indicative errors on these quantities are on the order of 20–30%.
The conversion factor of the CO luminosity at the surface of clouds
to the H$_2$ mass in their interior, $X_{\text{CO}}$, differs too. For mechanical
heating, $X_{\text{CO}}$ drops from nearly Galactic values at ambient regions
to ~4–5 times lower values along the jet trail, where dispersed
accelerated gas exists$^{3,4,32}$. For CR heating, $X_{\text{CO}}$ is nearly Galactic
everywhere, as the emission from dense layers adds to the emis-
sion from dispersed layers. This result is of extreme importance
for the mass measurement of jet-impacted clouds. A promis-
ing diagnostic of the most appropriate model is the flux ratio of
CO(7–6) to CI(2–1); that is, of two lines that can be simultane-
ously observed near 809 GHz. In several regions along the jet trail,
CO(7–6) is brighter than CI(2–1) for mechanical heating but
CI(2–1) is brighter than CO(7–6) for CR heating (Methods and
Extended Data Fig. 4). This inversion owes to the ionization of He
by CRs. He$^+$ then engages in the reaction He$^+$ + CO → O + C$^+$ + He,
efficiently destroying CO molecules and enhancing atomic
C. Finally, C$^+$ easily captures an electron, becoming neutral$^{33}$. It
is noteworthy that, in general, X-rays could be leading to similar
chemistry and contributing to $\Gamma$ (ref. $^{34}$). In IC 5063, diffuse X-ray
emission linked both to the weak AGN and to the jets was recently
discovered$^{35}$, but it is unlikely to play a dominant role in the gas exci-
tation. A gas heating rate of $10^{-20} \text{ erg s}^{-1} \text{ cm}^{-3}$ could only be attained
at the northwest radio lobe if all clouds are near (for example, within
several parsecs of) a single X-ray source of integrated luminosity
of ~$10^{39} \text{ erg s}^{-1}$ (ref. $^{36}$). A more realistic spatial distribution of the
molecular gas, for example, including clouds hundreds of pc away
from the X-ray emitters, would considerably decrease $\Gamma_{\text{Xrays}}$, demot-
ing or eliminating the role of X-ray heating. Still, we note that even if
X-rays were fully capable of providing the needed heating rate, they
would lead to similar spatial energy distribution as CRs$^{34}$. For these
reasons, we do not implement them further but instead focus on the
role of the abundant jet-related CRs.

The internal pressure $P_i$ of a typical line-of-sight cloud is nearly
identical for all excitation models. This is a very powerful result for
studies of cloud stability, which owes to lower $T_{\text{kin}}$ but higher $n_{\text{H}_2}$
values when the gas is CR heated instead of mechanically heated.
In three dimensions, $P_i$ is $5n_{\text{H}_2}kT_{\text{kin}}/2 + 3n_{\text{H}_2}m_{\text{H}}U_{\text{turb}}^2$, where $k$ is
the Boltzmann constant, $m_{\text{H}}$ is the mass of a hydrogen atom, and
$U_{\text{turb}}$ is the cloud microturbulent velocity, 0.6 km s$^{-1}$. $P_i$ displays clear
spatial variations along the jet trail (Fig. 3). Unaffected gas clouds
at ambient locations have $P_i/k = 8 \times 10^5 \text{ K cm}^{-3}$. Heavily impacted
clouds at the radio lobes have $P_i/k = 7 \times 10^6 \text{ K cm}^{-3}$, which is similar
to the values of giant molecular clouds in the distant Universe$^{37}$, and
indicates different conditions for future star formation. Under the
assumption that all clouds had similar initial properties, the spatial
gradients indicate how the internal pressure of individual clouds
changed after the jet passage. Unaffected clouds kept memory of

the initial properties, while impacted clouds changed properties.
Indeed, the ratio of the total hydrogen nucleus column density ($N_{\text{tot}}$)
to $n_{\text{H}_2}$, which is a measure of radius, is two times larger for clouds
at ambient regions (0.8 pc) than for clouds at the northwest radio
lobe (0.4 pc).

Along with internal pressure gradients, external pressure gradients determine the hydrodynamic fate of molecular clouds. The external medium pressure $P_e$ is assumed to be that of the ionized gas seen in VLT MUSE data. We used classical optical line diagnostics, [S II] 6,716 Å/[S II] 6,732 Å and ([N II] 6,548 Å + [N II] 6,584 Å)/[N II] 5,758 Å (ref. 38), to determine the ionized gas density and temperature $T$, respectively, under the assumption of thermal equilibrium for the electronic level populations of atoms (Methods and Extended Data Fig. 5). We computed $P_e$ as (3/2)$n_{\text{H}}kT$, where $n_{\text{H}}$ is the total H-nucleus number density, and found that $P_e/k$ can reach $2 \times 10^7 \text{ K cm}^{-3}$ (Fig. 3), peaking near an ionized gas bifurcation point seen in [O III] (ref. 9). It indicates a severe jet-cloud collision there. The drop happens along the axis of a known biconical [Fe II] and H$_2$ outflow$^9$. Coaxial to this outflow, there are regions with low ionized gas density (Fig. 4) that can be explained as pathways that the jet opened through the ISM once scattered at the impact point. Indeed, faint 3 cm emission from the jet is seen parallel to these regions, at a direction perpendicular to its initial path$^{10}$. With the exception of these regions, the ionized and molecular media have otherwise surprisingly comparable pressures, with $P_e$ being on average two times $P_i$. Pressure measurements of the gas from VLT SINFONI data yielded a similar range of values at the radio-lobe vicinity (Methods and Extended Data Fig. 6). From [Fe II] 1.533 and 1.644 µm, we deduced shocked ionized medium densities of 300–1,100 cm$^{-3}$ (refs. 38,39). From H$_2$ (0–0) S(1) and S(3), we obtained $T_{\text{kin}}$ of 1,500–2,000 K for the thermalized molecular gas, noting that this range could be somewhat elevated if there is a non-negligible contribution of fluorescence to the upper quantum level populations of many molecules. Combining these numbers, even though they reflect properties of different ISM phases that may not coincide or be in equilibrium, yields an estimated shocked gas pressure $P_{\text{shocked}}/k \approx 10^6-10^7 \text{ K cm}^{-3}$.

With $P_e$ frequently exceeding $P_i$ (Fig. 3), we are observing the
expansion of an over-pressurized ionized medium bubble, the
so-called cocoon of the jet, that engulfs and compresses clouds as
it propagates in all directions. Otherwise, $P_i$ would exceed $P_e$ and
the stability of molecular clouds would be achieved thanks to self
gravity, as in the normal ISM. Whether the compression leads to
a temporary enhancement of star formation is unknown, depend-
ing on other factors such as microturbulence. The over-pressurized
region would be more extended if the microturbulence of the ion-
ized gas exceeded that of the molecular gas. Simultaneously, under-
pressurized regions exist along the pathways that the jet opened
up. In clouds with kinetic energy exceeding the virial value, which
<figure><img src="image_5.png" /><figcaption>**Fig. 3** | Pressure maps. Internal pressure (left) of typical line-of-sight molecular clouds, deduced from CO and HCO+. Ionized medium pressure (middle), considered external to molecular clouds, from [S II] and [N II] (Fig. 4). Ratio of these two pressures (right), showing that due to the jet, $P_e$ can be either higher or lower than $P_i$. Masks, contours and crosses are as in Fig. 1.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 4** | Optical gas properties derived from VLT MUSE data. Temperature map based on the [N II] flux ratio (left). Density map based on the [S II] flux ratio (right). Masks and crosses are as in Fig. 1.</figcaption></figure>

happens in most regions with $P_i/P_e \gtrsim 2$ (Fig. 3), evaporation of cloud layers that become gravitationally unbound is possible. For the region ~350 pc west of the nucleus, the evaporation will be possible even if the ionized gas $v_{\text{turb}}$, the one-dimensional microturbulent velocity of the cloud, is as high as $10 \text{ km s}^{-1}$. It could be followed by mass loading of the outflow thanks to the aid of ram pressure from co-spatial ionized gas winds. Direct proof that cloud evaporation has occurred is the decrease of the optical depth along the jet trail$^*$ (Methods). Conditions of both suppression and enhancement of star formation are thus simultaneously plausible.

Following the modelling of the CO and HCO+ emission in ALMA data of IC 5063, we found that the jet passage leads to an increase of the internal pressure of molecular clouds (from $P_i/k = 8 \times 10^5 \text{ K cm}^{-3}$ to $P_i/k = 7 \times 10^6 \text{ K cm}^{-3}$). From the fluxes of [S II] and [N II] lines in VLT MUSE data, we found that the external pressure of molecular clouds increases enough in several regions to exceed $P_i$. This result led us to conclude that we are observing the expansion of an ionized over-pressurized cocoon that compresses molecular clouds and that could lead to their collapse. Contrarily, some jet-impacted clouds near pathways that the jet cleared have decreased $P_e$ compared with $P_i$. They are likely to undergo evaporation of their outer layers and could mass load the outflow. The observed pressure changes thus led us to conclude that both star formation enhancement and suppression could simultaneously occur. These changes owe to higher temperatures and densities of the line-emitting cloud layers, which can be caused by mechanical and CR heating. This conclusion could concern all outflows: jet-driven, AGN-driven and starburst-driven, as CRs can be found

in all of them because of Fermi acceleration. Given that CR heating is probably associated with momentum exchange, CRs could also help the launching of outflows, potentially diminishing the role of the adiabatic expansion of radiation-heated plasma bubbles to momentum boosting of the gas. Even in cosmological simulations with jet-related feedback, the inclusion of CRs along with mechanical heating could be essential for the proper derivation of galaxy masses.

## Methods

**The ALMA CO emission maps.** To create two-dimensional maps of the CO spectral line energy distribution, we used all archived ALMA data of IC 5063. For the CO(1–0) and CO(3–2) cube reconstruction, we used the observations for the programmes 2015.1.00420.S (PI F.C.) and 2015.1.00467.S (PI R. Morganti) with total on-source integration times of 2.5 h. For CO(2–1), we merged the data of the programmes 2012.1.00435.S (PI R. Morganti) and 2016.1.01279.S (PI S. Carniani) with a total on-source integration time of 1.5 h. For the second programme, we used the central pointing of the mapping observations. For CO(4–3), only data from the 2015.1.00420.S (PI F.C.) programme were available, with an on-source integration time of 40 min.

The reduction was carried out with the Common Astronomy Software Applications (CASA). To calibrate the data of the individual programmes, we used the respective CASA version and functions that were used by the ALMA staff. For all following tasks, we used CASA version 5.4.0-70. We ran the function 'uvcontsub' for each measurement set, employing zero-order polynomials to model and subtract the continuum emission at frequencies free of line emission. For programmes with spectral coverage overlap, we used the functions 'mtransform' and 'msconcat' to merge the measurement sets once we regridded them to the common velocity grid and resolution. For the CO(2–1) data, we also ran the 'statwt' function to recompute the visibility weights of the 2012 data before concatenating them with the 2016 data, as needed for CASA versions before 4.2.2.
To transform the merged visibilities into images, we used the 'tclean' function, running it iteratively to clean the images from the synthesized beam sidelobes. The primary beam response correction was applied. To ensure that the extended flux is well recovered, in particular for the low-frequency lines, we opted to use a Briggs robustness parameter of 2, corresponding to natural weighing of the interferometric u-v plane. This choice mattered for the CO(1-0) flux, but made no difference for the measured CO(4-3) flux. The image reconstruction with natural weighing led to synthesized beams of 0.72" × 0.65", 0.63" × 0.59", 0.66" × 0.56" and 0.36" × 0.33" for CO (1-0), (2-1), (3-2) and (4-3), respectively. To reach a common resolution for all lines, we convolved the cubes to the CO(1-0) beam using the function 'ia.convolve2d'. For the adopted cosmology (Λ cold dark matter model, where Λ is the cosmological constant, with: current value of the Hubble constant, $H_0 = 70 \text{ km s}^{-1} \text{Mpc}^{-1}$; matter density parameter $\Omega_M = 0.3$ and dark energy density parameter $\Omega_\Lambda = 0.7$) and for the measured source redshift (0.01128), the beam corresponds to a scale of 165 × 150 pc². All cubes were built at a common spectral resolution channel of 20 km s⁻¹. In the common spatial and spectral resolution, the sensitivity levels were 0.25, 0.23, 0.39 and 1.8 mJy per beam for CO (1-0), (2-1), (3-2) and (4-3), respectively. For the inner 5 kpc (22"), the respective line fluxes were 30.6(±2.4) Jy km s⁻¹, 77.0(±4.5) Jy km s⁻¹, 107(±8.0) Jy km s⁻¹ and 140(±20) Jy km s⁻¹.

Observations of HCO⁺(1-0) and HCO⁺(4-3), with respective exposure times of 38 min and 82 min, were acquired for the same programmes and were used to constrain the CO SLED excitation modelling. The HCO⁺ data reduction was carried out similarly to that of the CO. The emission was more compact than CO(4-3) allowing us to adopt a robustness parameter of 1 for the image reconstruction. We then convolved the cubes to the resolution of CO(1-0). The measured line fluxes were 5.6(±0.8) Jy km s⁻¹ for HCO⁺(1-0), 3.2(±0.2) Jy km s⁻¹ for HCO⁺(4-3). The noise, at the common resolution, was 0.4 and 0.5 mJy per beam, respectively.

We created moment maps of the cubes using the weighted-flux algorithm that is described in the CASA 'immoments' function, with different integration limits around each spectral line per pixel. While immoments uses a fixed velocity range for all pixels, our routine determines the integration limits in each pixel by performing an initial Gaussian line profile fit. This way, the resulting moment maps suffer less from noise contamination in regions with low dispersion (for example, in the spiral arms) and include high-velocity gas emission in regions with outflows. The derived line fluxes, recession velocity and velocity dispersion maps, averaged over all CO lines, are shown in Fig. 1. Information on the HCO⁺ line extent and ratio is also given in Fig. 1. All maps include both the disk and outflow emission, as a global solution is sought: the emission of the quiescent gas in the disk cannot be distinguished from that of the jet-heated or jet-accelerated gas perpendicular to the line of sight. Before fitting the data, we divided the observed line fluxes by their full-width at half-maximum in each pixel, so that all computations are performed for individual molecular clouds of intrinsic line width of 1 km s⁻¹ in the line of sight.

**Radiative transfer modelling of the CO SLED.** To model any SLED, RT for spectral lines is needed. The most traditional and simple RT approach is the modelling of one molecular species at a time, using fixed abundance ratios between species and assuming that the gas is isothermal (as, for example, in the one dimensional non-equilibrium radiative transfer code RADEX). Modelling of molecular line fluxes with RADEX⁴⁰ has been presented before for IC 5063³⁹, indicating that the plausible range of $T_{\text{kin}}$ and $n_{\text{H}_2}$ increases with proximity to the radio lobes. However, not all pre-selected gas column and volume density combinations correspond to dynamically meaningful cloud states (for example, some solutions correspond to collapsed or evaporated states). Likewise, pre-selected temperatures and abundances do not necessarily match their appropriate values for the available gas heating sources (and heating rate). Moreover, the isothermal assumption is unrealistic.

To improve the RT handling we opted to use a code that implements heating and cooling calculations until thermal balance is self-consistently reached, together with astrochemical calculations that provide the appropriate abundance ratios for the enhancement or destruction of species at each temperature. For this purpose, we used the code 3D-PDR²⁴. The astrochemical calculations used for the modelling of the observed CO and HCO⁺ line intensities depend on the implemented heating sources. We included far-ultraviolet radiation heating, CR heating and mechanical (shock-related) heating. Considering various cooling processes, 3D-PDR converges once the total heating and cooling functions are balanced. It uses a subset of the UMIST Database For Astrochemistry 2012 chemical network¹¹, comprising 33 species and 330 reactions. The metallicity is assumed solar, the dust-to-gas ratio is 100 and the initial gas-phase chemical abundances, normalized to hydrogen, are $[\text{C}^+/\text{H}] = 1 \times 10^{-4}$ (ref. ⁴²), $[\text{O}/\text{H}] = 3 \times 10^{-4}$ (ref. ⁴³) and $[\text{He}/\text{H}] = 0.1$. The thermochemical grid is built upon one-dimensional uniform-density distributions. Taking the level populations provided by the astrochemical calculations, the code performs RT calculations to compute the coolants' emission. It is assumed that the cloud is embedded in the cosmic microwave background radiation field and that the molecules get excited by collisions with up to seven partners (for example, hydrogen, helium). The collision rates are taken from the Leiden Atomic and Molecular Database (LAMBDA) database¹⁴. The calculations are performed under

the large velocity gradient approximation for an expanding sphere. The code provides the kinetic temperature $T_{\text{kin}}$, the excitation temperature $T_{\text{ex}}$, the radiation temperature $T_{\text{R}}$, the abundance and the volume density of each requested coolant. The coolants that we requested were several CO and HCO⁺ lines. Simultaneous solutions of the CO SLED and HCO⁺ (4-3)/(1-0) ratios were found.

To construct the grid of viable models, we set the input parameters that relate to the gas excitation as follows. We used Galaxy Evolution Explorer (GALEX) near- and far-ultraviolet flux measurements, $1.2 \times 10^{-4}$ and $5.2 \times 10^{-5}$ Jy, respectively, to deduce the power irradiated in the 6-13.6 eV energy range⁴⁵. Based on the GALEX images, about 5% of the flux is within ~500 pc from the nucleus of IC 5063. This leads to an average observed, attenuated ultraviolet flux $f_{\text{UV}}$ of order 1 Draine field (i.e., $0.936 \times 10^{-3}$ erg cm⁻² s⁻¹) at that distance. We thus assigned the intrinsic, unattenuated $f_{\text{UV}}$ values of 1, 5 or 10 Draine to allow for extinction, noting that this parameter mainly plays a role for the CO SLED for $n_{\text{H}_2} \le 10^3$ cm⁻³ (ref. ⁴⁶). For higher densities, CRs, which exist in the ISM and are abundant in jets, can play a considerable role instead. We allowed for CR ionization rates in the range $10^{-17} \le \zeta_{\text{CR}} \le 10^{-12}$ s⁻¹ (with a logarithmic step of 0.1 dex). This range covers $\zeta_{\text{CR}}$ values as low as those observed in the solar neighbourhood¹⁷ and as high as those estimated in extreme extragalactic objects such as the type-1 Seyfert Mrk231 (ref. ⁴⁹). Typical $\zeta_{\text{CR}}$ values are $10^{-16}$ s⁻¹ for Milky Way regions far from the Galactic Centre⁴⁰⁻⁵³ and $10^{-14}$ s⁻¹ near the Galactic Centre⁵⁴,⁵⁵. To encapsulate other mechanisms of heat deposition, including shocks, we inserted a mechanical heating function with free but constant rate $\Gamma_m$ (ref. ⁵⁶). To calculate the $\Gamma_m$ range, we used the observed molecular outflow kinetic luminosity $1/2 M_{\text{outflow}} V_{\text{outflow}}^3 / d$ (as in energy-conserving starburst-driven outflows)⁵⁷, distributing it in a homogeneous sphere with radius equal to the travelled distance $d$. For the outflow properties, we adopted a mass $M_{\text{outflow}} = 2 \times 10^6 M_{\odot}$ and a velocity $V_{\text{outflow}} = 300 \text{ km s}^{-1}$ (refs. ⁴,⁹). For $d$, we assumed it to be in the range 60–500 pc. The latter value, 500 pc, is similar to the distance between the northwest radio lobe and the nucleus. It is meaningful if the outflow started spreading from the nucleus. Alternatively, the outflow could have formed locally from the collision of the jet with a cloud ~500 pc away from the nucleus. For this scenario, we use $d=60$ pc, which is the extent of the jet–cloud interaction region as deduced from the spatial offset between the CO(1–0) peak and the CO(4–3) peak in the corresponding maps. This selection makes sense if CO(4–3) traces the densest sites of ongoing jet–cloud interaction while CO(1–0) traces sites closer to the nucleus, containing clouds rarefied by previous interactions with the jet. These assumptions lead to a $\Gamma_m$ range of $10^{-24}-10^{-20}$ erg s⁻¹ cm⁻³, which we spanned with a logarithmic step of 0.1 dex. To distinguish between the jet-related heating mechanisms, we created two additional grids: one implementing only mechanical heating ($\Gamma_m$ variable, $\zeta_{\text{CR}} = 10^{-16}$ s⁻¹, as in our Galaxy) and another implementing only CR enhancement ($\zeta_{\text{CR}}$ variable, $\Gamma_m = 0$).

We set the input parameters that relate to the spatial distribution and intrinsic properties of the clouds as follows. To ensure that the modelled CO lines met the observed ~1 km s⁻¹ line width, the gas internal microturbulent velocity was set to 0.6 km s⁻¹ (ref. ⁴⁰). The total H-nucleus number density, $n_{\text{H}_2}$, extended from $10^2$ to $10^4$ cm⁻³ with a step of 0.05 dex. The visual extinction $A_{\text{V}}$ which relates as $6.3 \times 10^{-23} \text{N}_{\text{H}_2} \text{mag cm}^2$ to the H-nucleus column density $N_{\text{H}_2}$ (refs. ⁵⁸,⁵⁹), had a minimum value of $10^{-3}$ mag, a maximum value of 10 mag and a step of 0.02 dex. Cloud dynamics further constrain the viable parameter range combinations of any grid. An $A_{\text{V}}-n_{\text{H}_2}$ relation that was established by hydrodynamical simulations of molecular clouds at parsec to kiloparsec scales⁶⁰⁻⁶³ has been used in the past to indicate parameter combinations that probe stable clouds⁶⁰. This relation is plotted with a cyan line in Extended Data Fig. 1 for $N_{\text{H}_2}/N_{\text{H}_1} = 2$ and typical $N_{\text{H}_2}/N_{\text{CO}}$ abundances⁶⁴, where $N_{\text{CO}}$ is the column density of CO molecules. However, it may not hold for jet-impacted clouds, which may or may not be in dynamical equilibrium. Instead, only grid limits ought to be set to reflect dynamical constraints. For this purpose, we provide an analytic expression for the stability of a molecular cloud as a function of three basic parameters of any RT code: the cloud column density, volume density and kinetic temperature. For H₂ molecules, the expression is

$$N_{\text{H}_2} = \frac{1}{m_{\text{H}_2}} \sqrt{\frac{20}{9\pi} \frac{n_{\text{H}_2}}{G} (5kT_{\text{kin}} + 3m_{\text{H}_2} v_{\text{turb}}^2)}, \quad (1)$$

where G is the gravitational constant, $v_{\text{turb}}$ is the one-dimensional microturbulent velocity of a cloud and $m_{\text{H}_2}$ is the mass of an H₂ molecule. The existence of this three-dimensional plane can be derived in a Jeans-like computation, starting from the assumption that a cloud is a virialized sphere of uniform density, supported against gravity by thermal and turbulent motions. From the virial theorem, the cloud virial mass is $5R (5kT + 3m_{\text{H}_2} v_{\text{turb}}^2)/(3Gm_{\text{H}_2})$. Equation (1) is obtained by equating the virial mass with $\pi R^2 N_{\text{H}_2} m_{\text{H}_2}$, where the cloud radius R is $3N_{\text{H}_2}/4n_{\text{H}_2}$, and solving for $N_{\text{H}_2}$. In this computation, we assumed that the cloud is fully molecular, that its mass is dominated by H₂ and that its heat capacity ratio is 7/5, appropriate for linear diatomic molecules. Limiting $T_{\text{kin}}$ and $v_{\text{turb}}$ reveals the viable parameter range in the column density–volume density plane. One limit is at $T_{\text{kin}} = 3$ K (the cosmic microwave background temperature) and $v_{\text{turb}} = 0$ km s⁻¹. The other limit is at $T_{\text{kin}} = 200$ K (that is, the maximum temperature indicated by an analytic solution for the CO(4-3)/CO(2-1) ratio in optically thin conditions and local thermodynamical equilibrium⁵¹) and $v_{\text{turb}} = 1$ km s⁻¹. The viable subgrids
(Extended Data Fig. 1) thus contain out-of-dynamical-equilibrium solutions, such as collapsing cores and expanding cloud shells, while encompassing the $A_V-n_{HII}$ relation.

To select the optimal model, we performed an error-weighed $\chi^2$ minimization between observed and modelled CO line fluxes in logarithmic space⁶⁵, which yields better results than a minimization in linear space for fluxes of different orders of magnitude. Before the minimization, we divided the observed fluxes by a beam filling factor $b_{\text{th}}$, which parameterizes the fraction of a telescope beam occupied by clouds, that converts observed beam-averaged brightness temperatures to intrinsic brightness temperatures and that is common for all CO lines. We allowed $b_{\text{th}}$ to vary between $1 \times 10^{-3}$ and 0.1, with a step of 0.5 dex. We ran the minimization algorithm for all CO spectral lines simultaneously requesting that the solution also satisfies the HCO⁺ (4-3)/(1-0) flux ratio, when available. Owing to the use of a line ratio, the HCO⁺ $b_{\text{th}}$ is irrelevant for our analysis. The fitting was performed for all pixels with signal-to-noise ratio (S/N) exceeding 3 in CO (1-0), (2-1) and (3-2). A common mask was created for these pixels. In CO(4-3), the detection is confined close to the jet trail (Fig. 1); to enable the study of quiescent regions, we used flux upper limits for all pixels between the common mask and the CO(4-3) S/N=3 contours (Fig. 1). The best-fit solutions for characteristic regions are shown in Extended Data Fig. 1.

For comparison with the results of the variable $\Gamma_m$ model (Extended Data Fig. 2), the results of the variable $\zeta_{CR}$ model are shown in Extended Data Fig. 3. For both models, the results are nearly identical for the pressure and similar for most other parameters, with CRs leading to emission from somewhat cooler and denser ISM. Ambient clouds in the spiral arm are the closest to their theoretical virial state, while jet-impacted clouds at the radio lobes are the furthest away from their virial state. The change in the cloud surface is reflected in the $X_{CO}$ maps, showing whether we see more or less dissipated cloud layers. For CRs, cold, dense gas layers and warm, shallow gas layers are likely to both contribute to the emission, leaving no impact to the overall $X_{CO}$. In all models, the typical values of the beam filling factor were between 0.01 and 0.1, with most clouds being (reasonably) located where spiral arms meet the central disk. Predictions of the CO(7-6)/CI(2-1) line ratio that can serve as diagnostic of the most appropriate model are shown in Extended Data Fig. 4.

**Optical and near-infrared data.** To determine the properties of the diffuse gas layers surrounding the CO-emitting dense molecular clouds, we analysed VLT MUSE and SINFONI data cubes in optical and near-infrared (NIR) wavelengths that are publicly available in the European Southern Observatory (ESO) archive. The MUSE data, taken for the programme 60.A-9339, were processed by the ESO MUSE pipeline version 1.6.1. The cube that we obtained is nearly identical to that produced when manually running the ESO Reflex routines⁶⁶. We fitted the profile of each line in the MUSE data with two Gaussian functions, one for the disk and one for the outflow, as required for the deblocking of blended lines. Therefore, four Gaussians were used for the fitting and proper flux measurement of [SII] 6,716 Å and 6,731 Å, which overlap in regions with outflows. An extinction correction was also applied for lines far in wavelength used in some line ratios; for example, for [NII] at 5,755 Å and for the [NII] doublet at 6,548 Å and 6,584 Å. The correction was derived from the Balmer ratio, assuming an intrinsic Hα/Hβ decrement of 3.1 (ref. ⁶⁷), a foreground extinction and a Milky Way extinction curve⁶⁸. Six Gaussians were, in total, used for the deblocking of the Hα and above-mentioned [NII] doublet complex, which was fitted simultaneously with (another two Gaussians for) Hβ. Similar kinematics were assumed for [NII] 6,548 Å and 6,584 Å on one hand and Hα and Hβ on the other. Both the disk and the outflow of these same-species lines were assumed to have the same central velocity and velocity dispersion. Before fitting all gas emission lines, we created a template of the typical spectrum emitted by all stellar populations in the centre of IC 5063 from gas-free regions of the cube (from 13-32" away from the jet trail, parallel to it). The template was scaled and removed from each pixel's spectrum to ensure that stellar absorption does not alter the emission line flux measurements.

In the SINFONI data, which were acquired for the programmes 075.B-0348A and 089.B-0971A and reduced according to the recipes in their initial publication⁷, the smaller field of view of the instrument did not allow for a similar removal of the stellar absorption. Instead, we removed the CO(7-4) absorption shortward of [FeII] 1.644 µm using a rescaled CO(6-3) bandhead profile. The use of CO(2-0) instead of CO(6-3) yielded undistinguishable results. Similarly, the [FeII] 1.533 µm line was contaminated by weak Brackett 18-4 emission. To remove it, we subtracted the Brackett 13-4 line profile from the data cube after scaling its peak to the Brackett 18-4 peak.

In the optical regime, the flux ratio [SII] 6,716 Å/[SII] 6,731 Å is a classical density probe, involving transitions of comparable energy that result in the same electronic quantum level⁶⁷. The code PyNeb⁶⁸ was used for the computation of this ratio as a function of density with updated atomic coefficients (Extended Data Fig. 5), providing results similar to those of other recent codes⁶⁹. The debtalled [SII] line ratio map from the MUSE data was used to produce the density map shown in Fig. 4, under the assumption that the optically emitting ionized gas is at 10,000 K. The results are, nonetheless, little dependent on this temperature. Temperatures were measured using the auroral [NII] line at 5,755 Å and the [NII] doublet at 6,548 Å, 6,584 Å, since the brighter oxygen lines cannot be used instead (as [OIII]).

at 4,363 Å is out of the MUSE spectral range). Again, the temperature was derived using PyNeb, from the [NII] 6,548 Å + [NII] 6,584 Å/[NII] 5,755 Å ratio and assuming an electron density $n_e = 100 \text{ cm}^{-3}$. Expected variations with density are negligible as shown in Extended Data Fig. 5. The results for the temperature are also shown in Fig. 4.

In the NIR, the H₂ rovibrational lines (1-0) S(1) at 2.1218 µm and (1-0) S(3) at 1.9576 µm were used for the analytic determination of $T_{ex}$ of the warmest molecular gas layers, under the assumption that this gas is optically thin and in local thermodynamic equilibrium⁹. $T_{\text{kin}}$ is then equal to $T_{ex}$, which is found by comparing the measured with the theoretical luminosity ratio of the two lines. The theoretical luminosity is the product of the Boltzmann-distribution-related number of photons emitted during a transition times the photon energy times the transition Einstein coefficient. The $T_{ex}$ map for the H₂-emitting gas indicates a temperature range of 1,500–4,000 K (Extended Data Fig. 6), with fluorescence adding uncertainty to $T_{ex}$ because of its potential contribution to the emission. Still, the average $T_{ex}$ lies in the lower part of this range in most regions⁹. Higher values are mainly seen in regions with high fractions of outflowing gas, for example, in the biconical H₂ outflow¹⁰. A density probe in the NIR that is effective for gas of several thousand Kelvin is the flux ratio [FeII] 1.533 µm/[FeII] 1.644 µm. These lines probe transitions from the same upper quantum level to lower quantum levels of comparable energy, so that the main factor that determines the number of ions that de-excite radiatively relates to the critical density. Their critical densities differ enough that their ratio can be used as a density measure in the range $10^2-10^6 \text{ cm}^{-3}$ (ref. ⁹⁹). We obtained an updated analytic function that connects the [FeII] 1.533 µm/[FeII] 1.644 µm flux ratio with the ISM density by running PyNeb (Extended Data Fig. 5). Our calculations indicated ionized gas densities in the range 300–30,000 cm⁻³ (Extended Data Fig. 6), for a typical temperature of 5,000 K. This temperature is well chosen for two reasons. First, the [FeII] ratio was limited to a maximum of 0.3, not requiring higher temperatures to be produced (Extended Data Fig. 5). Second, the [FeII] emission is co-spatial with the few-thousand-Kelvin H₂ gas in IC 5063¹⁰. The derived density appears to be lowest in the jet-impacted regions and to increase with distance from the jet trail (Extended Data Fig. 6). This overdensity, which needs to be confirmed with higher S/N infrared data, could be a nice tracer of the border of the cocoon, which shows a mild CO outflow¹¹.

## Data availability

All ALMA raw data can be retrieved from the ALMA archive using the project identifiers 2015.1.00420.S, 2015.1.00467.S, 2012.1.00435.S and 2016.1.01279.S. The line images and all other data presented in Fig. 1 can be obtained in .fits format from Fighshare at https://doi.org/10.6084/m9.fighshare.19742125. Fully reduced data cubes can be provided from the corresponding author upon reasonable request. The pipeline-processed MUSE data can be directly downloaded from the ESO data portal using the project identifier 60.A-9339.

## Code availability

The RT code 3D-PDR is a publicly available radiative transfer code that can be obtained at https://uclchem.github.io/3dpdr.

Received: 7 January 2022; Accepted: 27 May 2022; Published online: 21 July 2022

## References

1. Fluetsch, A. et al. Cold molecular outflows in the local Universe. *Mon. Not. R. Astron. Soc.* **483**, 4586–4614 (2019).

2. Veilleux, S., Maiolino, R., Bolatto, A. & Aalto, S. Cool outflows in galaxies and their implications. *Astron. Astrophys.* **28**, 2 (2020).

3. Bolatto, A. D., Wolfire, M. & Leroy, A. K. The CO-to-H₂ conversion factor. *Annu. Rev. Astron. Astrophys.* **51**, 207–268 (2013).

4. Dasyra, K. M. et al. ALMA reveals optically thin, highly excited CO gas in the jet-driven winds of the galaxy IC 5063. *Astron. Astrophys.* **595**, L7 (2016).

5. Aalto, S. et al. A precessing molecular jet signaling an obscured, growing supermassive black hole in NGC 1377? *Astron. Astrophys.* **590**, 73 (2016).

6. Croft, S. et al. Minkowski's object: a starburst triggered by a radio jet, revisited. *Astrophys. J.* **647**, 1040–1055 (2006).

7. Crockett, R. M. et al. Triggered star formation in the inner filament of Centaurus A. *Mon. Not. R. Astron. Soc.* **421**, 1603–1623 (2012).

8. Maiolino, R. et al. Star formation inside a galactic outflow. *Nature* **544**, 202–206 (2017).

9. Dasyra, K. M. et al. A radio jet drives a molecular and atomic gas outflow in multiple regions within one square kiloparsec of the nucleus of the nearby galaxy IC5063. *Astrophys. J.* **815**, 34 (2015).

10. Oosterloo, T. et al. A strong jet-cloud interaction in the Seyfert galaxy IC 5063: VLBI observations. *Astrophys. J.* **119**, 2085–2091 (2000).

11. Morganti, R., Oosterloo, T., Oonk, R., Frieswijk, W. & Tadhunter, C. The fast molecular outflow in the Seyfert galaxy IC 5063 as seen by ALMA. *Astron. Astrophys.* **580**, A1 (2015).
12. Fabian, A. C. et al. A very deep Chandra observation of the Perseus cluster: shocks, ripples and conduction. *Mon. Not. R. Astron. Soc.* **366**, 417–428 (2006).

13. Ponti, G. et al. An X-ray chimney extending hundreds of parsecs above and below the Galactic Centre. *Nature* **567**, 347–350 (2019).

14. Hopkins, P. F. & Elvis, M. Quasar feedback: more bang for your buck. *Mon. Not. R. Astron. Soc.* **401**, 7–14 (2010).

15. Zubovas, K. & King, A. Galaxy-wide outflows: cold gas and star formation at high speeds. *Mon. Not. R. Astron. Soc.* **439**, 400–406 (2014).

16. Dannen, R., Progra, D., Waters, T. & Dyda, S. Clumpy AGN outflows due to thermal instability. *Astrophys. J.* **893**, L34 (2020).

17. Wagner, A. Y., Bicknell, G. V., Umura, M., Sutherland, R. S. & Silk, J. Galaxy-scale AGN feedback—theory. *Astron. Nachr.* **337**, 167–174 (2016).

18. Ruszkowski, M., Karen Yang, H.-Y. & Zweibel, E. Global simulations of galactic winds including cosmic ray streaming. *Astrophys. J.* **834**, 2 (2017).

19. Klessen, R. S. & Glover, S. C. O. In Star Formation in Galaxy Evolution: Connecting Numerical Models to Reality SaaS-Fee Advanced Course book series Vol. 43 (eds Revaz, Y. et al.) 85–249 (Springer, 2016).

20. Papadopoulos, P. P. A cosmic-ray-dominated interstellar medium in ultra luminous infrared galaxies: new initial conditions for star formation. *Astrophys. J.* **720**, 226–232 (2010).

21. Padovani, M., Ivlev, A. V., Galli, D. & Caselli, P. Cosmic-ray ionisation in circumstellar discs. *Astrophys. J.* **614**, 111 (2018).

22. Gaches, B. A. L., Offner, S. S. R. & Bisbas, T. G. The astrochemical impact of cosmic rays in protoclusters. I. Molecular cloud chemistry. *Astrophys. J.* **878**, 105 (2019).

23. Girichidis, P. et al. Launching cosmic-ray-driven outflows from the magnetized interstellar medium. *Astrophys. J.* **816**, L19 (2016).

24. Bisbas, T. G. et al. 3D-PDR: a new three-dimensional astrochemistry code for treating photodissociation regions. *Mon. Not. R. Astron. Soc.* **427**, 2100–2118 (2012).

25. Dubois, Y., Devriendt, J., Slyz, A. & Teyssier, R. Self-regulated growth of supermassive black holes by a dual jet-heating active galactic nucleus feedback mechanism: methods, tests and implications for cosmological simulations. *Mon. Not. R. Astron. Soc.* **420**, 2662–2683 (2012).

26. Gaibler, V., Khochfar, S., Krause, M. & Silk, J. Jet-induced star formation in gas-rich galaxies. *Mon. Not. R. Astron. Soc.* **425**, 438–449 (2012).

27. Wagner, A. Y., Bicknell, G. V. & Umura, M. Driving outflows with relativistic jets and the dependence of active galactic nucleus feedback efficiency on interstellar medium inhomogeneity. *Astrophys. J.* **757**, 136 (2012).

28. Draine, B. T. Photoelectric heating of interstellar gas. *Astrophys. J. Suppl. Ser.* **36**, 595–619 (1978).

29. Oosterloo, T. et al. Properties of the molecular gas in the fast outflow in the Seyfert galaxy IC 5063. *Astron. Astrophys.* **608**, 38 (2017).

30. Kennicutt, R. C. Star formation in galaxies along the Hubble sequence. *Annu. Rev. Astron. Astrophys.* **36**, 189–231 (1998).

31. Schuppan, F., Becker, J. K., Black, J. H. & Casanova, S. Cosmic-ray-induced ionization in molecular clouds adjacent to supernova remnants. Tracing the hadronic origin of GeV gamma radiation. *Astron. Astrophys.* **541**, A126 (2012).

32. Richings, A. J. & Faucher-Giguère, C.-A. The origin of fast molecular outflows in quasars: molecule formation in AGN-driven galactic winds. *Mon. Not. R. Astron. Soc.* **474**, 3673–3699 (2018).

33. Bisbas, T. G., Papadopoulos, P. P. & Viti, S. Effective destruction of CO by cosmic rays: implications for tracing H₂ gas in the Universe. *Astrophys. J.* **803**, 37 (2015).

34. Meijerink, R. & Spaans, M. Diagnostics of irradiated gas in galaxy nuclei. I. A FUV and X-ray dominated region code. *Astron. Astrophys.* **436**, 397–409 (2005).

35. Travascio, A. et al. AGN-host interaction in IC 5063. I. Large-scale X-ray morphology and spectral analysis. *Astrophys. J.* **921**, 129 (2021).

36. Maloney, P. R., Hollenbach, D. J. & Tielens, A. G. G. M. X-ray irradiated molecular gas. I. Physical processes and general results. *Astrophys. J.* **466**, 561–584 (1996).

37. Dessauges-Zavadsky, M. et al. Molecular clouds in the Cosmic Snake normal star-forming galaxy 8 billion years ago. *Nat. Astron.* **3**, 1115–1121 (2019).

38. Luridiana, V., Morisset, C. & Shaw, R. PyNeb: a new tool for analyzing emission lines I. Code description and validation of results. *Astron. Astrophys.* **573**, A42–56 (2015).

39. Pradhan, A. K. & Zhang, H. L. New excitation rates and line ratios for [Fe ii]. *Astrophys. J.* **409**, L77–79 (1993).

40. van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J. & van Dishoeck, E. F. A computer program for fast non-LTE analysis of interstellar line spectra. With diagnostic plots to interpret observed line intensity ratios. *Astron. Astrophys.* **468**, 627–635 (2007).

41. McElroy, D. et al. The UMIST Database For Astrochemistry 2012. *Astron. Astrophys.* **550**, A36 (2013).

42. Cardelli, J. A., Meyer, D. M., Jura, M. & Savage, B. D. The abundance of interstellar carbon. *Astrophys. J.* **467**, 334–340 (1996).

43. Cartledge, S. I. B. et al. The homogeneity of interstellar oxygen in the Galactic Disk. *Astrophys. J.* **613**, 1037–1048 (2004).

44. Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F. & Black, J. H. An atomic and molecular database for analysis of submillimetre line observations. *Astron. Astrophys.* **432**, 369–379 (2005).

45. Marino, A. et al. Nearby early-type galaxies with ionized gas: the UV emission from GALEX observations. *Mon. Not. R. Astron. Soc.* **411**, 311–331 (2011).

46. Bisbas, T. G., Schruba, A. & van Dishoeck, E. F. Simulating the atomic and molecular content of molecular clouds using probability distributions of physical parameters. *Mon. Not. R. Astron. Soc.* **485**, 3097–3111 (2019).

47. Cummings, A. C. et al. Galactic cosmic rays in the local interstellar medium: Voyager 1 observations and model results. *Astrophys. J.* **831**, 18–39 (2016).

48. González-Alfonso, E. et al. Outflowing OH⁺ in Markarian 231: the ionization rate of the molecular gas. *Astrophys. J.* **857**, 66 (2018).

49. van Dishoeck, E. F. & Black, J. H. Comprehensive models of diffuse interstellar clouds: physical conditions and molecular abundances. *Astrophys. J. Suppl. Ser.* **62**, 109–145 (1986).

50. McCall, B. J. et al. An enhanced cosmic-ray flux towards ζ Persei inferred from a laboratory study of the HⅠ⁻e⁻ recombination rate. *Nature* **422**, 500–502 (2003).

51. Dalgarno, A. The galactic cosmic ray ionization rate. *Proc. Natl Acad. Sci. USA* **103**, 12269–12273 (2006).

52. Neufeld, D. A. et al. Herschel/HIFI observations of interstellar OH⁺ and H₂O⁺ towards W49N: a probe of diffuse clouds with a small molecular fraction. *Astron. Astrophys.* **521**, L10 (2010).

53. Indriolo, N. et al. Herschel survey of galactic OH⁺, H₂O⁺, and H₂O⁺: probing the molecular hydrogen fraction and cosmic-ray ionization rate. *Astrophys. J.* **800**, 40 (2015).

54. Oka, T. et al. Hot and diffuse clouds near the Galactic Center probed by metastable HⅠ⁺. *Astrophys. J.* **632**, 882–893 (2005).

55. Goto, M. et al. Absorption line survey of HⅠ⁺ toward the Galactic Center sources. II. Eight infrared sources within 30 pc of the Galactic Center. *Astrophys. J.* **688**, 306–319 (2008).

56. Loenen, A. F., Spaans, M., Baan, W. & Meijerink, R. Mechanical feedback in the molecular ISM of luminous IR galaxies. *Astron. Astrophys.* **488**, L5–L8 (2008).

57. Heckman, T. M., Armus, L. & Miley, G. K. On the nature and implications of starburst-driven galactic superwinds. *Astrophys. J.* **74**, 833–868 (1990).

58. Weingartner, J. C. & Draine, B. T. Dust grain-size distributions and extinction in the Milky Way, Large Magellanic Cloud, and Small Magellanic Cloud. *Astrophys. J.* **548**, 296–309 (2001).

59. Röllig, M. et al. A photon dominated region code comparison study. *Astron. Astrophys.* **467**, 187–206 (2007).

60. Glover, S. C. O. et al. Modelling CO formation in the turbulent interstellar medium. *Mon. Not. R. Astron. Soc.* **404**, 2–29 (2010).

61. Van Loo, S., Butler, M. J. & Tan, J. C. Kiloparsec-scale simulations of star formation in disk galaxies. I. The unmagnetized and zero-feedback limit. *Astrophys. J.* **764**, 36 (2013).

62. Safranek-Shrader, C. et al. Chemistry and radiative shielding in star-forming galactic discs. *Mon. Not. R. Astron. Soc.* **465**, 885–905 (2017).

63. Seifried, D. et al. SILCC-zoom: the dynamic and chemical evolution of molecular clouds. *Mon. Not. R. Astron. Soc.* **472**, 4797–4818 (2017).

64. Smith, R. J., Glover, S. C. O., Clark, P. C., Klessen, R. S. & Springel, V. CO-dark gas and molecular filaments in Milky Way-type galaxies. *Mon. Not. R. Astron. Soc.* **441**, 1628–1645 (2014).

65. Viti, S. et al. Molecular line emission in NGC 1068 imaged with ALMA. II. The chemistry of the dense molecular gas. *Astron. Astrophys.* **570**, A28 (2014).

66. Venturi, G. et al. MAGNUM survey: compact jets causing large turmoil in galaxies. *Astron. Astrophys.* **648**, 17 (2021).

67. Osterbrock, D. E., Ferland, G. J. Astrophysics of Gaseous Nebulae and Active Galactic Nuclei (Univ. Science Books, 2006).

68. Cardelli, J. A., Clayton, G. C. & Mathis, J. S. The relationship between infrared, optical, and ultraviolet extinction. *Astrophys. J.* **345**, 245–256 (1989).

69. Kewley, L. J. et al. Theoretical ISM pressure and electron density diagnostics for local and high-redshift galaxies. *Astrophys. J.* **880**, 16–40 (2019).

## Acknowledgements

K.M.D., G.F.P. and T.G.B. acknowledge financial support by the Hellenic Foundation for Research and Innovation (HFRI), under the first call for the creation of research groups by postdoctoral researchers that was launched by the General Secretariat For Research and Innovation (project number 1882, PI K.M.D.). G.F.P. also acknowledges support for this research by the International Max-Planck Research School for Astronomy and Astrophysics at the University of Bonn and Cologne. T.G.B. further acknowledges support from Deutsche Forschungsgemeinschaft (DFG grant No. 424563772, T.G.B.). Financial support by the Spanish Ministry of Science and Innovation and the European Union—NextGenerationEU through the Recovery and Resilience Facility project ICTS-MRR-2021-03-CEFCA is acknowledged for J.A.F.-O. This paper makes use of the ALMA data with identifiers 2012.1.00435.S, 2015.1.00420.S, 2015.1.00467.S and 2016.1.01279.S. ALMA is a partnership of the ESO (representing its member states),
National Science Foundation (United States) and National Institute of Natural Sciences (NINS; Japan), together with the National Research Council (NRC; Canada), National Science Council (NSC) and Academia Sinica Institute of Astronomy and Astrophysics (ASIA; Taiwan), and Korea Astronomy and Space Science Institute (KASI; Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, Associated Universities Inc. (AUI)/National Radio Astronomy Observatory (NRAO) and National Astronomical Observatory of Japan (NAOJ).

## Author contributions

K.M.D. led the project, prepared the HFRI grant proposal, reduced the ALMA and SINFONI data and led their interpretation and write-up. G.F.P. wrote and ran the codes for the spatially resolved CO and HCO$^+$ SLED modelling and benchmarked 3D-PDR and RADEX. T.G.B., an author of 3D-PDR, provided the appropriate grids following astrochemical and radiative transfer calculations for the different gas-excitation sources. F.C. led the ALMA proposal providing one of the datasets and contributed to the data interpretation. J.A.F.-O. assisted with the fitting of the optical data.

## Competing interests

The authors financed by HFRI grant 1882 (K.M.D., G.F.P., T.G.B.) for this project declare the existence of a financial/non-financial competing interest (radiative transfer modelling of some common datasets) with ERC grant 320745.

## Additional information

Extended data is available for this paper at https://doi.org/10.1038/s41550-022-01725-9.

Correspondence and requests for materials should be addressed to Kalliopi M. Dasyra.

Peer review information *Nature Astronomy* thanks Roberto Maiolino and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Reprints and permissions information is available at www.nature.com/reprints.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2022
<figure><img src="image_7.png" /></figure>

**Extended Data Fig. 1** | SLED models. Top left: CO SLED results for characteristic regions. The SLEDs are shown as velocity-integrated line fluxes as a function of the rotational number $J$ of the upper energy level. Mechanical plus CR heating models are shown in purple, mechanical heating models are shown in orange, and CR heating models are shown in blue. Top right panels: $T_{\text{kin}} - n_{\text{H}_2}$ grid and results for each model. Bottom left: Corresponding $N_{\text{tot}}$ and $n_{\text{H}_2}$ solutions plotted over the input grid $N_{\text{tot}} - n_{\text{H}_2}$ projection. The cyan curve depicts a relation for stable clouds$^{46}$ suggested by hydrodynamic simulations (see Methods for more information). Shaded areas indicate the areas occupied by virialized clouds (see equation (1)): the blue area is for $T_{\text{kin}}=20$ K and $0 < v_{\text{turb}} < 1$ km s$^{-1}$, the purple area is for $T_{\text{kin}}=100$ K and $0 < v_{\text{turb}} < 1$ km s$^{-1}$. Orange, blue, and purple points correspond to mechanical, CR, and combined heating model results, respectively. Other bottom panels: $N_{\text{CO}} - n_{\text{H}_2}$ grid projection for all models. The initial grids and the grids trimmed for dynamical considerations are shown in light and dark grey, respectively in all pertinent panels.
<figure><img src="image_8.png" /></figure>

**Extended Data Fig. 2** | Spatially-resolved SLED fitting results for the mechanical heating model. In this model, Γ<sub>m</sub> is variable and ζ<sub>CR</sub>=10<sup>-16</sup> s<sup>-1</sup>. The quantities shown are the heating rate (total, mechanical, and FUV), the molecular gas kinetic temperature and volume density (top panels), and the cloud internal pressure, the CO column density along the line of sight, the CO beam filling factor, the CO(1-0) optical depth, and X<sub>CO</sub> (bottom panels).
<figure><img src="image_9.png" /><figcaption>**Extended Data Fig. 3** | Spatially-resolved SLED fitting results for the CR heating model. In this model, $\zeta_{CR}$ is variable and $\Gamma_m=0$. The quantities shown are similar to those in Extended Data Fig. 2.</figcaption></figure>
<figure><img src="image_10.png" /></figure>

**Extended Data Fig. 4 | Predictions of the CO(7-6)/CI(2-1) flux ratio.** This ratio can be used as a diagnostic between mechanical and CR heating.
<figure><img src="image_11.png" /></figure>

**Extended Data Fig. 5** | Optical line diagnostics of ionized gas conditions. [Fe II] and [S II] flux ratios used as density diagnostics (**left and center**) and [N II] flux ratio used as temperature diagnostic (**right**) as obtained by PyNeb³⁸.
<figure><img src="image_12.png" /></figure>

**Extended Data Fig. 6** | Gas properties derived from NIR VLT SINFONI data. **Left:** Temperature based on the H₂ (1-0) S(1) and S(3) line fluxes. The highest *T*<sub>kin</sub> values are close to regions of high turbulence perpendicular to the jet, associated with it<sup>9,66</sup>. **Right:** Density based on the [Fe II] 1.533 and 1.644 µm line fluxes.