# Axial-Current Anomaly in Euler Fluid

A.G. Abanov  Department of Physics and Astronomy and Simons Center for
Geometry and Physics, Stony Brook University, Stony Brook, NY 11794, USA  P.B.
Wiegmann  Kadanoff Center for Theoretical Physics, University of Chicago, 5640
South Ellis Ave, Chicago, IL 60637, USA  Institute for Information
Transmission Problems, Moscow, Russia

(August 29, 2024; August 29, 2024)

###### Abstract

We argue that a close analog of the axial-current anomaly of quantum field
theories with fermions occurs in the classical Euler fluid. The conservation
of the axial current (closely related to the helicity of inviscid barotropic
flow) is anomalously broken by the external electromagnetic field as
$\partial_{\mu}j_{A}^{\mu}=2\,\bm{E}\\!\cdot\\!\bm{B}$ similar to that of the
axial current of a quantum field theory with Dirac fermions such as QED.

## Introduction.

Axial-current anomaly of QED asserts that while the electric (vector) current
of Dirac fermions $j^{\mu}=\bar{\psi}\gamma^{\mu}\psi$ is conserved, the axial
current $j^{\mu}_{A}=\bar{\psi}\gamma_{5}\gamma^{\mu}\psi$ is not

$\displaystyle\partial_{\mu}j^{\mu}=0,$ (1)
$\displaystyle\partial_{\mu}j_{A}^{\mu}=\tfrac{k}{4}{}^{\star}\\!FF\,,$ (2)

where
${}^{\star}\\!F^{\mu\nu}=\tfrac{1}{2}\epsilon^{\mu\nu\lambda\rho}F_{\lambda\rho}$
is the dual electromagnetic tensor. The constant $k$ is integer valued when
electromagnetic tensor $F$ is measured in units of the magnetic flux quantum
$\Phi_{0}=hc/e$. In QED it is $k=2$, the number of Weyl fermions in the Dirac
multiplet. In terms of electric and magnetic field the anomaly (2) reads

$\partial_{\mu}j_{A}^{\mu}=k\bm{E}\\!\cdot\\!\bm{B}\,.$ (3)

The term anomaly emphasizes that while simultaneous transformation
$\psi_{L,R}\to e^{i\alpha}\psi_{L,R}$ of left and right components of the
Dirac multiplet by virtue of the Noether theorem yields the conservation of
the electric charge $Q=\int j^{0}d\bm{x}$, the axial transformation

$\psi_{L,R}\to e^{\pm i\alpha_{A}}\psi_{L,R}$ (4)

does not warrant the conservation of the chirality $Q_{A}=\int
j^{0}_{A}d\bm{x}$ even though it leaves the classical Dirac equation
unchanged. It follows from (2) that

$\displaystyle\tfrac{d}{dt}Q=0,\quad\tfrac{d}{dt}Q_{A}=2\int\bm{E}\\!\cdot\\!\bm{B}d\bm{x}.$
(5)

Obtained in 1969 by Adler Adler  in QED and Bell and Jackiw Bell-Jackiw  in
the liner $\sigma$-model the axial-current anomaly (or PCAC, partial
conservation of axial current) is a fundamental non-perturbative result in
gauge field theories which goes well beyond QED, proven experimentally at
different scales of high energy. The most recent advances are taken place in
heavy-ion collision Kharzeev , the subject which initiated a search for
anomalies in relativistic hydrodynamics Son  (see, also Abanov ). Anomalies
also have important applications in quantum fluids, the most notably in the
superfluid $He^{3}$ Volovik .

## Main Results.

In this paper we show that axial-current anomaly is also a property of a
classical Euler’s hydrodynamics of the ordinary inviscid barotropic charged
fluid. Such fluid is described by the Euler equations with Lorentz force

$\displaystyle\dot{\rho}+\bm{\nabla}(\rho\bm{v})=0,$ (6)
$\displaystyle(\partial_{t}+\bm{v}\\!\cdot\\!\bm{\nabla})m\bm{v}+\bm{\nabla}\mu=e\bm{E}+(e/c)\bm{v}\times\bm{B},$
(7)

where $\mu$, a function of the density $\rho$, is the chemical potential
related to pressure as $dp=\rho d\mu$.

Like QED, the barotropic fluid possesses two locally conserved charges. One is
electric charge (the mass in units of $e$) $Q=\int\rho d\bm{x}$. Its current,
a 4-vector $j^{\mu}=(\rho,\ \rho\bm{v})$ is manifestly divergence-free as it
is stated by the continuity equation (6) and expressed by (1).

The axial charge is the fluid helicity defined in Moffat  as

$\displaystyle\mathcal{H}=(1/\Gamma^{2})\int\bm{v}\cdot(\nabla\\!\times\\!\bm{v})d\bm{x}\,.$
(8)

It is conserved in the absence of external fields. Here $\Gamma$ is a
normalization scale of the dimension of the velocity circulation converting
helicity to a dimensionless number.

Even though we deal with classical flows, it is convenient to normalize
helicity by the Onsager circulation quantum $\Gamma=h/m$ dependent on the
Planck constant $h$. In these units, helicity is an even integer, as in the
case of a superfluid. In the classical fluid, this number is very large but is
still even integer. We adopt the units where $h,e,c=1$, but keep $m$ to
distinguish between the fluid momentum and the velocity. That distinction has
a profound physical significance in the presence of external fields; when the
difference between momentum and velocity becomes important, the helicity is
defined through the density of canonical momentum

$\bm{\pi}=m\bm{v}+\bm{A}$ (9)

as

$\mathcal{H}=\frac{1}{h^{2}}\int\bm{\pi}\cdot(\bm{\nabla}\times\bm{\pi})d\bm{x},$
(10)

where $\bm{A}$ is the electromagnetic vector potential.

We comment that while the fluid momentum (9) is defined up to a gradient of a
function, the helicity (10) is uniquely defined for appropriate boundary
conditions (e.g., a closed manifold).

With or without external field helicity defined by (10) is conserved

$\tfrac{d}{dt}{\mathcal{H}}=0.$ (11)

We argue that the axial charge and its density should be identified with

$\displaystyle Q_{A}$ $\displaystyle=\int\rho_{A}d\bm{x}\,,$ (12)
$\displaystyle\rho_{A}$
$\displaystyle=m\bm{v}\\!\cdot\\!\left(\bm{\omega}+2\bm{B}\right)\,,\
\bm{\omega}=\bm{\nabla}\\!\times(m\bm{v})\,,$ (13)

where $\bm{\omega}$ (defined as a curl of the fluid momentum) is the vorticity
of the fluid. We refer to $Q_{A}$ as a fluid chirality bringing the
terminology closer to that of QED. Chirality is the sum of the fluid helicity
(the first term in (13)) and twice the cross-helicity (the second term).

The factor 2 in (13) eventually converts to the coefficient $k=2$ in the
anomaly equation (2). It is itself is a topological number Volovik . Later we
discuss the relationship between the anomaly and linking numbers and briefly
touch on the topological interpretation of this factor (see (16-18)).

We will show that the chirality (12) obeys the anomaly equation (5), while the
chirality density (13) obeys the local anomaly equation (2)

$\displaystyle\dot{\rho}_{A}+\bm{\nabla}\bm{j}_{A}=2\bm{E}\\!\cdot\\!\bm{B}\,,$
(14)

where the chirality flux $\bm{j}_{A}$ (defined modulo a curl) is explicitly
given by

$\displaystyle\\!\\!\\!\\!\\!\bm{j}_{A}\\!=\\!\rho_{A}\bm{v}\\!+\\!(\bm{\omega}\\!+\\!2\bm{B})\left(\\!\mu\\!-\\!\frac{m\bm{v}^{2}}{2}\\!\right)\\!-\\!m\bm{v}\\!\times\\!(\bm{E}\\!+\\!\bm{v}\\!\times\\!\bm{B}).$
(15)

If the flow is incompressible $(\bm{\nabla}\\!\cdot\\!\bm{v}=0)$ then $\mu$ in
(15) is expressed through velocity by the solution of the equation
$\bm{\nabla}^{2}\mu=\bm{\nabla}[\bm{E}+\bm{v}\\!\times\\!\bm{B}-(\bm{v}\\!\cdot\\!\bm{\nabla})m\bm{v}]$.

Later we will clarify the geometric meaning of the cumbersome looking formula
(15).

Eqs. (13,15) could be considered as a realization of the anomaly equations (2)
by hydrodynamics similar to the Wess-Zumino non-linear $\sigma$-model WZ .
Also, the ‘bosonization’ of Dirac fermions in one spatial dimension could be
seen as a precursor of an axial current anomaly in hydrodynamics in higher odd
spatial dimensions as Eqs. (13,15) demonstrate. Similarly, the axial anomaly
appears in non-inertial reference frames. For example, in the case of rotating
fluid subject to a potential force, say, gravity, one replaces the magnetic
field with the frequency of rotation in the equations above.

There is a simple physical picture behind these formulas. Consider a fluid in
electric and magnetic fields in a local reference frame moving and rotating
with the fluid. In this frame the chirality density (13) is locally
approximated by $\rho_{A}=2m\bm{v}\\!\cdot\\!\bm{B}$ and the chirality flux
$\bm{j}_{A}$ is divergent free. Electric field accelerates the fluid
$m\dot{\bm{v}}=\bm{E}$. Hence $\dot{\rho}_{A}=2(\bm{E}\\!\cdot\\!\bm{B})$ as
in (14). Going back to the laboratory frame the magnetic field transforms
$\bm{B}\to\bm{B}+\bm{\omega}/2$ (Larmor precession) and the formula for
chirality transforms as $\rho_{A}=2m\bm{v}\\!\cdot\\!\bm{B}\to
m\bm{v}\cdot(\bm{\omega}+2\bm{B})$. At the same time,
$2(\bm{E}\\!\cdot\\!\bm{B})$ being an invariant undergoes no change. This
yields the formulas (13) and (14). The term $2\bm{B}\mu$ in (15) and the
extension $2\bm{B}\mu\to(\bm{\omega}+2\bm{B})\mu$ is reminiscent the chiral
magnetic effect Kharzeev  and the chiral vortical effect of Ref.
kharzeev2011testing , although in these papers the term
$(\bm{\omega}+2\bm{B})$ appeared in the vector current, not in the axial
current as in (15).

To elucidate the global aspect of the axial anomaly in hydrodynamics we assume
that the vorticity and magnetic field are approximated by vortex and magnetic
filaments having the same circulation $\Gamma$ and the same magnetic flux
$\Phi_{0}$, respectively, as if the fluid were a superfluid (our calculations
do not rely on the assumptions of discreteness but are more of a kinematic
nature). The Helmholtz law warrants that once created, vortex lines and their
bundles can not be destroyed. If vortices and magnetic lines are discrete then
in units $\Gamma$ and $\Phi_{0}$ the fluid helicity
$\int(\bm{v}\cdot\bm{\omega})d\bm{x}$ (10), the magnetic helicity
$\int(\bm{A}\cdot\bm{B})d\bm{x}$, the cross-helicity
$\int(\bm{v}\cdot\bm{B})d\bm{x}$ are topological invariants. In the respective
order they are twice the linkages between vortex lines ${\rm
2\,Link}[\bm{\omega}]$ Moffat , magnetic flux lines ${\rm 2\,Link}[\bm{B}]$
MoffattRicca  and the cross-helicity is the cross-linkage between vortex and
magnetic lines ${\rm Link}[\bm{\omega},\bm{B}]$ Moffatt78 . Hence, the
chirality $Q_{A}$ (12) and helicity (10) are even integers written as sums of
the linkages

$\displaystyle Q_{A}$
$\displaystyle=2{\rm\,Link}[\bm{\omega}]+2{\rm\,Link}[\bm{\omega},\bm{B}]\,,$
(16) $\displaystyle\mathcal{H}\ $
$\displaystyle=2{\rm\,Link}[\bm{\omega}]+2{\rm\,Link}[\bm{\omega},\bm{B}]+{\rm
2\,Link}[\bm{B}]\,.$ (17)

Then, the relation between chirality and helicity is

$\displaystyle Q_{A}=\mathcal{H}-2{\rm\,Link}[\bm{B}].$ (18)

While the individual linkages in the RHS of (17) may change in the course of
the flow, the total helicity $\mathcal{H}$ does not (see (11)). Then the time
derivative of (18) gives the anomaly equation (5). If we treat the chirality
$Q_{A}$ as a difference between the number of left and right moving fermions
as in QED, we may say that an extra link to magnetic lines in accord to (18)
changes the chirality by 2, by flipping the chirality of a fermion from right
to left, and according to (17) changes the sum of linkages of the vortex lines
and the cross-linkage by 1.

As an example of the linkage changing evolution, consider an instantaneous
process of formation of a closed magnetic filament with magnetic flux
$\Phi_{0}$. A fast change of magnetic field triggers a strong electric field
that spins the fluid around the filament. As a result, the vorticity loop is
formed along the magnetic filament. The magnetic field and vorticity thus
created satisfy $\bm{\nabla}\times\bm{\pi}=\bm{\omega}+\bm{B}=0$. The magnetic
and vortical linkages are changed in that process while the total helicity
(10,17) does not.

Eqs.(13-15) and their global version (16) are the major results of this study.

Since Eqs.(13-15) follow from the Euler equation (7) one can validate them by
elementary means being equipped by no more than the vector calculus. Below are
the evolution equations for helicity and cross-helicity densities obtained
directly from the Euler equations (7)

$\displaystyle\begin{split}&\partial_{t}(m\bm{v}\bm{B})+\bm{\nabla}\left[\bm{v}(m\bm{v}\cdot\bm{B})+\bm{B}\left(\mu-\frac{m\bm{v}^{2}}{2}\right)-m\bm{v}\times(\bm{E}+\bm{v}\times\bm{B})\right]+\bm{\omega}(\bm{E}+\bm{v}\times\bm{B})=\bm{E}\bm{B}\,,\\\
&\partial_{t}(\bm{v}\bm{\omega})+\bm{\nabla}\left[\bm{v}(m\bm{v}\cdot\bm{\omega})+\bm{\omega}\left(\mu-\frac{m\bm{v}^{2}}{2}\right)+m\bm{v}\times(\bm{E}+\bm{v}\times\bm{B})\right]-2\bm{\omega}(\bm{E}+\bm{v}\times\bm{B})=0\,.\end{split}$
(19)

Combining we obtain the equation identical to (13-15) with the helicity flux
of the form equivalent to (15).

We emphasize that we do not discuss magnetohydrodynamics in this work. The
formulas (11) and (14) have been derived when the electromagnetic field is
treated as an external having no feedback from the charged fluid motion, nor
influence the fluid through Ohm’s law. In the limit of ideal
magnetohydrodynamics, it is additionally assumed that the plasma is infinitely
conducting so that the Ohm’s law gives $\bm{E}+\bm{v}\times\bm{B}=0$ and
consequently $\bm{E}\\!\cdot\\!\bm{B}=0$. In this case, the cross-helicity and
helicity densities conserve separately as it immediately follows from (19). In
plasmas with finite conductivity, only the total helicity is conserved. The
interconversion of various helicities as in (17) in poorly conducting
relativistic plasmas has been discussed in beckenstein1987helicity .

Behind the straightforward algebra yielding to (19) and (13,14,15) there are
deeper symmetry and geometry-based reasons (see Arnold  as a general reference
for geometric view on hydrodynamics). Here we only touch the surface leaving a
more comprehensive discussion to future publications. We start from the
derivation of the helicity conservation in the form that makes the
conservation of chirality an easy corollary.

## Vorticity transport and helicity conservation.

We will use the 4-dimensional space-time formalism. The formalism is standard
in relativistic hydrodynamics but is not common in studies of non-relativistic
flows. Still, we find that it is the most compact way to expose the geometric
nature of fundamental laws of the Euler flow: Helmholtz law for advection of
vorticity and the conservation of helicity with or without external fields and
in non-relativistic or relativistic hydrodynamics alike. A reason is that
these laws are expressed in terms of differential forms and, therefore, are
not sensitive to the space-time metric.

We start by writing the Euler equation (7) in terms of the fluid momentum

$\displaystyle\rho(\dot{\bm{\pi}}-\bm{\nabla}\pi_{0})-\rho\bm{v}\times(\bm{\nabla}\times\bm{\pi})=0.$
(20)

Here $\pi_{0}$ denotes the Bernoulli function

$\displaystyle\pi_{0}=\Phi+A_{0},\quad-\Phi=\mu+\tfrac{1}{2}mv^{2},$ (21)

where $A_{0}$ is the electrostatic potential.

Next we recognize the mass 4-current $j^{\mu}=(\rho,\rho v^{i})\,$ as a
4-vector field and the 4-momentum $\pi_{\mu}=(\pi_{0},\pi_{i})$ as a covector
field. In these terms the continuity equation has the form (1) and the Euler
equation (2) appears in a remarkably compact form. It follows from (20)

$\displaystyle\quad j^{\mu}\Omega_{\mu\nu}=0\,,$ (22)

where

$\displaystyle\Omega_{\mu\nu}=\partial_{\mu}\pi_{\nu}-\partial_{\nu}\pi_{\mu}$
(23)

the 4-vorticity antisymmetric tensor extended by electromagnetic field (also
referred as canonical symplectic form or Khalatnikov canonical vorticity
tensor see, e.g., G2 ; CarterKhalatnikov  and references therein).

Euler equations in various forms (7,20,22) are all equivalent. The advantage
of the form (22) is that it is insensitive to the space-time structure. For
example, it stays the same regardless of whether the space-time is Galilean or
Minkowski. The information of the space-time structure is delegated to the
relation between the fluid momentum $\pi$ and the vector-current $j$. For the
Galilean barotropic fluid these relations are (9,21). In the relativistic
context, the Euler equation in the form (22) is known as the Carter-
Lichnerowitz equation (see G  for an excellent review of the subject). Eq.(22)
is the basis of geometric interpretation of the Euler flow. It states that the
vorticity vector field spans a bundle of integral two-dimensional surfaces
normal to the space-time flow and that these surfaces form a foliation of the
space-time.

The spatial components of the 4-vorticity are:

$\displaystyle\Omega_{ij}=\epsilon_{ijk}(\omega^{k}+B^{k}).$ (24)

The Euler equation connects the space-time component
$\Omega_{0i}=\dot{\pi}_{i}-\partial_{i}\pi_{0}$ to the spatial components by

$\Omega_{0i}=v^{j}\Omega_{ji}.$ (25)

This relation states that 4-vorticity, similar to the density, is advected by
the flow (the Helmholtz law): vorticity can not be destructed or created;
instead, it moves with the flow. Equivalently, it implies that the Lie
derivative of vorticity vanishes along the flow with the mass current
$j^{\mu}$.

It is customary to present these equations in terms of differential forms. We
assemble the 4-momentum 1-form $\pi=\pi_{\mu}dx^{\mu}$ and 4-vorticity 2-form
$\Omega=\Omega_{\mu\nu}dx^{\mu}\wedge dx^{\nu}$, a closed 2-form equal to the
exterior derivative of the momentum

$\displaystyle\Omega=d\pi,\quad d\Omega=0.$ (26)

The Euler equation in the form of Carter-Lichnerowitz (22) is the statement
that the 1-form obtained by the interior product between the current vector
field and vorticity 2-form vanishes

$\displaystyle\iota_{j}\Omega=0.$ (27)

Now we turn to helicity. In the 4-dimensional formalism helicity is the 3-form

$\displaystyle h=\pi\wedge d\pi=\pi\wedge\Omega.$ (28)

The components of $h$ are the helicity density $h_{0}$ and the flux $\bm{h}$.
They read

$\displaystyle\begin{split}h_{0}&=\bm{\pi}\\!\cdot(\bm{\nabla}\\!\times\\!\bm{\pi}),\\\
\bm{h}&=\bm{\pi}\times(\dot{\bm{\pi}}\\!-\bm{\nabla}\pi_{0})\\!-\\!\pi_{0}(\bm{\nabla}\\!\times\\!\bm{\pi})\\\
&=h_{0}\bm{v}-(\bm{\nabla}\times\bm{\pi})(\bm{\pi}\cdot\bm{v}+\pi_{0}).\end{split}$
(29)

We note that counter to vorticity, helicity is not frozen into the flow as
$\bm{h}\neq h_{0}\bm{v}$.

We would like to show that the helicity 3-form (28) is closed

$\displaystyle dh=\Omega\wedge\Omega=0.$ (30)

Eq. (30) amounts to the conservation of helicity (11)

$\displaystyle\dot{h}_{0}+\bm{\nabla}\bm{h}=0.$ (31)

One can check (31) either by elementary algebra with the help of the Euler
equation or apply the following arguments. In four dimensions,
$\Omega\wedge\Omega=d\pi\wedge d\pi$ is a 4-form, hence it is proportional to
the 4-volume form. On the other hand, it follows from (27) that
$i_{j}(\Omega\wedge\Omega)=0$. Assuming that the fluid density does not vanish
and, therefore, $j\neq 0$, we conclude that the proportionality coefficient
between $\Omega\wedge\Omega$ and the volume form is zero. This implies the
conservation of helicity in the form (30).

To have the same argument in component notations we use the identity

$\displaystyle\epsilon^{\lambda\rho\nu\alpha}\partial_{\lambda}\pi_{\rho}\partial_{\mu}\pi_{\nu}=\delta^{\alpha}_{\mu}\epsilon^{\lambda\rho\nu\tau}\partial_{\lambda}\pi_{\rho}\partial_{\nu}\pi_{\tau}.$
(32)

The RHS of the identity is $(\Omega\wedge\Omega)\delta^{\alpha}_{\mu}$. If we
contract one free index, say $\mu$, of the LHS with the current $j^{\mu}$,
then the equation (22) prompts that the contraction vanishes and we get
$0=j^{\alpha}(\Omega\wedge\Omega)$, hence (30).

We comment that the derivation of the helicity conservation does not utilize
the relation between the current $j$ and the momentum $\pi$. It relies on the
Eq.(22) (or (27)) which states the geometric property of the flow: surfaces
spanned by the vorticity vector field orthogonal to the current form a
foliation of the space time.

## Axial anomaly.

When we invoke the relation between the current and the momentum (9) we
encounter a caveat common to gauge theories. The canonical momentum $\pi$,
and, therefore, helicity 3-form (28,29) are local in terms of the gauge
potential but cannot be locally expressed through $\bm{E}$ and $\bm{B}$. At
the same time without electromagnetic field the helicity form (28) is a local
functional of the Eulerian fields

$\displaystyle\begin{split}&h_{0}=m\bm{v}\\!\cdot\\!\bm{\omega},\quad\bm{h}=h_{0}\bm{v}-\bm{\omega}\left({m\bm{v}^{2}}+\Phi\right).\end{split}$
(33)

[This form of helicity are equivalent to the chirality (29) up to exact form
and identical to (13,15)].

It is desirable to extend the formulas (33) for non-vanishing gauge potential
in such a manner that they remain local in terms of electric and magnetic
fields. The 4-dimensional formalism yields the result in few lines. In this
framework the 4-chirality is the 3-form

$\displaystyle j_{A}=(\pi-A)\wedge(d\pi+dA),$ (34)

where $A=A_{\mu}dx^{\mu}$ is the gauge potential 1-form. Components of the
chirality form are defined by (13,15). More precisely, (34) yields the
expression for the flux

$\bm{j}_{A}=m\bm{v}\times(m\dot{\bm{v}}-\nabla\Phi-2\bm{E})-\Phi(\bm{\omega}+2\bm{B})\,,$
(35)

which under substitution $m\dot{\bm{v}}$ from the Euler equation appears in
the form which differs from (15) by a curl.

Helicity and chirality are related by the identity

$\displaystyle\pi\wedge d\pi=j_{A}+d(A\pi)-A\wedge dA\,.$ (36)

We take the exterior derivative of both parts of (36) and recall that the
helicity is a closed form (30). We obtain the axial anomaly (2) announced in
the introduction in terms of the chirality 3-form and the electromagnetic
2-form $F=dA=\tfrac{1}{2}F_{\mu\nu}dx^{\mu}\wedge dx^{\nu}$

$dj_{A}=F\wedge F\,.$ (37)

This simple procedure explains the origin of the anomaly in the context of
hydrodynamics: helicity $\mathcal{H}$ is conserved, but its current (28) is
not local in terms of Eulerian fields. At the same time chirality is comprised
locally by $\bm{E}$ and $\bm{B}$ but the Chern-Simons term $A\wedge dA$ in
(36) makes its conservation anomalous. The exterior derivative of the Chern-
Simons is $\bm{E}\cdot\bm{B}$. It contributes the anomaly term in the
divergence of the axial current.

## Relativistic hydrodynamics.

The derivation above remains unchanged for _relativistic_ barotropic flow.
Relativistic flow is characterized by the specific Gibbs energy
$\mu_{R}=mc^{2}+\mu$ which differs from the chemical potential $\mu$ by the
rest energy, 4-velocity $u_{\mu}$ normalized as $u^{\mu}u_{\mu}=-1$, and
canonical momentum 1-form $\pi=(c^{-1}\mu_{R}u_{\mu}+A_{\mu})dx^{\mu}$. In
this case the mass current is $j^{\mu}=nu^{\mu}$, where $n$ is particle number
and the 4-momentum is $\pi_{\mu}=\mu_{R}u_{\mu}$.

The formula for the chirality 3-form (34) remains intact and obeys the anomaly
equation (2). In components, the chirality current reads

$\displaystyle
j_{A}^{\mu}=\tfrac{1}{2}\epsilon^{\mu\nu\lambda\rho}\mu_{R}^{2}(u_{\nu}\partial_{\lambda}u_{\rho}+2u_{\nu}F_{\lambda\rho})\,.$
(38)

It is instructive to trace how does the compact relativistic expression
$(\ref{42})$ turns into the cumbersome non-relativistic formula (15). The
temporal component of (38) $j^{0}_{A}$ reproduces the chirality density (13)
as $v/c\to 0$ in a straightforward manner. But the spatial component of (38),
the flux, as it stands, does not have the non-relativistic limit. It diverges
as $j_{A}^{i}\to c(\omega^{i}+2B^{i})$ with $c\to\infty$. However, the
divergent term is a curl. We recall that the formula (38) as well as its non-
relativistic version (35) is defined up to a curl. Hence, the divergent term
can be removed (by an appropriate choice of the frame). As a result, the non-
relativistic limit is a mix of terms originated from the rest energy in
$\mu_{R}$ and the expansion in $v/c$.

## Summary.

Alike quantum field theories with Dirac fermions, Euler hydrodynamics, too,
possesses a conserved vector current and an axial current, helicity. Both are
conserved. The conservation of the vector current is explicitly imposed as the
continuity equation and is associated with the gauge symmetry. The origin of
helicity conservation is more subtle. It follows from the Euler equation and
could be traced to different but ultimately related roots: a topological
nature of helicity as a linkage of vortex lines, or to the degeneracy of
Poisson structure and its foliations, or to the relabeling symmetry acting in
the extended phase space (see, e.g., Morrison ; Zakharov  ). Our results
indicate that these properties are represented by the axial gauge
transformation (4). We defer these discussions pending further studies.

These results also suggest that there may be special hydrodynamic flows
representing fermions of quantum field theory if the analogy is to be followed
further.

From a hydrodynamics perspective, the origin of the anomaly may be summarized
as follows. Electromagnetic forces do not destroy the conservation of
helicity. However, the helicity current, although conserved (31) is expressed
through the canonical fluid momentum. That prevents treating the helicity
current as a local functional of $\bm{E}$ and $\bm{B}$. Similar to QED, the
“conflict” is resolved at the expense of the conservation of the axial
current. The chirality defined by (12,13), is identical to the helicity in the
absence of electromagnetic fields and is local in terms of $\bm{E}$ and
$\bm{B}$ but does not conserve, instead, obeying the equation (14) identical
to the axial anomaly equation (2) of the quantum field theory.

The appearance of the axial-current anomaly in 3D Euler flow is not
accidental. All formulas presented in this work have a straightforward
generalization to hydrodynamics in higher odd spatial dimensions matching the
axial-current anomaly of quantum field theories with Dirac fermions. Also, we
expect that not only perturbative anomalies but also global anomalies appear
in classical hydrodynamics.

Emergence of the axial current anomaly in Euler fluids helps to clarify and to
illustrate aspects of anomalies in quantum field theories also providing an
interesting angle at fluid dynamics.

The authors thank A. Cappelli and B. Khesin for helpful discussions. P.W.
acknowledges helpful discussions with V.P. Nair, A. Polyakov and G. Volovik.
P.W. thanks G. Volovik for pointing out the association of the axial current
anomaly with the fluid helicity discussed in V1  in the context of superfluid
$\operatorname{{}^{3}He-A}$. P.W. thanks Simons Center for Geometry and
Physics for the hospitality at various stages of this work. Both authors
acknowledge INFN Galileo Galilei Institute where the work had been completed.
A.G.A. thanks S. Vitouladitis for insightful comments. The work of P.W. was
supported by the NSF under the Grant NSF DMR-1949963 and by the Ministry of
Science and Higher Education of the Russian Federation, agreement No
075-15-2021-602. The work of A.G.A. was supported by NSF DMR-1606591.

## References

  * (1) S.L. Adler. Axial-vector vertex in spinor electrodynamics.  Phys. Rev. 177: 2426, 1969.
  * (2) J. S. Bell, R. Jackiw. A PCAC puzzle: $\pi_{0}\to\gamma\gamma$ in the $\sigma$ -model. Il Nuovo Cimento A. 60: 47, 1969.
  * (3) D. E. Kharzeev. The Chiral Magnetic Effect and anomaly-induced transport. Progress in Particle and Nuclear Physics 75: 133, 2014.
  * (4) D.T. Son and P. Surowka. Hydrodynamics with triangle anomalies. Phys. Rev. Lett., 103: 191601, 2009.
  * (5) G.M. Monteiro, A.G. Abanov, and V.P. Nair. Hydrodynamics with gauge anomaly: Variational principle and Hamiltonian formulation. Phys. Rev. D91:125033, 2015.
  * (6) G.E. Volovik The universe in a helium droplet: Oxford University, 2009.
  * (7) D.E. Kharzeev and D.T. Son, Testing the chiral magnetic and chiral vortical effects in heavy ion collisions. Phys. Rev. Lett., 106(6): 062301, 2011\.
  * (8) J. Wess and B. Zumino. Consequences of anomalous Ward identities. Phys. Lett. B 37: 95, 1971.
  * (9) H.K. Moffat. The degree of knottedness of tangled vortex lines. Fluid Mech. 35: 117, 1969.
  * (10) H.K. Moffatt and R.L. Ricca. Helicity and Calugareanu invariant. Proc. Roy. Soc. Ser. A 439: 411, 1992.
  * (11) H.K. Moffatt. Magnetic Field Generation in Electrically Conducting Fluids, Cambridge University Press ,1978.
  * (12) J.D. Bekenstein, Helicity conservation laws for fluids and plasmas, The Astrophys. Journal 319: 207, 1987.
  * (13) V. I. Arnold and B. A. Khesin. Topological methods in hydrodynamics. Springer, 2008.
  * (14) M. Charalampos, Kō. Uryū, E. Gourgoulhon, J.-P. Nicolas, N. Andersson, A. Pouri, and V. Witzany. Conservation laws and evolution schemes in geodesic, hydrodynamic, and magnetohydrodynamic flows.  Phys. Rev.  96: 064019 (2017).
  * (15) B. Carter and I.M. Khalatnikov. Canonically covariant formulation of Landau’s Newtonian superfluid dynamics Rev. Math. Phys. 6: 277 (1994).
  * (16) P.J. Morrison, Hamiltonian description of the ideal fluid. Rev. Mod. Phys., 70: 467, 1998.
  * (17) V.E. Zakharov, E.A. Kuznetsov. Hamiltonian Formalism for Nonlinear Waves. Phys.Uspekhi 40: 1087, 1997.
  * (18) E. Gourgoulhon. An introduction to relativistic hydrodynamic. EASPubl.Ser. 21:43, 2006.
  * (19) G.E. Volovik. Topological Invariants for Standard Model: From Semi-Metal to Topological Insulator. JETP Lett. 113(8):, 538, 2021.


11institutetext: Instituto de Astrofísica de Canarias, Calle Vía Láctea, s/n,
E-38205, La Laguna, Tenerife, Spain 11email<EMAIL_ADDRESS>22institutetext: Departamento de Astrofísica, Universidad de La Laguna,
E-38206, La Laguna, Tenerife, Spain 33institutetext: Observatorio Astronómico
Nacional (OAN-IGN)-Observatorio de Madrid, Alfonso XII, 3, 28014 Madrid, Spain
44institutetext: Observatoire de Paris, LERMA, Collège de France, PSL
University, CNRS, Sorbonne University, Paris, France 55institutetext:
Dipartimento di Fisica, Sezione di Astronomia, Universitá di Trieste, via
Tiepolo 11, 34143 Trieste, Italy 66institutetext: Inter-University Centre for
Astronomy and Astrophysics, Savitribai Phule Pune University Campus, Pune-
411007, India 77institutetext: Research School of Astronomy and Astrophysics,
The Australian National University, Canberra, ACT 2611, Australia
88institutetext: Center for Computational Sciences, University of Tsukuba,
Tennodai 1-1-1, 305-0006, Tsukuba, Ibaraki, Japan

# Jet-induced molecular gas excitation and turbulence in the Teacup

A. Audibert 1122 C. Ramos Almeida 1122 S. García-Burillo 33 F. Combes 44 M.
Bischetti 55 M. Meenakshi 66 D. Mukherjee 66 G. Bicknell 77 A. Y. Wagner 88

(Received January 20, 2023; accepted February 26, 2023)

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_{\rm L}$=388 Mpc) radio-quiet type-2 quasar (QSO2) with a compact radio
jet ($P_{\rm jet}\approx$1043 erg s-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 find this enhancement in the
velocity dispersion but also a higher brightness temperature ratio ()
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
outflow 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.

###### Key Words.:

galaxies: active – galaxies: Individual: Teacup – galaxies: kinematics and
dynamics – galaxies: jets – ISM: jets and outflows

##  1 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 sub-grid 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, are jets launched by the supermassive black hole (SMBH).
These jets can strongly impact (sub-)kpc 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 kpc 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
$\sim$0$\aas@@fstack{\prime\prime}$2 (400 pc) resolution of 7 radio-quiet
type-2 quasars (QSO2s, i.e., obscured quasars) at redshifts z$\sim$0.1, with
bolometric luminosities of L${}_{\rm bol}\approx$1045-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 ultra-luminous
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 $\sim$0.8 kpc, PA=60∘)
detected in Very Large Array (VLA) data at 0$\aas@@fstack{\prime\prime}$25
resolution (Harrison et al. 2015), makes this object a promising target to
study the jet-ISM interaction. It 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 of the Teacup stems from the 12 kpc ionized 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 Figure 1. The [Oiii] kinematics reveal a
$\sim$1 kpc scale nuclear outflow with maximum velocities of $\sim$-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$\alpha$ and [Si VI] emission lines, also extending
$\sim$1 kpc and oriented at PA$\sim$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 characterized in different gas phases
(ionized, warm, and cold molecular). This, together with its compact radio
jet, makes it an ideal target to advance in our understanding of how compact
jets impact the multi-phase ISM of radio-quiet AGN.

Figure 1: The Teacup as seen in [Oiii], radio continuum, and CO. The left
panel shows the large scale (15″$\times$15″) [Oiii]$\lambda$5007 Å HST image
in colour with the low-resolution (LR) VLA 6 GHz continuum contours overlaid
in blue. The middle panel shows the
2$\aas@@fstack{\prime\prime}$8$\times$2$\aas@@fstack{\prime\prime}$8 peak
intensity map of CO(3-2) with the high-resolution (HR) VLA 6 GHz continuum
contours in gray. The right panel corresponds to the CO(2-1) peak intensity
map with the ALMA 220 GHz continuum contours overlaid in magenta. The beam
sizes are shown in the bottom left corner of each panel with corresponding
colours.

Here we present new ALMA CO(3-2) observations at
$\sim$0$\aas@@fstack{\prime\prime}$6 (960 pc) resolution, in addition to the
CO(2-1) 0$\aas@@fstack{\prime\prime}$2 data from  RA22, to investigate the
influence of the jet on the molecular gas excitation and kinematics of the
Teacup. We adopt a flat $\Lambda$CDM cosmology with H0=70 km s-1Mpc-1,
$\Omega_{\rm M}$=0.3, and $\Omega_{\rm\Lambda}$=0.7.

##  2 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$\aas@@fstack{\prime\prime}$21$\times$0$\aas@@fstack{\prime\prime}$18 and
a root mean square (rms) noise of 0.39 mJy/beam 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 done 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 $\sim$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. As flux, bandpass, and phase
calibrator the sources J1550+0527, J1337-1257 and J1446+1721 were used. 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 0th 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$\aas@@fstack{\prime\prime}$60 $\times$0$\aas@@fstack{\prime\prime}$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$\aas@@fstack{\prime\prime}$60$\times$0$\aas@@fstack{\prime\prime}$54 (960
pc$\times$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$\aas@@fstack{\prime\prime}$25, HR) A-array and the low-resolution
(1$\aas@@fstack{\prime\prime}$0, LR) B-array maps. The data were presented and
analysed in full by Jarvis et al. (2019).

##  3 Results and analysis

The CO(3-2) peak intensity map, shown in the middle panel of Figure 1, shows
the double peaked morphology first detected in CO(2-1) and discussed in  RA22.
The distance between the two peaks is $\sim$0$\aas@@fstack{\prime\prime}$8
(1.3 kpc). The integrated intensity, intensity weighted velocity field, and
velocity dispersion ($\sigma$) maps of the CO(3-2) emission are shown in
Appendix A. Overall, the 1st and 2nd moment maps (velocity field and $\sigma$)
are similar to their CO(2-1) counterparts, although with less 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 2nd order moment map shows high values of
$\sigma$, reaching $\sim$120 km s-1 with PA$\sim$-30∘, orthogonal to the radio
jet orientation (PAjet=60∘), similar to the $\sigma$ map of the CO(2-1) data
in  RA22.

In order to investigate the possible influence of the jet on 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$\aas@@fstack{\prime\prime}$6, equivalent the common resolution of
the CO(3-2) and convolved CO(2-1) cubes. The PVDs are shown in Figure 2 and
the position and orientation of the slits are indicated as shaded rectangles
in the left panel of Figure 3. A broad gradient of velocities is detected both
along and perpendicular to the jet axis, spanning up to $\pm$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), specially in the case of the
blueshifted velocities. The PVDs perpendicular to the radio jet show an
inverted S-shape in the central $\sim$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 resemble 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 there considered.

Figure 2: PVDs along (PA=60∘) and perpendicular (PA=-30∘) to the radio jet
direction, extracted in a slit of 0$\aas@@fstack{\prime\prime}$6 width and
showing CO emission. The contours are 2, 4, 8, and 16$\sigma_{\rm rms}$, with
$\sigma_{\rm rms}$=1 mJy beam-1 (left panel) for CO(2-1) and $\sigma_{\rm
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 3DBAROLO model at 2, 4
and 6$\sigma_{\rm 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 ).

We analysed the CO(2-1) and CO(3-2) kinematics at
0$\aas@@fstack{\prime\prime}$6 resolution using the “3D-Based Analysis of
Rotating Objects from Line Observations” (3DBAROLO) software by Di Teodoro &
Fraternali (2015), following the methodology described in  RA22. For the
3DFIT, we applied a mask using a threshold of 3$\sigma_{\rm 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: PAdisc=4∘and $i$=38∘, and the centre
position was fixed to the continuum peak at 220GHz (RA=14h30m29.88s and
Dec=+13d39m11.93s). Then, we allowed 3DBAROLO to fit the rotation velocity and
velocity dispersion. From the model datacube created by 3DBAROLO, we produced
PVDs extracted in slits with the same width and orientations described above,
and overlaid them to the CO(2-1) and CO(3-2) PVDs (see red contours in Figure
2). The bulk of the CO(2-1) and CO(3-2) emission at
0$\aas@@fstack{\prime\prime}$6 resolution can be well reproduced by our
3DBAROLO 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 $\sim$100 km s-1. For
comparison, in the case of the CO(2-1) data at 0$\aas@@fstack{\prime\prime}$2
resolution, only 55% of the gas can be accounted for by rotation. This is
because a significant fraction of non-circular and tangential motions are
diluted due to the larger beam size of the CO(3-2) data. However, as shown in
Fig. 2, the molecular gas with high velocities ($\sim\pm$400 km s-1) cannot be
explained by rotation, as the model reproduces the rotation pattern up to
maximum values of around $\pm$300 km s-1.

The right panels of Figure 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 (T32 and T21, in units of K), considering
only CO emission above 2$\sigma_{\rm rms}$. A clear increase in the line ratio
is seen in regions of 0$\aas@@fstack{\prime\prime}$4 radius ($\sim$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 line
ratio of the integrated intensity (moment 0) maps, which is shown in the left
panel of Figure 3. This line ratio map shows the highest values, of $\sim$0.8,
along the direction perpendicular to the radio jet, while other regions of the
disc present ratios of 0.3$<T\rm_{32}/T\rm_{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 in the integrated line ratio map
shown in Figure 3 are lower than those in the PVDs shown in Figure 2 because
the former are average values along the LOS. In the right panel of Figure 3 we
also show the $\sigma$ map of the original CO(2-1) data (i.e., at
0$\aas@@fstack{\prime\prime}$2 resolution). The region with larger $\sigma$
values (i.e., highest turbulence) coincides with that having the highest
values of . As discussed in Section 4, enhanced 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. 2017, 2019), but to the best of our knowledge, this is
the first detection of increased along the direction perpendicular to the
radio jet in a radio-quiet AGN.

Figure 3: Brightness temperature ratio () and CO(2-1) velocity dispersion
maps. The left panel shows the map at 0$\aas@@fstack{\prime\prime}$6
resolution with the VLA 6 GHz HR contours overlaid in black (at 4, 8, 16, 32,
and 64$\sigma_{\rm rms}$, with $\sigma_{\rm rms}$=23 $\mu$Jy/beam). The
PAjet=60∘ is indicated as a black solid line. The orientation and width
(0$\aas@@fstack{\prime\prime}$6) of the slits used to produce the PVDs shown
in Figure 2 are indicated as grey shaded areas. The middle panel corresponds
to the $\sigma$ map of the CO(2-1) line emission at
0$\aas@@fstack{\prime\prime}$2 resolution. The cyan star in the left and
middle panels indicates the peak of the 220 GHz ALMA continuum. On the right
panel, we show the $\sigma$ map of the gas from Simulation E in Meenakshi et
al. (2022) for a jet power of $P_{\rm jet}=10^{45}$ erg s-1, inclined 20∘ from
the CO disc plane, and with an angle of $\sim$25∘ relative to the LOS.

Using the 1.5 GHz VLA flux and the spectral index of $\alpha$=-0.87 reported
by Jarvis et al. (2019), it is possible to estimate the jet power of the
Teacup from the relations of e.g., Bîrzan et al. (2008) and Cavagnolo et al.
(2010). If we do so, we obtain $P_{\rm jet}\simeq$1-3$\times$1043 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 into the galaxy ISM as it is the case
of the Teacup. As pointed out by Godfrey & Shabala (2016), the $P_{\rm
jet}$-${L_{\rm radio}}$ relations may be quite 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, we are
only interested in its order of magnitude, as we aim at comparing with the
simulations described in Section 4.2.

To perform this comparison, we also need a rough approximation 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 by using
the core to extended radio flux ratio parameter, $R_{\rm c}=\frac{S_{\rm
core}}{S_{\rm ext}}$. $S_{\rm core}$ is the core flux density of the
unresolved beamed nuclear jet, and $S_{\rm ext}$ is the flux density of the
unbeamed lobes. Statistically, $R_{\rm c}$ should be larger in more beamed
objects, but the extended flux is time-dependent and it 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, $R_{\rm
c}$=(2.3/0.6)=3.8, which would correspond to an angle of $\sim$20∘ relative to
the line-of-sight (LOS) if we use Eq. 1 in Orr & Browne (1982) with $\gamma$=5
and $R_{\rm T}$=0.024. Considering the inclination of the CO disc derived from
the modelling with 3DBAROLO ($i_{\rm 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$\sim$-30∘.

##  4 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 on the ISM of galaxies.

###  4.1 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_{\rm 31}\equiv T_{\rm 32}/T_{\rm
10}\sim$5 (i.e., involving the CO(3-2) and CO(1-0) emission lines) and $R_{\rm
21}\equiv T_{\rm 21}/T_{\rm 10}\gtrsim$4 (involving CO(2-1) and CO(1-0) lines)
correspond to hot dense gas (kinetic temperatures $T>$150 K and $\eta>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 multi-phase
outflow (Tadhunter et al. 2014; Morganti et al. 2015). The jet power estimated
for IC 5063 is $P_{\rm jet}\sim$1044 erg s-1 and the PVD along the jet axis
revealed values ranging from 1.0$-$1.5 in the outflowing regions, clearly
different from the gas following regular rotation ($R_{\rm 32}<0.7$). Under
local thermodynamic equilibrium (LTE) conditions, $R_{\rm 32}>$1 corresponds
to excitation temperatures ($T_{\rm ex}$) of $\sim$50 K (Oosterloo et al.
2017). A similar temperature was found for the fast outflow gas component of
IC 5063 using $R_{\rm 42}\equiv T{\rm{}_{43}}/T{\rm{}_{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 as those reported in the outflow region
of IC 5063 (i.e., $>$1; see Figure 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_{\rm ex}\sim$50 K, probably excited by the cocoon of shocked gas that is
driving the lateral outflow/turbulence that we observe in CO.

Regarding the gas kinematics, examples of enhanced $\sigma$ of the ionised gas
along the direction perpendicular to 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, 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 $\sigma$ values were interpreted as 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 has
revealed similar features in [Oiii] for some of the radio galaxies (Balmaverde
et al. 2019, 2022). In the Teacup, we measure the highest values of the
CO(2-1) and CO(3-2) $\sigma$, of up to 140 km s-1, perpendicular to the jet
direction. We note that the enhancement in is shifted by
$\sim$0$\aas@@fstack{\prime\prime}$2 to the east relative to the gas with the
highest $\sigma$ (see Figures 3 and 4). High angular resolution CO(3-2)
observations are required for exploring this shift in more detail.

###  4.2 Is the radio jet driving a lateral molecular outflow?

Simulations show that strong jet-induced feedback can occur even when $P_{\rm
jet}\sim$1043-44 erg s-1 (Mukherjee et al. 2016; Talbot et al. 2022). These
low-power jets are trapped in the galaxy’s ISM 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, since 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 Section 3, this seems to be the
case of the Teacup.

We can compare the PVDs shown in Figure 2 with simulation E in Meenakshi et
al. (2022). This simulation consists of a jet of $P_{\rm jet}\sim$1045 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 ($\sim$1043 erg s-1), but in addition to the
caveats mentioned in Section 3, $P_{\rm 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 as it is 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 does not only induce turbulence
due to the injection of energy within its immediate vicinity but also causes a
large broadening in the PVDs, with velocities reaching up $\sim$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 Figure 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 Figure 5) show high
velocity components reaching $\pm$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 reproduce qualitatively the main kinematic
features of the CO(3-2) and CO(2-1) observations.

###  4.3 Molecular outflow scenarios

Based on the CO(2-1) data at 0$\aas@@fstack{\prime\prime}$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
$\dot{M}_{\rm 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 Figure 2 show high velocities, of up to $\pm$400 km s-1, that
cannot be reproduced with rotation, and both and $\sigma$ (see Figure 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$\times$107 and 6$\times$108
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 3DBAROLO model,
is outflowing. Based on this model, we find that half of the molecular gas
mass is rotating, and the mass outflow rate is 44 M$\rm{}_{\odot}\,yr^{-1}$
using an average outflow velocity of 100 km s-1. Scenario ii assumes that only
the high-velocity gas (faster than $\pm$300 km s-1) participates in the
outflow, and for this we measure 41 M$\rm{}_{\odot}\,yr^{-1}$. Finally, in
Scenarios iii and iv we measure outflow rates of 15.1 and 6.7
M$\rm{}_{\odot}\,yr^{-1}$ by subtracting the rotation curve from the CO(2-1)
datacube and then considering high-velocity gas only (faster than $\pm$300 km
s-1 in Scenario iii and just the high-velocity wings shown in Figure 8 in
Scenario iv).

Based on the results presented in Section 3, aided by the simulations
discussed in Section 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 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$\times$107 M☉, and the
outflow mass rate within 15-41 M$\rm{}_{\odot}\,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,
Pjet/L${}_{\rm Edd}>$10-4 (Wagner et al. 2012). In the Teacup, this ratio is
$\sim$9$\times$10-4, and the jet power is enough to drive the molecular
outflow because it is about two or three orders of magnitude higher than the
outflow kinetic power (see Table 1). However, the fate of the molecular gas
during this feedback episode is a galactic fountain, since the escape velocity
of the galaxy v${}_{\rm esc}\gtrsim$530 km s-1, estimated from its dynamical
mass as in Feruglio et al. (2020), is larger than the fastest velocities in
the outflow ($\sim\pm$400 km s-1). Regardless, this will have a global effect
on the galaxy’s evolution by redistributing mass, metals, and delaying further
star-formation.

Based on the results presented here, together with the predictions from the
models described in Section 4.2, jet-induced turbulence and shocks are the
most likely mechanisms to explain the high values of $\sigma$ and reported
here for the cold molecular gas in the Teacup. Larger samples of radio-quiet
quasars showing low-power radio jets are required to better quantify the
impact of jet-ISM coupling in the kinematics of the molecular gas.

###### Acknowledgements.

We are very grateful to the referee for very helpful and detailed comments,
which improved the presentation of the paper. AA and CRA acknowledge the
projects “Quantifying the impact of quasar feedback on galaxy evolution”, with
reference EUR2020-112266, funded by MICINN-AEI/10.13039/501100011033 and the
European Union NextGenerationEU/PRTR, and the Consejería de Economía,
Conocimiento y Empleo del Gobierno de Canarias and the European Regional
Development Fund (ERDF) under grant “Quasar feedback and molecular gas
reservoirs”, with reference ProID2020010105, ACCISI/FEDER, UE. AA acknowledges
funding from the MICINN (Spain) through the Juan de la Cierva-Formación
program under contract FJC2020-046224-I. SGB acknowledges support from the
research project PID2019-106027GA-C44 of the Spanish Ministerio de Ciencia e
Innovación. This paper makes use of the following ALMA data:
ADS/JAO.ALMA#2018.1.00870.S, and ADS/JAO.ALMA#2016.1.01535.S. ALMA is a
partnership of ESO (representing its member states), NSF (USA) and NINS
(Japan), together with NRC(Canada) and NSC and ASIAA (Taiwan), in cooperation
with the Republic of Chile. The Joint ALMA Observatory is operated by ESO,
AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of
the National Science Foundation operated under cooperative agreement by
Associated Universities, Inc.

## References

  * Balmaverde et al. (2022) Balmaverde, B., Capetti, A., Baldi, R. D., et al. 2022, A&A, 662, A23
  * Balmaverde et al. (2019) Balmaverde, B., Capetti, A., Marconi, A., et al. 2019, A&A, 632, A124
  * Bîrzan et al. (2008) Bîrzan, L., McNamara, B. R., Nulsen, P. E. J., Carilli, C. L., & Wise, M. W. 2008, ApJ, 686, 859
  * Carilli & Walter (2013) Carilli, C. L. & Walter, F. 2013, ARA&A, 51, 105
  * Cavagnolo et al. (2010) Cavagnolo, K. W., McNamara, B. R., Nulsen, P. E. J., et al. 2010, ApJ, 720, 1066
  * Cazzoli et al. (2022) Cazzoli, S., Hermosa Muñoz, L., Márquez, I., et al. 2022, A&A, 664, A135
  * Cielo et al. (2018) Cielo, S., Bieri, R., Volonteri, M., Wagner, A. Y., & Dubois, Y. 2018, MNRAS, 477, 1336
  * Couto et al. (2017) Couto, G. S., Storchi-Bergmann, T., & Schnorr-Müller, A. 2017, MNRAS, 469, 1573
  * Croton et al. (2006) Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11
  * Dasyra et al. (2016) Dasyra, K. M., Combes, F., Oosterloo, T., et al. 2016, A&A, 595, L7
  * Di Teodoro & Fraternali (2015) Di Teodoro, E. M. & Fraternali, F. 2015, MNRAS, 451, 3021
  * Downes & Solomon (1998) Downes, D. & Solomon, P. M. 1998, ApJ, 507, 615
  * Dubois et al. (2016) Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948
  * Feruglio et al. (2020) Feruglio, C., Fabbiano, G., Bischetti, M., et al. 2020, ApJ, 890, 29
  * Fiore et al. (2017) Fiore, F., Feruglio, C., Shankar, F., et al. 2017, A&A, 601, A143
  * Fluetsch et al. (2019) Fluetsch, A., Maiolino, R., Carniani, S., et al. 2019, MNRAS, 483, 4586
  * García-Burillo et al. (2019) García-Burillo, S., Combes, F., Ramos Almeida, C., et al. 2019, A&A, 632, A61
  * Girdhar et al. (2022) Girdhar, A., Harrison, C. M., Mainieri, V., et al. 2022, MNRAS, 512, 1608
  * Godfrey & Shabala (2016) Godfrey, L. E. H. & Shabala, S. S. 2016, MNRAS, 456, 1172
  * Harrison et al. (2014) Harrison, C. M., Alexander, D. M., Mullaney, J. R., & Swinbank, A. M. 2014, MNRAS, 441, 3306
  * Harrison et al. (2015) Harrison, C. M., Thomson, A. P., Alexander, D. M., et al. 2015, ApJ, 800, 45
  * Heckman et al. (1984) Heckman, T. M., van Breugel, W. J. M., & Miley, G. K. 1984, ApJ, 286, 509
  * Jarvis et al. (2019) Jarvis, M. E., Harrison, C. M., Thomson, A. P., et al. 2019, MNRAS, 485, 2710
  * Keel et al. (2012) Keel, W. C., Chojnowski, S. D., Bennert, V. N., et al. 2012, MNRAS, 420, 878
  * Keel et al. (2015) Keel, W. C., Maksym, W. P., Bennert, V. N., et al. 2015, AJ, 149, 155
  * Leroy et al. (2009) Leroy, A. K., Walter, F., Bigiel, F., et al. 2009, AJ, 137, 4670
  * McMullin et al. (2007) McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127
  * Meenakshi et al. (2022) Meenakshi, M., Mukherjee, D., Wagner, A. Y., et al. 2022, MNRAS[arXiv:2203.10251]
  * Morganti et al. (2021) Morganti, R., Oosterloo, T., Murthy, S., & Tadhunter, C. 2021, Astronomische Nachrichten, 342, 1135
  * Morganti et al. (2015) Morganti, R., Oosterloo, T., Oonk, J. B. R., Frieswijk, W., & Tadhunter, C. 2015, A&A, 580, A1
  * Mukherjee et al. (2016) Mukherjee, D., Bicknell, G. V., Sutherland, R., & Wagner, A. 2016, MNRAS, 461, 967
  * Mukherjee et al. (2018a) Mukherjee, D., Bicknell, G. V., Wagner, A. Y., Sutherland, R. S., & Silk, J. 2018a, MNRAS, 479, 5544
  * Mukherjee et al. (2018b) Mukherjee, D., Wagner, A. Y., Bicknell, G. V., et al. 2018b, MNRAS, 476, 80
  * Murthy et al. (2022) Murthy, S., Morganti, R., Wagner, A. Y., et al. 2022, Nature Astronomy, 6, 488
  * Nelson et al. (2019) Nelson, D., Pillepich, A., Springel, V., et al. 2019, MNRAS, 490, 3234
  * Nyland et al. (2018) Nyland, K., Harwood, J. J., Mukherjee, D., et al. 2018, ApJ, 859, 23
  * Oosterloo et al. (2019) Oosterloo, T., Morganti, R., Tadhunter, C., et al. 2019, A&A, 632, A66
  * Oosterloo et al. (2017) Oosterloo, T., Raymond Oonk, J. B., Morganti, R., et al. 2017, A&A, 608, A38
  * Orr & Browne (1982) Orr, M. J. L. & Browne, I. W. A. 1982, MNRAS, 200, 1067
  * Ramos Almeida et al. (2022) Ramos Almeida, C., Bischetti, M., García-Burillo, S., et al. 2022, A&A, 658, A155
  * Ramos Almeida et al. (2017) Ramos Almeida, C., Piqueras López, J., Villar-Martín, M., & Bessiere, P. S. 2017, MNRAS, 470, 964
  * Rupke et al. (2005) Rupke, D. S., Veilleux, S., & Sanders, D. B. 2005, ApJS, 160, 115
  * Solomon & Vanden Bout (2005) Solomon, P. M. & Vanden Bout, P. A. 2005, ARA&A, 43, 677
  * Speranza et al. (2022) Speranza, G., Ramos Almeida, C., Acosta-Pulido, J. A., et al. 2022, A&A, 665, A55
  * Tadhunter et al. (2014) Tadhunter, C., Morganti, R., Rose, M., Oonk, J. B. R., & Oosterloo, T. 2014, Nature, 511, 440
  * Talbot et al. (2022) Talbot, R. Y., Sijacki, D., & Bourne, M. A. 2022, MNRAS, 514, 4535
  * Venturi et al. (2021) Venturi, G., Cresci, G., Marconi, A., et al. 2021, A&A, 648, A17
  * Villar Martín et al. (2014) Villar Martín, M., Emonts, B., Humphrey, A., Cabrera Lavers, A., & Binette, L. 2014, MNRAS, 440, 3202
  * Viti et al. (2014) Viti, S., García-Burillo, S., Fuente, A., et al. 2014, A&A, 570, A28
  * Wagner & Bicknell (2011) Wagner, A. Y. & Bicknell, G. V. 2011, ApJ, 728, 29
  * Wagner et al. (2012) Wagner, A. Y., Bicknell, G. V., & Umemura, M. 2012, ApJ, 757, 136

##  Appendix A Moment maps of the CO(3-2) emission

We computed the moment maps of the CO(3-2) emission at
$\sim$0$\aas@@fstack{\prime\prime}$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 4.

Figure 4: Moment maps of the CO(3-2) emission of the Teacup. Top, middle, and
bottom panels correspond to the integrated intensity (moment 0, in Jy beam-1
km s-1 units), intensity weighted velocity field (moment 1, in km s-1), and
$\sigma$ (moment 2, in km s-1), respectively. East is to the left and north to
the top. In the $\sigma$ map, the contours show the of Figure 3 at 0.5, 0.6,
and 0.7, and the red ellipse in the bottom left corner indicates the beam
size.

##  Appendix B Simulation of jet propagating in the molecular disc

Figure 5: Snapshots of the mean velocity field and PVDs from Sim E (top
panels, jet inclined 20∘ relative to the disc) and Sim D (bottom panels, jet
inclined at 45∘) for $P_{\rm jet}$=1045 erg s-1 at 1.24 Myr and 1 Myr,
respectively. Left panels: mean LOS velocity fields with the position of the
slits of 500 pc width along (dashed lines) and perpendicular (solid lines) to
the jet. Middle and right panels: PVDs along and perpendicular to the jet. The
green curves are the mean velocity curve with $\pm$2$\sigma$ deviation, with a
maximum deviations along and perpendicular to the jet reaching 408 (448) km
s-1 and 434 (412) km s-1, respectively, for Sim E (Sim D).

Figure 6: $\sigma$ maps from Sim E (top panels) and D (bottom panels). Left
panels: maps at a resolution of 100 pc. Right panels: simulated maps and jet-
tracer smoothed by convolving with a 2D Gaussian of width 100 pc and 60 pc,
respectively. The jet tracer is shown at 0.05 contour level.

In order to elucidate the nature of 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}\leavevmode\nobreak\ \mbox{erg s}^{-1}$ are launched at
angles of $45^{\circ}$ and $20^{\circ}$ 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 $\sim$2 kpc is
placed in the $X-Y$ plane. A pair of relativistic jets are launched from the
centre. The jet’s axis lie in the $X-Z$ plane, inclined from the disc at
different angles of choice. According to the simulation’s conventions,
$\theta_{\rm I}$ = 90∘, $\Phi_{\rm 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 $\theta_{\rm I}=52^{\circ}$ to match that of the Teacup. The
azimuthal angle $\phi_{\rm 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 $\mbox{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 5, the top left panel shows the projected velocity map for
Sim E at a $\theta_{\rm I}=52^{\circ}$ and $\phi_{\rm I}=340^{\circ}$,
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 $\sim\pm$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. 5, we show
the mean velocity field and PVDs for Sim D, which are produced for an image
plane oriented at $\theta_{\rm I}=52^{\circ}$ and $\phi_{\rm I}=30^{\circ}$.
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 5 also exhibit
velocities of up to $\sim\pm$400 km s-1 within the central 2 kpc.

The simulated $\sigma$ 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 6. The $\sigma$ 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 $\sigma$, shown in Figure 3. The $\sigma$
map from Sim D (see bottom panels of Fig. 6) 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 8 in Meenakshi et al. 2022). This indicates that the jet-gas
interaction in the nuclear regions of the disc can lead to high-velocity
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 outflow mass rates, from the least to the most
conservative. For this analysis, we only used the CO(2-1) data at the original
0$\aas@@fstack{\prime\prime}$2 resolution to avoid artefacts or wrong
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 just
integrated the emission the 3DBAROLO 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$\times$108 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\sim$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$\sigma_{\rm rms}$. Thus, there is some extended low-velocity gas that is not
accounted for by our rotating model, and some nuclear high-velocity gas that
is considered by 3DBAROLO as rotation while it might not be. This uncertainty
is not considered in the errors listed in Table 1. In this scenario, we
adopted an outflow velocity of $v_{\rm out}$=100 km s-1, estimated from the
mean velocity residuals from the 3DBAROLO fit in  RA22.

ii) In the second scenario, we consider that only the highest velocities,
faster than $\pm 300$ km s-1, correspond to outflowing gas. We chose this
value from the comparison of the PVDs shown in Figure 2 and the 3DBAROLO
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
$\pm$300 km s-1, shown in Figure 7. In this scenario, the mass corresponding
to the outflow is 1.6$\times$108 M⊙.

Figure 7: Integrated intensity of the high velocity components of the CO(2-1)
emission at 0$\aas@@fstack{\prime\prime}$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$=$\pm$300 km s-1 are shown as red
and blue contours at (0.35,0.5,0.65,0.8,0.95)$\times\sigma_{\rm max}$, with
$\sigma_{\rm max}^{\rm red}$= 0.28 Jy beam-1 km s-1 and $\sigma_{\rm max}^{\rm
blue}$=0.14 Jy beam-1 km s-1. The beam size is indicated with a red ellipse in
the bottom left corner.

iii) The third scenario is a mix of the first two. 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 re-shuffled 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 8. Since this method is suited to optimize
the signal-to-noise of the emission associated to non-circular motions, it can
reveal high-velocity wings that otherwise are too faint to be detected. Using
the residual profile and considering only gas with velocities above $\pm$300
km s-1, we derive a mass of 5.7$\times$107 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 8,
corresponds to turbulence unrelated to the outflow. This turbulence can be due
to the merger, or produced by the jet, but here consider that it does not
correspond to outflowing gas. The residual profile can be fitted with a single
Gaussian of FWHM$\sim$330 km s-1, as shown in Figure 8. If we subtract this
Gaussian from the total residual profile, we measure a mass of 2.6$\times$107
M⊙, which mainly corresponds to the high-velocity wings.

Figure 8: Integrated CO(2-1) residual profile after subtraction of the
rotation curve. High-velocity blue and red wings are detected after
subtraction. The cyan solid line corresponds to the fit with a single Gaussian
component of FWHM$\sim$330 km s-1 used in Scenario iv. The blue and red shaded
areas correspond to the high-velocity gas considered to compute the outflow
flux in Scenario iii.

We computed the outflow masses from the integrated flux measurements in the
outflow. First, the integrated fluxes are converted to CO luminosities
L${}_{\rm CO(2-1)}^{\prime}$, in units of ${\rm K\,km\,s^{-1}\,pc^{2}}$ using
Equation 3 of Solomon & Vanden Bout (2005) and then translated into masses
using the CO-to-H2 ($\alpha_{\rm CO}$) conversion factor M${}_{\rm
H_{2}}$=$\alpha_{\rm CO}$R12L${}^{\prime}_{\rm CO(2-1)}$. Under the assumption
that the gas is thermalised and optically thick, the brightness temperature
ratio $R_{\rm 12}={L_{\rm CO(1-0)}^{\prime}}/{L_{\rm CO(2-1)}^{\prime}}$=1. A
conservative $\alpha$=0.8 is usually adopted for deriving outflow masses (see
RA22 and references therein). The values for each scenario are listed in Table
1.

Table 1: Outflow measurements for the four scenarios proposed for the Teacup.

Scenario | $S_{\rm CO}$ | $M_{\rm out}$ | $\dot{M}_{\rm out}$ | $M_{\rm out}/M_{\rm H_{2}}$ | $\dot{E}_{\rm kin}/P_{\rm jet}$
---|---|---|---|---|---
|  (Jy km s-1) |  (107 M⊙) |  ( M$\rm{}_{\odot}\,yr^{-1}$) | (%) |
i |  7.6$\pm$0.9 |  51.5$\pm$48.2 |  44.0 $\pm$ 32.7† | 8.3 | 0.004
ii |  2.4$\pm$0.1 |  16.0$\pm$10.9 |  41.0$\pm$27.9 | 2.6 | 0.035
iii |  0.9$\pm$0.1 |  5.6$\pm$4.3 |  15.1$\pm$11.0 | 0.9 | 0.013
iv |  0.4$\pm$0.3 |  2.6$\pm$3.4 |  6.7$\pm$8.7 | 0.4 | 0.006

111The uncertainty in $\alpha_{\rm CO}=$0.8$\pm$0.5 M⊙(K km s-1 pc2)-1 (Downes
& Solomon 1998) is included in the mass error estimates. ${\dagger}$: computed
using $v_{\rm out}$=100 km s-1, while for the others scenarios we adopted
$v_{\rm out}$=300 km s-1. The outflow mass fractions were computed using the
total mass of $M_{\rm H_{2}}$=6.2$\times 10^{9}$ M⊙ from  RA22, estimated with
$\alpha_{\rm CO}$=4.35 M⊙(K km s-1 pc2)-1.

To compute outflow mass rates, we use an outflow radius of $r_{\rm
out}$=0$\aas@@fstack{\prime\prime}$75 (1.2 kpc), as constrained from the
extent of the high-velocity gas with high $\sigma$ and (see Figures 2 and 3).
The outflow velocities reach up to $\pm$400 km s-1, but here we adopt a more
conservative velocity of $v_{\rm out}$=300 km s-1, except for Scenario i. For
a time-averaged thin expelled shell geometry (Rupke et al. 2005),
$\dot{M}_{\rm out}$=$v_{\rm out}\left(M_{\rm out}/r_{\rm out}\right)$, which
corresponds to the outflow mass averaged over the flow timescale, $t_{\rm
flow}={\rm R_{\rm out}}/{\rm v_{\rm out}}$=3.9 Myr. Considering all the
previous, the outflow mass rates are in the range $\dot{M}_{\rm out}$=6.7$-$44
M$\rm{}_{\odot}\,yr^{-1}$. The ratio between the kinetic power of the
molecular outflow ($\dot{E}_{\rm kin}=0.5\dot{M}_{\rm out}{v_{\rm out}}^{2}$)
and the jet power varies from 0.004$<\dot{E}_{\rm kin}/P_{\rm jet}<$0.035 (see
Table 1), indicating that the jet has enough kinetic energy to drive the
molecular outflow.