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

# Resolved Neutral Outflow from a Lensed Dusty Star-forming Galaxy at $z=2.09$

Kirsty M. Butler<fnref n="1" />,*, Paul P. van der Werf<fnref n="1" />, Matus Rybak<fnref n="1" />,<fnref n="2" />,*, Tiago Costa<fnref n="3" />, Pierre Cox<fnref n="4" />, Axel Weiß<fnref n="5" />,*, Michał J. Michałowski<fnref n="6" />, Dominik A. Riechers<fnref n="7" />,*, Dimitra Rigopoulou<fnref n="8" />, Lucia Marchetti<fnref n="9" />,<fnref n="10" />, Stephen Eales<fnref n="11" />, and Ivan Valtchanov<fnref n="12" />,†

## Abstract

We report the detection of a massive neutral gas outflow in the $z=2.09$ gravitationally lensed dusty star-forming galaxy HATLAS J085358.9+015537 (G09v1.40), seen in absorption with the $OH^+(1_1-1_0)$ transition using spatially resolved ($0.′′5 \times 0.′′4$) Atacama Large Millimeter/submillimeter Array (ALMA) observations. The blueshifted $OH^+$ line is observed simultaneously with the CO(9–8) emission line and underlying dust continuum. These data are complemented by high-angular-resolution ($0.′′17 \times 0.′′13$) ALMA observations of $CH^+(1-0)$ and underlying dust continuum, and Keck $2.2 \mu$m imaging tracing the stellar emission. The neutral outflow, dust, dense molecular gas, and stars all show spatial offsets from each other. The total atomic gas mass of the observed outflow is $6.7 \times 10^9 M_\odot$, >25% as massive as the gas mass of the galaxy. We find that a conical outflow geometry best describes the $OH^+$ kinematics and morphology and derive deprojected outflow properties as functions of possible inclination ($0.°38–64°$). The neutral gas mass outflow rate is between 83 and $25,400 M_\odot \text{ yr}^{-1}$, exceeding the star formation rate ($788 \pm 300 M_\odot \text{ yr}^{-1}$) if the inclination is $>3.°6$ (mass-loading factor = 0.3–4.7). Kinetic energy and momentum fluxes span ($4.4–290) \times 10^9 L_\odot$ and ($0.1–3.7) \times 10^{37}$ dyne, respectively (energy-loading factor = 0.013–16), indicating that the feedback mechanisms required to drive the outflow depend on the inclination assumed. We derive a gas depletion time between 29 and 1 Myr, but find that the neutral outflow is likely to remain bound to the galaxy unless the inclination is small and may be reaccreted if additional feedback processes do not occur.

*Unified Astronomy Thesaurus concepts: Starburst galaxies (1570); Galaxy evolution (594); Galaxy processes (614); High-redshift galaxies (734); Interstellar absorption (831); Strong gravitational lensing (1643)*

## 1. Introduction

The formation and evolution of galaxies are intrinsically linked to the cosmic web. Dark matter halos accrete gas from the intergalactic medium (IGM) which, by the dissipation of energy, cools and condenses to form a central galaxy (Rees & Ostriker 1977; White & Rees 1978). Gas within galaxies may then collapse to form stars or accrete onto supermassive black holes, injecting energy back into the interstellar medium (ISM) via stellar winds, radiation pressure, supernova (SN) explosions, or through strong feedback associated with an active galactic nucleus (AGN), respectively. Mild forms of these feedback processes heat and disturb the surrounding ISM, prolonging its collapse into new stars while in their extremes, eject gas from the galaxy. Ejected gas is either recycled through the circumgalactic medium (CGM) where it can be reaccreted back onto the galaxy at a later time or be lost to the IGM, removing the fuel for star formation (SF) altogether (Bregman 1980; Bregman et al. 2013; Fluetsch et al. 2019; Nelson et al. 2019; Mitchell et al. 2020; Spilker et al. 2020a).

This model of self-regulated galaxy growth became evident in early cosmological simulation work, which failed to reproduce disky galaxy morphologies without invoking sufficiently strong

SN feedback processes capable of removing low-angular-momentum material from the centers of galaxies (later coined the angular momentum catastrophe; Governato et al. 2010). Today’s state-of-the-art theoretical models (e.g., EAGLE, Schaye et al. 2015; Illustris-TNG Nelson et al. 2019; Pillepich et al. 2019; L-GALAXIES Henriques et al. 2020) similarly rely on a series of feedback processes that tap into the energy released by stars and AGNs in order to regulate stellar mass growth (the overcooling problem, e.g., Somerville & Primack 1999; Cole et al. 2000; Benson et al. 2003; Kereš et al. 2009; Bower et al. 2012), reproduce metallicity gradients and the galaxy mass function, and to pollute the CGM/IGM with metals (Veilleux et al. 2005).

In fact, the identification of metals in the low-density CGM/IGM via absorption in QSO sight lines (Meyer & York 1987; Simcoe et al. 2004) provided the first observational evidence that some fraction of the enriched matter within galaxies must be ejected. Optical and X-ray observations led to the first direct evidence of outflowing material from galaxies in the form of ionized gas (Heckman et al. 1990; Strickland et al. 2004), which has since been complemented by observations of the molecular and atomic phases probed by IR/submillimeter wavelengths (e.g., Walter et al. 2002; Sturm et al. 2011;

<footnote n="1">¹Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands; kirstymaybutler@gmail.com</footnote>

<footnote n="2">²THz Sensing Group, Faculty of Electrical Engineering, Mathematics and Computer Science, TU Delft, The Netherlands</footnote>

<footnote n="3">³Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, D-85748 Garching b. München, Germany</footnote>

<footnote n="4">⁴Sorbonne Université, UPMC Université Paris 6 & CNRS, UMR 7095, Institut d'Astrophysique de Paris, 98b boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="5">⁵Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69 D-53121 Bonn, Germany</footnote>

<footnote n="6">⁶Astronomical Observatory Institute, Faculty of Physics, Adam Mickiewicz University, ul. Słoneczna 36, 60-286, Poznań, Poland</footnote>

<footnote n="7">⁷Department of Astronomy, Cornell University, Space Sciences Building, Ithaca, NY 14853, USA</footnote>

<footnote n="8">⁸Astrophysics, Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="9">⁹Department of Astronomy, University of Cape Town, Private Bag X3, 7701 Rondebosch, Cape Town, South Africa</footnote>

<footnote n="10">¹⁰INAF—Institute for Radio Astronomy, Via Gobetti 101, I-40129, Bologna, Italy</footnote>

<footnote n="11">¹¹School of Physics and Astronomy, Cardiff University, The Parade, Cardiff CF24 3AA, UK</footnote>

<footnote n="12">¹²Telespazio UK for ESA, European Space Astronomy Centre, Operations Department, E-28691 Villanueva de la Cañada, Spain</footnote>

<footnote>Received 2021 February 26; revised 2021 May 31; accepted 2021 June 17; published 2021 September 16</footnote>
Bolatto et al. 2013), corresponding to an enormous range in temperatures (10–10⁸ K, Veilleux et al. 2005) and densities (~10–10⁵ cm⁻³, Shopbell & Bland-Hawthorn 1998; Aalto et al. 2015).

The bright emission and absorption lines associated with the ionized gas phase have led the majority of galaxy outflow observations to focus on this component and have successfully shown that galaxy outflows are ubiquitous in the local universe from dwarf galaxies to luminous infrared galaxies (LIRGs; $L_{IR} > 10^{11} L_{\odot}$). Over this range, outflow velocities have been found to correlate with star formation rate (SFR), stellar mass ($M_*$), and SFR surface density (e.g., Lehnert & Heckman 1996; Rupke et al. 2002; Martin 2005; Westmoquette et al. 2012; Rubin et al. 2014; Chisholm et al. 2016; Heckman & Borthakur 2016), suggesting a close connection between galaxy outflows and the ongoing evolution of their host galaxies. The hot phase, however, only dominates the thermal and kinetic energy of the outflow, while the cooler, denser molecular, and neutral phases are believed to dominate the mass and momentum budget (Walter et al. 2002; Rupke et al. 2005, 2017; Feruglio et al. 2010; Alatalo et al. 2011; Rupke & Veilleux 2013; Herrera-Camus et al. 2020; Fluetsch et al. 2021).

The first detections of molecular outflows from local ultraluminous infrared galaxies (ULIRGs; $L_{IR} > 10^{12} L_{\odot}$) were achieved almost simultaneously using very deep ground-based CO(1–0) spectra (Feruglio et al. 2010) and spectra of OH lines obtained with the Herschel satellite (Fischer et al. 2010) of the ULIRG/AGN Mrk 231. Now, with the addition of new facilities such as the Atacama Large Millimeter Array (ALMA) and the Northern Extended Millimeter Array (NOEMA), outflows have been detected in a large number of local LIRGs and ULIRGs, using both CO and OH rotational lines (Sturm et al. 2011; Spoon et al. 2013; Veilleux et al. 2013; Cicone et al. 2014). In all cases, mass outflow rates of the order or even significantly larger than the star formation rate in the galaxy were derived, suggesting that galactic winds regulate star formation in these systems.

At redshifts $z=1–3$ where the cosmic star formation and black hole accretion peak (Madau & Dickinson 2014), outflows are expected to be ubiquitous. However, observing molecular gas outflows using CO observations in high-$z$ galaxies is extremely challenging. In local (U)LIRGs, the CO emission from the outflowing gas typically represents only a few percent of the total CO emission of the galaxy and requires high signal-to-noise ratio (S/N) observations at high spatial resolution, in order to observationally separate the outflowing gas from the bulk CO emission (Cicone et al. 2014; García-Burillo et al. 2014, 2015; Pereira-Santaella et al. 2018, 2020). Such observations at high-$z$ have provided mostly tentative results and only in galaxies hosting an AGN (see, e.g., Weiß et al. 2012; Carniani et al. 2017; Feruglio et al. 2017; Vayner et al. 2017; Brusa et al. 2018; Fan et al. 2018; Herrera-Camus et al. 2019) and are generally beyond present observational capabilities for other classes of galaxies, even with ALMA.

High-velocity wings in [C II] 158 µm spectra have provided strong evidence of outflowing gas in one main-sequence star-forming galaxy at $z \sim 5$ (Herrera-Camus et al. 2021) and in a handful of particularly extreme high-redshift QSO spectra (Maiolino et al. 2012; Cicone et al. 2015). In the stacked [C II] spectra of somewhat less extreme systems, the evidence of high-velocity wings range from suggestive (Gallerani et al. 2018) to undetected (Decarli et al. 2018) in QSOs, and strong

(Ginolfi et al. 2020) to undetected (even considering only galaxies with known molecular outflows, Spilker et al. 2020b) in star-forming galaxies. It is further uncertain to what extent [C II] traces ionized, neutral, and molecular gas in these outflows and at $z \sim 2$, the high-frequency observations needed to observe this fine-structure atomic line further complicate its use.

High-excitation water transitions present a promising probe of the dense warm molecular outflowing component but have so far only been observed in one starburst galaxy at $z = 5.656$ (Jones et al. 2019). Alternatively, one can utilize the blueshifted absorption features of outflowing gas situated between the observer and the host galaxy. The OH 119 µm doublet absorption line provides promising strength based on low-(e.g., Sturm et al. 2011; Spoon et al. 2013; Veilleux et al. 2013; Calderón et al. 2016; Stone et al. 2016; González-Alfonso et al. 2017) and high-redshift investigations (Zhang et al. 2018) but is only observable with ground-based facilities at redshifts $z > 4$ (e.g., Spilker et al. 2018, 2020a, 2020b).

Herschel SPIRE spectra of OH⁺ in local (U)LIRGs also reveal blueshifted absorption lines or (in a minority of cases) even complete P Cygni profiles (van der Werf et al. 2010; Rangwala et al. 2011; González-Alfonso et al. 2018). OH⁺ spectral lines lie at much lower frequencies than the OH 119 µm line, allowing us to probe them with ALMA Bands 3–7 at $z > 1.75$ or NOEMA Band 3 at $z > 2.74$, and at even lower $z$ with higher-frequency weather-sensitive bands. Additionally, the OH⁺($1_1-1_0$) line at 1033.1 GHz lies close to the CO(9–8) line at 1036.9 GHz, which traces warm dense gas in the host galaxy disk. Both OH⁺ and CO(9–8) can be observed with a single ALMA tuning, providing simultaneous velocity measurements of the host galaxy disk (in CO(9–8)) and any outflowing gas (through OH⁺ absorption, if blueshifted). Early observations of OH⁺ at high redshift have detected the line in absorption toward the massive starburst galaxy HFLS3 at $z = 6.34$ (Riechers et al. 2013). More recent observations of OH⁺ and H₂O⁺ in two $z \sim 2.3$ lensed SMGs SMM J2135–0102 and SDP 17b have been used to constrain the cosmic-ray ionization rate within these galaxies (Indriolo et al. 2018), finding rates much lower than predicted for the star-forming regions in these galaxies. The consequence of this finding is that OH⁺ likely traces the diffuse, turbulent, and predominantly neutral gas halos, also seen in CH⁺ (Falgarone et al. 2017), surrounding high-$z$ galaxies.

Recent observations of OH⁺ have further demonstrated the importance of this molecular ion to trace fueling and feedback in high-$z$ galaxies with the detection of inflowing gas via the three ground-state transitions of OH⁺ (together with CH⁺) in the star-forming galaxy HerBS-89a at $z = 2.95$ (Berta et al. 2021) and a powerful outflow traced by a P Cygni profile in the hyperluminous ($>10^{13} L_{\odot}$) starbursting merger ADFS-23 at $z = 5.655$ (Riechers et al. 2021a). As discussed in Berta et al. (2021) and shown by the results of Riechers et al. (2021b) who observed a sample of 18 starburst galaxies at $z = 2–6$, most OH⁺ measurements to date have revealed cases of outflow with only a few examples displaying clear detections of infall activity. In all 18 galaxies studied by Riechers et al. (2021b), OH⁺ is detected in either absorption (14), emission (10), or both (8).

Even in the brightest high-redshift sources, however, detecting and spatially resolving OH⁺ outflows remain observationally expensive (on the order of days) and thus limit analysis to
unresolved studies. Fortunately, observing time can be sig-
nificantly reduced with the aid of strong gravitational lensing,
which magnifies the light emitted by the background source,
increasing its on-sky size while maintaining surface brightness
and thereby improving the source plane resolution and total
observed flux.

This paper presents spatially resolved ALMA band 6
observations of the z = 2.0924 gravitationally lensed dusty
star-forming galaxy (DSFG; a galaxy selected at infrared or
submillimeter wavelengths, see Casey et al. 2014) HATLAS
J085358.9+015537 (hereafter G09v1.40), revealing a large-
scale neutral outflow traced by OH⁺ (and CH⁺) in addition to
the warm dense gas component and dust in the host galaxy
traced by CO(9–8) and underlying dust continuum emission,
respectively. Identified in the Herschel Astrophysical Terahertz
Large Area Survey, H-ATLAS (Negrello et al. 2010, 2017),
G09v1.40, and its lens have previously been modeled in the
submillimeter (Bussmann et al. 2013; Enia et al. 2018) and
near-infrared (NIR) (Calanog et al. 2014), providing excellent
input parameters for gravitational modeling with our high-
resolution ALMA observations. Additional studies of
G09v1.40 include accurate redshift determination and CO
SLED modeling using multiple CO transitions (Yang et al.
2017), analysis of the turbulent halo of diffuse gas surrounding
the galaxy seen in CH⁺ absorption (Falgarone et al. 2017), and
rest-frame optical spectral energy distribution modeling
including Spitzer/IRAC imaging at 3.6 and 4.5 µm Ma et al.
(2015).

Throughout our work we assume a flat ΛCDM cosmology
with Ω_m = 0.307 and H₀ = 67.7 km s⁻¹ Mpc⁻¹ (Planck
Collaboration et al. 2016). At the redshift of G09v1.40,
z = 2.0924, 1" corresponds to 8.53 kpc.

**2. Observations and Data Reduction**

**2.1. ALMA Band 6 Observations and Reduction**

In this paper, we present ALMA Band 6 data of the DSFG G09v1.40. The ALMA Band 6 observations of G09v1.40 were taken in ALMA Cycle 3 as part of project 2015.1.01042.S (PI: P. van der Werf). The observations were taken on 2016 April 22 using 36 antennas of the 12 m array with baseline lengths spanning 15–462 m. The uv-plane coverage provides sensitivity down to spatial scales of 3.5 kpc (and lower in the source plane) at z = 2.0924. The average precipitable water vapor level was 1.00 mm and the average system temperature of 168.8 K. G09v1.40 was observed for a duration of 21.60 minutes, with an additional 35.82 minutes allocated to phase calibration (J0909+0121), atmosphere and water vapor radiometry calibration (J0854+2006, J0909+0121, and J085358.9+015537), and bandpass, flux, and pointing calibration (J0854+2006).

The Band 6 receivers were tuned to observe the OH⁺(1ᵢ⁻₁₀) and CO(9–8) lines simultaneously in two slightly overlapping spectral windows to ensure continuous coverage of both lines without a decrease in sensitivity through the intermediate frequencies. This provided a 3.24 GHz bandwidth equating to ~2890 km s⁻¹. One more 2 GHz wide spectral window was placed to detect the rest-frame 1034 GHz continuum at high sensitivity. All spectral windows were configured with a channel resolution of 15.625 MHz.

The data were reduced with the ALMA Cycle 3 pipeline
using Common Astronomy Software Applications (CASA;

McMullin et al. 2007) version 4.53. We use Briggs weighting with a robust parameter 0, resulting in a beam with dimensions of 0.′′52 × 0.′′41 and position angle −78°2. This choice of weighting provides the optimal combination between sidelobe suppression and surface brightness sensitivity. The rest-frame 1034 GHz continuum map was created using 112 channels resulting in an rms = 0.12 mJy beam⁻¹. The OH⁺ (1₀−1₁) and CO(9−8) lines were separated into two data cubes with velocity resolutions of 14.5 km s⁻¹ and rms = 0.56 and 0.58 mJy beam⁻¹, respectively.

**2.2. Ancillary NIR Data**

We include ancillary NIR imaging of G09v1.40, captured by
the Keck II Near-Infrared Camera 2 (NIRC2) using the K
(λ = 2.2 µm) filter with laser guide star adaptive optics. The
observation and reduction of these data are presented in detail
by Calanog et al. (2014).

**2.3. Ancillary CH<sup>+</sup> Data**

To supplement our analysis of G09v1.40, we incorporate
two ancillary CH⁺(1−0) and rest-frame 836 GHz continuum
data sets observed with ALMA during Cycles 2 and 4. The
Cycle 2 observations (ALMA project 2013.1.00164.S; P.I. E.
Falgarone) have been discussed and analyzed by Falgarone
et al. (2017) with a focus on the CH⁺ line. In this study, we
simply utilize the pipeline product available from the ALMA
archive, reduced in CASA version 4.2.2 and imaged with a
Briggs weighting of 0.5, providing a beam size of
0.′′59 × 0.′′46.

The Cycle 4 data was observed as part of the 2016.1.00282.S ALMA program (P.I. E. Falgerone), aiming to measure the CH⁺(1−0) and dust continuum emission at a higher spatial resolution. We reduce the data using CASA version 4.7.0-1 with a Briggs weighting of 0, resulting in a beam size of 0.′′17 × 0.′′13. The CH⁺(1−0) line is detected in both absorption and emission, and the detected continuum emission provides the highest-spatial-resolution data of the dust profile in our analysis.

We refer to these data sets as the “low-” and “high-”spatial-resolution CH⁺(1−0) line and underlying rest-frame 836 GHz continuum data, respectively, throughout the paper.

**3. Results**

The ALMA Band 6 imaging successfully resolves and
detects the 1034 GHz dust continuum, CO(9–8) line emission
and the OH⁺(1ᵢ⁻₁₀) line seen in absorption. The dust
continuum reveals a complete Einstein ring and two bright
images, one to the east and a brighter one to the west
(Figure 1). This configuration is indicative of a single extended
source with a bright central dust region that is lensed by a
single, almost perfectly aligned foreground galaxy. The
brighter western image indicates that the source sits just to
the west of the foreground lens.

The CO(9–8) similarly displays a complete Einstein ring; however, in contrast to the continuum, the CO(9–8) emission reveals three distinct peaks. The southwestern, and brightest, peak is separated into two emission peaks by the deconvolved model provided by the cleaning procedure (see Figure 20 in Appendix A for more discussion). The spectra observed across all peaks (outer panels of Figure 1) do not display significant variation in shape or central velocity, suggesting that we are
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Middle row—far-left panel: rest-frame 1034 GHz dust continuum emission observed with ALMA. The peak of the dust emission is lensed into a prominent double-image configuration with a brighter western image, accompanied by a complete Einstein ring. The synthesized beam (0''52 × 0''41) is shown by the faded orange ellipse in the bottom left. The gray contours indicate continuum levels of 20%, 40%, 60%, 80%, and 100% that of the peak and are overlaid in the following three panels. Middle row—center-left panel: intensity map of the rest-frame 1036.9 GHz CO(9–8) emission line observed with ALMA, produced via a single Gaussian fitting procedure. The CO(9–8) emission reveals itself in three distinctly different peaks compared to the continuum but again accompanied by a complete Einstein ring. Middle row—center-right panel: integrated optical depth map of the blueshifted component of the OH⁺(1₁–1₀) absorption line observed with ALMA, produced via a double Gaussian fitting procedure. The peak OH⁺ optical depth is approximately spatial with the peaks in the continuum emission albeit with an elongated morphology stretching from east to west. Middle row—far-right panel: Keck *Ks* (λ = 2.2 µm) stellar intensity map with the lens subtracted (Calanog et al. 2014). As in the continuum and CO(9–8) emission, a distinct Einstein ring is observed, but with gaps at the positions of peak FIR continuum emission. Outer four panels: Each of the four outer panels corresponds to a region with a matching label (E: east, NW: northwest, SW: southwest, W: west) indicated by the solid white boxes in the central two panels and are placed over peaks in the continuum and CO(9–8) emission. The outer panels display continuum-subtracted CO(9–8) and OH⁺ spectra of the individual spaxels (pixel area = 0''05 × 0''05) of each region in transparent blue and green, respectively. The average spectra of each region are shown by the solid line. The CO(9–8) and OH⁺ spectra do not display significant variation in spectral shape or central velocity between the four regions.</figcaption></figure>

seeing the same CO(9–8) component lensed into four images, of which two are blended together in the southwest. This requires the peak of the CO(9–8) component to be lying directly over a section of the inner caustic, where the magnification of the source goes to infinity (diamond in bottom two rows of Figure 6).

The two peaks in OH⁺ optical depth (and hence column density) are approximately spatial with the two continuum images, albeit with elongated morphologies stretching from east to west and dropping off in intensity to the north and south. As for CO(9–8), the OH⁺ spectra do not display significant variation in spectral shape or central velocity across its images. Unlike both the dust continuum and CO(9–8), there is no discernible Einstein ring in the OH⁺ optical depth, indicating already in the image plane that the alignment of the OH⁺ component in the source plane does not lie, even partially, over the inner caustic and must therefore be located fully to the west of it.

For comparison, we also include the lens-subtracted Keck near-IR stellar intensity map from Calanog et al. (2014). There are no distinct peaks but a clear Einstein ring indicates that the stellar component is lying directly over the inner caustic.

The close spectral proximity of the OH⁺(1₁–1₀) and CO(9–8) lines allows us to capture both transitions

simultaneously, with a single ALMA tuning. This is highlighted in Figure 2, displaying a single, typical spaxel in our data set. Remarkably, at every position, the OH⁺ absorption peak is blueshifted by ~200 km s⁻¹, with respect to the central velocity of the CO(9–8) emission line. Because the redshift determined from the CO(9–8) emission is consistent with the multiline redshift, z = 2.0924 ± 0.0001, precisely determined by Yang et al. (2017)), and because we detect the OH⁺ line in absorption and are therefore tracing gas located in front of the dust continuum, the blueshift of the OH⁺ absorption indicates that it is tracing gas outflowing from the host galaxy toward us.

To create the CO(9–8) emission and OH⁺(1₁–1₀) optical depth maps presented in Figure 1, we feed the CO(9–8) and OH⁺ data cubes through a robust Gaussian spectral line fitting algorithm whereby each spaxel is individually fitted. We find, by examining the residuals produced by this fitting process, that a Gaussian, or combination of Gaussian profiles, provides a good fit to the observed spectra. The CO(9–8) and OH⁺(1₁–1₀) spectral lines do not overlap in frequency at any source location and so this process is performed separately for each of the spectral line data cubes. The fitted intensity (*S*), velocity (*V*), and velocity dispersion (σ_v) of each spaxel then form a 2D data map in the R.A.–decl. plane. This technique is
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Continuum-subtracted OH⁺(1₁⁻₁₀) absorption (green) and CO(9⁻⁸) emission (blue) spectra in a single ALMA spaxel (pixel area = 0''⁰⁵ × 0''⁰⁵) observed simultaneously in the DSFG G09v1.40 Top: spectra plotted as a function of observed frequency. Bottom: spectra plotted as a function of velocity with respect to the galaxy’s systemic velocity determined from the z = 2.0924 redshift measurement by Yang et al. (2017). We fit the OH⁺ (1₀⁻₁₁) absorption and CO(9⁻⁸) emission lines with double and single Gaussian functions, respectively; the shaded regions of the same color indicate the individual Gaussian components of each spectral line, while the dashed green line presents the full OH⁺ fit. The main OH⁺ component is blueshifted ~200 km s⁻¹ with respect to the bulk molecular gas as traced by the CO(9⁻⁸) emission, revealing a neutral outflow at these velocities. Note also the larger width of the OH⁺ line.</figcaption></figure>

favored over the use of moment maps as spectral fitting is capable of cleanly disentangling separate overlapping velocity components, as discussed below.

We attempted to fit one, two, and three Gaussian functions to the CO(9⁻⁸) spectra in order to fit the more complex spectral shape, as seen in Figure 1. The multi-Gaussian fits, however, did not produce smooth velocity or intensity fields, indicating that this method was not capable of extracting separate kinematic CO components. We therefore apply single Gaussian fits over the full source, finding this to best trace the bulk gas component. The intensity, velocity, and velocity dispersion of CO(9⁻⁸) were left as free parameters in this process.

The OH⁺ spectrum is relatively consistent across the entire source, with absorption typically peaking at velocities ~200–300 km s⁻¹ blueshifted with respect to the source’s systemic velocity (Figure 2). Additional absorption at systemic velocities is responsible for the skewed spectra at all source locations (Figure 1) and is cleanly disentangled from the outflowing component by our spectral fitting described below.

We first fit a double Gaussian function to the OH⁺ absorption spectra where one Gaussian is fixed at the systemic velocity and the velocity of the blueshifted component is left free (intensity and velocity dispersion are also left free for

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Top left: image plane velocity map of the CO(9⁻⁸) emission produced via a single Gaussian fitting procedure described in the text. Top right: image plane velocity map of the OH⁺(1₁⁻₁₀) blueshifted absorption component, produced via a double, Gaussian fitting procedure described in the text. Bottom left: image plane velocity dispersion map of the CO(9⁻⁸) emission produced via a single Gaussian fitting procedure described in the text. Bottom right: image plane velocity dispersion map of the OH⁺(1₁⁻₁₀) blueshifted absorption component produced via a double Gaussian fitting procedure described in the text. The beam (0''⁰⁵ × 0''⁰⁴¹) is shown by the faded orange ellipse in the bottom-left corner of the upper-left panel. Contours indicate flux levels at 20%, 40%, 60%, and 80% of the peak flux and optical depth of the CO(9⁻⁸) and OH⁺ spectral lines, respectively.</figcaption></figure>

both lines). This provides the central velocity and velocity dispersions of the OH⁺ in each spaxel. We then convert the OH⁺ absorption data cube into an optical depth cube, τ, via

$$ \tau_{i,j}(\nu) = -\ln\left(\frac{S_{\text{abs},i,j}(\nu)}{S_{\text{cont},i,j}}\right), \qquad (1) $$

where S<sub>abs,i,j</sub>(ν) is the absolute value of the OH⁺ flux density in pixel [i, j] at frequency ν and S<sub>cont,i,j</sub> is the continuum flux in the same pixel. A double Gaussian function is then fitted to the OH⁺ optical depth profiles where again, one Gaussian is fixed at the systemic velocity and the other is fixed at the central velocity found in the previous fit to the absorption spectra. This second fitting step provides us with integrated optical depths for each spaxel.

We discard CO(9⁻⁸) and OH⁺ spaxel fits that return flux amplitudes below the noise level of each cube (with integrated fluxes reaching S/Ns much higher than 1), velocities outside the observed bandwidth, or velocity dispersion narrower than two channels or wider than the full velocity bandwidth. OH⁺ spaxels are further rejected if the continuum flux is below 5σ. The remaining CO(9⁻⁸) and OH⁺ spectral fits (e.g., shaded Gaussians in the bottom panel of Figure 2) are then used to make clean CO(9⁻⁸) intensity and blueshifted OH⁺ optical depth maps as presented in Figure 1 along with their respective velocity and velocity dispersion maps, shown in Figure 3. Note that we only present and analyze the blueshifted component of the OH⁺ because we aim to study the outflow.
The maximum blueshift in the OH⁺ line is found along the peak in OH⁺ optical depth (indicated by the contours in Figure 3). The velocity dispersion similarly peaks along the main outflow component, with both the projected outflow velocity and velocity dispersion increasing mildly toward the far western and eastern edges, as seen in the optical depth. These regions are likely where the outflow is at its widest; however, it must be noted that the background continuum in these extreme regions is faint and care should be taken when interpreting the OH⁺ absorption.

# 4. Gravitational Lens Modeling and Source Plane Reconstruction

## 4.1. Lens Modeling: VISILENS

Because we are interested in the intrinsic properties of G09v1.40 and its outflow, we must first model and remove the effects of gravitational lensing. ALMA, as an interferometer, observes the Fourier transform of the sky intensity distribution over a range of two-dimensional spatial frequencies in the uv plane. Noise properties and resolution effects are much better understood in the uv plane than in the inverted images where uncertainties become correlated and may bias further measurements. To avoid such biases affecting our lens modeling, we have chosen to employ the parametric reconstruction code visilens (Hezaveh et al. 2013; Spilker et al. 2016), which fits a lens model directly to visibility measurements. We invoke the `modelcal` option available in visilens, which corrects for calibration errors caused by, e.g., uncertain antenna positions and atmospheric conditions, allowing for multiplicative amplitude rescaling and astrometric drift. If present and not corrected for, these calibration errors can result in a shift of the model parameters away from their intrinsic values. See Hezaveh et al. (2013) and Spilker et al. (2016) for more in-depth discussions of the code.

To improve the accuracy of our lens model further, we take advantage of the underlying 836 GHz dust continuum emission in the “high-”spatial-resolution CH⁺(1–0) data. The higher spatial resolution of this data provides visilens more information on the dust distribution in G09v1.40 from which a more accurate model can be derived (Figure 4).

We use a single Sérsic profile to represent the continuum emission of G09v1.40, characterized by a Sérsic index $n_s$, half-light radii $a_s$ with axis ratio $b_s/a_s$, position angle east of north $\phi_s$, flux density $F_s$, and position $\Delta x_s, \Delta y_s$ with respect to the lens. We model the mass profile of the lens with a singular isothermal ellipsoid (SIE), fitting for the lens position $x_L, y_L$ with respect to the ALMA phase center, mass $M_L$ (and corresponding Einstein radius, $\theta_{EL}$, within which the mass is parameterized), ellipticity $\epsilon_L$, and position angle, east of north $\phi_L$. We do not invoke any external shear in our model.

visilens begins by creating a 2D source plane parametric model of the source dust emission and lenses this model into the image plane for a given lens model. The two-dimensional lensed emission is then Fourier transformed into the uv plane where it is directly compared to our observed interferometric data. We initiate this procedure with values taken from Bussmann et al. (2013), and a Markov Chain Monte Carlo (MCMC) sampling algorithm explores the model parameter space of both the source emission and lens mass profiles, using the `emcee` code (Foreman-Mackey et al. 2013). For each point in the parameter space, a source plane, image plane, and

<figure><img src="image_5.png" /><figcaption>**Figure 4.** visilens modeling of the high-spatial-resolution CH⁺(1–0) continuum emission in G09v1.40. Top left: dirty image of the data with beam (0′′17 × 0′′13) shown as the orange ellipse in the bottom-left corner. Top right: dirty image of the visilens model, recovering 85% of the flux in the dirty image. Bottom: residuals, with the rms beam⁻¹ shown in the bottom-left corner.</figcaption></figure>

uv-plane visibilities are generated and checked for quality of fit to the data using a $\chi^2$ metric, where the best-fit parameters minimize the $\chi^2$ value. Note that the `emcee` routine is known to underestimate uncertainties in some circumstances and may be the cause of the extremely small uncertainty provided for our source position.

From the high-resolution data, we extract 280 channels of continuum over the three available spectral windows to model, for which visilens finds the best-fit lens and source parameters presented in Table 1. The model recovers the bulk (85%) of the dirty image flux, as shown in Figure 4. The residuals exhibit structure that are not consistent with the rms, suggesting that the dust continuum of the source contains more complex structure than can be captured by a single Sérsic profile. We attempted to model the source with two Sérsic profiles but did not find significant changes to the lens model. Because we are primarily interested in the lens model and do not use the parametric models of the source in our following analysis, we opt to use the single-source model.

## 4.2. Source Reconstruction: LENSTOOL

From the best-fit lens model obtained in Section 4.1 we can derive a lens operator that maps our data from the image plane into the source plane. For this task we employ the pixelated reconstruction code `LENSTOOL` (Kneib et al. 1996; Jullo et al. 2007; Jullo & Kneib 2009). A pixelated reconstruction has the advantage of capturing details that cannot be easily modeled by a parametric model. This is of particular interest in the case of an outflowing component, which does not need to follow a Sérsic profile and which is implied by the unusual OH⁺ optical depth morphology.

The reconstructed dust continuum emission, from the “high-”spatial-resolution CH⁺(1–0) data set used to model the gravitational lens, reveals a single elongated and compact dust
<caption>Table 1</caption>

visilens Parameters of the Best-fit SIE Lens and Sérsic Source Models, Using the Underlying Dust Continuum of the High-spatial-resolution CH$^{+}$(1-0) Observations

<table><thead><tr><th></th><th>Lens (SIE)</th><th></th></tr></thead><tbody><tr><td>x<sub>L</sub></td><td>["]</td><td>2.824 ± 0.677</td></tr><tr><td>y<sub>L</sub></td><td>["]</td><td>2.194 ± 0.042</td></tr><tr><td>M<sub>L</sub></td><td>[10<sup>11</sup>M<sub>☉</sub>]</td><td>0.459 ± 0.119</td></tr><tr><td>e<sub>L</sub></td><td></td><td>0.037 ± 0.042</td></tr><tr><td>θ<sub>L</sub></td><td>[deg] CCW from E</td><td>162.933 ± 23.426</td></tr></tbody><thead><tr><th></th><th>Source (Sérsic)</th><th></th></tr></thead><tbody><tr><td>Δx<sub>S</sub></td><td>["]</td><td>0.087 ± 0.029</td></tr><tr><td>Δy<sub>S</sub></td><td>["]</td><td>0.001 ± 0.001</td></tr><tr><td>F<sub>S</sub></td><td>[mJy]</td><td>2.971 ± 1.097</td></tr><tr><td>a<sub>S</sub></td><td>["]</td><td>0.058 ± 0.009</td></tr><tr><td>b<sub>S</sub>/a<sub>S</sub></td><td>["]</td><td>0.679 ± 0.120</td></tr><tr><td>n<sub>S</sub></td><td></td><td>4.805 ± 1.638</td></tr><tr><td>φ<sub>S</sub></td><td>[deg] CCW from E</td><td>164.264 ± 52.123</td></tr></tbody></table>

**Note.** Lens positions are given with respect to the ALMA phase center: (J2000) 08:53:58.68+01:55:35.45 and source positions with respect to the lens position. Parameter descriptions are as follows: $x_L$, lens position in R.A.; $y_L$, lens position in decl.; $M_L$, lens mass inside the Einstein radius; $e_L$, lens ellipticity; $θ_L$, lens position angle; $Δx_S$, source position in R.A.; $Δy_S$, source position in decl.; $F_S$, source flux; $a_S$, source major axis; $b_S/a_S$, source axis ratio; $n_S$, source Sérsic index; and $φ_S$, source position angle.

continuum profile (Figure 5). The weak feature directly to the east of the inner caustic is a result of beam smearing across the outer caustic in the image plane. When reconstructed, the flux that has been smeared to the opposite side of the outer caustic is traced to the wrong side of the inner caustic in the source plane. This effect becomes more severe with larger beam sizes (Figure 19) and is more difficult to interpret when the galaxy component is lying directly over the caustic line (Figure 20; see Appendix A for a detailed discussion).

We use *LENSTOOL* to reconstruct the OH$^{+}$, CO(9–8), underlying 1034 GHz dust continuum and ancillary Keck NIR data using the same lens model (Figure 6). As in the high-resolution continuum data, the rest-frame 1034 GHz dust continuum emission reveals a single compact and elongated source. The eastern artifact seen in the high-resolution continuum reconstruction is also present here but at higher significance (39% of the source plane flux as opposed to 9% for the high-resolution data). The more prominent eastern artifact in the 1034 GHz dust continuum reconstruction is due to the larger beam size and thus enhanced beam smearing (see Appendix A). We examine the model intensity maps of the 1034 GHz dust continuum produced during cleaning, which when reconstructed reveals only a single western source (no eastern artifact). This is consistent with our initial interpretation of the image plane dust continuum emission as a single, double-imaged extended source. The brighter eastern image indicates that the peak in the dust continuum lies outside and to the west of the inner caustic, with the eastern edge of the extended emission crossing the caustic and forming the faint Einstein ring. Because the continuum peaks outside of the inner caustic, the effects of beam smearing on the source plane reconstruction appear to be straightforward and can be rectified by simply masking the eastern artifact (maps shown in Figure 6 display maps with the eastern artifact removed). We do this for

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

**Figure 5.** Source plane reconstruction of dust continuum emission from the “high”-spatial-resolution CH$^{+}$(1-0) data set using *LENSTOOL*. The left panel displays the observed image plane dust continuum intensity with the beam (0′′17 × 0′′13) shown in the lower left by the shaded orange ellipse. The right panel is a zoom in of the white dashed region in the left panel, displaying the source plane reconstruction of the dust continuum. The lens caustics are indicated by the solid white lines and a physical scale of 1 kpc is given in the lower-right corner. A single compact and elongated source is revealed directly to the west of the inner caustic. The weak eastern feature is an artifact of beam smearing over the outer caustic in the image plane.

all the dust continuum maps (including those from other data sets) in the rest of our analysis.

In agreement with the distinct CO(9–8) image plane morphology, a distinct source plane morphology is found. The CO(9–8) emission forms a “three-pronged” morphology connecting over the inner caustic. The image plane quadruply imaged CO(9–8) emission indicated that the CO(9–8) component crosses the inner caustic, which is indeed seen in the reconstruction. The striking source plane morphology, however, is not already evident in the image plane and is likely a manifestation of beam smearing of the image plane. We investigate the effects of beam smearing on the source plane reconstruction by taking the model produced by the cleaning procedure and convolving it with different beam sizes (see Appendix A, Figure 20). In particular, the northwestern “prong” of the source plane morphology is deemed to be an artifact of beam smearing, but overall the effects of the beam on the small-scale source plane structure are complex. For this reason, we do not attempt any source plane masking of the reconstructed CO(9–8) emission. The more extended distribution and spatial offset of the CO(9–8) emission with respect to the dust continuum are already evident in the image plane and can be interpreted as intrinsic differences between the CO(9–8) and dust components. A more detailed analysis of the source structure would require higher-spatial-resolution observations.

By construction the OH$^{+}$ absorption is only observed where there is background continuum; however, the reconstructed OH$^{+}$ optical depth (and thus column density) exhibits a distinctively different source plane morphology from that of all the other components. We find an elongated triangular OH$^{+}$ distribution, with its central axis running from the east and flaring out toward the west. Unlike the CO(9–8) emission, the elongated morphology of the OH$^{+}$ optical depth is already evident in the image plane distribution. In particular, the steep falloff in optical depth to the north and south, compared to the background dust continuum, is indicative of a sharp physical transition and is interpreted as a true feature. Despite the background continuum displaying a full Einstein ring, this feature is absent in the OH$^{+}$ optical depth, indicating that,
<figure><img src="image_7.png" /><figcaption>**Figure 6.** ALMA Cycle 3 observations of the gravitationally lensed galaxy HATLAS J085358.9+015537 at redshift $z = 2.09$. Row 1: image plane intensity maps of the dust continuum, CO(9–8) emission tracing the warm dense gas in the disk, OH⁺(1₁–1₀) optical depth tracing neutral outflowing gas, and Keck Ks-band stellar emission. The beam (0.′′52 × 0.′′41) is shown in the bottom-left corner of the left-hand panel. Row 2: reconstructed source plane intensity maps of the components shown in row 1. Lensing caustics are shown in white. Row 3 Left panel: false-color image of all the source plane intensity maps of the dust (red), CO(9–8) (green), OH⁺ (white), and stars (blue). A single contour of each distribution is drawn at 60% of their respective peak values. Row 3: reconstructed source plane velocity maps. The second panel shows the velocity map of the CO(9–8) emission, and the third panel shows the OH⁺ velocity field. The first panel is the OH⁺ velocity map minus the CO(9–8) velocity map, which, if the outflow was a flat sheet of gas lifting off the galaxy disk, should be uniform. This is clearly not the case and this scenario is therefore ruled out (see Section 5.1). Pixel area (0.′′05 × 0.′′05) is the same in both the image and source planes and a source plane physical scale of 1 kpc is given in the lower-right corner of the last panel.</figcaption></figure>

unlike the dust, CO(9–8), and stellar components, the OH⁺ distribution does not cross the inner caustic. This is further evidence that the eastern artifact seen in the continuum reconstruction, and which is also found in the OH⁺ reconstructions, is an artifact caused by beam smearing in the image plane. We mask this feature in the source plane OH⁺ maps in the same fashion as for the continuum. Lastly, the OH⁺ optical depth appears to increase toward and peak in the west. Higher sensitivity and spatial-resolution observations are required to show if this trend is indeed true, particularly in the extreme western edge where the background continuum, and therefore OH⁺ spectral S/N, decreases significantly. It is likely, however, that the OH⁺ distribution extends past the background continuum in this direction, causing a sharp cutoff in the OH⁺ distribution visible to us via absorption.

The Keck NIR emission, which shows a similar image plane morphology to the CO(9–8) emission but with a more prominent Einstein ring, displays a similarly extended source plane component, partially crossing the inner caustic. The peaks in stellar and CO(9–8) emission are not cospatial, with the stellar peak lying to the southwest of the inner caustic, below that of the peak in dust continuum and OH⁺ optical depth. We do not apply any source plane mask to the NIR emission as beam smearing is not an issue with this data set.

In the lower-left panel of Figure 6 we display a false-color image of all the reconstructed intensity maps. We highlight the bulk offsets of each distribution, in particular the dip in stellar light at the position of the peak dust continuum. Spatial offsets between the gas, dust, and stellar components of galaxies at high redshift have previously been observed in many cases
<table><caption>**Table 2**
Intrinsic Source Properties</caption><tr><td>L<sub>FIR</sub></td><td>[10<sup>11</sup>L<sub>☉</sub>]</td><td>45.5 ± 25</td></tr><tr><td>L'<sub>CO(9–8)</sub></td><td>[10<sup>9</sup> K km s<sup>-1</sup> pc<sup>2</sup>]</td><td>3.1 ± 1.2</td></tr><tr><td>L<sub>NIR</sub></td><td>[μJy]</td><td>1.7 ± 0.7</td></tr><tr><td>μ<sub>cont. model</sub></td><td></td><td>11.6 ± 4.5</td></tr><tr><td>μ<sup>high res cont</sup></td><td></td><td>14.5 ± 5.6</td></tr><tr><td>μ<sub>836 GHz</sub></td><td></td><td>11.1 ± 4.3</td></tr><tr><td>μ<sub>CO(9–8)</sub></td><td></td><td>7.7 ± 3.0</td></tr><tr><td>μ<sup>OH+</sup></td><td></td><td>9.3 ± 3.6</td></tr><tr><td>μ<sub>NIR</sub></td><td></td><td>11.4 ± 4.4</td></tr><tr><td>μ<sup>CH+</sup></td><td></td><td>10.0 ± 3.9</td></tr><tr><td>r<sub>eff,cont.model</sub></td><td>[pc]</td><td>408 ± 7.0</td></tr><tr><td>r<sub>eff,NIR</sub></td><td>[pc]</td><td>1100 ± 2.0</td></tr><tr><td>SFR</td><td>[M<sub>☉</sub> yr<sup>-1</sup>]</td><td>788 ± 300</td></tr><tr><td>Σ<sub>SFR</sub></td><td>[M<sub>☉</sub> yr<sup>-1</sup> kpc<sup>-2</sup>]</td><td>753 ± 290</td></tr></table>

Note. Parameter descriptions are as follows: $L_{FIR}$, lens-corrected total infrared luminosity (8–1000 µm) using the $\mu_{L_{FIR}}$ from Bussmann et al. (2013) and the magnification factor of the high-spatial-resolution CH$^{+}$(1–0) dust continuum. $L_{CO(9–8)}$, delensed CO(9–8) line luminosity. $M_{NIR}$, delensed NIR magnitude using the apparent NIR luminosity from Calanog et al. (2014) and NIR magnification factor derived in this work. $\mu_{cont, model}$, magnification factor provided by the *visilens* model of the OH$^{+}$ underlying continuum. $\mu_{high res cont}$, magnification factor of the spatial resolution of the CH$^{+}$(1–0) dust continuum. $\mu_{836 GHz}$, magnification factor of the OH$^{+}$ underlying dust continuum. $\mu_{CO(9–8)}$, magnification factor of the CO(9–8) component. $\mu_{OH+}$, magnification factor of the masked, outflowing OH$^{+}$ component. $\mu_{NIR}$, magnification factor of the NIR stellar component. $\mu_{CH+}$, magnification factor of the low-spatial-resolution CH$^{+}$(1–0) component. $r_{eff,cont, model}$, effective radius of the dust continuum derived from the *visilens* model parameters, where $r_{eff} = a_S \sqrt{b_S/a_S} \cdot r_{eff,cont, model}$, the NIR effective radius from Calanog et al. (2014). SFR, lens-corrected SFR derived from the $L_{IR}$ using the Kennicutt (1998) calibration, $SFR = 1.73 \times 10^{-10} L_{IR} M_{☉} yr^{-1}$ and assuming a Salpeter IMF, $\Sigma_{SFR}$ lens-corrected SFR surface density.

(e.g., Riechers et al. 2010; Hodge et al. 2015; Rybak et al. 2015; Chen et al. 2017; Fujimoto et al. 2017; Simpson et al. 2017; Calistro Rivera et al. 2018; Cochrane et al. 2021). The most likely explanation for the offset between stellar and dust components is that the optical component does indeed extend into the dusty regions but experiences high extinction. This scenario is consistent with the dip in stellar intensity observed in G09v1.40 (evident already in the image plane distributions). Additional observations would be needed, however, to concretely justify this claim, and we cannot reject the possible scenario of true physical misalignments between gas, dust, and stars or the possible scenario of two merging systems, one with extreme dust extinction and another optically bright galaxy.

For each component we measure the image plane and source plane luminosities, providing a value for the magnification specific to each component (Table 2). Source plane luminosities are measured with the eastern artifacts removed for the continuum and OH$^{+}$ maps. We additionally provide the SFR and SFR surface density derived from the $L_{IR}$ measured by Bussmann et al. (2013) and corrected using our new dust continuum magnification factor and source size. Because the source size is the same for both our model and that of Bussmann et al. (2013), this results in a change of 6% for both values.

We reconstruct the CO(9–8) and blueshifted OH$^{+}$ velocity and velocity dispersion maps. The OH$^{+}$ line is significantly blueshifted with respect to the CO(9–8) at all locations across the source and displays an opposite velocity gradient across the

northeast to southwest axis as expected from the image plane velocity maps. Further insight into the intrinsic velocity structure of the CO(9–8) emission should not be read from the reconstructed CO(9–8) velocity map due to beam smearing, blending components of the CO(9–8) emission on opposite sides of the inner caustic together. Indeed, the source plane CO(9–8) velocity field appears comparatively more chaotic than what would be expected from the smooth image plane velocity field (Figure 3). This may in turn disguise signatures of rotation in the host galaxy, if present. Fortuitously, this is not an issue in the OH$^{+}$ maps because the entire OH$^{+}$ component lies to the west of the inner caustic, and the negative gradient to the west, and positive gradients to the north and south are interpreted as real kinematic features.

### 4.3. Comparison with Previous Lens Models

Gravitational lens models of G09v1.40 have been previously derived by Bussmann et al. (2013) and Calanog et al. (2014) using 880 µm Submillimeter Array (SMA) and Keck II Near-Infrared Camera 2 (NIRC2) 2.2 µm observations, respectively. Our lens parameters are consistent with those measured by Bussmann et al. (2013) and Calanog et al. (2014), summarized in Table 3. The contrast in image plane morphology between the NIR and submillimeter SMA imaging (Einstein ring and double image, respectively) was interpreted by Calanog et al. (2014) as a consequence of poor spatial resolution in the submillimeter data compared to the Keck AO. With our high-spatial-resolution analysis of the dust continuum in this paper, it is now clear that contrasts in image plane morphology are due to intrinsic differences in the source plane morphology of these components. The almost perfect alignment of G09v1.40 with the lensing galaxy means that small variations and offsets of the stellar, dust, and gas components in the source plane produce strikingly different image plane morphologies and magnifications.

## 5. Outflow Geometry

With the source plane reconstructed maps in hand, we now investigate possible geometries of the outflowing gas. Given the limitations in spatial resolution of our data, we compare the suitability of three simple outflow geometries: a flat sheet lifting off a star-forming disk, a spherical outflow originating and expanding from a single location in the galaxy, and a conical outflow.

### 5.1. Sheet

For a galactic disk with extended star formation, it is easy to imagine a flat sheet of outflowing gas lifting off perpendicularly from the disk. In this scenario the velocity signature of the disk can be imprinted onto that of the outflow. Both velocity maps will therefore exhibit the same velocity gradients, albeit offset in the R.A. and decl. plane depending on the inclination of the disk and height of the outflow. It is obvious directly from the image plane CO(9–8) and OH$^{+}$(1$_{1}$–1$_{0}$) velocity maps (Figure 3) that the molecular gas in the host galaxy and outflowing neutral gas display opposite velocity gradients across their 2D projections. The uncertainty in the velocity fields is likely less than the velocity resolution (<14.5 km s$^{-1}$) of the spectra and much less than the observed velocity gradients, given the high integrated S/N of each spaxel. The observed velocity gradients seen in the CO(9–8) and OH$^{+}$ are
<table><caption>**Table 3**

Best-fit Parameters from Previous Gravitational Lens Models by Bussmann et al. (2013) and Calanog et al. (2014)</caption><thead><tr><th>θ<sub>E</sub><br>(")</th><th>e<sub>L</sub></th><th>φ<sub>L</sub><br>(deg) E of N</th><th>n<sub>S</sub></th><th>a<sub>S</sub><br>(")</th><th>e<sub>S</sub><br>(deg)</th><th>r<sub>eff</sub><br>(kpc)</th><th>φ<sub>S</sub><br>(deg) E of N</th><th>μ<sub>source</sub></th></tr></thead><tbody><tr><td colspan="9">Bussmann et al. (2013; SMA 880 μm source and lens model)</td></tr><tr><td>0.553 ± 0.004</td><td>0.06 ± 0.02</td><td>70 ± 12</td><td>2 ± 0.7</td><td>0.06 ± 0.01</td><td>0.33 ± 0.14</td><td>0.41 ± 0.08</td><td>83 ± 17</td><td>15.3 ± 3.5</td></tr><tr><td colspan="9">Calanog et al. (2014; NIR 2.2 μm source and lens model)</td></tr><tr><td>0.56<sup>+0.01</sup><sub>-0.02</sub></td><td>0.0<sup>+0.1</sup><sub>-0.2</sub></td><td>-57<sup>+4</sup><sub>-1</sub></td><td>0.51<sup>+0.02</sup><sub>-0.04</sub></td><td>0.18<sup>+0.01</sup><sub>-0.01</sub></td><td>0.49<sup>+0.02</sup><sub>-0.06</sub></td><td>1.1 ± 0.002</td><td>87<sup>+6</sup><sub>-4</sub></td><td>11.4<sup>+0.9</sup><sub>-1</sub></td></tr><tr><td colspan="9">Calanog et al. (2014; NIR source model using SMA lens model)</td></tr><tr><td></td><td></td><td></td><td></td><td>0.18<sup>+0.01</sup><sub>-0.01</sub></td><td>0.51<sup>+0.03</sup><sub>-0.1</sub></td><td>1.1 ± 0.002</td><td></td><td>10<sup>+1</sup><sub>-1</sub></td></tr></tbody></table>

**Note.** Parameter descriptions are as follows: $θ_E$, Einstein radius of the lens. $e_L$, lens ellipticity. $θ_L$, lens position angle. $a_S$, source major axis. $e_S$, source ellipticity (where $e_S = 1 - b_S/a_S$ and $b_S/a_S$ is the axis ratio). $n_S$, source Sérsic index. $φ_S$, source position angle. $μ_S$, magnification of the source.

therefore believed to be truly disparate, and we dismiss this outflow geometry.

## 5.2. Spherical

In this scenario, we explore a spherical shell expanding uniformly from a single location, such as an AGN or central star-forming region. To explore this geometry we have created a simple toy model that converts an expanding 3D spherical shell into a 2D projected velocity map. As we are using an absorption line in our analysis and therefore only probe gas in front of the galaxy, we consider only the front-facing hemisphere in our model. The process of creating and comparing the model outflow with the observed $OH^+$ outflow is as follows.

To build a hemispherical outflow with radius $R$ and finite thickness $dR$, we first create a 3D box with dimensions $2(R+1) \times 2(R+1) \times (R+1)$ and grid it with 10 times the spatial resolution of our observed source plane data maps ($0''05$ or $\sim 0.43$ kpc). In this model the third, shorter axis is parallel with the line of sight (LOS), and the center of the outflow is placed at the farthest distance along this axis from the observer and in the center of the other two axes such that the full outflow fits within the box. Grid elements of the box that lie within $R$ and $RdR$ are then assigned values equal to their LOS positions. The box is then collapsed and averaged along the LOS axis to create a 2D map of the average LOS position for each location across the face of the outflow. Because we assume the outflow is moving radially outwards, this 2D map has an identical gradient to the average deprojected outflow velocity and can thus be used to compare to the 2D deprojected velocity map of the observed outflow (see the bottom-right panel of Figure 7).

Because in reality, we can only observe the outflowing gas situated in front of the dust continuum we must then select small regions of the model outflow velocity field to compare with our data. This is done systematically: moving pixel by pixel across the face of the projected model velocity field, we cut regions with matching pixel dimensions of the observed reconstructed $OH^+$ velocity map. This allows us to determine where the observed outflow may lie with respect to the ejection point. Once cropped, the model velocity field is spatially averaged to match the resolution of the observed field and both maps are normalized such that only the gradient of the fields are compared. A residual map and $\chi^2$ value are derived for each region selected from the model velocity map, allowing us to find the most probable position of our observed outflow with respect to the ejection location (see Figure 7 for an example).

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Example spherical outflow model compared to the observed $OH^+(1_1-1_0)$ velocity map. The model suggests the outflow is ejected from outside the galaxy and is therefore ruled out. Top left: reconstructed and normalized velocity map of the observed outflow in $OH^+$. Top right: spatially averaged model velocity map from the region in the dashed box in the panel below. The yellow marker indicates the ejection location predicted by this model. Bottom left: residual velocity map of the normalized model velocity map and normalized $OH^+$ velocity map. Bottom right: full 2D LOS velocity map of the spherical outflow model with a dashed box outlining the edges of the panel above.</figcaption></figure>

We find that increasing the radius of the toy model monotonically reduces the $\chi^2$ value of the best-fit velocity map. This in turn predicts a larger and larger distance of the outflow ejection point with respect to the observed outflow, placing it far outside the host galaxy (Figure 7). Adopting a spherically symmetric model would therefore imply accepting the unphysical situation that the origin of the outflow lies far outside the host galaxy and the model is therefore rejected. We further note that a spherical outflow with an ejection point directly at the peak of the dust continuum would produce a deprojected velocity field with the peak in projected outflow velocity at the same position. This is not the case for the observed $OH^+$ velocity field (see Figure 6) and is obvious already in the image plane velocity fields (Figure 3).
### 5.3. Conical

Following the failure of the sheet and spherical geometries in characterizing the observed outflow morphology in G09v1.40, we consider a conical outflow geometry. In this scenario, gas is ejected from a localized region in the host galaxy and expands as it flows radially away from the galaxy. We therefore expect to observe the vertex of the conical outflow cospatial with signatures of the ejection mechanism, i.e., the peak in the dust continuum where star formation is assumed to be at a maximum.

As discussed in Section 4.2, the absence of an Einstein ring in the image plane $OH^+$ optical depth map indicates that the outflow does not extend over the inner caustic toward the east. Similarly, the $OH^+$ optical depth decreases toward the north and south (evident in both the image and source planes) away from the peak of the dust continuum intensity, indicating that we are observing the true edge of the $OH^+$ distribution to the east, north, and south. This cannot be said for the $OH^+$ distribution in the western direction where the $OH^+$ component likely extends farther out than the dust component but becomes invisible to us without the background continuum to absorb.

With these points in mind, the elongated triangular morphology revealed in the source plane $OH^+$ optical depth map can be interpreted as the 2D projection of a 3D conical structure viewed from outside the opening angle (i.e., not observed "down the barrel"). The vertex of this rough isosceles triangle sits cospatial with the peak of the dust continuum (Figure 6) and flares out and away from the dust continuum peak toward the west. Additionally, there is a mild negative gradient in LOS velocity observed along the east to west axis, and an increase in velocity dispersion. This may be indicative of an outflow driven over an extended period of time, as opposed to a single ejection event.

To illustrate this model we employ our simple toy model introduced in Section 5.2 to construct multiple 3D cone models over a range of inclinations, such that they all have the same 2D projected radius $R_{obs}$, opening angle $\Delta\phi_{obs}$, and position angle $\alpha_{obs}$ observed in the source plane $OH^+$ optical depth map (Figure 8). We note that the radius of the outflow measured directly from the source plane $OH^+$ optical depth map is larger than its true value due to the effects of beam smearing and is used in this toy model for illustrative purposes only.

We find that all conical models produce similar projected morphologies (Figure 8) and reproduce the general characteristics observed in the $OH^+$ optical depth map. The conical outflow geometry is therefore the most suitable geometry to describe the observed $OH^+$ outflow in G09v1.40 and is therefore adopted throughout the rest of our analysis and expanded upon in Section 7.

An obvious difference between the data and model projected outflow morphologies is the offset in peak optical depth, with the optical depth in the data peaking toward the "edge" of the cone. If real, this can be explained by an outflow extending past the background continuum, causing the fainter far edge of the cone to become invisible to the observer and thus bringing the observed edge closer to the peak in observed optical depth. We stress however that such detailed interpretation of the source plane structure should be taken with extreme care, given the effects of the beam. Higher-resolution observations are needed to investigate the structure within the outflow.

While we do not attempt to model the velocity structure of the conical outflow, we note that there is a slight negative radial

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Example conical outflow model compared to the observed $OH^+$ optical depth map. Top left: source plane $OH^+$ optical depth map with the measured 2D projected cone geometry overlaid in white. The projected radius is denoted by $R_{obs}$ and 2D projected opening angle by $\Delta\phi$. Top right: 2D projection of a model conical outflow with constant radial density. Bottom left: source plane $OH^+$ line-of-sight velocity map with measured 2D projected cone geometry overlaid in white. Outflow velocity increases with radius. Bottom right: source plane $OH^+$ velocity dispersion map with measured 2D projected cone geometry overlaid in white. Outflow velocity dispersion increases with radius.</figcaption></figure>

gradient in LOS velocity, indicating an increase in outflow speed at larger radii, which is coupled with an increase in velocity dispersion. Assuming the cone geometry is the correct choice for this outflow, the positive trend of outflow velocity and velocity dispersion with radius is indicative of an outflow driven over an extended period of time, as opposed to a single ejection event. This interpretation however should be taken simply as a suggestion given the velocity difference from one end of the outflow to the other is a mere $\sim 40 \text{ km s}^{-1}$. Additional observations and analyses are required to confirm this scenario.

## 6. Chemical Properties of the Gas

In addition to dynamical and morphological information, molecular spectral lines offer insight into the physical state of their media through an understanding of the chemical and physical requirements to form them. We first address the state of the ISM in the host galaxy as traced by the CO(9–8) emission line observed at systemic velocities in Section 6.1. We then discuss the formation processes and necessary environmental conditions required to produce the $OH^+$ ($1_0-1_1$) and $CH^+$($1-0$) transitions separately (Sections 6.2 and 6.3), before commenting on the physical state of the outflowing gas as traced by both of these transitions.

### 6.1. CO(9–8)

High-redshift DSFGs contain large molecular gas reservoirs (Carilli & Walter 2013), which fuel their rapid ongoing star formation. CO rotational lines can be used to constrain the kinetic temperature and gas density when observed over multiple-$J$ lines.
Nonlocal thermodynamic equilibrium (non-LTE) models of the CO spectral line energy distributions (SLEDs) in high-redshift heavily star-forming galaxies suggest there are likely two (or more) excitation components dominating the CO emission (e.g., Ivison et al. 2010; Danielson et al. 2011; Yang et al. 2017) in these galaxies. The low-excitation component corresponds to a cooler extended molecular gas reservoir, dominating the global CO SLED at low-J transitions. The high-excitation component needed to explain the mid-/high-J transitions, on the other hand, is produced by a warmer, denser, and more compact molecular gas reservoir believed to closely trace the ongoing star formation.

Fifteen $z \sim 2$ SMG CO SLEDs, including that of G09v1.40, were analyzed by Yang et al. (2017) using a large velocity gradient (LVG) statistical equilibrium method. Fitting the CO(2-1), (4-3), (6-5), and (7-6) transition lines, G09v1.40 required a two-component model, indicating that the emission in CO(9-8) traces the warm, dense, and more compact molecular gas dominating the high-excitation component likely associated with ongoing star formation.

High CO excitation is also observed in galaxies harboring powerful AGNs capable of dominating the IR luminosity. Significant boosting of the high CO($J > 10$) transition lines is believed to be caused by AGN heating within the central hundreds of pc around the nucleus (van der Werf et al. 2010) and may contribute significantly to the excitation of CO(9-8). From our data, we derive an intrinsic CO(9-8) luminosity in G09v1.40 of $L'_{CO(9-8)} = 3.1 \times 10^9$ K km s$^{-1}$ pc$^2$, which is consistent with the two-component CO SLED presented by Yang et al. (2017) and does not indicate obvious signs of an AGN boost. Furthermore, in a scenario where the AGN is significantly contributing to the thermal dust continuum and high-J CO transition lines, we would expect these two components to be cospatial. This is not the case for G09v1.40 where the CO(9-8) emission is both offset and more extended than that of the dust continuum (Figure 6). We therefore maintain our assumption that the CO(9-8) emission observed in G09v1.40 is excited predominantly via mechanisms associated with ongoing star formation.

## 6.2. $OH^+$

The observed $OH^+$ absorption line is clearly blueshifted with respect to the systemic velocity of G09v1.40 (Figures 2, 3, and 6) and therefore must trace outflowing gas in this system. It is important however to further constrain what phase of the outflowing gas is traced by $OH^+$ in order to accurately derive properties and interpret the multiphase outflow as a whole. We do this by addressing the chemistry required to produce $OH^+$ absorption.

In the cool diffuse neutral ISM, neutral-neutral reactions advance slowly, allowing ion-neutral reactions to dominate the chemistry when in the presence of an external ionization field. In these conditions, chemical species with a first-ionization potential less than that of neutral hydrogen (13.6 eV) will be predominantly ionized by the incident far-ultraviolet (UV) radiation, and species with a first-ionization potential > 13.6 eV, such as oxygen, O, are shielded by the abundant atomic hydrogen, H. Many of the reaction networks of the latter are therefore kicked off by $H^+$ and $H_3^+$, predominantly formed via cosmic-ray (CR) ionization.

Indeed, the dominant formation pathway of $OH^+$ in the cool diffuse ISM begins with the ionization of neutral H by a CR:

$$H + CR \rightarrow H^+ + e^- + CR'. \quad (2)$$

$O^+$ can then be formed via an endothermic charge transfer between O and $H^+$,

$$H^+ + O + \Delta E \leftrightarrow O^+ + H, \quad (3)$$

which proceeds backwards, uninhibited, in an exothermic reaction. The neutralization of $H^+$ via $e^-$ capture and charge transfer with polycyclic aromatic hydrocarbons (PAHs) also counteract the production of $O^+$.

The remaining $O^+$ can react with $H_2$,

$$O^+ + H_2 \rightarrow OH^+ + H, \quad (4)$$

to produce $OH^+$, which can then be destroyed by dissociative recombination, photodissociation or proceed further along the oxygen chemistry network through abstraction reactions with $H_2$,

$$OH^+ + H_2 \leftrightarrow H_2O^+ + H. \quad (5)$$

The rapidity with which this abstraction process proceeds creates a sensitive relationship between the $OH^+$ and $H_2$ abundances.

$OH^+$ can alternatively form via

$$O + H_3^+ \rightarrow OH^+ + H_2. \quad (6)$$

However, this reaction requires a significantly higher molecular hydrogen fraction in combination with low $e^-$ abundances (fuller discussions on oxygen chemistry in the ISM can be found in Hollenbach et al. 2012; Indriolo et al. 2015).

As the first oxygen-bearing ion to form after the ionization of H, $OH^+$ is a key ingredient in constraining the physical and chemical properties of the ISM, including the cosmic-ray density and molecular hydrogen fraction. Analyses of Milky Way sight lines have shown that $OH^+$ predominantly traces the cool diffuse gas in the ISM where hydrogen is primarily neutral (Gerin et al. 2016). It is therefore expected that the bulk of $OH^+$ present in the ISM forms via Equation (4), with formation via Equation (6) only dominating within the opaque and predominantly molecular interiors of molecular clouds (Hollenbach et al. 2012; Indriolo et al. 2018).

Bialy et al. (2019) further investigated the large scatter in $OH^+$-to-neutral hydrogen column density ratios, $N(OH^+)/N(H)$, measured in Milky Way sight lines, in the context of a turbulent medium. The abstraction of hydrogen in Equation (5) means $OH^+$ is highly sensitive to the abundance of $H_2$ which in turn is sensitive to density fluctuations in the underlying turbulent medium (Bialy et al. 2017). Using magnetohydrodynamic (MHD) simulations they modeled increasingly turbulent density fields and postprocessed them with chemical models to obtain probability density functions of the abundances. The model that best reproduced the observations required high levels of turbulence suggesting that turbulence is an important factor in the production of $OH^+$ in the cool diffuse ISM.

Observational studies analyzing $OH^+$ (and $H_2O^+$) absorption in the high-redshift galaxies SMM J2135-0102 and SDP 17b ($z \sim 2.3$, Indriolo et al. 2018) similarly conclude that the bulk of the $OH^+$ resides in cool diffuse gas, surrounding the galaxies in massive extended halos.
We therefore interpret the blueshifted $OH^+$ absorption measured in our observations of G09v1.40 as a tracer of the predominantly atomic gas-phase component in the outflow.

### 6.3. CH$^+$

To further constrain our analysis of the neutral outflowing gas as traced by $OH^+$, we include an analysis of the two ancillary observations of the $CH^+(1-0)$ transition in G09v1.40. Before comparing the two light hydride data sets directly, we first explore the chemistry required to form $CH^+$ and the physical conditions needed to produce its spectral lines.

$CH^+$ can form via the endothermic reaction between ionized C and molecular hydrogen,

$$C^+ + H_2 + h\nu \rightarrow CH^+ + H, \quad (7)$$

but requires temperatures $\gtrsim 10^3$ K.

In the diffuse ISM, such temperatures can be reached locally via the intermittent dissipation of turbulent energy (Godard et al. 2009). Due to the high critical density of the $J=1-0$ transition, $n_{crit} \sim 10^7$ cm$^{-3}$, most of the $CH^+$ in this diffuse environment will be in the ground state, causing high $J=1-0$ line opacities. $CH^+$ absorption in the MW has been observed in several hundreds of sight lines, with abundances of the molecule scaling positively with the turbulent energy transfer rate of the diffuse molecular gas, supporting the scenario where $CH^+$ is formed predominantly via the reaction shown in Equation (7) (Godard et al. 2014). Local extragalactic observations similarly find $CH^+$ absorption in turbulent environments, such as the medium surrounding SN 2014J, in M82 (Ritchey et al. 2015), and the ISM of the starburst galaxy Arp 220 (Rangwala et al. 2011).

In regions of dense gas ($n_H > 10^5$ cm$^{-3}$), sufficiently illuminated by ultraviolet radiation and undergoing suprathermal heating, $CH^+$ can be observed in emission (Godard & Cernicharo 2013). Star-forming regions such as DR21 are well modeled by a C-shock scenario where the approaching magnetic field causes sudden heating of the upstream neutral gas via ion-neutral friction, resulting in very wide emission signatures (Falgarone et al. 2010). $CH^+$ emission observed in the nearby ULIRG Mrk 231 (van der Werf et al. 2010) is also likely due to the very strong UV field and photon-dense regions present in this source.

Falgarone et al. (2017) presented results finding both $CH^+$ absorption and emission in G09v1.40, as part of a sample of six $z \sim 2$ DSFGs. They conclude that the broad ($>1000$ km s$^{-1}$) $CH^+$ emission arises in shocked gas associated with galactic winds driven by the central starbursts, while the narrower $CH^+$ absorption lines must trace gas outside the galaxy. $CH^+$ absorption traces the cool diffuse gas sitting in turbulent halos around these galaxies, which are mechanically fueled by the outflowing gas from the central starburst. In the sample of six DSFGs, there are in total four blueshifted $CH^+$ absorption lines (only three are reported by Falgarone et al. (2017) as the spectrum of G09v1.40 was wrongly displayed due to an incorrect redshift).

### 6.4. Comparison of $OH^+$ and $CH^+$

Given the similarities between $OH^+$ and $CH^+$ chemistry, we perform a simple analysis of the “low-” and “high-”spatial-resolution $CH^+(1-0)$ observations.

The low-resolution observations were previously studied by Falgarone et al. (2017) and included self-calibration in their reduction of the data. Falgarone et al. (2017) report overlapping

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Continuum-subtracted $CH^+(1-0)$ spectra of a single ALMA spaxel (pixel area = 0.0"05 × 0.0"05) in the low-spatial-resolution $CH^+$ data set, compared with the $OH^+(1_1-1_0)$ and CO(9-8) spectra previously presented in Figure 2, at the same physical location in G09v1.40. Top: $CH^+$ spectrum plotted as a function of observed frequency. Bottom: $CH^+$ spectrum (yellow solid line) plotted over the $OH^+$ (green) and CO(9-8) (blue) lines as a function of velocity with respect to the galaxy’s systemic velocity (Yang et al. 2017). The $CH^+$ spectrum is fitted with a single Gaussian function shown in shaded yellow ($OH^+$ and CO(9-8) spectral fits are as indicated as described in Figure 2). The $CH^+$ line is blueshifted $\sim 250$ km s$^{-1}$ with respect to the bulk molecular gas, tracing the same kinematic component as the $OH^+$ absorption.</figcaption></figure>

$CH^+$ absorption and emission lines at central velocities of $111 \pm 7$ km s$^{-1}$ and $28 \pm 34$ km s$^{-1}$ and velocity dispersions of $1124 \pm 87$ km s$^{-1}$ and $361 \pm 24$ km s$^{-1}$, respectively, using an incorrect redshift of 2.0894 (quoted in Bussmann et al. 2013). In our reanalysis of the data, we make use of the automatic ALMA pipeline products and adopt the precise redshift of $z = 2.0924 \pm 0.0001$ derived by Yang et al. (2017) from multi-J CO spectra, which agrees within errors with the redshift derived from our CO(9-8) of $z = 2.093$. With this correction, the $CH^+$ absorption has a central velocity of $\sim -250$ km s$^{-1}$, blueshifted with respect to the bulk molecular gas and closely follows the $OH^+$ absorption, as shown in Figure 9. We find $\sim 4 \times$ weaker $CH^+$ emission compared to that reported by Falgarone et al. (2017) when we stack the spectra over the source.

To investigate this further, we analyze the high-resolution $CH^+(1-0)$ observations. We image this data set using a robust weighting of 0 and select spectral channels greater than 1500 km s$^{-1}$ away from the line center for the continuum modeling and subtraction so as to avoid any contamination of the emission. We then stack all spaxels with a continuum S/N $> 3$ and again find only a weak signature of the wide $CH^+$ emission line, approximately $3 \times$ lower than that reported by Falgarone et al. (2017) (see Figure 10).
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Stacked CH⁺(1−0) spectra of all spaxels with continuum signal to noise > 3, in the high-spatial-resolution data set. The spectra are fitted with two Gaussian components, shown separately with the dashed yellow lines and combined with the shaded yellow region. The blueshifted absorption and emission lines reveal outflowing neutral and shocked gas with line of sight velocities of V<sub>abs</sub> = −221 km s<sup>−1</sup> and V<sub>emi</sub> = −409 km s<sup>−1</sup>, respectively. The strength of the CH⁺ emission is much weaker, and its central velocity is much more blueshifted than previously reported by Falgarone et al. (2017) using the low-spatial-resolution data set.</figcaption></figure>

To analyze the spatial distributions of the CH⁺(1−0) absorption and emission in G09v1.40 we return to the low-resolution data set as these observations provide higher S/N in the spectra of the CH⁺ line in each spaxel and a more comparable beam size to the OH⁺ observations.

The same spectral fitting routine introduced in Section 3 is applied to CH⁺(1−0) absorption line in each spaxel of the low-resolution data set. We do not attempt to fit simultaneously for the emission and absorption due to the low S/N of the emission line. From each spaxel, we subtract the best-fit Gaussian to the CH⁺ absorption and sum the residuals. We interpret a positive summed residual as excess CH⁺(1−0) emission and provide an approximate distribution of the CH⁺(1−0) in both absorption and emission. The weak CH⁺ emission appears to be compact and cospatial with the dust continuum while the absorption covers a more extended area, cospatial with that of the OH⁺ absorption in our data (Figure 11).

The kinematic and spatial coincidence of the CH⁺ and OH⁺ absorption lines indicate these molecules are tracing the same diffuse, turbulent, and predominantly atomic gas reservoir. The compact and central spatial distribution of the CH⁺ excess emission (to be interpreted cautiously) is consistent with the scenario of CH⁺ emission tracing shocked regions of dense gas in close proximity to a strong ultraviolet radiation source (i.e., the central starburst region).

## 7. Derived Outflow Properties

In this section we derive key outflow properties (total neutral gas mass, mass outflow rate (M<sub>OF</sub>), kinetic energy flux E, and momentum flux p: Table 4) to further investigate its ejection mechanism and to quantify the impact of the outflow on the evolution of its host galaxy.

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Top row: image plane maps of the CH⁺(1−0) absorption line optical depth and excess CH⁺(1−0) emission in the low-resolution CH⁺ data set. The beam (0′′59 × 0′′46) is displayed in the bottom-left corner of the left-hand panel. Bottom row: source plane reconstructions of the above panels. While the CH⁺ optical depth map follows the OH⁺ optical depth map morphology, the CH⁺ emission follows more closely the dust continuum intensity distribution. Pixel area (0′′1 × 0′′1) is the same in both the image and source planes.</figcaption></figure>

### 7.1. Outflow Mass

The OH⁺(1<sub>1</sub>−1<sub>0</sub>) column density N<sub>OH⁺</sub> and integrated optical depth are related through

$$ \int \tau_{\nu} dv = \frac{\lambda^3}{8\pi} AN_{OH^+}, \qquad (8) $$

where A = 2.11 × 10⁻² s⁻¹ is the Einstein coefficient and λ ≈ 0.029 cm is the wavelength of the transition (see Appendix B for the derivation). This gives an average N<sub>OH⁺</sub> value in the source plane of 2.4 × 10¹⁴ cm⁻².

To derive a total neutral H column density, we must now adopt an OH⁺ abundance, X<sub>OH⁺</sub> = N<sub>OH⁺</sub>/N<sub>H</sub>. The formation of OH⁺ is sensitive to the CR ionization rate of atomic hydrogen (see Equation (2)), allowing its abundance relative to that of the H₂O⁺ ion to be used as a constraint on the radiation field (Gerin et al. 2010; Neufeld et al. 2010; Porras et al. 2014; Indriolo et al. 2015; Zhao et al. 2015). Observational studies utilizing this method have uncovered a wide range in ionization rates from as low as 10⁻¹⁷ s⁻¹ in dense gas of the Milky Way, to rates of 10⁻¹³ s⁻¹ in the nuclear regions of ULIRGs (González-Alfonso et al. 2013, 2018). Naively we may expect observations such as the latter or those targeting the disks of nearby starburst galaxies (van der Tak 2016) to match the conditions of high-z DSFGs such as G09v1.40, as is indeed reflected in the ionization rate estimate of the star-forming region in SMM J2135−0102 (10⁻¹³−10⁻¹¹ s⁻¹, Danielson et al. 2013). Observations of OH⁺ and H₂O⁺ absorption in the same galaxy (and in that of SDP 17b), however, show ionization
<table><caption>**Table 4**

Derived Properties of the Neutral Outflow in G09v1.40 Using the Mean OH⁺ Abundance, log(XOH⁺), Modeled by Bialy et al. (2019) and the Highest and Lowest Observed Abundances in the Milky Way Sight Lines Used as Comparison in Their Analysis</caption><thead><tr><th>X<sub>OH+</sub></th><th>M<sub>neut</sub><br>(10<sup>9</sup>M<sub>☉</sub>)</th><th>ṁ<sub>OF</sub><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>Ė<sub>OF</sub><br>(10<sup>9</sup>L<sub>☉</sub>)</th><th>Ė<sub>OF</sub>/Ė<sub>SF</sub></th><th>ṁ<sub>OF</sub><br>(10<sup>37</sup> dyne)</th><th>ṁ<sub>OF</sub>/ṁ<sub>ej</sub></th><th>ṁ<sub>OF</sub>/(L<sub>IR</sub>/c)</th><th>τ<sub>OF</sub><br>(Myr)</th></tr></thead><tbody><tr><td>Model: Mean</td><td>6.7</td><td>83–25400</td><td>1.6–2100</td><td>0.013–16</td><td>0.040–27</td><td>0.13–91</td><td>0.067–47</td><td>300–0.98</td></tr><tr><td>Observed: High</td><td>16</td><td>200–6200</td><td>4.0–5200</td><td>0.031–40</td><td>0.096–66</td><td>0.32–220</td><td>0.16–110</td><td>250–0.83</td></tr><tr><td>Observed: Low</td><td>2.5</td><td>31–9550</td><td>0.62–790</td><td>0.0047–6.1</td><td>0.015–10</td><td>0.049–34</td><td>0.025–18</td><td>340–1.1</td></tr></tbody></table>

rates several orders of magnitude lower (Indriolo et al. 2018), reflecting ionization levels and locations in $N(\text{OH}^+)/N(\text{H})$ versus $N_{\text{OH}^+}/N_{\text{H}_2\text{O}^+}$ parameter space seen in diffuse clouds of the Milky Way (Indriolo et al. 2015; Neufeld & Wolfire 2017). This supports the scenario where OH⁺ absorption traces the same extended, turbulent halos of neutral gas surrounding high-z DSFGs traced by CH⁺ absorption (Falgarone et al. 2017), where the large distances from the central starburst region supplying the CR flux results in a dramatic decrease in ionization rate.

Following the evidence provided above, we therefore adopt an OH⁺ abundance of $\log_{10}(X_{\text{OH}^+}) = \log_{10}(N_{\text{OH}^+}/N_{\text{H}}) \approx -7.8 \pm_{0.075}^{0.05}$ from Bialy et al. (2019). This is the mean value derived in their best-fitting and most turbulent isothermal magnetohydrodynamic (MHD) simulation of the diffuse neutral medium which they compare with values observed in Milky Way sight lines. They stress that the dispersions derived in their best-fit models do not reflect the scatter measured in the Milky Way observations, and so we do not adopt the formal uncertainties of their models in our derivations. We instead consider the highest and lowest observed abundances ($\log(X_{\text{OH}^+}) \approx -7.4$ and $-8.2$) from the comparison sample of Milky Way sight lines as extreme cases. Outflow properties derived using the extreme observed abundances are listed in Table 4 and shown as dashed curves in Figures 12, 14, and 15. We further note that the average H column densities found in G09v1.40 are $\sim 2 \times$ higher than in the models of Bialy et al. (2019) and caution that all abundances derived above from both models and observations do not contain nonequilibrium chemistry and may not match other physical properties of outflowing gas and gas at high redshift.

With $X_{\text{OH}^+}$ in hand, the total neutral gas mass of the outflow, by summing over all the pixels $(i, j)$ in the source plane, is given by

$$M_{\text{neut}} = 1.36 m_{\text{H}1} \sum_{i,j} N_{\text{OH}^+} A_{\text{pix}} / X_{\text{OH}^+}, \quad (9)$$

where the 1.36 factor is the correction for the helium abundance, $m_{\text{H}1}$ is the mass of a hydrogen atom, and $A_{\text{pix}}$ is the area of a single pixel in cm².

Following this method we measure a total neutral gas mass of the outflow in G09v1.40 of $M_{\text{neut}} = 6.7 \times 10^9 M_{\odot}$. This is more than 25% as massive as the molecular gas mass (corrected for the helium abundance) in the host galaxy, $2.5 \times 10^{10} M_{\odot}$, derived via non-LTE radiative transfer modeling of multi-J CO lines (Yang et al. 2017). While the uncertainty associated with the OH⁺ abundance likely dominates, we also note that this measurement excludes any part of the outflow that extends past, or lies behind the dust continuum but cannot be observed in absorption.

As mentioned earlier, the diffuse turbulent and predominantly atomic gas component traced by OH⁺ absorption is

believed to be the same component seen in CH⁺ absorption (Falgarone et al. 2017). Using a turbulent framework to analyze the global CH⁺ absorption spectra, Falgarone et al. (2017) derived a radius of the full turbulent reservoir around G09v1.40 of $r_{\text{TR}} = 12$ kpc. They then extrapolate the column densities observed over the dust continuum (0.41 kpc) to both sides of the galaxy and out to 12 kpc, finding a total neutral gas mass of the full turbulent reservoir of $1.1 \times 10^{10} M_{\odot}$. Thus to compare with the mass derived from OH⁺, we scale the mass full reservoir mass by a factor of $0.41 \text{ kpc}/(2 \times 12 \text{ kpc})$, giving a value of $0.19 \times 10^9 M_{\odot}$, 3% that of our derived mass. We note however that the turbulent framework used by Falgarone et al. (2017) to convert the observed CH⁺ absorption into a total neutral gas mass is a very different approach from our own.

For a more direct comparison with our work, we additionally convert the CH⁺ optical depth map (Figure 11) into a total neutral outflow gas mass in the same way that we have done for the OH⁺. We first convert the source plane CH⁺ optical depth map into CH⁺ column density following the equation presented by Falgarone et al. (2017) in their supplementary methods. We then convert this into a neutral H column density using an average observed CH⁺ abundance of $7.6 \times 10^{-9}$ (Godard et al. 2014). Summing over the entire outflow, we find a total neutral outflow gas mass of $4.2 \times 10^9 M_{\odot}$, comparable with that derived from OH⁺.

## 7.2. Mass Outflow Rate

Methods of estimating the mass outflow rate, $\dot{M}_{\text{OF}}$, vary among the literature and depend on the assumed geometry of the outflow (Veilleux et al. 2005, 2020). As concluded in Section 5, the most likely basic geometry of the neutral outflow observed in G09v1.40 is conical. The $\dot{M}_{\text{OF}}$ is therefore given by

$$\dot{M}_{\text{OF}} = 3 \frac{MV}{R}, \quad (10)$$

where M is the total mass, V the velocity, and R the maximum radius of the cone. The factor of 3 accounts for two assumptions: (1) the cone is filled and not just a thin shell; (2) the density of the cone is constant with radius. These assumptions are based on the observed morphology of the OH⁺ optical depth presented in Figure 8, which is best reproduced by a toy model implementing these assumptions.

As previously discussed, we only observe and measure the 2D projected values of these parameters and require the inclination (inc) of the outflow, with respect to the observer’s LOS, to correct for this effect. As our data is not of high enough quality to measure the inclination of the outflow with confidence, we instead opt to place sensible bounds on the inclination and derive a range of possible deprojected geometrical parameters and $\dot{M}_{\text{OFs}}$.

We begin by measuring the 2D projected opening angle, $\psi_{\text{obs}} \approx 60^{\circ}$ rad, of the outflow using the OH⁺ optical depth map
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Deprojected outflow parameters as a function of inclination where inc = 90° corresponds to an outflow central axis perpendicular to the line of sight. Top: radius (solid line) and 2D opening angle (dashed line). Middle: velocity. Bottom: $\dot{M}_{\text{OF}}$ (solid line) with the $1\sigma$ spread of the modeled $OH^+$ abundance from Bialy et al. (2019) shown by the shaded region. The upper and lower dashed lines indicate the $\dot{M}_{\text{OF}}$ if the highest and lowest observed $OH^+$ abundances from Milky Way sight lines are assumed, respectively (Bialy et al. 2019). The right-hand axis displays the $\dot{M}_{\text{OF}}$ normalized by the SFR indicating that if the inclination of the outflow is >3°6, then the neutral gas $\dot{M}_{\text{OF}}$ exceeds the SFR in G09v1.40. Hatched regions indicate regions outside of the considered inclination range.</figcaption></figure>

(Figure 8). If the outflow is flowing perpendicular to the LOS, inc = 90°, the deprojected 2D opening angle, $\psi$, is indeed $\psi = \psi_{\text{obs}}$, however if inc < 90° then $\psi < \psi_{\text{obs}}$. If we then assume that all the outflowing gas is situated in front of the galaxy (i.e., no part of the cone may have an inclination larger than 90°), then the maximum possible inclination of the cone's central axis is given by $\psi/2 = 90° - \text{inc} \approx 64°$.

We cannot determine a lower limit to the inclination based only on the measured projected parameters. Instead, we assume that the outflow does not extend farther than the halo virial radius $r_h$, thus putting a limit on the radius of the outflow which increases rapidly at small inclinations (Figure 12). Taking the stellar radius–halo radius value (SRHR $\approx$ 0.018, defined as the ratio of galaxy radius to halo virial radius) measured by Somerville et al. (2018) and the effective half-light radius of the reconstructed NIR distribution in G09v1.40 (Table 3), we estimate a halo radius of $r_h = 0.018 r_{\text{eff,NIR}} = 61$ kpc. The minimum inclination of the outflow possible is then given by $\text{inc}_{\text{min}} = \ar sin(R_{\text{obs}}/R_h) \approx 0.38$.

Thus, we take an inclination range of inc = 0.4–64° that we believe the observed conical outflow may have. This corresponds to a range of possible outflow radii of $R = R_{\text{obs}}/\sin(\text{inc}) = 0.45–61$ kpc, 2D opening angles of $\psi = \arctan(\tan(\psi_{\text{min}})\sin(\text{inc})) = 0.2–24°$, where $\psi_{\text{min}} = 2(90°-\text{inc}_{\text{max}})$, and outflow velocities of $V = V_{\text{obs}}/\cos(\text{inc}) = 250–570$ km s$^{-1}$, where $V_{\text{obs}}$ is the maximum LOS velocity in the source plane $OH^+$ velocity map. The neutral gas $\dot{M}_{\text{OF}}$ in G09v1.40 may then have a value within the range of $\dot{M}_{\text{OF}} = 83–25400 M_{\odot}$ yr$^{-1}$, which exceeds the SFR = 788 $\pm$ 300 $M_{\odot}$ yr$^{-1}$ if the inclination is above 3°6 (see Figure 12). This corresponds to mass-loading factors of the neutral gas outflow between 0.11 and 32.

In Figure 13 we compare the neutral gas $\dot{M}_{\text{OF}}$ of G09v1.40 with local neutral and molecular outflow samples from both active and purely star-forming galaxies (Cicone et al. 2014; Cazzoli et al. 2016; Fluetsch et al. 2019), along with the small sample of measured molecular and neutral outflows at high redshift (Feruglio et al. 2017; Brusa et al. 2018; Herrera-Camus et al. 2019, 2021; Jones et al. 2019; Spilker et al. 2020a) as a function of SFR. For similar SFRs, it is evident in the low-redshift samples that molecular outflows and outflows driven by AGNs have higher $\dot{M}_{\text{OF}}$s than neutral outflows and outflows driven by star formation. For G09v1.40, the range in deprojected $\dot{M}_{\text{OF}}$ spans more than two orders of magnitude, comparable to the scatter seen in the full comparison sample.

Because it is difficult to compare the range of $\dot{M}_{\text{OF}}$ derived for G09v1.40 with those in the literature, where deprojection of outflow parameters is either not attempted or derived for a single assumed inclination, we consider three additional derivations of $\dot{M}_{\text{OF}}$: (1) using only projected values, (2) assuming an axially symmetric disk and perpendicular outflow, and (3) assuming an alternative maximum outflow radius equal to the radius of the diffuse turbulent halo surrounding G09v1.40.

To derive $\dot{M}_{\text{OF}}$ using the projected outflow parameters we take the LOS outflow velocity, $V_{\text{obs}}$, and the effective radius of the dust continuum, $r_{\text{eff,cont.model}}$, provided by the *visilens* model. This is a comparable method to that used by, e.g., Spilker et al. (2020a) and gives a $\dot{M}_{\text{OF}}$ of 12500 $M_{\odot}$ yr$^{-1}$,
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Mass outflow rate as a function of star formation rate. The shaded light orange region indicates the parameter space of the neutral gas outflow in G09v1.40 over the possible inclination range discussed in Section 7. We additionally indicate the $\dot{M}_{\text{OF}}$ estimated for an axially symmetric galaxy with a perpendicular outflow (upper dashed line), using only projected outflow parameters (middle dashed line), and for an outflow with a radius equal to that of the turbulent reservoir surrounding G09v1.40 (Falgarone et al. 2017) (lower dashed line). For comparison we plot molecular and neutral gas outflows measured in local galaxy samples with gray and black symbols, respectively (Cicone et al. 2014; Cazzoli et al. 2016; Fluetsch et al. 2019). Open and filled symbols indicate sources with and without AGN activity, respectively. With orange symbols we show high-redshift molecular outflows measured in the following sources, listing their galaxy type, outflow tracer and redshift in brackets: XID2028 (QSO, CO, $z=1.593$; Brusa et al. 2018), zC400528 (AGN, CO, $z=2.387$; Herrera-Camus et al. 2019), APM08279+5255 (QSO, CO, $z=3.912$; Feruglio et al. 2017), SPT 0346-52 (DSFG, H$_2$O, $z=5.656$; Jones et al. 2019), and various SPT sources (DSFG, OH, $z>4$; Spilker et al. 2020a). Also in orange, we include the neutral gas outflow in HZ4 (main-sequence star-forming galaxy, [C II], $z \sim 5$ Herrera-Camus et al. 2021).</figcaption></figure>

which we mark in Figure 13, placing it at the most extreme end of observed outflows at all SFRs.

Assuming an axially symmetric disk model, we estimate an inclination of 49° for G09v1.40 using $\sin(\text{inc}) = (1 - (b/a)^2)/(1 - q_0)$, where the axial ratio, $(b/a) = 0.679$, is provided by the best-fit *visilens* model, and a typical intrinsic thickness of $q_0 = 0.2$ is assumed (Förster Schreiber & Wuyts 2020). This is the method used by Herrera-Camus et al. (2021) for the neutral outflow in HZ4 and corresponds to a deprojected velocity of 380 km s$^{-1}$, an outflow radius of 540 pc, and an $\dot{M}_{\text{OF}}$ of 14100 $M_{\odot}$ yr$^{-1}$, comparable to that derived for the projected case.

As an alternative maximum outflow radius, we consider the radius of the turbulent halo of diffuse neutral gas seen in CH$^+$ (Falgarone et al. 2017), believed to be the same reservoir containing OH$^+$ (Indriolo et al. 2018). Falgarone et al. (2017) analyzed the CH$^+$ halo surrounding G09v1.40 using a turbulent framework, estimating a radius of $r_{\text{TR}} = 12$ kpc. This corresponds to lower limits on the outflow velocity and

$\dot{M}_{\text{OF}}$ of 250 km s$^{-1}$ and 420 $M_{\odot}$ yr$^{-1}$, respectively, comparable with other high-$z$ molecular outflows.

Despite the large range in $\dot{M}_{\text{OF}}$ that can be derived for G09v1.40, it appears that the neutral outflow is at least comparable to, if not considerably more extreme than, that of molecular outflows observed at high redshift. In a study of eight nearby AGN and star-formation-driven outflows observed across their ionized, atomic, and molecular phases, the molecular gas was found, on average, to dominate the total mass and $\dot{M}_{\text{OF}}$ of their outflows (Fluetsch et al. 2021). The neutral gas dominates in only two purely star-forming galaxies. The authors suggest that the more powerfully driven AGN outflows are compacted more by the ambient CGM, leading to the observed higher gas densities and thus higher molecular gas fractions in these outflows. This, however, is certainly not the rule and comparable or greater neutral gas fractions have also been found in sources hosting an AGN, such as in the low-redshift Seyfert systems Mrk 231, Mrk 273, and IRAS F08572+3915 (see compilation by Herrera-Camus et al. 2020). It is
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Deprojected outflow kinetic energy flux (kinetic luminosity) of the outflow $\dot{E}_{OF}$ as a function of inclination (left axis) and normalized by the kinetic energy flux injected by SNe, $\dot{E}_{SF}$ (a.k.a. energy loading, right axis). The bottom gray region indicates the energy transferred to the ISM assuming an energy coupling up to 40%, which can occur if SNe are clustered. Higher $\dot{E}_{OF}$ require either unusually high-energy coupling or an AGN contribution. The upper and lower dashed lines indicate the $\dot{E}_{OF}$ if the highest and lowest observed $OH^+$ abundances from Milky Way sight lines are assumed, respectively (Bialy et al. 2019). An outflow flowing perpendicular to the LOS has an inc = 90°. Hatched regions indicate regions outside of the considered inclination range.</figcaption></figure>

therefore perhaps not surprising that the $\dot{M}_{OF}$ of the neutral outflowing gas component found in G09v1.40 is comparable to or greater than the average molecular outflows observed in similar galaxies.

In the case of a highly inclined outflow in G09v1.40, the extreme $\dot{M}_{OF}$ derived can be reconciled theoretically with a highly obscured QSO scenario (Costa et al. 2018). When multiscattering radiation pressure from IR radiation on dust grains is taken into account, outflows of predominantly cool gas, and the peak $\dot{M}_{OF}$ on the order of $10^3-10^{-4} M_{\odot} \text{yr}^{-1}$, can be produced. This phase is short lived (<10 Myr) and requires the QSO to be heavily enshrouded in dense gas, a possible scenario for a compact DSFG like G09v1.40. We venture further into the required driving mechanisms in the following subsections.

### 7.3. Outflow Energetics

In this section, we consider the possible sources and mechanisms required to drive energy and momentum flux of the observed neutral outflow in G09v1.40. We again approach this problem considering the range of possible inclinations (see Section 7), as both the momentum and energy flux of the outflow depend strongly on the inclination assumed and therefore presents multiple possible feedback and outflow scenarios.

We begin by deriving the kinetic energy flux (a.k.a. kinetic luminosity/ kinetic power) of the outflowing neutral gas $\dot{E}_{OF}$:

$$ \dot{E}_{OF} = \frac{1}{2} \dot{M}_{OF} (V^2 + 3\sigma^2), \qquad (11) $$

using the minimum value in the source plane velocity dispersion map, $\sigma = 77 \text{ km s}^{-1}$. We find possible values ranging between $\dot{E}_{OF} = 1.6 - 2100 \times 10^9 L_{\odot}$ which we compare to the

expected kinetic energy flux injected by SNe, $\dot{E}_{SF} = 7 \times 10^{41} (\text{SFR}/M_{\odot} \text{yr}^{-1}) \text{ erg s}^{-1} = 130 \times 10^9 L_{\odot}$ (Veilleux et al. 2005) (Figure 14) providing ratios between $\dot{E}_{OF}/\dot{E}_{SF} = 0.012 - 16$. The fraction of $\dot{E}_{SF}$ that is ultimately coupled to the ISM, and therefore used in driving the outflow, depends strongly on the clustering of SNe in the galaxy and gas-phase metallicity and structure of the ISM. If clustering is strong, SN-driven superbubbles can retain as much as 40% of the input energy (e.g., Sharma et al. 2014; Fielding et al. 2018), which would provide enough energy to drive the neutral outflow in G09v1.40 if the inclination is $\le 12^\circ$. If inc > 12°, then the outflow requires either an extremely high and unusual energy coupling efficiency or an additional source of energy flux, i.e., a past, low-luminosity, or obscured AGN, to be driven. If inc > 24°, a coupling efficiency of > 100% is required and an AGN contribution must certainly be playing a role to preserve energy conservation.

Kinetic coupling efficiencies of AGN-driven outflows can be up to ~0.1 (see Figure 2 in a compilation by Harrison et al. 2018), which would imply an AGN with the luminosity of at least ~8 $\times 10^{46} \text{ erg s}^{-1}$ in G09v1.40 for the most extreme scenario. This is well within the range of observed AGN luminosities and does not exclude any of the high-inclination outflow scenarios.

Next, we consider whether the outflow is consistent with a momentum or energy-driven scenario. Momentum-driven outflows occur when the thermal energy of the shocked SN ejecta is efficiently radiated away. The observed momentum flux of the outflow must then be supplied by the momentum flux deposited directly by SNe ejecta or by radiation pressure from young stars on dust grains in the outflow.

First, we derive the momentum flux of the outflowing neutral gas $\dot{P}_{OF}$:

$$ \dot{P}_{OF} = \dot{M}_{OF} V, \qquad (12) $$

finding a possible range of $\dot{P}_{OF} = 0.040 - 27 \times 10^{37} \text{ dyne}$, which we display on left axes in Figure 15. We then estimate the momentum flux expected to be deposited by SN ejecta as the product of the SN rate (~15 SNe yr$^{-1}$) and the momentum associated with the ejecta of a single SN. For an ejecta mass of $10 M_{\odot}$ and launch velocity $v = 3000 \text{ km s}^{-1}$ (see, e.g., Section 2.2 in Murray et al. 2005), the total injected momentum flux is $\dot{P}_{ej} = 3.0 \times 10^{36} \text{ dyne}$, giving $\dot{P}_{OF}/\dot{P}_{ej} = 0.13 - 91$. This is sufficient to drive the neutral outflow through a momentum-driven phase if inc $\le 2^\circ 9$.

In addition to the momentum flux deposited by SN ejecta, momentum flux supplied by UV radiation pressure from young stars onto dust grains in the outflow may also be contributing. The effectiveness of this mechanism depends on the optical depth of the outflow:

$$ \dot{P}_{\text{rad}} = \tilde{\tau} L_{\text{bol}}/c, \qquad (13) $$

where $\tilde{\tau} = (1 - e^{-\tau_{\text{single}}})(1 + \tau_{\text{eff,IR}})$, which includes both single- and multiple-scattering events (Hopkins et al. 2014, 2020). $\tilde{\tau}$ therefore ranges from $\tau_{\text{single}} = \tau_{\text{UV/optical}} \ll 1$ when the outflow is optically thin, to $\sim (1 + \tau_{\text{eff,IR}})$ when optically thick, where $\tau_{\text{UV/optical}}$ and $\tau_{\text{eff,IR}}$ are the optical depths in the UV/optical and IR regimes, respectively (Murray et al. 2005). If we make the conservative assumption that $L_{\text{bol}} \approx L_{\text{IR}}$ (i.e., all
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Deprojected momentum flux of the outflow $\dot{p}_{\text{OF}}$ (left axes) as a function of inclination. Top: comparing $\dot{p}_{\text{OF}}$ to the momentum flux injected by SN ejecta $\dot{p}_{\text{ej}}$ (right axis), where the gray region indicates the momentum flux available from the estimated SNe rate $\sim 15$ SNe yr$^{-1}$. Bottom: comparing $\dot{p}_{\text{OF}}$ to the possible momentum flux contributed by radiation pressure from young stars on outflow dust grains. Values on the right axes indicate the required effective IR optical depth for a given $\dot{p}_{\text{OF}}$, where the optically thin regime is indicated by the gray region. The $1\sigma$ spread in the $OH^+$ abundance modeled by Bialy et al. (2019) is shown by the dark gray shading, and the upper and lower dashed lines indicate the $\dot{E}_{\text{OF}}$ if the highest and lowest observed $OH^+$ abundances from Milky Way sight lines are assumed, respectively (Bialy et al. 2019). An outflow flowing perpendicular to the LOS has an inc $= 90^\circ$. Hatched regions indicate regions outside of the considered inclination range.</figcaption></figure>

the UV stellar radiation is absorbed and reradiated in the IR), then $\dot{p}_{\text{OF}}/(L_{\text{IR}}/c) = 0.067 - 47$ and radiation pressure could deposit a momentum flux of the order $3.9 \times 10^{36}$ dyne for an optically thick outflow ($\tau \approx 1$). In addition to the momentum flux deposited by SN ejecta, this is sufficient to drive the neutral gas as a momentum-driven outflow if inc $\le 8^\circ 4$. In the case of an optically thin outflow, single-scattering radiation could provide a maximum momentum of $(1 - e^{-\tau_{\text{single}}})L_{\text{IR}}/c$.

If, however, the thermal energy of the shocked SN ejecta is not efficiently radiated away, it may be used up in doing work against the ambient medium, driving an energy-driven outflow. This provides a boost in the momentum flux of the outflow in addition to the momentum flux deposited by the ejecta and thus drives a stronger outflow. If radiation pressure is negligible and

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Possible outflow-driving mechanisms as a function of outflow inclination. If the inclination is $< 2^\circ 9$ the momentum provided by SN ejecta is sufficient to drive the outflow through a momentum-driven phase (blue bar). If inc $> 2^\circ 9$, a momentum-driven phase is still possible if radiation pressure from young stars on dust grains in the outflow is nonnegligible. This requires $\tau_{\text{eff,IR}} > 0$ and $\tau_{\text{eff,IR}} > 1$ if inc $> 8^\circ 4$. In the absence of radiation pressure, the outflow requires a momentum boost if inclinations are $> 2^\circ 9$, which can be provided by an energy-driven phase of the shocked SN ejecta if thermal energy is not immediately and completely radiated away. The energy coupling required to drive the outflow via SN feedback is indicated below the yellow bar. If the inclination is $> 12^\circ$ the outflow requires either an unusually high-energy coupling of $> 40\%$, which is higher than that expected from a clustered SN scenario, or a contribution from a low-luminosity, obscured, or fossil AGN (red bar). In reality, all driving mechanisms may be contributing simultaneously and if the full multiphase and double-sided outflow were to be taken into account the inclination ranges shown here would be shifted to smaller inclinations.</figcaption></figure>

the inclination of the outflow is $> 2^\circ 9$, it is possible that we are observing an energy-driven outflow. The effectiveness of this mechanism depends on the coupling efficiency of energy to the ISM (see Figure 16), which is unlikely to exceed 40% in the case of highly clustered SNe. Thus, if the inclination is $> 12^\circ$, an unusually high-energy coupling efficiency or past, low-luminosity, or obscured AGN activity is needed to explain both the energy and momentum flux of the outflow.

We summarize the possible driving scenarios of the neutral gas outflow in G09v1.40 in Figure 16 but again note that if the full multiphase and double-sided outflow is taken into account, the ranges presented here would be shifted to smaller inclinations. We also note that our results are sensitive to the $OH^+$ abundance assumed and provide alternative curves (dashed lines in Figures 14/15/17) using the extreme observational $OH^+$ abundance taken from Bialy et al. (2019), which would again significantly shift the inclination ranges of the driving mechanism scenarios summarized in Figure 16.

### 7.4. Impact on the Host Galaxy and Fate of the Outflowing Neutral Gas

Cool gas outflows remove the direct fuel for star formation and therefore must have an impact on the future growth and activity of the host galaxy. Disregarding the possibility of gas accretion or a change in SFR and/or $\dot{M}_{\text{OF}}$ in G09v1.40, we estimate a depletion time of the host galaxy's gas reservoir due to the observed neutral gas mass outflow rate using $\tau_{\text{OF}} = M_{\text{gas}}/\dot{M}_{\text{OF}}$. We derive this timescale over the range of possible inclinations, finding a $\tau_{\text{OF}}$ between 300 Myr for the lowest possible inclination (i.e., the star-formation-driven end of the spectrum), and 0.98 Myr at the highest possible inclination (i.e., the AGN-driven end of the spectrum). If depletion due to star formation is also taken into account, $\tau_{\text{OF}+\text{SFR}}$, these timescales reduce to 29 and 0.95 Myr, respectively (Figure 17), where the depletion time due to star
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Depletion time of the host galaxy’s gas reservoir given the observed neutral gas mass outflow rate (black lines), and in combination with the SFR (gray lines), as a function of possible outflow inclination. The upper and lower dashed lines (black and gray) indicate the depletion time if the highest and lowest observed OH⁺ abundances from Milky Way sight lines are assumed, respectively (Bialy et al. 2019). The outflow-driving mechanisms from Figure 16 are shown at the top. An outflow flowing perpendicular to the LOC has an inc = 90°. Hatched regions indicate regions outside of the considered inclination range.</figcaption></figure>

formation alone is $\tau_{SFR} = 32$ Myr. Thus, if the inclination is low, star formation likely plays a major role in the depletion of the host galaxy gas reservoir.

Typical depletion times derived from SFRs in compact star-forming galaxies and quasars at $z \sim 2-4$ are on the order of $\sim 50$ and $\sim 100$ Myr, respectively (Spilker et al. 2016; Stacey et al. 2021), consistent with the timescales derived for the low-inclination scenarios in G09v1.40. The $\sim 1$ Myr depletion times derived for high inclinations are instead consistent with timescales predicted for DSFGs to transition into unobscured gas-poor QSOs, via a far-infrared bright QSO phase (Simpson et al. 2012; Costa et al. 2018). Thus, if a high-inclination scenario in G09v1.40 is assumed, this would suggest the galaxy is currently in an evolutionary stage just prior to or at the beginning of a highly obscured QSO phase.

Ejected gas may, however, be reaccreted back onto the galaxy at a later time, replenishing the galaxy’s gas reservoir and prolonging $\tau_{OF}$. We therefore investigate the likelihood of the neutral gas outflow in G09v1.40 escaping the galaxy’s potential well by considering the mass required, $M_{req}$, to gravitationally bind an outflow with velocity $V$ and radius $R$,

$$M_{\text{req}} = \frac{V^2 R}{2G}, \qquad (14)$$

where $G$ is the gravitational constant.

For an outflow with a Gaussian velocity distribution, this equation will provide the $M_{req}$ capable of containing half the outflowing material if the central velocity is used. We therefore derive required binding masses, over the range of possible inclinations, using deprojected velocities at the 50th, 60th, 70th, 80th, and 90th percentiles (again assuming a velocity

<figure><img src="image_19.png" /><figcaption>**Figure 18.** The mass required to gravitationally bind the outflowing neutral gas, as a function of possible inclination. The solid black curve is derived using the central velocity of the blueshifted OH⁺ line (i.e., indicating the mass required to bind 50% of the outflowing neutral gas at a given inclination). The dashed curves are derived using an outflow velocity at the 60th, 70th, 80th, and 90th percentiles (i.e., indicating the mass required to bind 60%, 70%, 80%, and 90% of the outflowing neutral gas), using a velocity dispersion of 77 km s⁻¹. The dotted horizontal lines indicate the gas mass $M_{gas}$ (Yang et al. 2017), stellar mass $M^*$ (Ma et al. 2015), and total mass $M_{gas} + M^*$ of the host galaxy. The outflow-driving mechanisms from Figure 16 are shown at the top. An outflow flowing perpendicular to the LOS has an inc = 90°. Hatched regions indicate regions outside of the considered inclination range.</figcaption></figure>

dispersion of 77 km s⁻¹), corresponding to an $M_{req}$ capable of containing 50%, 60%, 70%, 80%, and 90% of the outflow (Figure 18). The results do not significantly change if the maximum velocity dispersion $\sigma_{v,max} = 130$ km s⁻¹ of the outflowing neutral gas is taken instead (indicated by the arrow in Figure 18).

$M_{req}$ is largest ($8.6 \times 10^{11} M_{\odot}$) if the outflow inclination is small, due to the very large deprojected radii in this regime. If the inclination is $> 14°$, the gas mass of the galaxy alone is capable of containing 90% or more of the neutral outflow. If we include the stellar mass $M^* = 0.8 \pm 0.1 \times 10^{11} M_{\odot}$ (Ma et al. 2015), G09v1.40 is capable of retaining this fraction of the outflow down to an inclination of 3°, implying that for the majority of possible outflow scenarios, most or all of the outflowing neutral gas will remain bound to the galaxy. This material is then available to be reaccreted by the galaxy at a later time unless additional feedback processes, such as thermal feedback from an unobscured QSO phase (Costa et al. 2018), cause the gas to remain in the circumgalactic medium. For an inclination larger than 11°, where an AGN is required to drive the outflow and heating of the circumgalactic medium is likely, the galaxy may still be expected to quench on timescales of $\tau_{OF} \sim 1$ Myr.

## 8. Conclusions

We have presented resolved ($0''52 \times 0''41$) ALMA Band 6 observations of a massive $z = 2.09$ neutral outflow in the gravitationally lensed DSFG G09v1.40 (HATLAS J085358.9 +015537). We detect the outflow in absorption with the
1033 GHz OH⁺(1₁–1₀) transition, exploiting its close proximity to the CO(9–8) transition to observe both lines and the underlying 1034 GHz dust continuum, simultaneously with a single ALMA tuning. We obtain a spatially and spectrally resolved view of the cool neutral gas in the outflow as traced by the OH⁺, blueshifted with respect to the warm dense gas at systemic velocities as traced by CO(9–8). We perform spectral fitting on all spaxels to obtain clean intensity and velocity maps of the outflowing OH⁺ absorption and systemic CO(9–8) emission. In addition, we use ancillary data from ALMA tracing the CH⁺(1–0) absorption and underlying continuum at both low and high angular resolution, and from the Keck K-band at 2.2 µm tracing the stellar emission of the background galaxy.

The CO(9–8) displays a strikingly different image plane morphology to that of the dust continuum, following more closely that of the stellar distribution (Calanog et al. 2014). The image plane optical depth distribution of the OH⁺ absorption follows the continuum by construction but displays a comparatively more elongated morphology, falling off dramatically to the north and south, and does not display an Einstein ring, as seen in the dust, CO(9–8) emission and stars.

We obtain a lens model, exploiting the high-resolution (0.′′17 × 0.′′13) dust continuum observations, with *visilens* and reconstruct all our 2D data maps into the source plane using the pixelated reconstruction code *LENSTOOL*. The dust continuum reveals itself as a compact ellipse with the CO(9–8) and stellar components offset to the east and displaying more extended distributions. The blueshifted OH⁺ forms an extended triangular morphology flaring out toward the west. There is a dip in the stellar emission at the position of the peak dust continuum which we believe is most likely due to extreme extinction in this region.

Three simple outflow geometries (a sheet, spherical bubble, and cone) are considered and compared with the observed and reconstructed OH⁺ morphology and kinematics. We find that a conical outflow geometry, where outflowing gas is ejected from the central dusty star-forming region toward the west, is the most suitable choice.

The physical conditions necessary for forming OH⁺ in the ISM suggests that OH⁺ absorption traces the diffuse and predominantly atomic gas component of a turbulent outflow. Comparing the absorption of OH⁺ with that of the CH⁺, which similarly probes diffuse atomic gas (Falgarone et al. 2017), we find that both absorptions lines are spatially in G09v1.40, tracing the same blueshifted kinematic component, confirming these lines trace the same gas phase. We therefore adopt an OH⁺ to H I abundance from Bialy et al. (2019) of log₁₀(n_HO⁺/n_H) ≈ -7.8±0.05, finding a total atomic gas mass of the outflow of M Neut = 6.7 × 10⁹ M⊙, which is more than 25% as massive as the molecular gas component in the host galaxy (Yang et al. 2017).

We consider a range of possible 2D projections of the conical outflow, deriving possible inclinations of the central axis with respect to the observer's line of sight between inc = 0.°4–64.°. Over this inclination range, we derive possible deprojected outflow radii between R = 0.45 and 61 kpc, 2D opening angle between inc = 0.°2–24.°, and velocity between V = 250–570 km s⁻¹.

Physical properties of the conical outflow are also derived as functions of possible inclination. The total neutral gas $\dot{M}_{\text{OF}}$ is between 83 and 25,400 $M_{\odot}$ yr⁻¹, which exceeds the SFR of

788 ± 300 $M_{\odot}$ yr⁻¹ if the inclination is greater than 3.°6. We find ranges in the kinetic and momentum fluxes of $\dot{E}_{\text{OF}} = 1.6 – 2100 \times 10^9 L_{\odot}$ and $\dot{p}_{\text{OF}} = 0.040 – 27 \times 10^{36}$ dyne, respectively.

We compare these values to the kinetic energy ($\dot{E}_{\text{SF}} = 130 \times 10^9 L_{\odot}$) and momentum flux injected by SNe ($\dot{p}_{\text{ej}} = 3.0 \times 10^{36}$ dyne) and radiation from young stars to determine the likely driving mechanism of the outflow, finding that this depends strongly on the inclination assumed. If the inclination is ≤2.9, the outflow may be momentum-driven by SN ejecta. If inc > 2.9, the outflow may still be momentum driven, provided radiation pressure from young stars onto dust grains in the outflow is taken into account. For inc > 8.4, this requires the outflow to be optically thick.

In the case where thermal energy deposited by SNe into the ISM is not efficiently radiated away, it may be used to do work on the ambient medium, providing a momentum boost for an energy-driven outflow. If radiation pressure is negligible, then the outflow may be energy driven if the inclination is inc > 2.9 up to a maximum inclination of 12°, where a coupling efficiency of the thermal energy to the ISM of 40% is needed. If the inc > 12°, either an extremely high coupling efficiency or an additional driving mechanism is needed, e.g., an AGN.

Depletion times of the host galaxy gas reservoir, due to the SFR and neutral outflow range from 29 Myr in the regime of a stellar-driven outflow, down to 0.95 Myr at the extreme end of the AGN-driven regime. This is consistent with timescales derived for other intensely star-forming galaxies at the same redshift (Stacey et al. 2021) and with timescales predicted for DSFGs to transition into unobscured gas-poor QSOs via a far-infrared bright QSO (Simpson et al. 2012; Costa et al. 2018). In the latter case, this would imply that G09v1.40 is in a phase just prior to a highly obscured QSO phase.

Most or all of the gas in the neutral outflow, however, is likely to remain bound to the galaxy in all but the least inclined scenarios, where the deprojected radii are large. This gas may then be reaccreted by the galaxy at a later time, replenishing the gas reservoir, unless additional feedback, such as thermal feedback from a previously obscured QSO, causes the gas to remain in the CGM.

While the current observations provide sufficient information to determine global properties and offsets between the dust, gas, and stellar components, analysis of the detailed morphological and kinematic structures will require new observations at the higher spatial resolution, including a determination of the true inclination of the neutral outflow in G09v1.40.

Finally, we note that our analysis of the outflow in G09v1.40 using OH⁺/CH⁺ only probes the diffuse neutral component, on one side of the galaxy. The full multiphase, double-sided outflow will carry even more mass, momentum, and energy, likely shifting the conclusions of this paper to more extreme scenarios. Future observations, targeting other phases of the outflow (e.g., molecular and ionized) in both absorption and emission lines, are needed to fully constrain the impact of this outflow on the evolution of G09v1.40.

The authors thank the referee for their many, and very appreciated, suggestions, questions, and corrections. The authors thank Jae Calanog for providing the reduced and lens-subtracted Keck NIR images introduced in Section 2.2. D.R. acknowledges support from the National Science
Foundation under grant Nos. AST-1614213 and AST-1910107. D.R. also acknowledges support from the Alexander von Humboldt Foundation through a Humboldt Research Fellowship for Experienced Researchers. M.J.M. acknowledges the support of the National Science Centre, Poland through the SONATA BIS grant 2018/30/E/ST9/00208. M.R. acknowledges the support of the Veni research program with project number 202.225 and the Vidi research program with project number 639.042.611, which are (partly) financed by the Dutch Research Council (NWO). This paper makes use of the following ALMA data: ADS/JAO. ALMA#2015.1.01042. S, ADS/JAO. ALMA#2013.1.00164. S, and ADS/JAO. ALMA2016.1.00282. 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.

*Facilities:* ALMA, Keck.

*Software:* CASA (v4.53; McMullin et al. 2007), *visilens* (Hezaveh et al. 2013; Spilker et al. 2016) and *LENSTOOL* (Kneib et al. 1996; Jullo et al. 2007; Jullo & Kneib 2009).

## Appendix A

### Beam-smearing Effects on Source Reconstruction of Gravitational Lenses

When we observe the sky, the spatial distribution of the sky emission is convolved with the shape of the beam. This has the effect of smearing the light emitted from structures smaller than the beam over its point-spread function. The consequence of this in our observations is that our sources appear larger and fuzzier than they intrinsically are (see upper rows of Figures 19 and 20). For a gravitationally lensed source sitting on or very close to the inner caustic of the gravitational lens, the light from

one side of the caustic will be smeared over and onto the other side of the caustic. When mapping this light back into the source plane, it will not be reconstructed in the correct position.

This is most obviously demonstrated by the dust continuum in G09v1.40, whose simple double-image configuration and faint Einstein ring in the image plane (Figure 1) imply an intrinsic source plane geometry of a single extended source, partially overlapping the inner caustic but with its peak situated just to the west. However, when the image plane intensity map is reconstructed, two peaks appear in the source plane, one to the west of the inner caustic as expected and another weaker mirror image of the peak on the opposite side of the caustic.

To investigate the role of beam smearing in this process, we reconstruct the model of the dust continuum intensity produced by the cleaning procedure during data reduction. This provides an indication of the continuum intensity distribution without beam convolution. We also reconstruct maps of the model map convolved with artificial beams with axes 25%, 50%, 75%, and 100% as long as the beam in our observations. The source plane reconstruction of the fully deconvolved model reveals a single source to the west of the inner caustic with no mirroring image on the east. As the beam size is increased, the artifact to the east of the caustic appears and grows (Figure 19).

We repeat this experiment on the moment 0 map of the model CO(9-8) emission. Again we note that the model is only an indication of the true deconvolved intensity distribution. Again, the fully deconvolved map produces a single source when reconstructed in the source plane. The CO(9-8) emission is more extended and elongated than the dust continuum and sits directly over the inner caustic and extending to the northeast and southwest of the caustic. A source overlapping the inner caustic should form a quadruply imaged lens configuration in the image plane, which is not obvious in our observations but is evident in the middle-upper panel of Figure 20. The southwest image is in fact two merged

<figure><img src="image_20.png" /><figcaption>**Figure 19.** Illustrating the effect of beam smearing on our observations and reconstruction of the dust continuum in G09v1.40. Top row: image plane model of the dust continuum, derived from the cleaning method during data reduction, with no beam convolution (far-left panel) and convolved with a beam FWHM at 25%, 50%, 75%, and 100% of the beam in our observations (consecutive panels to the right). Beam sizes are shown by the shaded orange ellipse in the upper left of each panel and the dash white region indicates the region enlarged in the lower panels. Bottom row: source plane reconstructions of the panels above with lensing caustics shown in white. As the data are convolved with increasingly larger beam sizes, more of the flux in the image plane is smeared past the Einstein ring and is then reconstructed on the wrong side of the inner caustic in the source plane. This causes the artifact to the east of the inner caustic in the source plane to become more severe.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 20.** Illustrating the effect of beam smearing on our observations and reconstruction of the CO(9–8) emission in G09v1.40. Top row: image plane model of the CO(9–8) emission, derived from the cleaning method during data reduction, with no beam convolution (far-left panel) and convolved with a beam FWHM 25%, 50%, 75%, and 100% that of the beam in our observations (consecutive panels to the right). Beam sizes are shown by the shaded orange ellipse in the upper left of each panel, and the dashed white region indicates the region enlarged in the lower panels. Bottom row: source plane reconstructions of the panels above with the inner lensing caustics shown in black and Einstein ring shown in white. As the data are convolved with increasingly larger beam sizes, more of the flux in the image plane is smeared across the Einstein ring and is then reconstructed on the wrong side of the inner caustic in the source plane. This causes flux from locations in the source plane on opposite sides of the caustic (northeast and southwest) to be blended together after reconstruction and results in the spurious northwestern artifact which becomes more severe with larger beam sizes.</figcaption></figure>

images of the source, which are not separately resolved by our beam. As the model is convolved with larger beam sizes, flux emitted from positions within the source on either side of the caustic are blended and a spurious artifact appears and grows to the northwest of the caustic (see Figure 20).

## Appendix B
### OH⁺ 1₀⁻¹₁ Optical Depth to Column Density

We trace the neutral gas outflow in G09v1.40 with the 1033 GHz OH⁺ 1₀⁻¹₁ line in absorption. Both levels are split into two by hyperfine structure, resulting in four energy levels, which we label 0, 1, 2, and 3, in order of increasing energy. The line that we detect is thus the sum of four different absorption lines: 0 → 2, 0 → 3, 1 → 2, and 1 → 2. The Einstein A coefficients, level energies, and quantum numbers of each transition can be found at Splatalogue (https://splatalogue.online/).

For a single-absorption line (ignoring stimulated emission) the optical depth is given by

$$ \tau_{\nu}^{\text{single}} = N_{l} \frac{\lambda^{3}}{8\pi} \frac{g_{u}}{g_{l}} \frac{A_{ul}}{\sqrt{2\pi\sigma}} e^{-\nu^{2}/2\sigma^{2}}, \quad (B1) $$

where $g_u$ and $g_l$ are the statistical weights of the upper and lower levels respectively, $A_{ul}$ the Einstein coefficient of the transition from the upper to lower energy state, $N_l$ the column density of the lower level, and $\lambda_{ul}$ is the wavelength of the line, which is assumed to be Gaussian centered on $\nu$ with velocity dispersion $\sigma$. Assuming all the OH⁺ molecules are in the ground state ($N_l \gg N_u$) the quantity of interest is $N_l$, the column

density in the lower, absorbing level, which in this case is the combination of the two hyperfine ground-state levels $N_1 + N_0$.

Integrating the optical depth over velocity we obtain for a single line

$$ \int \tau_{\nu}^{\text{single}} dv = N_{l} \frac{\lambda^{3}}{8\pi} \frac{g_{u}}{g_{l}} \frac{A_{ul}}{\sqrt{2\pi\sigma}} \int e^{-\nu^{2}/2\sigma^{2}} dv \quad (B2) $$

$$ = N_{l} \frac{\lambda^{3}}{8\pi} \frac{g_{u}}{g_{l}} A_{ul}. \quad (B3) $$

Summing all four lines we therefore obtain

$$ \int \tau_{\nu} dv = \frac{1}{8\pi} \left( \frac{g_3}{g_0} A_{30} \lambda_{30}^3 N_0 + \frac{g_3}{g_1} A_{31} \lambda_{31}^3 N_1 \right. \\ \left. + \frac{g_2}{g_0} A_{20} \lambda_{20}^3 N_0 + \frac{g_2}{g_1} A_{21} \lambda_{21}^3 N_1 \right). \quad (B4) $$

From Splatalogue: $A_{20} = 7.03 \times 10^{-3} \text{s}^{-1}$, $A_{30} = 1.76 \times 10^{-2} \text{s}^{-1}$, $A_{21} = 1.41 \times 10^{-2} \text{s}^{-1}$, and $A_{31} = 3.53 \times 10^{-3} \text{s}^{-1}$; $g_1 = g_2 = 2$ and $g_0 = g_3 = 4$; and $\lambda_{30} \approx \lambda_{31} \approx \lambda_{20} \approx \lambda_{21} \approx 0.29 \text{ cm}$. This leads to

$$ \int \tau_{\nu} dv = \frac{\lambda^3}{8\pi} \left( \frac{1}{2} A_{20} N_0 + A_{30} N_0 + A_{21} N_1 + 2 A_{31} N_1 \right) \\ = \frac{\lambda^3}{8\pi} (2.11 \times 10^{-2} N_0 + 2.11 \times 10^{-2} N_1), \quad (B5) $$

which can simply be written as

$$ \int \tau_{\nu} dv = \frac{\lambda^3}{8\pi} A_{\text{OH}^+}, \quad (B6) $$
where the Einstein coefficient is given by $A = 2.11 \times 10^{-2} \text{ s}^{-1}$.

## ORCID iDs

Kirsty M. Butler https://orcid.org/0000-0001-7387-0558

Matus Rybak https://orcid.org/0000-0002-1383-0746

Axel Weiß https://orcid.org/0000-0003-4678-3939

Dominik A. Riechers https://orcid.org/0000-0001-9585-1462

Lucia Marchetti https://orcid.org/0000-0003-3948-7621

Ivan Valtchanov https://orcid.org/0000-0001-9930-7886

## References

Aalto, S., Garcia-Burillo, S., Muller, S., et al. 2015, A&A, 574, A85

Alatalo, K., Blitz, L., Young, L. M., et al. 2011, ApJ, 735, 88

Benson, A. J., Bower, R. G., Frenk, C. S., et al. 2003, ApJ, 599, 38

Berta, S., Young, A. J., Cox, P., et al. 2021, A&A, 646, A122

Bialy, S., Burkhart, B., & Sternberg, A. 2017, ApJ, 843, 92

Bialy, S., Neufeld, D., Wolfire, M., et al. 2019, ApJ, 885, 109

Bolatto, A. D., Warren, S. R., Leroy, A. K., et al. 2013, Natur, 499, 450

Bower, R. G., Benson, A. J., & Crain, R. A. 2012, MNRAS, 422, 2816

Bregman, J. N. 1980, ApJ, 236, 577

Bregman, J. N., Miller, E. D., Seitzer, P., et al. 2013, ApJ, 766, 57

Brusa, M., Cresci, G., Daddi, E., et al. 2018, A&A, 612, A29

Bussmann, R. S., Pérez-Fournon, I., Amber, S., et al. 2013, ApJ, 779, 25

Calanog, J. A., Fu, H., Cooray, A., et al. 2014, ApJ, 797, 138

Calderón, D., Bauer, F. E., Veilleux, S., et al. 2016, MNRAS, 460, 3052

Calistro Rivera, G., Hodge, J. A., Smail, I., et al. 2018, ApJ, 863, 56

Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105

Carniani, S., Marconi, A., Maiolino, R., et al. 2017, A&A, 605, A105

Casey, C. M., Narayanan, D., & Cooray, A. 2014, PhR, 541, 45

Cazzoli, S., Arribas, S., Maiolino, R., et al. 2016, A&A, 590, A125

Chen, C.-C., Hodge, J. A., Smail, I., et al. 2017, ApJ, 846, 108

Chisholm, J., Tremonti, C. A., Leitherer, C., et al. 2016, MNRAS, 457, 3133

Cicone, C., Maiolino, R., Gallerani, S., et al. 2015, A&A, 574, A14

Cicone, C., Maiolino, R., Sturm, E., et al. 2014, A&A, 562, 21

Cochrane, R. K., Best, P. N., Smail, I., et al. 2021, MNRAS, 503, 2622

Cole, S. M., Baugh, C., Frenk, C., et al. 2000, RSPTA, 358, 2093

Costa, T., Rosdahl, J., Sijacki, D., et al. 2018, MNRAS, 479, 2079

Danielson, A. L. R., Swinbank, A. M., Smail, I., et al. 2011, MNRAS, 410, 1687

Danielson, A. L. R., Swinbank, A. M., Smail, I., et al. 2013, MNRAS, 436, 2793

Decarli, R., Walter, F., Venemans, B. P., et al. 2018, ApJ, 854, 97

Enia, A., Negrello, M., Gurwell, M., et al. 2018, MNRAS, 475, 3467

Falgarone, E., Ossenkopf, V., Gerin, M., et al. 2010, A&A, 518, L118

Falgarone, E., Zwaan, M. A., Godard, B., et al. 2017, Natur, 548, 430

Fan, L., Knudsen, K. K., Fogasy, J., et al. 2018, ApJL, 856, L5

Feruglio, C., Ferrara, A., Bischetti, M., et al. 2017, A&A, 608, A30

Feruglio, C., Maiolino, R., Piconcelli, E., et al. 2010, A&A, 518, L155

Fielding, D., Quataert, E., & Martizzi, D. 2018, MNRAS, 481, 3325

Fischer, J., Sturm, E., González-Alfonso, E., et al. 2010, A&A, 518, L41

Fluetsch, A., Maiolino, R., Carniani, S., et al. 2019, MNRAS, 483, 4586

Fluetsch, A., Maiolino, R., Carniani, S., et al. 2021, MNRAS, 505, 5753

Foreman-Mackey, D., Hogg, D. W., Lang, D., et al. 2013, PASP, 125, 306

Förster Schreiber, N. M., & Wuyts, S. 2020, ARA&A, 58, 661

Fujimoto, S., Ouchi, M., Shibuya, T., et al. 2017, ApJ, 850, 83

Gallerani, S., Pallottini, A., Feruglio, C., et al. 2018, MNRAS, 473, 1909

García-Burillo, S., Combes, F., Usero, A., et al. 2014, A&A, 567, A125

García-Burillo, S., Combes, F., Usero, A., et al. 2015, A&A, 580, A35

Gerin, M., De Luca, M., Black, J., et al. 2010, A&A, 518, L110

Gerin, M., Neufeld, D. A., & Goicoechea, J. R. 2016, ARA&A, 54, 181

Ginolfi, M., Jones, G. C., Béthermin, M., et al. 2020, A&A, 633, A90

Godard, B., & Cernicharo, J. 2013, A&A, 550, A8

Godard, B., Falgarone, E., & Pineau Des Forêts, G. 2009, A&A, 495, 847

Godard, B., Falgarone, E., & Pineau des Forêts, G. 2014, A&A, 570, A27

González-Alfonso, E., Fischer, J., Bruderer, S., et al. 2013, A&A, 550, A25

González-Alfonso, E., Fischer, J., Bruderer, S., et al. 2018, ApJ, 857, 66

González-Alfonso, E., Fischer, J., Spoon, H. W. W., et al. 2017, ApJ, 836, 11

Governato, F., Brook, C., Mayer, L., et al. 2010, Natur, 463, 203

Harrison, C. M., Costa, T., Tadhunter, C. N., et al. 2018, NatAs, 2, 198

Heckman, T. M., Armus, L., & Miley, G. K. 1990, ApJS, 74, 833

Heckman, T. M., & Borthakur, S. 2016, ApJ, 822, 9

Henriques, B. M. B., Yates, R. M., Fu, J., et al. 2020, MNRAS, 491, 5795

Herrera-Camus, R., Förster Schreiber, N., Genzel, R., et al. 2021, A&A, 649, A31

Herrera-Camus, R., Janssen, A., Sturm, E., et al. 2020, A&A, 635, A47

Herrera-Camus, R., Tacconi, L., Genzel, R., et al. 2019, ApJ, 871, 37

Hezaveh, Y. D., Marrone, D. P., Fassnacht, C. D., et al. 2013, ApJ, 767, 132

Hodge, J. A., Riechers, D., Decarli, R., et al. 2015, ApJL, 798, L18

Hollenbach, D., Kaufman, M. J., Neufeld, D., et al. 2012, ApJ, 754, 105

Hopkins, P. F., Grudić, M. Y., Wetzel, A., et al. 2020, MNRAS, 491, 3702

Hopkins, P. F., Kereš, D., Oñorbe, J., et al. 2014, MNRAS, 445, 581

Indriolo, N., Bergin, E. A., Falgarone, E., et al. 2018, ApJ, 865, 127

Indriolo, N., Neufeld, D. A., Gerin, M., et al. 2015, ApJ, 800, 40

Ivison, R. J., Smail, I., Papadopoulos, P. P., et al. 2010, MNRAS, 404, 198

Jones, G. C., Maiolino, R., Caselli, P., et al. 2019, A&A, 632, L7

Jullo, E., & Kneib, J.-P. 2009, MNRAS, 395, 1319

Jullo, E., Kneib, J.-P., Limousin, M., et al. 2007, NJPh, 9, 447

Kennicutt, R. C. 1998, ARA&A, 36, 189

Kereš, D., Katz, N., Davé, R., et al. 2009, MNRAS, 396, 2332

Kneib, J.-P., Ellis, R. S., Smail, I., et al. 1996, ApJ, 471, 643

Lehnert, M. D., & Heckman, T. M. 1996, ApJ, 462, 651

Ma, B., Cooray, A., Calanog, J. A., et al. 2015, ApJ, 814, 17

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Maiolino, R., Gallerani, S., Neri, R., et al. 2012, MNRAS, 425, L66

Martin, C. L. 2005, ApJ, 621, 227

McMullin, J. P., Waters, B., Schiebel, D., et al. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw et al. (San Francisco, CA: ASP), 127

Meyer, D. M., & York, D. G. 1987, ApJL, 315, L5

Mitchell, P. D., Schaye, J., & Bower, R. G. 2020, MNRAS, 497, 4495

Murray, N., Quataert, E., & Thompson, T. A. 2005, ApJ, 618, 569

Negrello, M., Amber, S., Amvrosiadis, A., et al. 2017, MNRAS, 465, 3558

Negrello, M., Hopwood, R., De Zotti, G., et al. 2010, Sci, 330, 800

Nelson, D., Pillepich, A., Springel, V., et al. 2019, MNRAS, 490, 3234

Neufeld, D. A., Goicoechea, J. R., Sonnentrucker, P., et al. 2010, A&A, 521, L10

Neufeld, D. A., & Wolfire, M. G. 2017, ApJ, 845, 163

Pereira-Santaella, M., Colina, L., García-Burillo, S., et al. 2018, A&A, 616, A171

Pereira-Santaella, M., Colina, L., García-Burillo, S., et al. 2020, A&A, 643, A89

Pillepich, A., Nelson, D., Springel, V., et al. 2019, MNRAS, 490, 3196

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13

Porras, A. J., Federman, S. R., Welty, D. E., & Ritchey, A. M. 2014, ApJL, 781, L8

Rangwala, N., Maloney, P. R., Glenn, J., et al. 2011, ApJ, 743, 94

Rees, M. J., & Ostriker, J. P. 1977, MNRAS, 179, 541

Riechers, D. A., Bradford, C. M., Clements, D. L., et al. 2013, Natur, 496, 329

Riechers, D. A., Capak, P. L., Carilli, C. L., et al. 2010, ApJL, 720, L131

Riechers, D. A., Cooray, A., Perez-Fournon, I., et al. 2021a, ApJ, 913, 141

Riechers, D. A., Nayyeri, H., Burgarella, D., et al. 2021b, ApJ, 907, 62

Ritchey, A. M., Welty, D. E., Dahlstrom, J. A., et al. 2015, ApJ, 799, 197

Rubin, K. H. R., Prochaska, J. X., Koo, D. C., et al. 2014, ApJ, 794, 156

Rupke, D. S., Veilleux, S., & Sanders, D. B. 2002, ApJ, 570, 588

Rupke, D. S., Veilleux, S., & Sanders, D. B. 2005, ApJS, 160, 115

Rupke, D. S. N., Gültekin, K., & Veilleux, S. 2017, ApJ, 850, 40

Rupke, D. S. N., & Veilleux, S. 2013, ApJ, 768, 75

Rybak, M., Vegetti, S., McKean, J. P., et al. 2015, MNRAS, 453, L26

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Sharma, P., Roy, A., Nath, B. B., et al. 2014, MNRAS, 443, 3463

Shopbell, P. L., & Bland-Hawthorn, J. 1998, ApJ, 493, 129

Simcoe, R. A., Sargent, W. L. W., & Rauch, M. 2004, ApJ, 606, 92

Simpson, J. M., Smail, I., Swinbank, A. M., et al. 2012, MNRAS, 426, 3201

Simpson, J. M., Smail, I., Wang, W.-H., et al. 2017, ApJL, 844, L10

Somerville, R. S., Behroozi, P., Pandya, V., et al. 2018, MNRAS, 473, 2714

Somerville, R. S., & Primack, J. R. 1999, MNRAS, 310, 1087

Spilker, J. S., Aravena, M., Béthermin, M., et al. 2018, Sci, 361, 1016

Spilker, J. S., Aravena, M., Phadke, K. A., et al. 2020a, ApJ, 905, 86

Spilker, J. S., Marrone, D. P., Aravena, M., et al. 2016, ApJ, 826, 112

Spilker, J. S., Phadke, K. A., Aravena, M., et al. 2020b, ApJ, 905, 85

Spoon, H. W. W., Farrah, D., Lebouteiller, V., et al. 2013, ApJ, 775, 127

Stacey, H. R., McKean, J. P., Powell, D. M., et al. 2021, MNRAS, 500, 3667

Stone, M., Veilleux, S., Meléndez, M., et al. 2016, ApJ, 826, 111

Strickland, D. K., Heckman, T. M., Colbert, E. J. M., et al. 2004, ApJS, 151, 193

Sturm, E., González-Alfonso, E., Veilleux, S., et al. 2011, ApJL, 733, L16
van der Tak, F. F. S., Weiß, A., Liu, L., et al. 2016, A&A, 593, A43

van der Werf, P. P., Isaak, K. G., Meijerink, R., et al. 2010, A&A, 518, L42

Vayner, A., Wright, S. A., Murray, N., et al. 2017, ApJ, 851, 126

Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, ARA&A, 43, 769

Veilleux, S., Maiolino, R., Bolatto, A. D., et al. 2020, A&ARv, 28, 2

Veilleux, S., Meléndez, M., Sturm, E., et al. 2013, ApJ, 776, 27

Walter, F., Weiss, A., & Scoville, N. 2002, ApJL, 580, L21

Weiß, A., Walter, F., Downes, D., et al. 2012, ApJ, 753, 102

Westmoquette, M. S., Clements, D. L., Bendo, G. J., et al. 2012, MNRAS, 424, 416

White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341

Yang, C., Omont, A., Beelen, A., et al. 2017, A&A, 608, A144

Zhang, Z.-Y., Ivison, R. J., George, R. D., et al. 2018, MNRAS, 481, 59

Zhao, D., Galazutdinov, G. A., Linnartz, H., & Krawłowski, J. 2015, ApJL, 805, L12