# HARQ Retransmissions in C-V2X: A BSM Latency Analysis

Abdurrahman Fouda†, Randall Berry† and Ivan Vukovic‡ This project was
supported in part by the Ford-Northwestern University Alliance. †Department of
Electrical and Computer Engineering, Northwestern University, Evanston, IL
‡Ford Motor Company, Dearbon, MI
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Cellular vehicular-to-everything (C-V2X) systems offer the potential for
improving road safety, in part through the exchange of periodic basic safety
messages (BSMs) between nearby vehicles. The reliability and latency of these
messages is a key metric. Hybrid automatic repeat request (HARQ)
retransmissions are one technique used to this end. However, HARQ may come at
the expense of consuming the limited available wireless resources, especially
in highly congested scenarios. This paper studies BSM transmission latency and
reliability when HARQ retransmissions are used with the semi-persistent
scheduling (SPS) in C-V2X transmission mode 4. We do so through extensive
system-level simulations that closely follow the SPS process. Furthermore, we
provide an analytical model for the tail behavior of the BSM latency
distribution with HARQ retransmissions that is a good approximation to the
simulation results. Our study reveals the impact of several deployment
settings (e.g., bandwidth configurations and vehicle density).

##  I Introduction

Cellular vehicle-to-everything (C-V2X) was introduced by 3GPP to support
communication between vehicular users (VUEs) and between VUEs and
infrastructure [1]. The necessary LTE enhancements to support C-V2X were a
study item in Rel-15 [2] and continued in Rel-16 for new radio (NR)
enhancements [3]. 3GPP adopted two transmission modes (mode 3 and mode 4) to
support C-V2X (also known as LTE-V2X) with and without assistance from base
stations (eNBs). This paper focuses on mode 4, which uses a decentralized
approach for managing radio resources without assistance from eNBs called
semi-persistent scheduling (SPS) [4]. Application layer standards for C-V2X
have been defined by SAE, including the definition of basic safety messages
(BSM) for periodic updates exchanged between nearby VUEs [5, 6].

C-V2X supports hybrid automatic repeat request (HARQ) retransmission to
improve the packet reception ratio (PRR) and and Inter-packet gap (IPG) of BSM
transmissions [7]. However, HARQ retransmissions consume more resources, which
may impact performance in high VUE density scenerios. Our goal in this paper
is to study the impact of HARQ in terms of PRR and IPG.

The IPG performance of SPS scheduling has been studied in [8] and [9], where
the role of one-shot transmissions has also been highlighted. The IPG and PRR
performance of C-V2X in dense highway scenarios was compared with that of the
DSRC (another vehicular communication protocol) in [10]. An alternative method
to decrease the probability of losing consecutive BSMs has been discussed in
[11]. The performance of a blind retransmission method has been evaluated in
[12] for C-V2X with varying traffic speeds and for NR V2X in [13]. However, no
previous literature has studied the impact of HARQ retransmissions on the IPG
distribution in C-V2X as we do here.

This paper considers studying the tail distribution of IPG CCDFs when HARQ
retransmissions are used with the sensing-based SPS and one-shot reselections
for BSM transmission in C-V2X networks. IPG and PRR performances are evaluated
through extensive simulation campaigns using a C++ system-level simulator that
closely follows the C-V2X transmission mode 4 process. Our analysis reveals
that the best IPG tail improvement is achieved when HARQ retransmissions and
one-shot reselections are used at small vehicle-to-vehicle (V2V) distances
with high VUE densities and at any V2V distance with moderate and low
densities. Also, it is shown that using HARQ retransmissions and one-shot
reselections yields the best PRR performance when the system is noise-limited.
Furthermore, we propose an analytical model to characterize the the IPG tail
when HARQ retransmissions are used. This builds on our prior work in [4] which
did not include HARQ retransmissions.

##  II SPS with one-shot reselection and HARQ

We consider a setting where SPS interleaved with one-shot reselections is used
to allocate virtual radio resource blocks (VRBs) for BSM transmissions. VRBs
are autonomously (re)selected from a candidate list which is defined using a
pre-configured resource pool of time-frequency resources.111More precisely, we
consider a setting where each VRB consists of 10 physical resource blocks
(PRBs) in LTE which occupy a 1 msec subframe. The candidate VRBs that can be
selected depend in part on the measured signal strength in those VRBs. The VRB
reselection process is controlled by the SPS and the one-shot reselection
counters $C_{\mathrm{s}}$ and $C_{\mathrm{o}}$, which are chosen uniformly at
random between $\left[\alpha_{\mathrm{s}},\,\beta_{\mathrm{s}}\right]$ and
$\left[\alpha_{\mathrm{o}},\,\beta_{\mathrm{o}}\right]$, respectively. These
counters decrease by one for each transmissions and VRBs are only reselected
when the counter reaches zero. Namely, when $C_{\mathrm{s}}$ expires, VRBs are
reselected with a reselection probability $p_{\mathrm{r}}$ (see [4] for
further discussion). When $C_{\mathrm{o}}$ expires, a VUE reselects a VRB for
a one-shot transmission, and then returns to its previous VRB.

C-V2X mode 4 supports up to one redundant retransmission using HARQ to improve
the signal-to-interference-plus-noise ratio received (SINR) and decrease
packet losses [7]. VRBs for HARQ transmissions are selected from the same
candidate list such that the subframe of HARQ VRBs lies within 15 subframes
from the subframe of the initially selected set of VRBs. Whenever a vehicle
reselects a VRB (due to either $C_{\mathrm{s}}$ or $C_{\mathrm{0}}$ expiring),
it also reselects the VRB used for HARQ transmissions.

Given that BSMs are transmitted in a broadcast scenario, the HARQ
retransmission of a BSM is transmitted regardless of whether its first attempt
is received. Similar to the first transmission, the first two physical
resource blocks (PRBs) of a HARQ retransmission are reserved for the sidelink
control information (SCI) and the remaining PRBs (within two contiguous VRBs)
are reserved for the data packets. For each two transmitting attempts of the
BSM, the receiving VUE implements HARQ combining for the received SCIs and
data packets. Note that a BSM is successfully received if the SINR received
from either of the two transmission attempts or the combined is above the SINR
thresholds shown in Table I.

We use the SAE congestion control specified in [6], [5] to reduce the
probability of packet collision in highly congested scenarios. This is done by
adjusting the BSM generation rate based on the estimated number of nearby
transmitting vehicles. When congestion control is not active, BSMs are
generated once every 100 msec, which is also the underlying frame size from
which VRBs are selected. When the estimated number of vehicles within 100
meters exceeds 25, the BSM rate linearly reduces until it reaches a minimum
value of one BSM every 600 msec. In the simulations presented, congestion
control is not active when the VUE density is 125 VUE/km and is active for
higher densities, with the time between BSM having it maximum value of 600
msec when the density is 800 VUE/km.

##  III Simulation Results

In this section, we study the IPG distribution and PRR performance when HARQ
retransmissions are used in C-V2X systems. IPG is measured as the time elapsed
between each successive pair of successfully received BSMs between a pair of
vehicles. The value of the IPG CCDF at $i$ milliseconds and a V2V distance of
$d$ meters represents the fraction of IPG instances between vehicles separated
by $d$ meters that exceed $i$ milliseconds. PRR at a V2V distance of $d$
meters quantifies the transmission reliability by calculating the ratio
between the total number of BSM transmitted and successfully received between
any pairs of VUEs that are separated by $d$ meters.

IPG and PRR are analyzed using different deployment settings, including
different VUE densities, one-shot settings, V2V distances, and bandwidth
configurations. We use a high-fidelity system-level simulator that closely
follows the sensing-based SPS for BSM transmission using C-V2X transmission
mode 4. A single-lane highway scenario is used to carry out Monte Carlo
simulations. Numerical statistics are collected from vehicles located only in
the middle third of the highway to minimize the boundary effects. The
simulation settings are summarized in Table I.

TABLE I: Simulation parameters. Parameter | Value
---|---
Scenario | Layout: 5 km single-lane Highway
| Density: 125, 400 and 800 VUE/km
Channel model | Path loss: ITU-R UHF urban canyon [14]
| In-band emission: 3GPP TR 36.885 [2]
| Fast fading: Nakagami-$m$ distribution
Antenna settings |  $n_{t}=1$, $n_{r}=2$, MRC receiver
|  $C_{\mathrm{s}}\in{[5,\,15]}$, $p_{\mathrm{k}}=0.8$
Resource pool settings | Carrier bandwidth: 10, 20 MHz at 5.9 GHz
Power settings | Tx power: 20 dBm, Noise figure: 6 dB
| PSCCH power boost: 3dB
SINR threshold | PSSCH: 3 dB, PSCCH: 0 dB
HARQ settings | Maximum retransmissions: 1
| HARQ selection window: [-15, 15] ms
Simulation time | 500 seconds, warm-up: first 10 seconds.

(a)

(b)

(c)

Figure 1: 10 MHz, V2V distance = 300 m, with and without HARQ retransmissions.

###  III-A IPG Performance

Fig. 1 (a) compares the IPG tail distribution of the three configurations of
$C_{\mathrm{o}}$, OFF, $[2,\,6]$, and $[5,\,15]$ with and without HARQ
retransmissions at a 10 MHz bandwidth for a high VUE density of 800 VUE/km.
Figs. 1 (b) and (c) show similar comparisons for moderate and low VUE
densities of 400 and 125 VUE/km, respectively. As shown in Fig. 1 (a), the IPG
tails suffer substantial performance degradation when HARQ retransmissions are
used. This applies to all three configurations of $C_{\mathrm{o}}$.
Specifically, using the one-shot method (either $[2,\,6]$ or $[5,\,15]$
configuration) with the HARQ retransmissions does not improve the IPG tail
performance. In fact, deactivating the HARQ retransmissions yields a slightly
better tail performance even without using the one-shot method. This is due to
the scarcity of VRBs available for BSM transmissions. Essentially, with such a
highly congested scenario (800 VUEs/km) and low operating bandwidth (10 MHz),
it is expected that the number of available VRBs in the SPS (re)selection will
be significantly low. As shown, using the one-shot setting of $[2,\,6]$
without HARQ retransmissions produces the best IPG tail performance.

At lower VUE densities, as shown in Figs. 1 (b) and (c), using the HARQ
retransmissions improves the IPG tail distribution in all configurations of
$C_{\mathrm{o}}$. Furthermore, the improvement is better with the lower VUE
densities (i.e. Fig. 1 (c) vs. (b)). The statistics shown in Fig. 1 are
generated for a V2V distance of 300 m. Similar trends were observed for other
V2V distances. This suggests that at high VUE densities that HARQ
retransmissions may not be productive.

(a)

(b)

(c)

Figure 2: 20 MHz, V2V distance = 200 m, with and without HARQ retransmissions.

However, using HARQ retransmissions with higher bandwidth configurations
(e.g., 20 MHz) yields significantly better IPG tail performance, as shown in
Fig. 2. It can also be seen that the improvement is greater at lower VUE
densities. Fig. 2 (a) shows that using the one-shot configuration of $[2,\,6]$
without HARQ retransmission significantly improves the IPG tail performance
compared to the OFF and $[5,\,15]$ configurations of $C_{\mathrm{o}}$ with
HARQ retransmissions. However, using the same configuration (i.e. $[2,\,6]$)
with HARQ yields a slightly better tail. At lower densities (Fig. 2 (b) and
(c)), the gains are even more apparent.

Figure 3: IPG at different V2V densities, 800 VUE/km.

Fig. 3 (a) and (b) depict the effects of using HARQ retransmissions with 800
VUE/km at larger V2V distances of 300 and 500 m, respectively. This shows that
even with 20Mhz bandwidth, using the HARQ retransmissions may worsen the IPG
tail at larger V2V separation. As shown, the one-shot configuration of
$[2,\,6]$ without HARQ retransmissions yields better IPG tail performance than
any other configuration of the $C_{\mathrm{o}}$ counter when HARQ
retransmissions are used. This shows that the advantages of HARQ depend on
both the distance, the bandwidth configuration and the underlying vehicle
density.

###  III-B PRR Performance

(a)

(b)

(c)

Figure 4: PRR performance at 20 MHz.

The PRR performance versus distance of different configurations of
$C_{\mathrm{o}}$ and HARQ using a bandwidth of 20 MHz is shown in Fig. 4. As
shown, HARQ improves the PRR at large Tx-Rx separations. The improvement is
more pronounced at higher VUE densities. This happens because the system is
noise-limited at high Tx-Rx separations. On the other hand, Fig. 4 (a) reveals
that using HARQ retransmissions degrades the PRR performance at moderate Tx-Rx
separations, where the system is interference-limited.

Essentially, BSM transmissions suffer from a higher probability of packet
collision at high VUE densities because of the lack of available VRB
resources. Degradation is less pronounced at moderate and light VUE density in
Figs. 4 (b) and (c), respectively. It is worth noting that the use of HARQ
retransmissions in fact improved the PRR performance at small Tx-Rx
separations. In other words, the best PRR performance gains from HARQ are with
large Tx-Rx separations at fairly high VUE densities and with small Tx-Rx
separations at low VUE densities.

Similar trends were observed in a 10-MHz bandwidth configuration as shown in
Fig. 5. Figs. 5 (b) and (c) reveal that the PRR performance degrades at
moderate Tx-Rx separations with moderate and low VUE densities due to the
limited number of available VRB resources compared to the 20 MHz case. Fig. 5
(a) also shows that using HARQ retransmissions with the 10 MHz bandwidth at
high VUE densities barely improves the PRR performance at large Tx-Rx
separations. This is because, in the 10MHz case, the system is interference-
limited at small and large Tx-Rx separations.

(a)

(b)

(c)

Figure 5: PRR performance at 10 MHz.

##  IV Analytical Model for IPG Tail Distribution

In this section, we present an analytical model to approximate the IPG tail
distribution with HARQ for C-V2X communications. Our approach is based on the
model in [4] which does not include HARQ.

Unless mentioned otherwise, a BSM transmission means two BSM transmission
attempts representing the regular transmission and the HARQ retransmission.
For a given transmitter-receiver pair, conditioned on the fact that the first
BSM transmission is not received, let $T$ be a random variable that indicates
the length of IPG in terms of the required number of BSM transmission
opportunities until a successful transmission. The failure of a BSM
transmission means that the receiver did not receive the regular transmission
of the BSM and its HARQ retransmission.

Here, we assume that if a BSM is not received, it is due to interference from
a first interferer that is using overlapping VRB(s) with the VRB(s) of the
first transmission attempt and from a second interferer that is using
overlapping VRB(s) with the VRB(s) of the HARQ retransmission. In other words,
interference caused by two dominant interferer (one for each transmission) is
the main reason for packet losses. Other reasons for packet losses (e.g.,
half-duplex transmissions) can be fairly ignored on the basis of numerical
data collected from our simulations. In this case, under the SPS, the
interference will persist until either one of the interferers or the
transmitter reselects another resource or does a one-shot transmission.

We give a non-uniform geometric model for the SPS and the one-shot reselection
processes, in which, the probability of a resource reselection during the
$k^{\mathrm{th}}$ BSM transmission opportunity is given by
$q^{(k)}(\rho,\sigma,p)$, where
$q_{\mathrm{s}}^{(k)}=q^{(k)}\left(\alpha_{\mathrm{s}},\,\beta_{\mathrm{s}},\,p_{\mathrm{r}}\right)$
and
$q_{\mathrm{o}}^{(k)}=q^{(k)}\left(\alpha_{\mathrm{o}},\,\beta_{\mathrm{o}},\,1\right)$
for the SPS-triggered and one-shot-triggered reselection, respectively. Given
that there cannot be reselections in the first $\rho$ BSM transmission
opportunities, we set $q^{(k)}=0$ for $k\leq\rho$, and for $k>\rho$, we set
$q^{(k)}$ so that for large $k$, the expected number of reselections after $k$
opportunities is equal to $k/E(T_{r})$ where $E(T_{r})$ denotes the expected
time between reselections as:

$E(T_{r})=\left(\rho+\frac{\sigma-\rho}{2}\right)\frac{1}{p_{s}}.$ (1)

In this case, $q^{(k)}$ is selected to satisfy
$q^{(k)}(k-\rho)=\frac{2kp_{s}}{\sigma+\rho}$. Hence,
$q^{(k)}=\left(2kp_{s}\right)/\left((\sigma+\rho)(k-\rho)\right)$. The
reselection probability of the SPS interleaved with one-shot transmissions is
then:

$q_{\mathrm{\mathrm{i}}}^{(k)}=q_{\mathrm{o}}^{(k)}\left(1-q_{\mathrm{s}}^{(k)}\right)+q_{\mathrm{s}}^{(k)}\left(1-q_{\mathrm{o}}^{(k)}\right)+q_{\mathrm{s}}^{(k)}q_{\mathrm{o}}^{(k)},$
(2)

assuming that the two processes are independent.

Next, the probability of a successful BSM transmission (due to a reselection)
after $k$ transmission opportunities, where each $k^{\mathrm{th}}$ failure
represents the failure of the regular BSM transmission and its HARQ
retransmission, is calculated as:

$p_{\mathrm{u}}^{(k)}=P_{\mathrm{a}}P_{\mathrm{d}}+P_{\mathrm{n}}P_{\mathrm{r}}+P_{\mathrm{a}}P_{\mathrm{r}}.$
(3)

Here, $P_{\mathrm{a}}$ denotes the probability that at least one of the two
interferers reselects, $P_{\mathrm{d}}$ is the probability that the
transmitter does not reselect, $P_{\mathrm{n}}$ indicates the probability that
no interferer reselects new VRBs, and $P_{\mathrm{r}}$ is the probability that
a resource reselection is triggered at the transmitter.

In that, $P_{\mathrm{a}}$ can be calculated as:

$P_{\mathrm{a}}=1-\left(1-q_{\mathrm{i}}^{(k)}\right)^{2},$ (4)

where $\left(1-q_{\mathrm{i}}^{(k)}\right)$ represents the probability that a
single interferer does not reselect new VRBs when the reselection counters
expire. Therefore, $\left(1-q_{\mathrm{i}}^{(k)}\right)^{2}$ calculates the
probability that neither of the two dominant interferers reselects new VRBs.
Similarly, $P_{\mathrm{d}}$ is calculated as
$\left(1-q_{\mathrm{i}}^{(k)}\right)$ and $P_{\mathrm{r}}$ as
$\left(q_{\mathrm{i}}^{(k)}P_{\mathrm{f}}\right)$, where $P_{\mathrm{f}}$
represents the probability that the transmitter’s reselection results in a
successful transmission (i.e. it does not reselect into the VRB(s) as another
strong interferer). Note that the probability of successful transmission due
to one-shot reselection is the same as the average success probability of any
transmission (see [4] for further discussion). Finally, $P_{\mathrm{n}}$ can
be given by $P_{\mathrm{n}}=\left(1-q_{\mathrm{i}}^{(k)}\right)^{2}$.

The probability that the IPG exceeds $k$ BSM transmissions (including their
HARQ attempts) is then written as:

$P(T>k)=\prod^{k}_{i=1}\left(1-p_{\mathrm{u}}^{(i)}\right).$ (5)

The analytical model proposed for the tail distribution of the IPG CCDF can
then be validated by comparing its approximate slope (which is calculated
using (5)) with that obtained through simulations. Here, simulations are
carried out for $\approx\mathrm{10^{8}}$ samples to obtain numerical results
with a good level of confidence. Since $T$ is defined as conditioned on one
BSM not arriving, there will be a constant offset between $P(T>t)$ and the
probability that the IPG exceeds the same amount, which corresponds to the
probability that the first BSM does not arrive. Hence, (5) captures the slope
of the IPG CCDF but not the actual value.

Figs. 6 and 7 show the validation of the analytical model for 125 VUE/km using
10 MHz and 20 MHz bandwidths, respectively. Validation is done for the three
configurations of $C_{\mathrm{o}}$, OFF, $[2,\,6]$, and $[5,\,15]$. The dotted
and solid lines represent the slopes of the analytical and simulation results,
respectively. Note that $P_{\mathrm{f}}$ is calculated using Fig. 4 based on
the operating scenario and the desired V2V distance bin. Unless otherwise
mentioned, the validation presented in this section focuses on the V2V
separation of 200 m. Similar trends were observed at other V2V distances. As
shown, the IPG tail derived using (5) gives a reasonable approximation for the
actual (i.e., simulated) IPG slopes of different reselection modes at fairly
large IPG values (which is the main purpose of the proposed analytical model).

Figure 6: 125 VUE/km, 10 MHz. Figure 7: 125 VUE/km, 20 MHz.

For a high vehicle density of 800 VUE/km, using (4) leads to a noticeable gap
between the approximated and the simulated IPG tail distributions as shown in
Fig. 8 for a 10 MHz bandwidth (a similar trend was observed at the 20 MHz).
This is because, in highly congested scenarios and with HARQ retransmissions,
it is likely that each BSM transmission attempt suffers interference from more
than a single dominant interferer. When considering two dominant interferers
for each transmission attempt, $P_{\mathrm{a}}$ (for simplicity of
presentation, $q_{\mathrm{i}}^{(k)}$ is denoted as $q_{\mathrm{i},k}$ in the
rest of this paper) can be updated as:

$P_{\mathrm{a}}=1-\left(1-q_{\mathrm{i},k}^{2}\right)^{2}.$ (6)

Figure 8: 800 VUE/km, 10 MHz using (4) Figure 9: 800 VUE/km, 10 MHz. Figure
10: 800 VUE/km, 20 MHz.

The probability of a successful BSM transmission after $k$ transmission
opportunities can then be written as:

$\begin{split}p_{\mathrm{u}}^{(k)}=&\left(1-\left(1-q_{\mathrm{i},k}^{2}\right)^{2}\right)\left(1-q_{\mathrm{i},k}\right)+\\\
&\left(1-q_{\mathrm{i},k}^{2}\right)^{2}q_{\mathrm{i},k}P_{\mathrm{f}}+\left(1-\left(1-q_{\mathrm{i},k}^{2}\right)^{2}\right)q_{\mathrm{i},k}P_{\mathrm{f}},\end{split}$
(7)

where the first term represents the event that at least one interferer
reselects new VRBs and the transmitter continues using the same VRBs. The
second term denotes the event where neither one of the first and second
dominant interferers reselects new VRB(s) and the transmitter reselects with
probability $q_{\mathrm{i},k}$ and that reselection results in a successful
transmission with probability $P_{\mathrm{f}}$. The third term represents the
event where at least one of the dominant interferers and the transmitter
reselect new VRBs. By substituting (7) into (5) the approximated IPG slopes
now coincide much better with the simulated IPG tail distributions, as shown
in Figs. 9 and 10.222For VUE densities that result in a packet generation rate
that is not a multiple of 100 ms due to congestion control, then the
approximated slope in (5) needs to be updated to account for ”slippage
events”, where a transmitter may miss transmitting in a VRB due to the packet
generation rate not being matched to the VRB frequency. This is discussed in
[4] without HARQ and a similar approach can be adopted here.

##  V Conclusion

A high-fidelity system-level simulator is used to evaluate the latency of BSM
transmission (in terms of IPG) and reliability (in terms of PRR) when HARQ
retransmission are used with the SAE-based one-shot method for the sensing-
based SPS in C-V2X transmission mode 4. In addition, an analytical model is
presented to give a good approximation for the tail distribution of the IPG
CCDFs. The study reveals that the impact of HARQ retransmissions on the IPG
tail behavior and PRR performance are influenced by several deployment
settings, including, but not limited to, bandwidth configuration, vehicle
density, and Tx-Rx separation.

## References

  * [1] TSG SA, “Study on enhancement of 3GPP support for 5G V2X services,” 3GPP TR22.886 v16.2.0, Dec. 2021.
  * [2] TSG RAN, “Study on LTE-based V2X services,” 3GPP TR36.885 v14.0.0, Jul. 2016.
  * [3] ——, “Study on NR vehicle-to-everything (V2X),” 3GPP TR 38.885 v16.0.0, Mar. 2019.
  * [4] A. Fouda, R. Berry, and I. Vukovic, “Study of BSM inter-packet gap tails in C-V2X networks,” _submitted to IEEE Trans. Wireless Commun_ , May 2023. [Online]. Available: ResearchGateDOI:10.13140/RG.2.2.28607.69288
  * [5] SAE International, “On-board system requirements for LTE-V2X V2V safety communications,” SAE-J3161/1 DRAFT, Jun. 2021.
  * [6] ——, “On-board system requirements for V2V safety communications,” SAE-J2945/1, Mar. 2016.
  * [7] TSG RAN, “LTE; E-UTRA; medium control (MAC) protocol specification,” 3GPP TS 36.321 v16.3.0, Dec. 2020.
  * [8] A. Fouda, R. Berry, and I. Vukovic, “Interleaved one-shot semi-persistent scheduling for BSM transmissions in C-V2X networks,” in _Proc. IEEE Vehic, Netw. Conf. (VNC)_ , Nov. 2021, pp. 143–150.
  * [9] R. Gholmieh and H. I. Abbasi, “C-V2X (LTE-V2X) performance enhancement through SAE J3161/1 probabilistic one-shot transmissions,” in _Proc. IEEE Global Commun. Conf. (GLOBECOM)_ , Dec. 2021, pp. 1–6.
  * [10] H. I. Abbasi, R. Gholmieh, T. V. Nguyen, S. Patil, and J. Misener, “LTE-V2X (C-V2X) performance in congested highway scenarios,” in _Proc. IEEE Int. Conf. on Commun. (ICC)_ , May 2022, pp. 303–308.
  * [11] A. Bazzi, C. Campolo, A. Molinaro, A. O. Berthet, B. M. Masini, and A. Zanella, “On wireless blind spots in the C-V2X sidelink,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 8, pp. 9239–9243, Aug. 2020.
  * [12] D. Wang, R. R. Sattiraju, A. Qiu, S. Partani, and H. D. Schotten, “Effect of retransmissions on the performance of C-V2X communication for 5G,” in _Proc. IEEE Vehic. Technol.Conf. (VTC-Fall)_ , Feb. 2020.
  * [13] D. Wang, P. B. Mohite, Q. Zhou, A. Qiu, and H. D. Schotten, “Evaluating the impact of numerology and retransmission on 5G NR V2X communication at a system-level simulation,” _ArXiv e-prints_ , Jul. 2023. [Online]. Available: https://arxiv.org/abs/2307.14158.
  * [14] T. Shimizu, B. Cheng, H. Lu, and J. Kenney, “Comparative analysis of DSRC and LTE-V2X PC5 mode 4 with SAE congestion control,” in _Proc. IEEE Vehic, Netw. Conf. (VNC)_ , Dec. 2020, pp. 1–8.


# ALMA reveals a compact and massive molecular outflow driven by the young AGN
in a nearby ULIRG

Luke R. Holden1, Clive N. Tadhunter1, Anelise Audibert2,3, Tom Oosterloo4,5,
Cristina Ramos Almeida2,3, Raffaella Morganti4,5, Miguel Pereira-Santaella6,
and Isabella Lamperti7
1Department of Physics $\&$ Astronomy, University of Sheffield, S6 3TG
Sheffield, UK.
2Instituto de Astrofísica de Canarias, Calle Vía Láctea, s/n, 38205 La Laguna,
Tenerife, Spain.
3Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna,
Tenerife, Spain.
4ASTRON, the Netherlands Institute for Radio Astronomy, Oude Hoogeveensedijk
4, 7991 PD, Dwingeloo, The Netherlands.
5Kapteyn Astronomical Institute, University of Groningen, Postbus 800, 9700 AV
Groningen, The Netherlands.
6Instituto de Física Fundamental, CSIC, Calle Serrano 123, 28006 Madrid, Spain
7Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Ajalvir Km. 4, 28850
Torrejón de Ardoz, Madrid, Spain E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The ultra luminous infrared galaxy (ULIRG) F13451+1232 is an excellent example
of a galaxy merger in the early stages of active galactic nucleus (AGN)
activity, a phase in which AGN-driven outflows are expected to be particularly
important. However, previous observations have determined that the mass
outflow rates of the warm ionised and neutral gas phases in F13451+1232 are
relatively modest, and there has been no robust detection of molecular
outflows. Using high spatial resolution ALMA CO(1–0) observations, we detect a
kiloparsec-scale circumnuclear disk, as well as extended ($r\sim 440$ pc),
intermediate-velocity (300 < $|v|$ < 400 km s-1) cold molecular gas emission
that cannot be explained by rotational disk motions. If interpreted as AGN-
driven outflows, the mass outflow rates associated with this intermediate-
velocity gas are relatively modest ($\dot{M}_{\mathrm{out}}=22$–$27$ M⊙ yr-1);
however, we also detect a compact ($r_{\mathrm{out}}$ < 120 pc), high velocity
(400 < $v$ < 680 km s-1) cold molecular outflow near the primary nucleus of
F13451+1232, which carries an order of magnitude more mass
($\dot{M}_{\mathrm{out}}$$\sim$230 M⊙ yr-1) than (and several times the
kinetic power of) the previously-detected warmer phases. Moreover, the similar
spatial scales of this compact outflow and the radio structure indicate that
it is likely accelerated by the small-scale ($r$$\sim$130 pc) AGN jet in the
primary nucleus of F13451+1232. Considering the compactness of the nuclear
outflow and intermediate-velocity non-rotating gas that we detect, we argue
that high spatial-resolution observations are necessary to properly quantify
the properties of AGN-driven outflows and their impacts on host galaxies.

###### keywords:

galaxies: active – galaxies: evolution – galaxies: individual: F13451+1232 —
ISM: jets and outflows — galaxies: quasars: general — galaxies: interactions

††pubyear: 2023††pagerange: ALMA reveals a compact and massive molecular
outflow driven by the young AGN in a nearby ULIRG–B

##  1 Introduction

Although often required by models of galaxy evolution to explain the observed
properties of the galaxy population (e.g. Magorrian et al. 1998; Silk & Rees
1998; Di Matteo et al. 2005; Davé et al. 2019), considerable uncertainties
remain regarding the impact of gas outflows driven by active galactic nuclei
(AGN). Principally, this is because outflow masses and kinetic powers are
often not precisely determined (see Harrison et al. 2018), and their dominant
driving mechanisms (radiation pressure e.g. Hopkins & Elvis 2010; Meena et al.
2023, or jet-ISM interactions e.g. Mukherjee et al. 2016; Audibert et al.
2023) in different object types are unclear.

An excellent object to investigate the impact and launching mechanism(s) of
AGN-driven outflows is the ultraluminous infrared galaxy (ULIRG;
$L_{\mathrm{IR}}$ > $10^{12}$ L⊙: Sanders & Mirabel 1996) F13451+1232 — also
known as 4C12.50 — which is a merger (Gilmore & Shaw 1986; Heckman et al.
1986) in a late pre-coalescence phase (nuclear separation $\sim$4.5 kpc:
Tadhunter et al. 2018). ULIRGs such as F13451+1232 are among the most rapidly
evolving galaxies in the local universe due to merger-induced gas inflows
feeding the central AGN and triggering star formation — the exact situation
modelled in some hydrodynamic simulations of galaxy evolution that include
AGN-driven outflows as a feedback mechanism (e.g. Di Matteo et al. 2005;
Hopkins et al. 2008).

As a luminous radio source (L${}_{\mathrm{1.4\;GHz}}=1.9\times 10^{26}$ W
Hz-1), F13451+1232 contains radio emission on two distinct scales: small-scale
emission within a radius of $\sim$130 pc of the primary nucleus along position
angle PA=151∘ (Stanghellini et al., 1997; Lister et al., 2003; Morganti et
al., 2013), and larger-scale, diffuse emission extending to a radial distance
of $\sim$77 kpc from the primary nucleus along PA$\sim$180∘ (Stanghellini et
al., 2005). It has been proposed that this radio structure could be the result
of two epochs of jet activity (Stanghellini et al., 2005), with the smaller-
scale radio structure comprising a young, recently-restarted jet, and the
larger-scale emission originating from a previous jet cycle.

Because of the properties of its small-scale radio emission, F13451+1232 is
also classified as a gigahertz-peaked-spectrum (GPS) source (see O’Dea &
Saikia 2021 for a review); such sources are associated with powerful jet-ISM
interactions that may accelerate gas outflows in multiple phases (Mukherjee et
al., 2016; Holt et al., 2006; Santoro et al., 2018, 2020; Kukreti et al.,
2023). In addition, given the object’s high optical emission-line luminosity
(L${}_{\mathrm{[OIII]}}=1.45\times 10^{43}$ erg s-1: Rose et al. 2018) —
leading to its additional classification as a type-2 quasar
(QSO2)111F13451+1232 is included in the QUADROS sample of ULIRGs (Rose et al.,
2018), and the QSOFEED sample of type-2 quasars (Ramos Almeida et al., 2022).
— it might also be expected to accelerate outflows via radiation-pressure-
driven winds.

Despite the fact that outflows in F13451+1232 have been confirmed in the
coronal (Villar Martín et al., 2023), warm ionised (Holt et al., 2003; Holt et
al., 2011; Rose et al., 2018; Villar Martín et al., 2023), and neutral atomic
(HI + NaID: Morganti et al. 2005; Rupke et al. 2005; Morganti et al. 2013) gas
phases, there has so far not been a robust detection of molecular outflows in
this object (Fotopoulou et al. 2019; Lamperti et al. 2022; see discussion in
Villar Martín et al. 2023). In nearby AGN where molecular outflows have been
detected, it is often found that they contain more mass and kinetic power than
the warmer phases (e.g. Fiore et al. 2017; Ramos Almeida et al. 2019; Holden
et al. 2023; Holden & Tadhunter 2023; Speranza et al. 2023). In this context,
it is important to note that the warm ionised and neutral atomic outflows in
F13451+1232 have relatively modest mass outflow rates ($\sim$6–$12$ M⊙ yr-1:
Morganti et al. 2013; Rose et al. 2018); therefore, if molecular outflows are
present, they could potentially carry enough mass to change the interpretation
AGN feedback in this important object.

A cold molecular outflow in a similar source, PKS 1549-79, was found to be
compact ($r$ < 120 pc: Oosterloo et al. 2019), as were the previously-detected
neutral HI ($r_{\mathrm{HI}}$< 100 pc: Morganti et al. 2013) and warm ionised
($r_{\mathrm{[OIII]}}\sim$69 pc: Tadhunter et al. 2018) outflows in
F13451+1232. Therefore, here we present high-spatial resolution ($0.113\times
0.091$ arcsecond or $247\times 119$ pc beam-size), high sensitivity ALMA
observations of the inner few kiloparsecs of the primary nucleus of
F13451+1232 to search for (and characterise) cold molecular outflows in
CO(1–0) emission.

In Section 2, we describe the ALMA CO(1–0) observations and our data reduction
process; in Section 3 we detail our analysis of the cold molecular gas in the
central kiloparsecs of F13451+1232; in Section 4 we discuss our findings and
their wider implications, and in Section 5 we give our conclusions.

Throughout this work, we assume a cosmology of $H_{0}=70$ km s-1 Mpc-1,
$\Omega_{\mathrm{m}}=0.3$, and $\Omega_{\mathrm{\lambda}}=0.7$. This
corresponds222Calculated using Ned Wright’s Javascript Cosmology Calculator
(Wright, 2006). to an arcsecond-to-kpc spatial conversion factor of 2.189
kpc/arcsec and a luminosity distance of $D_{\mathrm{L}}=570$ Mpc at the
redshift of F13451+1232 ($z=0.121680$: Lamperti et al. 2022).

##  2 Observations and data reduction

F13451+1232 was observed in Band 3 of the Atacama Large
Millimeter/submillimeter Array (ALMA) on two nights using two different
configurations of the 12 m array (Project Code 2019.1.01757.S): the array
configuration used on the 10th September 2021 (baselines $=180$–16200 m)
resulted in a beam size of $0.083\times 0.069$ arcseconds, while the
configuration used on the 28th September 2021 (baselines $=70$–14400 m)
resulted in a beam size 0.166$\times$0.123 arcseconds. The total time on-
source was 5588 seconds ($\sim$1.5 hours). Four spectral windows were used for
the observations: one centred on 102.790 GHz with a bandwidth of 1.875 GHz
(covering the CO(1–0) line), and the remaining three centred on 104.581 GHz,
90.977 GHz, and 92.477 GHz with bandwidths of 2.000 GHz (covering the
continuum). J1337-1257 and J1353+1435 were observed to provide bandpass/flux
and phase calibration, respectively; phase self-calibration was performed
using the Miriad data reduction package (Sault et al., 1995).

The CO(1–0) transition was targeted for two reasons. First, the ALMA observing
band that covers the CO(1–0) line is wider than those that cover higher-order
transitions such as CO(2–1) or CO(3–2), and therefore facilitates more
accurate continuum subtraction. Second, the CO(1–0) transition is less
susceptible to gas excitation conditions than the higher-order CO transitions,
and hence may be a more reliable probe of molecular gas masses in ULIRGs.

In order to check for instrumental errors such as bandpass and spectral
calibration issues, we produced three CO(1-0) line datacubes: two produced
using the visibilities of the two different array configurations separately,
and one produced by combining the visibilities from both observations (with a
final beam size of $0.113\times 0.091$ arcseconds, PA$=-51.8^{\circ}$). To
optimise the column density sensitivity, and to ensure that the dirty beam was
not significantly non-Gaussian, each CO(1–0) line cube was made using Briggs
weighting with robust $=1$. When inspecting these cubes individually, we did
not find any indications of bandpass or spectral calibration issues.

The data products produced by the automated ALMA pipeline are potentially
affected by issues related to continuum subtraction, which are a consequence
of the bright and steep continuum emission produced by the primary nucleus of
F13451+1232. Therefore, for all three datacubes, we started with the
respective visibilities as calibrated by the ALMA pipeline and applied a
standard self-calibration using continuum images made from the line-free
channels. We then created the line datacubes, and removed the continuum in
each by subtracting a linear function that was fitted to the continuum on
either side of (but excluding) the CO(1–0) emission. A continuum image was
also generated from the combined visibilities of both observations using the
continuum spectral windows. We used least-squares fitting with the AstroPy
(Astropy Collaboration et al., 2013, 2018) Python module to fit a two-
dimensional Gaussian to the point source in this continuum image, and took the
centroid position (13:47:33.36 +12 17 24.23) to be the location of the
nucleus. The continuum image was made using uniform weighting and has a
resolution of $0.049\times 0.033$ arcseconds (PA$=-46.8^{\circ}$).

To ensure the highest sensitivity and signal-to-noise, we henceforth only
consider the combined CO(1–0) datacube in our main analysis, which was
produced using the combined visibilities of both observations; we discuss
results from the cubes produced using the two observations separately in
Appendix B. The combined cube has a root-mean-square (RMS) noise of $0.148$
mJy beam-1 for a velocity resolution of 28.8 km s-1.

##  3 Analysis and results

###  3.1 Moment maps

Figure 1: Integrated flux (moment 0), flux-weighted velocity (moment 1), and
velocity dispersion (moment 2) maps for the central $\sim$6 kiloparsecs of the
primary nucleus of F13451+1232. The white cross marks the continuum centre,
which we take to be the position of the AGN nucleus. The beam size of the
observations is shown as a pink ellipse in the bottom left corner, and a 1 kpc
scale bar is shown in black.

To investigate the flux and kinematic structures of the CO(1–0) emission
around the primary nucleus of F13451+1232, we first created moment maps of the
combined datacube using the SoFiA-2 source-finding pipeline (Serra et al.,
2015; Westmeier et al., 2021). First, SoFiA-2 applied different combinations
of spatial and velocity smoothing to the datacube: Gaussian filters with FWHMs
of 0, 3, 6, 12, and 15 pixels were used for the spatial smoothing, and boxcar
filters of widths 0, 3, 7, and 15 channels were used for the velocity
smoothing. For each combination of spatial and velocity smoothing, emission
above $4\sigma_{\mathrm{rms}}$ was added to a total mask. This mask was then
applied to the original, non-smoothed datacube, from which the maps in Figure
1 and 2 were produced.

In the $6\times 6$ kpc moment maps (Figure 1), cold molecular gas is detected
to radial distances of $\sim$2.5 arcseconds ($\sim$5.5 kpc) south of the
primary nucleus, and $\sim$1 arcsecond ($\sim$2 kpc) to the west. This large-
scale ($r$ > 1 kpc) emission, which has velocities in redshift and blueshift
of up to $v$$\sim$|500| km s-1, is indicative of gas that has been disturbed
by the merger.

Figure 2: Integrated flux (moment 0), flux-weighted velocity (moment 1), and
velocity dispersion (moment 2) maps for the inner kiloparsec region of the
primary nucleus of F13451+1232; the marker and symbol scheme is the same as in
Figure 1.

Within the central kiloparsec of the primary nucleus (Figure 2), there is
clear evidence for a disk with radius $r_{\mathrm{disk}}$$\sim$0.5 kpc centred
on the nucleus, blueshifted to the east and redshifted to the west, with a
projected maximum flux-weighted velocity of 200–250 km s-1. In addition, there
is an extended region (up to $r$$\sim$0.3 arcseconds or 660 pc from the
nucleus along PA$=120^{\circ}$) of emission with higher velocity dispersions
($\sigma$ > 150 km s-1) than that of the disk. Furthermore, a particularly
prominent redshifted component with a velocity dispersion of $\sim$250 km s-1
is seen within $\sim$0.1 arcseconds ($\sim$220 pc) of the continuum centre.

###  3.2 Properties of the kpc-scale disk

To separate non-gravitational kinematics from those expected from a rotating
gas disk, we modelled the disk using the 3DFIT task from the BBarolo tool (Di
Teodoro & Fraternali, 2015); the fitting process and results are detailed in
Appendix A. From the BBarolo model, we find that the disk has an average
inclination of $i\approx 54^{\circ}$ (where $i=0^{\circ}$ corresponds to the
disk being face-on) and that its major axis lies along PA$\approx
248^{\circ}$. The deprojected rotational velocity increases from
$v_{\mathrm{rot}}=246$ km s-1 to $v_{\mathrm{rot}}=307$ km s-1 over the radius
range 28 pc < $r_{\mathrm{disk}}$ < 560 pc. We highlight that the purpose of
this modelling was solely to give a reasonable estimate of the maximum
rotational velocity of the molecular gas, which was needed to robustly
identify non-rotational motions in our analysis. Considering this, it is
notable that the high velocity-dispersion emission ($\sigma$ > 150 km s-1)
seen near the disk in the moment maps (Figure 2) cannot be accounted for by
the BBarolo model (see Figure 6), and thus cannot be explained as being part
of the rotating gas disk.

###  3.3 Kinematics of the non-rotating gas

To investigate the non-rotational kinematics further, we binned the velocity
channels of the original cube by a factor of three and applied Hanning
smoothing using the CASA software suite (CASA Team et al., 2022) to improve
the signal-to-noise ratio of the emission. This resulted in a cube of velocity
resolution of 86.4 km s-1, with a RMS noise of $0.083$ mJy beam-1. Using the
pvextractor Python module333https://pvextractor.readthedocs.io/en/latest/ we
extracted $1.00\times 0.05$ arcsecond slices, centred on the nucleus, from
both this datacube and the BBarolo disk-model datacube along the major axis of
the disk (PA=248∘) and in the direction of the high-velocity dispersion gas
near the nucleus (PA=120∘). From these slices, we produced position-velocity
(PV) diagrams, which we present in Figure 3. In the PV diagrams, there is
clear evidence for non-rotational kinematics, with detected velocities being
well above the projected maximum rotational disk velocity given by the BBarolo
model. There is intermediate velocity emission (300 < $|v|$ < 400 km s-1) seen
in both blueshift and redshift along PA$=120^{\circ}$, in the ranges 0.1 < $r$
< 0.3 arcseconds and 0.0 < $r$ < 0.2 arcseconds, respectively, to the
southeast (SE) of the nucleus. Moreover, the compact ($r$ < 0.1 arcseconds)
redshifted feature visible in the moment maps (Figure 2) can be seen in the PV
diagrams between 400 < $v$ < 680 km s-1, located close to the nucleus.

Figure 3: Position-velocity (PV) diagrams for $1.00\times 0.05$ arcsecond
slices with PA=248∘ (left; the major axis of the disk) and PA=120∘ (right; the
direction of the high velocity dispersion emission seen in Figure 2), centred
on the continuum centre. CO(1–0) emission is shown in yellow/orange, with
black contours at the $-3\sigma_{\mathrm{rms}}$, $-1.5\sigma_{\mathrm{rms}}$,
1.5$\sigma_{\mathrm{rms}}$, 3$\sigma_{\mathrm{rms}}$,
6$\sigma_{\mathrm{rms}}$, and 12$\sigma_{\mathrm{rms}}$ levels (dashed:
negative; solid: positive); solid blue contours show the flux from the BBarolo
disk model at the same $\sigma_{\mathrm{rms}}$ levels. Extended, intermediate-
velocity emission is seen in both redshift and blueshift in the PA=120∘
diagram, and a compact, higher-velocity feature (400 < $v$ < 680 km s-1)
centred approximately on the nucleus (0.0 arcsec) can clearly be seen in both
diagrams.

###  3.4 Comparison to the small-scale radio structure

Due to the velocities at which the compact, high-velocity nuclear feature is
detected, it is unambiguously distinct from the rotating gas disk: it may be
interpreted as either an inflow moving towards the nucleus on the observer’s
side of the disk, or an outflow originating from the nucleus on the far side
of the disk. Spatial comparison of the high-velocity flow to the small-scale
radio structure may provide an indication of the nature of this flowing gas
and, if it is an outflow, the driving mechanism. In Figure 4a, we show the
flux map integrated over the range 400 < $v$ < 680 km s-1 — in which we detect
peak emission above the 6$\sigma_{\mathrm{rms}}$ level — along with VLBI 1266
MHz continuum imaging (from global VLBI experiment GM62B; data presented and
detailed in Morganti et al. 2013) of the central $0.20\times 0.20$ arcseconds
($440\times 440$ parsecs). The beam size of the VLBI observations was
$8.01\times 4.89$ milliarcseconds (PA$=-20.2^{\circ}$). Due to the small
uncertainties in the ALMA pointing, the VLBA image was offset from our
continuum centre by $\sim$0.1 arcseconds. We accounted for this pointing error
by aligning the compact core seen in the VLBA image with the continuum centre,
as it is expected that the same source is being imaged at different spatial
resolutions.

By fitting a 2D Gaussian profile to the integrated high-velocity CO(1–0)
emission in Figure 4a using the emcee Python module (Foreman-Mackey et al.,
2013), we find it to be offset to the southeast (SE) of the continuum centroid
by (3.2$\pm$0.2)$\times 10^{-2}$ arcseconds (70$\pm$4 pc) along PA$\sim$155∘.
We highlight that this is within the spatial extent of, and has a similar PA
to, the small-scale radio structure ($r\sim 0.07$ arcsec;
$\mathrm{PA}=151^{\circ}$)

Figure 4: a) Integrated flux map (orange–brown) of the velocity channels
between 400 < $v$ < 680 km s-1. The dashed black ellipse shows the aperture
extracted from the datacube to derive the mass of the outflow; the
corresponding line profile is shown in Figure 4b. The small-scale radio
structure, as seen in VLBI 1266 MHz continuum imaging (Morganti et al., 2013),
is shown as dark blue solid contours at the $0.3\times$(-1, 1, 2, 4, 8, 16,
32, 64, 128, 256, 512, 1024) mJy/beam levels, whereas the white cross marks
the position of the continuum centre seen in our observations (Section 2). b)
CO(1–0) line profile extracted from the aperture that covers the outflow; the
line profile shown was extracted from the 84 km s-1 velocity resolution
combined datacube. The part of the line profile that we take to represent
in-/out-flowing gas (400 < $v$ < 680 km s-1) is shaded in yellow.

###  3.5 Channel maps

To determine the level at which the high-velocity, nuclear CO(1–0) emission is
detected in different velocity channels, we created channel maps of the
central kiloparsec of the primary nucleus from the 86.4 km s-1 velocity
resolution cube. The high-velocity (400 < $v$ < 600 km s-1) channel maps —
presented in Figure 5 — show emission above the 3$\sigma_{\mathrm{rms}}$ level
(and up to 4$\sigma_{\mathrm{rms}}$) in all four channels on scales similar to
that of the beam size (0.11 arcseconds or 240 pc), indicating that the feature
is spatially-unresolved in the combined datacube.

Figure 5: Channel maps of the central $\sim$1 kpc of the primary nucleus of
F13451+1232 (centred on the nucleus: Section 2) between 448 < $v$ < 576 km
s-1, taken from the velocity-binned and smoothed cube with channel widths of
43.7 km s-1 and a velocity resolution of 86.4 km s-1. Black contours are shown
at the $-3\sigma_{\mathrm{rms}}$, $-1.5\sigma_{\mathrm{rms}}$,
1.5$\sigma_{\mathrm{rms}}$, 3$\sigma_{\mathrm{rms}}$ levels (dashed: negative;
solid: positive), where $\sigma_{\mathrm{rms}}$ is the RMS of the binned-and-
smoothed cube. The continuum centre is marked as a white cross in each
channel; the beam size is shown as a pink ellipse in the 448 km s-1 channel
map (top left panel),

###  3.6 Potential effects of continuum subtraction

Due to the bright, steep-spectrum continuum around the CO(1–0) line —
originating from primary nucleus of F13451+1232 (as seen in the cores of other
CSS/GPS sources, e.g. Oosterloo et al. 2019) — it is possible that
inaccuracies in the continuum subtraction process have contributed flux to the
compact, redshifted, high-velocity feature that is detected at the position of
the nucleus. However, we argue that this feature represents real CO(1–0)
emission, for several reasons. Firstly, the flux is detected above the
$6\sigma$ level in the flux map integrated over the velocity range 400 < $v$ <
680 km s-1 (Figure 4a), and in multiple velocity channels at the
3–4$\sigma_{\mathrm{rms}}$ levels (Figure 5). Secondly, it has a significantly
higher flux than the negative or positive residual features detected in the
integrated spectrum (Figure 4b). Finally, and most importantly, the feature is
detected at the same velocity — and at higher flux than any residual features
— in both of the two observations that were used to produce the combined
datacube, despite their differing spatial resolutions (see Appendix B). While
we cannot entirely rule out some contribution from instrumental effects, we
conclude that this feature is real.

It is notable that previous CO(2–1) observations of the primary nucleus of
F13451+1232 — that were of similar sensitivity to our CO(1–0) observations —
did not detect evidence for nuclear outflows (Lamperti et al., 2022). However,
in those observations, the other spectral windows were also used when fitting
a linear slope to the continuum because of the narrower observing band (a
factor of two narrower in velocity) of CO(2–1) observations relative to
CO(1–0) — this may have led to over-subtraction of the high-velocity outflow
component. Therefore, it is possible that outflow emission in CO(2–1) was not
detected due to uncertainties in the continuum subtraction process.

###  3.7 Outflow energetics

In order to provide an estimate of the mass inflow/outflow rates of the high-
velocity, spatially-unresolved redshifted feature located at the position of
the nucleus, we extracted an elliptical aperture matching the beam size
($0.113\times 0.091$ arcsec, PA$=-51.8^{\circ}$; shown in Figure 4a) from the
datacube. We then summed all pixels within the aperture to create a line
profile, which we present in Figure 4b, and selected emission between 400 <
$v$ < 680 km s-1 to isolate the non-rotating gas.

For the extended, intermediate-velocity (300 < $|v|$ < 400 km s-1) emission
seen in blueshift and redshift, we extracted two sub-apertures from the
$1.00\times 0.05$ arcsecond slice along PA$=120^{\circ}$ (used to create the
right PV diagram in Figure 3): one covering radii 0.0 < $r$ < 0.2 arcseconds
(for the extended redshifted emission), and the other covering the range 0.1 <
$r$ < 0.3 arcseconds (for the extended blueshifted emission). We then selected
emission between 300 < $v$ < 400 km s-1 and $-400$ < $v$ <$-300$ km s-1
respectively in these slices to isolate the non-rotating gas.

For every channel within the range of the selected velocities in the extracted
apertures, we calculated the CO(1–0) fluxes ($S_{\mathrm{CO}}\Delta V$) and
converted them to luminosities ($L^{\prime}_{\mathrm{CO}}$) using the relation
presented by Solomon & Vanden Bout (2005):

$L^{\prime}_{\mathrm{CO}}=3.25\times 10^{7}S_{\mathrm{CO}}\Delta
V{\nu}^{-2}_{\mathrm{obs}}D^{2}_{\mathrm{L}}(1+z)^{-3},$ (1)

where ${\nu}_{\mathrm{obs}}$ is the observed frequency of the CO(1–0) line
(102.789 GHz) and $D_{\mathrm{L}}$ is the luminosity distance to F13451+1232.
We then summed the CO(1–0) luminosities across all selected velocity channels,
resulting in total luminosities for each of the two intermediate-velocity
features and the high-velocity feature. Assuming a CO-to-H2 conversion factor
of $\alpha_{\mathrm{CO}}=0.8$ M⊙ (K km s-1 pc2)-1 (typical for ULIRGs: Downes
& Solomon 1998), we converted these total luminosities into molecular gas
masses. We also applied this method to the rotating kpc-scale CO disk, for
which we estimate a total mass of $1.1\times 10^{10}$ M⊙.

The corresponding mass flow rates of the non-rotating molecular gas were
estimated using

$\dot{M}_{\mathrm{flow}}=\frac{M_{\mathrm{flow}}v_{\mathrm{w}}}{\Delta
R}\cdot\mathrm{tan\theta},$ (2)

where $v_{\mathrm{w}}$ is the flux-weighted velocity, $\Delta R$ is the flow
radius, and $\theta$ is the angle of the outflow relative to our line of
sight. The $\mathrm{tan\theta}$ term is required to correct for projection
effects in both radius and velocity. To estimate this projection factor, we
assumed a solid-angle-weighted mean inclination for random orientations of the
line of sight to the outflow direction, giving $\theta=57.3^{\circ}$. Finally,
we estimated kinetic powers with

$\dot{E}_{\mathrm{kin}}=\frac{1}{2}\dot{M}_{\mathrm{flow}}\Bigl{(}\frac{v_{\mathrm{w}}}{\mathrm{cos\theta}}\Bigl{)}^{2}.$
(3)

We present the resulting CO(1–0) luminosities, masses, flux-weighted
velocities, flow radii, mass flow rates, kinetic powers, and coupling
efficiencies
($\epsilon_{\mathrm{f}}=\dot{E}_{\mathrm{kin}}$/$L_{\mathrm{bol}}$) in Table
1. For the high-velocity nuclear gas, we derive a mass flow rate of
$\dot{M}_{\mathrm{out}}\sim$230 M⊙ yr-1, which corresponds to a kinetic power
that is $\sim$1.4 per cent of the bolometric luminosity of the AGN
($L_{\mathrm{bol}}=4.8\times 10^{45}$ erg s-1: Rose et al. 2018). In contrast,
we find much lower mass flow rates ($\dot{M}_{\mathrm{out}}=22$–$27$ M⊙ yr-1)
and kinetic powers ((6.0–7.1)$\times 10^{-2}$ per cent of $L_{\mathrm{bol}}$)
for the extended, intermediate-velocity emission.

Component | $L^{\prime}_{\mathrm{CO}}$ | $M_{\mathrm{flow}}$ | $v_{\mathrm{w}}$ |  ${\Delta}R$ 1 | $\dot{M}_{\mathrm{flow,proj}}$ | $\dot{M}_{\mathrm{flow}}$ | $\dot{E}_{\mathrm{kin}}$ |  $\epsilon_{\mathrm{f}}$ 2
---|---|---|---|---|---|---|---|---
[K km s-1 pc2] | [M⊙] | [km s-1] | [arcsec] | [M⊙ yr-1] | [M⊙ yr-1] | [erg s-1] | [per cent]
Extended blueshifted | $2.7\times 10^{7}$ | $2.2\times 10^{7}$ | $-339$ | 0.2 | 17 | 27 | $3.4\times 10^{42}$ | $7.1\times 10^{-2}$
Extended redshifted | $2.2\times 10^{7}$ | $1.8\times 10^{7}$ | 350 | 0.2 | 14 | 22 | $2.9\times 10^{42}$ | $6.0\times 10^{-2}$
Nuclear redshifted | $4.2\times 10^{7}$ | $3.4\times 10^{7}$ | 514 | 0.0545 | 148 | 230 | $6.6\times 10^{43}$ | 1.4

1For the extended emission, we take ${\Delta}R$ to be the radial distance over
which the extended emission is seen; for the spatially unresolved, high-
velocity redshift feature, we take ${\Delta}R$ to be half of the beam major
axis.
2Calculated assuming $L_{\mathrm{bol}}=4.8\times 10^{45}$ erg s-1 (Rose et
al., 2018).

Table 1: CO(1–0) luminosities, masses, projected flux-weighted outflow
velocities, aperture radii, projected mass outflow rates
($\dot{M}_{\mathrm{flow,proj}}$), deprojected mass outflow rates
($\dot{M}_{\mathrm{flow}}$), and deprojected energetics for the extended
emission in both blueshift and redshift, and the nuclear redshifted feature
seen in the PV diagrams (Figure 3). When deprojecting, we assumed
$\theta=57.3^{\circ}$ (see Section 3.7).

##  4 Discussion

In addition to a kpc-scale disk, we have detected non-rotating cold molecular
gas in CO(1–0) emission near the primary nucleus of F13451+1232. This non-
rotating gas consists of spatially-extended (0.0 < $r$ < 0.3 arcseconds; 0 <
$r$ < 660 pc), intermediate-velocity (300 < $|v|$ < 400 km s-1) components
seen in both blueshift and redshift, and a compact ($r$ < 0.06 arcseconds; $r$
< 120 pc), spatially-unresolved, high-velocity (400 < $v$ < 680 km s-1)
nuclear component seen in redshift.

###  4.1 Interpreting the non-rotating gas as inflowing

Given that the AGN in F13451+1232 may have recently restarted and thus be
young (Stanghellini et al., 2005; Morganti et al., 2013), if the non-rotating
gas seen in intermediate and high-velocity CO(1–0) emission is inflowing, it
is possible that we are directly observing the gas flows responsible for its
triggering. The total gas flow rates (27–230 M⊙ yr-1: Table 1) are orders of
magnitude higher than those needed to sustain the AGN (0.8 M⊙yr-1, assuming
$\eta=0.1$ and $L_{\mathrm{bol}}=4.8\times 10^{45}$ erg s-1: Rose et al.
2018). Therefore, this is a plausible feeding mechanism even if only a small
fraction of the inflowing gas is accreted onto the central supermassive black
hole (SMBH). However, we argue below that the compact, high-velocity component
(and potentially also the intermediate velocity emission) is in reality
outflowing, rather than inflowing.

###  4.2 The spatially-extended, intermediate-velocity gas

The nature of the intermediate velocity emission (300 < $|v|$ < 400 km s-1)
that is seen in both redshift and blueshift in the PV diagram along
PA$=120^{\circ}$ (Figure 3) — and as a region of higher velocity dispersion to
the SE of the nucleus in the moment 2 maps (Figures 1 and 2) — is not clear.
On one hand, it may represent spatially-extended (220 < $r$ < 440 pc) AGN-
driven outflows. However, since the emission extends well beyond the compact
radio source, then these outflows would not be driven solely by the jet — an
additional outflow driving mechanism would be required. Based on the similar
PAs of the extended emission and the small-scale radio structure, one
possibility is that this gas is accelerated by radiatively-driven winds that
are collimated by a circumnuclear torus with an axis similar to that of the
compact radio source. Alternatively, due to the velocities not being much
above those expected for rotation (as given by the BBarolo model: Section
3.2), this extended emission may represent material from the merger settling
onto the disk. This latter interpretation is supported by the larger-scale
($r$ >1 kpc) filament-like structures seen to the south and west of the
primary nucleus (Figure 1), which likely represent gas disturbed by the
merger: these large-scale structures have similar velocity shifts and
dispersions to the intermediate-velocity gas seen near the nucleus.

If the intermediate-velocity emission truly did represent AGN-driven outflows,
then its mass outflow rates (22–27 M⊙ yr-1) are higher than those of the
previously detected warm ionised (11.3 M⊙ yr-1: Rose et al. 2018)444The warm
ionised mass outflow rates and kinetic powers were calculated following the
second method of Rose et al. (2018), which involves using the
$v_{\mathrm{05}}$ parameter as the outflow velocity to account for projection
effects; these estimates are likely to be upper limits. and neutral HI
($\sim$6 $M_{\odot}$ yr-1)555The neutral HI mass outflow rate has been
recalculated using the methodology in Appendix C1 of Holden et al. (2023) with
the values presented by Morganti et al. (2013); the kinetic power and coupling
efficiency have been recalculated by using these values with Equation 3 and
assuming $L_{\mathrm{bol}}=4.8\times 10^{45}$ erg s-1 (Rose et al., 2018).
phases. Furthermore, they are similar to those recently reported for the cold
molecular outflows in a sample of nearby QSO2s (8–16 M⊙ yr-1: Ramos Almeida et
al. 2022), which were identified using a similar method to that used here for
the intermediate velocity gas. When applying a factor of three to account for
assumed spherical outflow geometry, as was done by Fiore et al. (2017), we
find (in agreement with Ramos Almeida et al. 2022) that these mass outflow
rates (24–81 M⊙ yr-1) fall significantly below the correlation between mass
outflow rate and AGN bolometric luminosity given by Fiore et al. (2017).
However, it is important to note that Lamperti et al. (2022) reported a wide
range of mass outflow rates for ULIRGs with similar bolometric luminosities (6
< $\dot{M}_{\mathrm{flow}}$ < 300 M⊙ yr-1), indicating that the high mass
outflow rates found by Fiore et al. (2017) may represent the most extreme
cases, and therefore represent the upper envelope to any relationship between
mass outflow rate and AGN bolometric luminosity (see also Speranza et al.
2023).

###  4.3 The compact, high-velocity nuclear outflow

While it is possible that the intermediate-velocity emission does not
represent outflowing gas, we highlight that the observed velocities of the
redshifted nuclear feature (400 < $v$ < 680 km s-1) are much higher than what
would be expected due to settled gravitational motions ($v$ < 300 km s-1: the
projected velocities seen in the disk). Moreover, the emission of this
component is compact ($r$ < 120 pc), unresolved in the combined datacube, and
on a similar radial scale to both the compact radio structure ($r$ < 0.06
arcseconds; $r$ < 130 pc: Figure 4a) and the other outflow phases ($r$ < 100
pc: Morganti et al. 2013; Tadhunter et al. 2018). Therefore, it is much more
likely that we are seeing a nuclear outflow, rather than an inflow or settling
gas from the merger. Furthermore, the fact that the high-velocity emission is
spatially offset from the nucleus (by (3.2$\pm$0.2)$\times 10^{-2}$
arcseconds) along the direction of the small-scale radio structure
(PA=$151^{\circ}$) suggests a relation between the two. This is consistent
with acceleration by shocks, which may be induced by either a jet or a
broader, radiatively-driven wind that is collimated by the circumnuclear
torus. Given that the compact CO(1–0) outflow is seen to be offset along the
direction of the radio structure (which would not necessarily be expected in
the radiatively-driven scenario), and that the neutral atomic outflow
component (detected in HI absorption) is seen to be coincident with the radio
hotspot (Morganti et al., 2013), we favour the jet-driven interpretation. In
this context, it is interesting to note the southern bend in the small-scale
radio structure, which may be due to the jet colliding with dense molecular
gas, being deflected, and thereby launching the nuclear outflow.

The cold molecular mass outflow rate for this nuclear component
($\dot{M}_{\mathrm{flow}}$$\sim$230 M⊙yr-1) is more than an order of magnitude
higher than those previously derived for the other outflow phases in
F13451+1232, lies within the range of values for molecular outflows previously
detected in ULIRGs (a few to thousands of solar masses per year: Cicone et al.
2014; Pereira-Santaella et al. 2018; Fluetsch et al. 2019; Herrera-Camus et
al. 2020; Lamperti et al. 2022), and is close to (within 1$\sigma$) of the
correlation between cold molecular mass outflow rate and AGN bolometric
luminosity found by Fiore et al. (2017) when corrected for assumed spherical
outflow geometry.

Overall, the kinetic power of the cold molecular phase in F13451+1232
($6.6\times 10^{43}$ ergs-1) accounts for $\sim$1.4 per cent of Lbol —
approximately a factor of three higher than the combined kinetic power of the
neutral HI ($\sim$0.04 per cent of $L_{\mathrm{bol}}$: Morganti et al. 2013)
and warm ionised (0.49 per cent of $L_{\mathrm{bol}}$: Rose et al. 2018)
outflow phases. Therefore, for all phases, the total outflow kinetic power is
$\sim$2 per cent of the AGN bolometric luminosity. While this is consistent
with the requirements of simulations of outflows launched by AGN in galaxy
mergers (0.5 < $\epsilon_{\mathrm{f}}$ < 10 per cent: e.g. Di Matteo et al.
2005; Hopkins & Elvis 2010) — and thus may be interpreted as the outflows
having an impact on the evolution of the galaxy — we highlight that such
comparisons must be done with care, and that it is unclear if interpretations
can be made robustly (see Harrison et al. 2018).

Finally, we emphasize that even if we assume a CO-to-H2 conversion factor that
is at the lower end of the range of values used in the literature (i.e.
$\alpha_{\mathrm{CO}}=0.3$ M⊙ (K km s-1 pc2)-1 for optically thin gas:
Oosterloo et al. 2017, 2019), the derived cold molecular, nuclear mass outflow
rate is still much higher than those of the other phases ($\sim$86 M⊙yr-1).
Conversely, if a conversion factor that is characteristic of molecular clouds
in the Milky Way is assumed ($\alpha_{\mathrm{CO}}=4.3$ M⊙ (K km s-1 pc2)-1:
Bolatto et al. 2013), the mass and kinetic power of this outflow increase
significantly ($\dot{M}_{\mathrm{out}}\sim 1240$ M⊙ yr-1;
$\epsilon_{\mathrm{f}}\sim 7$ per cent).

##  5 Conclusions

Through the analysis of high-spatial-resolution CO(1–0) observations of the
primary nucleus of the ULIRG F13451+1232, we have detected and characterised a
kpc-scale circumnuclear disk and spatially-extended (0.0 < $r$ < 0.3
arcseconds; 0 < $r$ < 660 pc), intermediate-velocity emission (300 < $|v|$ <
400 km s-1). This intermediate velocity molecular gas may represent AGN
radiation-driven outflows, or merger material settling onto the disk.
Furthermore, we have also presented evidence for a compact ($r$ < 120 pc),
cold molecular outflow (400 < $v$ < 680 km s-1) at the position of the nucleus
with a mass outflow rate of $\dot{M}_{\mathrm{out}}\sim$230 M⊙ yr-1 and a
kinetic power that is $\sim$1.4 per cent of the AGN bolometric luminosity.

Overall, our detection of cold molecular outflow(s) in F13451+1232 changes the
scenario of the AGN-driven outflows in this object: previous studies found
that the neutral atomic HI and warm ionised phases had relatively modest mass
outflow rates, whereas here we find that the cold molecular outflow is
carrying over an order of magnitude more mass (and several times the kinetic
power) than the warm ionised and neutral phases. Similar results have been
found for other types of active galaxy, where the colder phases have been seen
to carry the majority of the outflowing mass (Ramos Almeida et al., 2019;
Holden et al., 2023; Holden & Tadhunter, 2023; Speranza et al., 2023). This
further demonstrates the need for multi-wavelength observations to fully
quantify the impact of AGN-driven outflows. Moreover, a similarly compact ($r$
< 120 pc) but more massive ($\dot{M}_{\mathrm{out}}\sim$650 M⊙yr-1) outflow
was found in the ULIRG PKS 1549-79 (Oosterloo et al., 2019), which is also
both a merger and a luminous radio source. Our results thus provide important
evidence that such powerful radio sources accelerate compact molecular
outflows, and that these outflows may carry more mass and power than other gas
phases.

Crucially, in demonstrating that cold molecular AGN-driven outflows can be
compact ($r$ < $120$ pc), our study emphasises the importance of high-spatial
resolution observations when robustly quantifying the impact of AGN-driven
outflows on their host galaxies.

## Acknowledgements

We thank the anonymous referee for their feedback, which helped to improve the
quality of this manuscript. LRH and CNT acknowledge support from STFC. This
paper makes use of the following ALMA data: ADS/JAO.ALMA#2019.1.01757.S. ALMA
is a partnership of ESO (representing its member states), NSF (USA) and NINS
(Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI
(Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA
Observatory is operated by ESO, AUI/NRAO and NAOJ. The European VLBI Network
is a joint facility of independent European, African, Asian, and North
American radio astronomy institutes. Scientific results from data presented in
this publication are derived from the following EVN project code(s): GM062. AA
and CRA acknowledge the projects “Feeding and feedback in active galaxies”,
with reference PID2019-106027GB-C42, funded by MICINN-
AEI/10.13039/501100011033, and “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. MPS
acknowledges funding support from the Ramón y Cajal program of the Spanish
Ministerio de Ciencia e Innovación (RYC2021-033094-I). IL acknowledges support
from PID2022-140483NB-C22 funded by AEI 10.13039/501100011033 and BDC 20221289
funded by MCIN by the Recovery, Transformation and Resilience Plan from the
Spanish State, and by NextGenerationEU from the European Union through the
Recovery and Resilience Facility.

## Data Availability

The ALMA CO(1–0) data used in this work is publicly available from the ALMA
Science Archive (https://almascience.eso.org/aq/) with Project Code
2019.1.01757.S. The data used to produce the VLBI 1266 MHz continuum image
presented in Figure 4a is publicly available from the EVN archive for
experiment GM062B (http://archive.jive.nl/scripts/arch.php?exp=GM062B_060605).

## References

  * Alonso-Herrero et al. (2018) Alonso-Herrero A., et al., 2018, ApJ, 859, 144
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, ApJ, 156, 123
  * Audibert et al. (2023) Audibert A., et al., 2023, A&A, 671, L12
  * Bolatto et al. (2013) Bolatto A. D., Wolfire M., Leroy A. K., 2013, ARA&A, 51, 207
  * CASA Team et al. (2022) CASA Team et al., 2022, PASP, 134, 114501
  * Cicone et al. (2014) Cicone C., et al., 2014, A&A, 562, A21
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827
  * Di Matteo et al. (2005) Di Matteo T., Springel V., Hernquist L., 2005, Nat, 433, 604
  * Di Teodoro & Fraternali (2015) Di Teodoro E. M., Fraternali F., 2015, MNRAS, 451, 3021
  * Domínguez-Fernández et al. (2020) Domínguez-Fernández A. J., et al., 2020, A&A, 643, A127
  * Downes & Solomon (1998) Downes D., Solomon P. M., 1998, ApJ, 507, 615
  * Fiore et al. (2017) Fiore F., et al., 2017, A&A, 601, A143
  * Fluetsch et al. (2019) Fluetsch A., et al., 2019, MNRAS, 483, 4586
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Fotopoulou et al. (2019) Fotopoulou C. M., Dasyra K. M., Combes F., Salomé P., Papachristou M., 2019, A&A, 629, A30
  * Gilmore & Shaw (1986) Gilmore G., Shaw M. A., 1986, Nature, 321, 750
  * Harrison et al. (2018) Harrison C. M., Costa T., Tadhunter C. N., Flütsch A., Kakkad D., Perna M., Vietri G., 2018, Nat. Astron, 2, 198
  * Heckman et al. (1986) Heckman T. M., Smith E. P., Baum S. A., van Breugel W. J. M., Miley G. K., Illingworth G. D., Bothun G. D., Balick B., 1986, ApJ, 311, 526
  * Herrera-Camus et al. (2020) Herrera-Camus R., et al., 2020, A&A, 635, A47
  * Holden & Tadhunter (2023) Holden L. R., Tadhunter C. N., 2023, MNRAS,
  * Holden et al. (2023) Holden L. R., Tadhunter C. N., Morganti R., Oosterloo T., 2023, MNRAS, 520, 1848
  * Holt et al. (2003) Holt J., Tadhunter C. N., Morganti R., 2003, MNRAS, 342, 227
  * Holt et al. (2006) Holt J., Tadhunter C., Morganti R., Bellamy M., González Delgado R. M., Tzioumis A., Inskip K. J., 2006, MNRAS, 370, 1633
  * Holt et al. (2011) Holt J., Tadhunter C. N., Morganti R., Emonts B. H. C., 2011, MNRAS, 410, 1527
  * Hopkins & Elvis (2010) Hopkins P. F., Elvis M., 2010, MNRAS, 401, 7
  * Hopkins et al. (2008) Hopkins P. F., Hernquist L., Cox T. J., Kereš D., 2008, ApJS, 175, 356
  * Kukreti et al. (2023) Kukreti P., Morganti R., Tadhunter C., Santoro F., 2023, A&A, 674, A198
  * Lamperti et al. (2022) Lamperti I., et al., 2022, A&A, 668, A45
  * Lister et al. (2003) Lister M. L., Kellermann K. I., Vermeulen R. C., Cohen M. H., Zensus J. A., Ros E., 2003, ApJ, 584, 135
  * Magorrian et al. (1998) Magorrian J., et al., 1998, ApJ, 115, 2285
  * Meena et al. (2023) Meena B., et al., 2023, ApJ, 943, 98
  * Morganti et al. (2005) Morganti R., Tadhunter C. N., Oosterloo T. A., 2005, A&A, 444, L9
  * Morganti et al. (2013) Morganti R., Fogasy J., Paragi Z., Oosterloo T., Orienti M., 2013, Science, 341, 1082
  * Mukherjee et al. (2016) Mukherjee D., Bicknell G. V., Sutherland R., Wagner A., 2016, MNRAS, 461, 967
  * O’Dea & Saikia (2021) O’Dea C. P., Saikia D. J., 2021, A&ARv, 29, 3
  * Oosterloo et al. (2017) Oosterloo T., Raymond Oonk J. B., Morganti R., Combes F., Dasyra K., Salomé P., Vlahakis N., Tadhunter C., 2017, A&A, 608, A38
  * Oosterloo et al. (2019) Oosterloo T., Morganti R., Tadhunter C., Raymond Oonk J. B., Bignall H. E., Tzioumis T., Reynolds C., 2019, A&A, 632, A66
  * Pereira-Santaella et al. (2018) Pereira-Santaella M., et al., 2018, A&A, 616, A171
  * Ramos Almeida et al. (2019) Ramos Almeida C., Acosta-Pulido J. A., Tadhunter C. N., González-Fernández C., Cicone C., Fernández-Torreiro M., 2019, MNRAS, 487, L18
  * Ramos Almeida et al. (2022) Ramos Almeida C., et al., 2022, A&A, 658, A155
  * Rose et al. (2018) Rose M., Tadhunter C., Ramos Almeida C., Rodríguez Zaurín J., Santoro F., Spence R., 2018, MNRAS, 474, 128
  * Rupke et al. (2005) Rupke D. S., Veilleux S., Sanders D. B., 2005, ApJ, 632, 751
  * Sanders & Mirabel (1996) Sanders D. B., Mirabel I. F., 1996, ARA&A, 34, 749
  * Santoro et al. (2018) Santoro F., Rose M., Morganti R., Tadhunter C., Oosterloo T. A., Holt J., 2018, A&A, 617, A139
  * Santoro et al. (2020) Santoro F., Tadhunter C., Baron D., Morganti R., Holt J., 2020, A&A, 644, A54
  * Sault et al. (1995) Sault R. J., Teuben P. J., Wright M. C. H., 1995, in Shaw R. A., Payne H. E., Hayes J. J. E., eds, Astronomical Society of the Pacific Conference Series Vol. 77, Astronomical Data Analysis Software and Systems IV. p. 433 (arXiv:astro-ph/0612759), doi:10.48550/arXiv.astro-ph/0612759
  * Serra et al. (2015) Serra P., et al., 2015, MNRAS, 448, 1922
  * Silk & Rees (1998) Silk J., Rees M. J., 1998, A&A, 331, L1
  * Solomon & Vanden Bout (2005) Solomon P. M., Vanden Bout P. A., 2005, ARA&A, 43, 677
  * Speranza et al. (2023) Speranza G., et al., 2023, arXiv e-prints, p. arXiv:2311.10132
  * Stanghellini et al. (1997) Stanghellini C., O’Dea C. P., Baum S. A., Dallacasa D., Fanti R., Fanti C., 1997, A&A, 325, 943
  * Stanghellini et al. (2005) Stanghellini C., O’Dea C. P., Dallacasa D., Cassaro P., Baum S. A., Fanti R., Fanti C., 2005, A&A, 443, 891
  * Tadhunter et al. (2018) Tadhunter C., et al., 2018, MNRAS, 478, 1558
  * Villar Martín et al. (2023) Villar Martín M., et al., 2023, A&A, 673, A25
  * Westmeier et al. (2021) Westmeier T., et al., 2021, MNRAS, 506, 3962
  * Wright (2006) Wright E. L., 2006, PASP, 118, 1711

##  Appendix A BBarolo modelling of the gas disk

In order to identify non-rotational motions in the central few kiloparsecs of
the primary nucleus of F13451+1232, we fit the disk-like kinematics seen in
our CO(1–0) data (Section 3; Figure 2) with a simple disk model using the
3DFIT task from the BBarolo tool (Di Teodoro & Fraternali, 2015). The 3DFIT
task fits a series of concentric rings to the data in three spatial dimensions
and three velocity dimensions (the rotational velocity, $v_{\mathrm{rot}}$;
the radial velocity, $v_{\mathrm{rad}}$, and the velocity dispersion
$v_{\mathrm{disp}}$).

Our fitting procedure followed the methodology outlined in Alonso-Herrero et
al. (2018), Domínguez-Fernández et al. (2020), and Ramos Almeida et al.
(2022). Throughout, the centre of the disk model was fixed to be the position
of the continuum centre (as measured from a two-dimensional Gaussian fit to
the continuum image: Section 2). We first performed an initial fit using the
3DFIT task in which the radial velocities of the rings were fixed to be 0 km
s-1 and the scale height of the disk ($z_{\mathrm{0}}$) was fixed to 0 pc. The
remaining parameters were allowed to vary; initial values for the rotational
velocity and velocity dispersion were based on the flux-weighted velocity
shift and velocity dispersion maps (Figure 2), the initial value of the
inclination ($i_{\mathrm{initial}}=38^{\circ}$) was based on the measured
ratio of the projected disk major and minor axes while assuming a circular
disk, and the initial value for the PA was set to that previously estimated
from two-dimensional spectroastrometry of CO(2–1) observations of the disk by
Lamperti et al. (2022). However, limits were imposed on the variation of these
free parameters — between successive rings, the rotational velocity was not
allowed to change by more than $\pm$50 km s-1, nor were the inclination and PA
allowed to change by more than $\pm 20^{\circ}$. The fits used uniform
weighting and local normalisation. It should be noted that, since the physical
size of the modelled rings is much less than the beam size of the
observations, the properties of adjacent rings are correlated. However, this
does not have a significant effect on the maximum rotational velocity of the
disk model (nor does changing the values of the initial parameters), and
therefore does not affect the interpretations made in our analysis.

In order to ensure that the overall radial extent and kinematics of the disk
were accurate — and to alleviate potential problems with fit degeneracy
arising from having a high number of free parameters — we then used the
results from the initial fit as a basis for a second fit. For this second run,
we fixed the velocity dispersion, inclination, and PA to the values determined
from the first run, so that only the rotational velocity was allowed to vary.
We present the final model parameters for each ring in Table 2; in Figure 6,
we present the moment maps and residuals of the final disk model. High flux-
weighted velocities ($\sim$200 km s-1) and velocity dispersions
($\sigma\sim$200 km s-1) due to non-rotational kinematics can be seen in the
residual maps near the centre of the disk.

By taking an aperture around the disk in our CO(1–0) datacube and integrating
between $-310$ < $v$ < $310$ km s-1, we find a disk mass of $1.1\times
10^{10}$ M⊙ (using Equation 1 and $\alpha_{\mathrm{CO}}=0.8$ M⊙ (K km s-1
pc2)-1). Lamperti et al. (2022) also detected this nuclear disk using CO(2-1)
spectroastrometry, and found a similar mass and overall properties to those
reported here.

Finally, we highlight that the purpose of this modelling was solely to
estimate the maximum CO(1–0) velocities that are explainable as rotational
motion, and hence identify outflow kinematics; the detailed parameters of the
disk model presented here are not intended to be a robust representation of
the disk’s properties.

$r_{\mathrm{ring}}$ (arcseconds) |  $r_{\mathrm{ring}}$ (pc) | vrot (km s-1) | vdisp (km s-1) |  $i$ $(^{\circ})$ | PA $(^{\circ})$
---|---|---|---|---|---
0.013 | 28.5 | 247.969 | 5.524 | 62.218 | 242.005
0.040 | 87.6 | 249.627 | 24.639 | 58.627 | 236.155
0.067 | 146.7 | 211.864 | 91.100 | 55.868 | 235.867
0.094 | 205.8 | 249.245 | 72.507 | 53.826 | 239.504
0.121 | 264.9 | 267.474 | 65.771 | 52.384 | 245.429
0.148 | 324.0 | 280.596 | 62.837 | 51.425 | 252.007
0.175 | 383.1 | 293.392 | 61.210 | 50.834 | 257.600
0.202 | 442.2 | 311.494 | 68.331 | 50.494 | 260.572
0.229 | 501.0 | 301.110 | 61.533 | 50.289 | 259.286
0.256 | 560.4 | 314.645 | 58.314 | 50.102 | 252.106
Table 2: Final parameters for each ring of the BBarolo 3DFit disk model. The
distance of each ring to the fixed disk centre ($r_{\mathrm{ring}}$) is given
in both arcseconds and parsecs. Here, vrot and vdisp are the rotational
velocities and velocity dispersions of each ring, respectively. Figure 6:
Integrated flux (moment 0; top), flux-weighted velocity (moment 1; middle),
and velocity dispersion (moment 2; bottom) maps of the CO(1–0) data (left),
BBarolo disk model (middle) and residuals (data - model; right). The ticks
have separations of 0.2 arcseconds. The black dashed line shows the PA of the
modelled disk, whereas the black cross shows the disk centre (which was set to
be the continuum centre position; see Section 2). The green lines in the
middle row panels are isovelocity contours at 0 km s-1.

##  Appendix B The CO(1–0) nuclear outflow observed at different spatial
resolutions

Figure 7: CO(1–0) line profiles of the primary nucleus of F13451+1232
extracted from the two sets of observations of differing array configurations
(hence beam sizes) that comprise the datacube presented in our analysis (see
Section 2). These profiles were extracted from the respective cubes using
apertures matched to the beam sizes. It can be seen that the flux in the high-
velocity redshifted component (shaded in grey) is less in the cube with the
smaller beam size.

The CO(1–0) data that we present in this work was produced from the combined
visibilities of two observations, which were taken with different array
configurations and thus have different beam sizes ($0.083\times 0.069$
arcseconds and 0.166$\times$0.123 arcseconds; see Section 2). Here, we discuss
the detection of the high-velocity, redshifted, compact feature seen at the
position of the nucleus (Section 3; Figures 2, 3, 5, and 4) in the datacubes
produced from the visibilities of these observations separately.

In Figure 7, we present CO(1-0) line spectra extracted from these cubes; the
size of the extraction aperture in each case was matched to the beam size. The
high-velocity redshifted feature is well detected in both datacubes in the
velocity range 400 < $v$ < 680 km s-1, while any residual features of
comparable flux are not present in both cubes at the same velocity. Moreover,
its absolute flux is less in the smaller-beam cube, indicating that the beam
is not recovering all of the flux and thus implying that the feature is
(partially) spatially resolved. Considering the arguments presented in Section
3.6, we conclude that this feature represents real emission from compact,
high-velocity molecular gas near the nucleus.

Since we argue that this molecular gas is a nuclear outflow (Section 4.3), the
fact that it is partially spatially resolved by the smaller beam observations
allows us to constrain the outflow radius to be in the range 88 < $\Delta{R}$
< 182 parsecs (0.042 < $\Delta{R}$ < 0.083 arcseconds; i.e. half of the beam
major axes of the two observations). Considering that the feature does not
appear to be spatially resolved in the channel maps (Section 3.5; Figure 5)
produced from the combined datacube (beam size: $0.113\times 0.091$
arcseconds; see Section 2), this indicates that it is partially spatially-
resolved our observations — we thus take 120 pc (0.0545 arcseconds) as an
upper limit on the outflow radius.