Jet-induced molecular gas excitation and turbulence in the Teacup

In order to investigate the impact of radio jets on the interstellar medium (ISM) of galaxies hosting active galactic nuclei (AGN), we present subarcsecond-resolution Atacama Large Millimeter / submillimeter Array (ALMA) CO(2-1) and CO(3-2) observations of the Teacup galaxy. This is a nearby ( D L = 388Mpc) radio-quiet type-2 quasar (QSO2) with a compact radio jet ( P jet ≈ 10 43 ergs − 1 ) that subtends a small angle from the molecular gas disc. Enhanced emission line widths perpendicular to the jet orientation have been reported for several nearby AGN for the ionised gas. For the molecular gas in the Teacup, not only do we ﬁnd this enhancement in the velocity dispersion but also a higher brightness temperature ratio ( T 32 / T 21 ) perpendicular to the radio jet compared to the ratios found in the galaxy disc. Our results and the comparison with simulations suggest that the radio jet is compressing and accelerating the molecular gas, and driving a lateral outﬂow that shows enhanced velocity dispersion and higher gas excitation. These results provide further evidence that the coupling between the jet and the ISM is relevant to AGN feedback even in the case of radio-quiet galaxies.


Introduction
The idea that outflows are almost ubiquitous in galaxies hosting active galactic nuclei (AGN) is supported by observational efforts in the past decades, as well as by the development of subgrid physics in cosmological models (Cielo et al. 2018;Nelson et al. 2019). The inclusion of AGN feedback in the recipes of simulations is necessary to explain the observed bright end of the galaxy luminosity function, as otherwise galaxies would grow too big and massive (Croton et al. 2006;Dubois et al. 2016).
One of the potential drivers of multi-phase outflows, even in the case of radio-quiet AGN, is the jets launched by the supermassive black hole (SMBH). These jets can strongly impact (sub)kiloparsec scales by altering the properties of the interstellar medium (ISM) of the host galaxies. Observational evidence for outflows driven by jets has been reported in ionised (e.g., Heckman et al. 1984;Jarvis et al. 2019;Cazzoli et al. 2022;Girdhar et al. 2022;Speranza et al. 2022) and molecular gas (e.g., Morganti et al. 2015;Oosterloo et al. 2019;García-Burillo et al. 2019;Murthy et al. 2022). Supporting the observational evidence, hydrodynamic simulations (Wagner & Bicknell 2011;Wagner et al. 2012;Mukherjee et al. 2018a) are able to reproduce the jet-ISM interaction and the impact of feedback induced by relativistic AGN jets inside the central kiloparsec of gas-rich radio galaxies. The models trace the evolution of a relativistic jet propagating in a two-phase ISM, gradually dispersing the clouds through the effect of the ram pressure and internal energy of the non-thermal plasma, and creating a cocoon of shocked material that drives multi-phase outflows as the bubble expands.
In recent work, Ramos Almeida et al. (2022, hereafter RA22) reported Atacama Large Millimeter/submillimeter Array (ALMA) CO(2-1) observations at ∼0 . 2 (400 pc) resolution of seven radio-quiet type-2 quasars (QSO2s; i.e., obscured quasars) at redshifts z ∼ 0.1, with bolometric luminosities of L bol ≈ 10 45−46 erg s −1 . These QSO2s are part of the Quasar Feedback (QSOFEED) sample. Cold molecular outflows with intermediate properties between those of Seyfert galaxies and ultraluminous infrared galaxies (ULIRGs) were detected in the five QSO2s with CO(2-1) detections. Their molecular mass outflow rates are lower than those expected from their AGN luminosities (Fiore et al. 2017;Fluetsch et al. 2019), suggesting that other factors, such as the jet power, the spatial distribution of the dense gas, and the coupling between jets and the dense gas, might also be relevant.
Among the QSO2s studied in RA22, the Teacup (SDSS J143029.88+133912.0; J1430+1339) revealed a peculiar CO morphology and disturbed kinematics. This, in addition to the compact radio jet (extent of ∼0.8 kpc, PA = 60 • ) detected in Very Large Array (VLA) data at 0 . 25 resolution (Harrison et al. 2015), makes this object a promising target to study the jet-ISM interaction. The Teacup is hosted in a bulge-dominated galaxy showing signatures of a recent merger, including stellar shell-like features seen in the Hubble Space Telescope (HST) images (Keel et al. 2015). The nickname 'Teacup' stems from the 12 kpc ionised gas bubble (Keel et al. 2012) that is also detected in the radio with the VLA (Harrison et al. 2015), both shown in the left panel of Fig. 1. The [Oiii] kinematics reveal a ∼1 kpc scale nuclear outflow with maximum velocities of ∼−750 km s −1 (Villar Martín et al. 2014;Harrison et al. 2014). In the near infrared (NIR), blueshifted broad components with maximum velocities of −1100 km s −1 were measured from the Paα and [Si VI] emission lines, also extending ∼1 kpc and oriented at PA ∼ 70 • (Ramos Almeida et al. 2017). The projected orientation and extension of the nuclear ionised outflow almost coincide with the compact jet, suggesting that the latter could be driving the outflow (Harrison et al. 2015). Therefore, the Teacup is one of the few nearby luminous AGN whose outflows have been characterised in different gas phases (ionised, warm, and cold molecular). This, together with its compact radio jet, makes it an ideal target to advance our understanding of how compact jets impact the multi-phase ISM of radio-quiet AGN.

Observations and data reduction
We present ALMA observations of the CO(2-1) and CO(3-2) emission lines observed in bands 6 and 7. The CO(2-1) observations and data reduction are described in detail in RA22. The CO(2-1) datacube has a synthesised beam size of 0 . 21×0 . 18 and a root mean square (rms) noise of 0.39 mJy beam −1 per channel of 10 km s −1 . The CO(3-2) data (2016.1.01535.S; PI: G. Lansbury) were retrieved from the ALMA archive. The observations were made in April 2017 in the C40-3 configuration with an on-source integration time of 30.3 min and in a single pointing, covering a field of view (FoV) of 18 . The spectral window of 1.875 GHz total bandwidth was centred at the CO(3-2) line, with a channel spacing of 7.8 MHz, corresponding to ∼7.2 km s −1 after Hanning smoothing.
The CO(3-2) data were calibrated using the casa software version 4.7.2 (McMullin et al. 2007) in the pipeline mode. The sources J1550+0527, J1337−1257 and J1446+1721 were used as flux, bandpass, and phase calibrators, respectively. Once calibrated, the imaging and cleaning were performed with the task tclean. The spectral line map was obtained after subtracting the continuum in the uv-plane using the tasks uvcontsub and a zeroth-order polynomial in the channels free from emission lines. The CO(3-2) data cube was produced with a spectral resolution of 10 km s −1 (10.6 MHz) and using Briggs weighting mode and a robust parameter set to 0.5 in order to achieve the best compromise between resolution and sensitivity. Finally, the datacube was corrected for primary beam attenuation, resulting in a synthesised beam of 0 . 60×0 . 54 at PA = 59 • and a rms noise of 1.9 mJy beam −1 per 10 km s −1 channel. In order to study the line ratios and compare the CO(2-1) to the new CO(3-2) datacube, we regridded the CO(2-1) data to the same pixel scale of the CO(3-2) and then convolved it to the same common beam size of 0 . 60×0 . 54 (960 pc×860 pc).
The auxiliary radio data used here are the 6 GHz VLA observations of the Teacup in C-band in two different configurations: the high-resolution (0 . 25, HR) A-array and the low-resolution (1 . 0, LR) B-array maps. The data were presented and analysed in full by Jarvis et al. (2019).

Results and analysis
The CO(3-2) peak intensity map shown in the middle panel of Fig. 1 shows the double-peaked morphology first detected in CO(2-1) and discussed in RA22. The distance between the two peaks is ∼0 . 8 (1.3 kpc). The integrated intensity, intensityweighted velocity field, and velocity dispersion (σ) maps of the CO(3-2) emission are shown in Appendix A. Overall, the first-and second-moment maps (velocity field and σ) are similar to their CO(2-1) counterparts, although with a lower level of detail due to the coarser angular resolution. For example, the velocity field shows the disc-rotation pattern, but the disturbed kinematics seen in CO(2-1) are not detected in CO(3-2). The second-order moment map shows high values of σ, reaching ∼120 km s −1 with PA ∼ −30 • , orthogonal to the radio jet orientation (PA jet = 60 • ), similar to the σ map of the CO(2-1) data in RA22.
In order to investigate the possible influence of the jet on the kinematics of the molecular gas, we created position-velocity diagrams (PVDs) along and perpendicular to the direction of the jet axis. We extracted them using a slit width of 0 . 6, equivalent the common resolution of the CO(3-2) and convolved CO(2-1) cubes. The PVDs are shown in Fig. 2 and the position and orientation of the slits are indicated as shaded rectangles in the left panel of Fig. 3. A broad gradient of velocities is detected both along and perpendicular to the jet axis, spanning up to ±400 km s −1 for both the CO(2-1) and CO(3-2) emission. We detect gas at higher velocities in the PVDs of CO(3-2) than in those of CO(2-1), especially in the case of the blueshifted velocities. The PVDs perpendicular to the radio jet show an inverted S-shape in the central ∼2 (3.2 kpc), which includes the bulk of the emission. There is also a contribution from slow velocities (<200 km s −1 ) at r > 1 that resembles a rotation pattern, unlike the central part of the PVDs. These features were not reported in RA22 because of the narrower slit width and different orientations considered there.
We analysed the CO(2-1) and CO(3-2) kinematics at 0 . 6 resolution using the '3D-Based Analysis of Rotating Objects from Line Observations' ( 3D BAROLO) software by Di Teodoro & Fraternali (2015), following the methodology described in RA22. For the 3DFIT, we applied a mask using a threshold of 3σ rms . We fixed the PA and the inclination of the disc to the values derived using the high-resolution CO(2-1) data in RA22: PA disc = 4 • and i = 38 • , and the centre position was fixed to the continuum peak at 220 GHz (RA = 14h30m29.88s and Dec = +13d39m11.93s). Then, we allowed 3D BAROLO to fit the rotation velocity and velocity dispersion. From the model datacube created by 3D BAROLO, we produced PVDs extracted in slits with the same width and orientations described above, and overlaid them on the CO(2-1) and CO(3-2) PVDs (see red contours in Fig. 2). The bulk of the CO(2-1) and CO(3-2) emission at 0 . 6 resolution can be well reproduced by our 3D BAROLO model. We find 77% and 90% of gas in rotation in CO(2-1) and CO(3-2), respectively, with a central velocity dispersion of ∼100 km s −1 . For comparison, in the case of the CO(2-1) data at 0 . 2 resolution, only 55% of the gas can be accounted for by rotation. This is because a significant fraction of noncircular and tangential motions are diluted because of the larger L12, page 2 of 9  The contours are 2, 4, 8, and 16σ rms , with σ rms = 1 mJy beam −1 (left panel) for CO(2-1) and σ rms = 1.9 mJy beam −1 (middle panel) for CO(3-2). The red contours shown in the left and middle panels correspond to the corresponding 3D BAROLO model at 2, 4, and 6σ rms , and the grey dots are the derived rotation velocities projected to the respective PA. Right panels are the PVDs of the CO(3-2)/CO(2-1) line ratios. The colour bars indicate the scales (flux density in mJy for the CO(3-2) and CO(2-1) PVDs and dimensionless in brightness temperature ratio T 32 /T 21 ). beam size of the CO(3-2) data. However, as shown in Fig. 2, the molecular gas with high velocities (∼±400 km s −1 ) cannot be explained by rotation, as the model reproduces the rotation pattern up to maximum values of around ±300 km s −1 .
The right panels of Fig. 2 show the PVDs of the CO(3-2)/CO(2-1) line ratio. To produce these PVDs, the units of the CO(3-2) and CO(2-1) PVDs were first converted to brightness temperature (T 32 and T 21 , in units of K), considering only CO emission above 2σ rms . A clear increase in the line ratio is seen in regions of 0 . 4 radius (∼600 pc) along and perpendicular to the radio jet, up to values of 1.2. In order to study this line ratio in a spatially resolved manner, we created the T 32 /T 21 line ratio of the integrated intensity (moment 0) maps, which is shown in the left panel of Fig. 3. This line ratio map shows the highest values, namely of ∼0.8, along the direction perpendicular to the radio jet, while other regions of the disc present ratios of 0.3 < T 32 /T 21 < 0.5, similar to typical values found in settled discs or in the Milky Way (Leroy et al. 2009;Carilli & Walter 2013). We note that the values of T 32 /T 21 in the integrated lineratio map shown in Fig. 3 are lower than those in the PVDs shown in Fig. 2 because the former are average values along the line of sight (LOS). In the right panel of Fig. 3 we also show the L12, page 3 of 9 σ map of the original CO(2-1) data (i.e., at 0 . 2 resolution). The region with larger σ values (i.e., highest turbulence) coincides with that having the highest values of T 32 /T 21 . As discussed in Sect. 4, enhanced T 32 /T 21 has been reported for a few nearby radio-loud AGN along the direction of the jet (e.g., Dasyra et al. 2016;Oosterloo et al. 2017Oosterloo et al. , 2019, but to the best of our knowledge, this is the first detection of increased T 32 /T 21 along the direction perpendicular to the radio jet in a radio-quiet AGN. Using the 1.5 GHz VLA flux and the spectral index of α = −0.87 reported by Jarvis et al. (2019), it is possible to estimate the jet power of the Teacup from the relations of Bîrzan et al. (2008) and Cavagnolo et al. (2010), for example. If we do so, we obtain P jet 1−3 × 10 43 erg s −1 . However, these scaling relations have large scatter (>0.78 dex) and are drawn from samples of galaxies dominated by more evolved jets. Therefore, they may not apply to compact jets confined within the galaxy ISM, as is the case of the Teacup. As pointed out by Godfrey & Shabala (2016), the P jet -L radio relations may be unreliable because, among other things, they neglect the effect of distance, morphology, and environment. Here, we do not attempt to derive an accurate measurement of the jet power, but are only interested in its order of magnitude, as we aim to compare our findings with the simulations described in Sect. 4.2.
To perform this comparison, we also need an approximate estimate of the jet inclination relative to the CO disc. Although it is not possible to infer a precise angle using the radio data available, we obtain an estimate using the core-to-extended-radioflux ratio parameter, R c = S core S ext , where S core is the core flux density of the unresolved beamed nuclear jet, and S ext is the flux density of the unbeamed lobes. Statistically, R c should be larger in more highly beamed objects, but the extended flux is time dependent and also depends on the frequency of observation. In the case of the Teacup, using the flux densities at 5.2 GHz reported by Jarvis et al. (2019) for the HR VLA data, we obtain R c = (2.3/0.6) = 3.8, which would correspond to an angle of ∼20 • relative to the LOS if we use Eq. (1) in Orr & Browne (1982) with γ = 5 and R T = 0.024. Considering the inclination of the CO disc derived from the modelling with 3D BAROLO (i CO =38 • ; i.e. 52 • relative to the LOS; RA22), we conclude that the jet likely subtends a small angle relative to the CO disc. The results presented here suggest that the radio jet is strongly coupled to the CO disc and because of that, it is driving turbulence and exciting the molecular gas along PA∼−30 • .

Discussion
In this Letter we present robust observational evidence that the radio jet is interacting with the cold ISM of the Teacup and causing the disturbed kinematics, large turbulence, and peculiar excitation conditions of the gas. In the following, we discuss the results and compare them with high-resolution simulations of jets propagating in the ISM of galaxies.

Physical conditions and kinematics of the molecular gas
Spatially resolved CO emission-line ratio analyses have been reported for a few nearby Seyferts and radio galaxies. Viti et al. (2014) modelled the excitation conditions of the molecular gas in the circumnuclear disc of NGC 1068, and claimed that elevated values of R 31 ≡ T 32 /T 10 ∼ 5 (i.e., involving the CO(3-2) and CO(1-0) emission lines) and R 21 ≡ T 21 /T 10 4 (involving CO(2-1) and CO(1-0) lines) correspond to hot dense gas (kinetic temperatures T > 150 K and η > 10 5 cm −3 ) excited by the AGN. Higher values of different emission-line ratios have been reported for the molecular gas interacting with the jet in the radio galaxy PKS 1549-79 (Oosterloo et al. 2019), and also in the Seyfert galaxy IC 5063 (Dasyra et al. 2016;Oosterloo et al. 2017). This indicates that jet-impacted regions have different gas-excitation conditions due to shocks and/or reduced optical thickness (Morganti et al. 2021). IC 5063 is a clear-cut example of jet coplanar with the galaxy disc, which is driving a multiphase outflow (Tadhunter et al. 2014;Morganti et al. 2015). The jet power estimated for IC 5063 is P jet ∼ 10 44 erg s −1 and the PVD along the jet axis revealed T 32 /T 21 values ranging from 1.0 to 1.5 in the outflowing regions, clearly different from the gas following regular rotation (R 32 < 0.7). Under local thermodynamic equilibrium (LTE) conditions, R 32 > 1 corresponds to excitation temperatures (T ex ) of ∼50 K (Oosterloo et al. 2017).
A similar temperature was found for the fast outflow gas component of IC 5063 using R 42 ≡ T 43 /T 21 > 1 (i.e., involving the CO(4-3) and CO(2-1) lines), suggesting that the gas in the outflow is optically thin (Dasyra et al. 2016).
The Teacup presents similar values of T 32 /T 21 to those reported in the outflow region of IC 5063 (i.e., T 32 /T 21 > 1; see Fig. 2). The difference is that in the Teacup, the highest values are found perpendicular and not along the jet direction, as shown in Fig. 3. This is indicative of the presence of hot dense gas, with T ex ∼ 50 K, probably excited by the cocoon of shocked gas that is driving the lateral outflow and/or turbulence that we observe in CO.
Regarding the gas kinematics, examples of enhanced σ of the ionised gas along the direction perpendicular to the jet axis have already been reported in the literature for nearby AGN. One example is the radio galaxy 3C 33, for which Couto et al. (2017) proposed a scenario in which the ISM is laterally expanding due to the passage of the jet, supported by the optical emission line ratios that they measured in the same region, which are typical of gas excitation induced by shocks. Venturi et al. (2021) reported high [Oiii] velocity widths in four nearby Seyfert galaxies from the MAGNUM survey hosting low-power radio jets. These high σ values were interpreted as being due to the passage of a low-inclination jet through the galaxy discs ISM, which produces shocks and induces turbulence. The MURALES survey of nearby 3CR radio sources revealed similar features in [Oiii] for some of the radio galaxies (Balmaverde et al. 2019(Balmaverde et al. , 2022. In the Teacup, we measure the highest values of the CO(2-1) and CO(3-2) σ perpendicular to the jet direction, namely of up to 140 km s −1 . We note that the enhancement in T 32 /T 21 is shifted by ∼0 . 2 to the east relative to the gas with the highest σ (see Figs. 3 and A.1). High-angular-resolution CO(3-2) observations are required in order to explore this shift in more detail.

Is the radio jet driving a lateral molecular outflow?
Simulations show that strong jet-induced feedback can occur even when P jet ∼ 10 43−44 erg s −1 (Mukherjee et al. 2016;Talbot et al. 2022). These low-power jets are trapped in the ISM of the galaxy for longer times and therefore they are able to disrupt the ISM over a larger volume than more powerful jets (Nyland et al. 2018). Another crucial element for having efficient feedback is the jet-ISM coupling. Simulations show that stronger coupling occurs when jets have low inclination relative to the gas disc, because jets are confined within the galaxy ISM, injecting more kinetic energy into the surrounding gas and inducing local turbulence and shocks (Mukherjee et al. 2018a;Meenakshi et al. 2022). According to the estimate described in Sect. 3, this seems to be the case for the Teacup.
We can compare the PVDs shown in Fig. 2 with simulation E in Meenakshi et al. (2022). This simulation consists of a jet of P jet ∼ 10 45 erg s −1 propagating in a gas disc with a central gas density of 200 cm −3 , and subtending an angle of 20 • (i.e., almost coplanar). This jet power is higher than our estimate for the Teacup (∼10 43 erg s −1 ), but in addition to the caveats mentioned in Sect. 3, P jet values estimated from the Bîrzan et al. (2008) and Cavagnolo et al. (2010) relations are usually lower limits for kpc-scale confined jets, such as in the case of the Teacup. For instance, in the case of IC 5063, the jet power required to model the jet-induced molecular gas dispersion was a factor of 3 to 10 higher than that inferred from scaling relations (Mukherjee et al. 2018b). The jet not only induces turbulence due to the injection of energy within its immediate vicinity but also causes a large broadening in the PVDs, with velocities reaching ∼400 km s −1 in the regions both along and perpendicular to the jet path (see Figs. 8 and 10 in Meenakshi et al. 2022).
Motivated by the resemblance between the ionised gas features in Sim E and the molecular gas features in the Teacup, we performed a comparison between our observations and the simulations in Meenakshi et al. (2022). Although the impact of jets on the dense molecular gas is not directly considered in the simulations, the emissivity of the CO(2-1) gas can be approximated with functions that depend on gas density, temperature, and cloud tracers, all extracted from the simulations. This method is similar to that adopted in Mukherjee et al. (2018b) to reproduce the molecular gas features observed in IC 5063. The details of the simulations are presented in Appendix B. As shown in the right panel of Fig. 3, the morphology of the gas velocity dispersion predicted by Sim E is similar to the Teacup, showing high values close to perpendicular to the jet direction. This happens because the main jet stream, as it propagates through the clumpy ISM, is split and deflected multiple times, carrying momentum in all directions including directions perpendicular to the jet axis and accelerating gas out of the disc plane. The simulated CO(2-1) PVDs along and perpendicular to the jet (see Fig. B.1) show high-velocity components reaching ±400 km s −1 within the central 2 kpc of the galaxy. Even though these simulations were not tailored to model the potential and radio power of the Teacup, they qualitatively reproduce the main kinematic features of the CO(3-2) and CO(2-1) observations.

Molecular outflow scenarios
Based on the CO(2-1) data at 0 . 2 resolution, RA22 reported evidence for a molecular outflow coplanar with the CO disc, with a velocity of 185 km s −1 , radius of 0.5 kpc, and deprojected outflow rate ofṀ out = 15.8 M yr −1 . This outflow rate is most likely a lower limit in the case of the Teacup because of the complex CO kinematics, and because RA22 only considered non-circular gas motions along the CO disc minor axis to compute the outflow mass. However, the PVDs along and perpendicular to the jet axis in Fig. 2 show high velocities, of up to ±400 km s −1 , that cannot be reproduced with rotation, and both T 32 /T 21 and σ (see Fig. 3) show higher values in the direction perpendicular to the jet (i.e., closer to the CO major axis). Therefore, it seems plausible that the higher gas excitation and increased turbulence that we find along PA = −30 • are signatures of a jet-driven molecular outflow. We then considered different scenarios to compute the molecular outflow mass, which varies between 3×10 7 and 6×10 8 M . The details are given in Appendix C.
Scenario i, the least conservative, consists of assuming that all the molecular gas that is not rotating, at least according to our 3D BAROLO model, is outflowing. Based on this model, we find that half of the molecular gas mass is rotating, and the outflow rate is 44 M yr −1 using an average outflow velocity of 100 km s −1 . Scenario ii assumes that only the high-velocity gas (faster than ±300 km s −1 ) participates in the outflow, and for this we measure 41 M yr −1 . Finally, in Scenarios iii and iv we measure outflow rates of 15.1 and 6.7 M yr −1 by subtracting the rotation curve from the CO(2-1) datacube and then considering high-velocity gas only (faster than ±300 km s −1 in Scenario iii and just the high-velocity wings shown in Fig. C.2 in Scenario iv).
Based on the results presented in Sect. 3, aided by the simulations discussed in Sect. 4.2, we favour a scenario where the jet is pushing the molecular gas out of the disc plane and producing a lateral expansion that compresses the gas, increasing turbulence and gas excitation. Therefore, in principle we favour Scenarios L12, page 5 of 9 ii and iii, which consider high-velocity gas only, without specifying any particular direction for the outflow. Considering this, the outflow mass would be in the range 5.6-16×10 7 M , and the mass outflow rate within 15-41 M yr −1 .
Jets produce efficient feedback by increasing the turbulence of the dense gas when the ratio between the radio power and the Eddington luminosity, P jet /L Edd > 10 −4 (Wagner et al. 2012). In the Teacup, this ratio is ∼9×10 −4 , and the jet power is sufficient to drive the molecular outflow because it is about two or three orders of magnitude higher than the outflow kinetic power (see Table C.1). However, the fate of the molecular gas during this feedback episode is a galactic fountain, because the escape velocity of the galaxy v esc 530 km s −1 , estimated from its dynamical mass as in Feruglio et al. (2020), is larger than the fastest velocities in the outflow (∼±400 km s −1 ). Regardless, this will have a global effect on the evolution of the galaxy, redistributing mass, metals, and delaying further star formation.
Based on the results presented here, together with the predictions from the models described in Sect. 4.2, jet-induced turbulence and shocks are the most likely mechanisms to explain the high values of σ and T 32 /T 21 reported here for the cold molecular gas in the Teacup. Larger samples of radio-quiet quasars showing low-power radio jets are required in order to better quantify the impact of jet-ISM coupling in the kinematics of the molecular gas.

Appendix A: Moment maps of the CO(3-2) emission
We computed the moment maps of the CO(3-2) emission at ∼0 . 6 resolution as described in Section 3. The integrated intensity, intensity-weighted velocity field, and velocity dispersion maps (moments 0, 1, and 2) are shown in Figure A.1.

Appendix B: Simulation of jet propagating in the molecular disc
In order to elucidate the nature of the jet-ISM interaction in the Teacup, we present here qualitative comparisons with the kinematics of the dense gas in simulations D and E of Mukherjee et al. (2018a, hereafter Sim D and E) using the diagnostics developed in Meenakshi et al. (2022). In these simulations, jets of power 10 45 erg s −1 are launched at angles of 45 • and 20 • from the disc plane for Sim D and E respectively. These inclinations resemble the angle subtended by the jet and the CO disc in the Teacup. It should be noted that although the simulations are not tailored to mimic the Teacup, we performed a qualitative comparison to try to understand the jet-gas interaction in this galaxy, as in Murthy et al. (2022).
In these simulations, a turbulent dense gas disc of radius ∼2 kpc is placed in the X − Y plane. A pair of relativistic jets are launched from the centre. The jet's axis lies in the X − Z plane, inclined from the disc at different angles of choice. According to the simulation's conventions, θ I = 90 • , Φ I = 360 • corresponds to the edge-on view, with the observer facing the X − Z plane. For the analysis below, we have used the value θ I = 52 • to match that of the Teacup. The azimuthal angle φ I is varied to match the observed relative orientation of the projected jet's axis to the disc. The image plane is inverted to the observed orientation of Teacup to facilitate the visual comparison of the simulated maps.
For the analysis, we constructed mean LOS velocities, and the PVDs weighted by synthetic CO luminosities, following the method outlined in Section 5.2 of Mukherjee et al. (2018b) and adopting the same parameters as in their work. Dense gas regions with temperatures lower than 5000 K and densities higher than 10 cm −3 were used for this analysis. The PVDs are then produced using slits of 500 pc in width oriented along and perpendicular to the jet. In Figure B.1, the top left panel shows the projected velocity map for Sim E at a θ I = 52 • and φ I = 340 • , overlaid by the jet tracer at a value of 0.05. At this orientation, the jet subtends an angle of 25 • from the LOS. The PVDs along and perpendicular to the jet are shown in the top right panels. The jet induces high velocities along both slits, reaching values of up to ∼ ±400 km s −1 , and with the largest dispersion appearing in the central 1 kpc regions. The morphology of the PVD for the perpendicular slit is also very similar to that obtained for the Teacup in Figure 2, which appears to be sharply curved due to contribution from the rotation pattern of the disc. In the bottom panels of Fig. B.1, we show the mean velocity field and PVDs for Sim D, which are produced for an image plane oriented at θ I = 52 • and φ I = 30 • . This corresponds to an angle of 23.5 • between the jet and the LOS. Similar to Sim E, the PVDs along and perpendicular to the jet in Figure B.1 also exhibit velocities of up to ∼ ±400 km s −1 within the central 2 kpc.
The simulated σ maps at the resolution of 100 pc and smoothed in the 2D domain by convolving with a Gaussian of similar kernel width are shown in Figure B.2. The σ distribution from Sim E (top panels) displays the highest values along the direction perpendicular to the jet path, similar to the observed spatial distribution of σ, shown in Figure 3. The σ map from Sim D (see bottom panels of Fig. B.2) also shows high values in the regions perpendicular to jet, although more compact than in Sim E, with the highest values along the jet. This extended velocity dispersion is caused by the jet-driven outflows along the minor axis, as well as along other directions throughout the disc, with the former being more prominent (e.g. see Figs. 7 and  dispersion in the molecular gas distribution that can extend in directions perpendicular to the jet. Such an enhanced dispersion can arise both from direct uplifting of the gas from the disc in the form of outflows (e.g. see Fig. 14 of Meenakshi et al. 2022), as well as turbulence introduced into the gas disc by the jet-driven bubble.

Appendix C: Outflow scenarios
Here we describe the four scenarios that we considered for estimating the outflow masses and mass rates, from the least to the most conservative. For this analysis, we only used the CO(2-1) data at the original 0 . 2 resolution to avoid artefacts or incorrect interpretations of the results due to beam smearing of the larger CO(3-2) beam size.
i) In the first scenario, we assume that all the molecular gas that cannot be reproduced with our rotating disc model is outflowing. Thus, we simply integrated the emission the 3D BAROLO model and subtracted it from the CO(2-1) datacube. We find that only 55% of the CO(2-1) emission can be reproduced with rotation, and the remaining 45% would correspond to an outflow mass of 5.1×10 8 M . However, as mentioned in Section 3, the PVDs perpendicular to the jet (PA=-30 • ; bottom panels in Figure 2) show an inverted S-shape in the central r ∼1 (1.6 kpc), which includes the bulk of the emission, and also low-velocity extended emission at r>1 that is only detected at 2σ rms . Thus, there is some extended low-velocity gas that is not accounted for by our rotating model, and some nuclear highvelocity gas that is considered by 3D BAROLO as rotation while it might not be. This uncertainty is not considered in the errors listed in Table C.1. In this scenario, we adopted an outflow velocity of v out =100 km s −1 , estimated from the mean velocity residuals from the 3D BAROLO fit in RA22. Notes. The uncertainty in α CO =0.8±0.5 M (K km s −1 pc 2 ) −1 (Downes & Solomon 1998) is included in the mass error estimates. †: computed using v out =100 km s −1 , while for the others scenarios we adopted v out =300 km s −1 . The outflow mass fractions were computed using the total mass of M H 2 =6.2×10 9 M from RA22, estimated with α CO =4.35 M (K km s −1 pc 2 ) −1 . Integrated intensity of the high-velocity components of the CO(2-1) emission at 0 . 2 resolution. The colour map corresponds to the moment 0 of the total CO(2-1) emission, and the contributions from velocities faster than v=±300 km s −1 are shown as red and blue contours at (0.35,0.5,0.65,0.8,0.95)×σ max , with σ red max = 0.28 Jy beam −1 km s −1 and σ blue max =0.14 Jy beam −1 km s −1 . The beam size is indicated with a red ellipse in the bottom left corner.
ii) In the second scenario, we consider that only the highest velocities, faster than ±300 km s −1 , correspond to outflowing gas. We chose this value from the comparison of the PVDs shown in Figure 2 and the 3D BAROLO rotation curve, shown in red the same figure. To calculate the flux of the high velocity gas, we created moment 0 maps by selecting only channels above ±300 km s −1 , shown in Figure C.1. In this scenario, the mass corresponding to the outflow is 1.6×10 8 M .
iii) The third scenario is amalgam of scenarios i and ii.
We assume that all the high-velocity gas is outflowing, but the contribution from rotation is subtracted following a similar approach as in García-Burillo et al. (2019). We deprojected the one-dimensional rotation curve to the corresponding velocity field on a pixel-by-pixel basis, and reshuffled the channels in order to remove the rotation component from the velocity axis of the CO(2-1) datacube. The result is a narrow residual CO profile, i.e. without the rotation component, which is shown in Figure C.2. As this method is suited to optimising the signal-tonoise ratio of the emission associated with non-circular motions, it can reveal high-velocity wings that are otherwise too faint to be detected. Using the residual profile and considering only gas with velocities above ±300 km s −1 , we derive a mass of 5.7×10 7 M for the outflow.
iv) The fourth scenario is based on the method described in scenario iii, but it assumes that the bulk of the residual CO profile shown in Figure C.2 corresponds to turbulence unrelated to the outflow. This turbulence can be due to the merger, or produced by the jet, but here we consider that it does not correspond to outflowing gas. The residual profile can be fitted with a single Gaussian of FWHM∼330 km s −1 , as shown in Figure C.2. If we subtract this Gaussian from the total residual profile, we measure a mass of 2.6×10 7 M , which mainly corresponds to the highvelocity wings.
We computed the outflow masses from the integrated flux measurements in the outflow. First, the integrated fluxes are converted to CO luminosities L CO(2−1) , in units of K km s −1 pc 2 using Equation 3 of Solomon & Vanden Bout (2005) and then translated into masses using the CO-to-H 2 (α CO ) conversion factor M H 2 =α CO R 12 L CO(2−1) . Under the assumption that the gas is thermalised and optically thick, the brightness temperature ratio R 12 = L CO(1−0) /L CO(2−1) =1. A conservative α=0.8 is usually adopted for deriving outflow masses (see RA22 and references therein). The values for each scenario are listed in Table C.1.
To compute mass outflow rates, we use an outflow radius of r out =0 . 75 (1.2 kpc), as constrained from the extent of the high-velocity gas with high σ and T 32 /T 21 (see Figures 2  and 3). The outflow velocities reach up to ±400 km s −1 , but here we adopt a more conservative velocity of v out =300 km s −1 , except for Scenario i. For a time-averaged, thin, expelled shell geometry (Rupke et al. 2005),Ṁ out =v out (M out /r out ), which corresponds to the outflow mass averaged over the flow timescale, t flow = R out /v out =3.9 Myr. Considering all the previous, the mass outflow rates are in the rangeṀ out =6.7−44 M yr −1 . The ratio between the kinetic power of the molecular outflow (Ė kin = 0.5Ṁ out v out 2 ) and the jet power varies from 0.004<Ė kin /P jet <0.035 (see Table C.1), indicating that the jet has sufficient kinetic energy to drive the molecular outflow.