Molecular Outflow in the Reionization-epoch Quasar J2054-0005 Revealed by OH 119 μm Observations

Molecular outflows are expected to play a key role in galaxy evolution at high redshift. To study the impact of outflows on star formation at the epoch of reionization, we performed sensitive Atacama Large Millimeter/submillimeter Array observations of OH 119 μm toward J2054-0005, a luminous quasar at z = 6.04. The OH line is detected and exhibits a P-Cygni profile that can be fitted with a broad blueshifted absorption component, providing unambiguous evidence of an outflow, and an emission component at near-systemic velocity. The mean and terminal outflow velocities are estimated to be v out ≈ 670 and 1500 km s−1, respectively, making the molecular outflow in this quasar one of the fastest at the epoch of reionization. The OH line is marginally spatially resolved for the first time in a quasar at z > 6, revealing that the outflow extends over the central 2 kpc region. The mass outflow rate is comparable to the star formation rate ( Ṁout/SFR∼2 ), indicating rapid (∼107 yr) quenching of star formation. The mass outflow rate in a sample of star-forming galaxies and quasars at 4 < z < 6.4 exhibits a positive correlation with the total infrared luminosity, although the scatter is large. Owing to the high outflow velocity, a large fraction (up to ∼50%) of the outflowing molecular gas may be able to escape from the host galaxy into the intergalactic medium.


INTRODUCTION
Quasar feedback is one of the fundamental processes that regulate galaxy evolution.Galaxies acquire gas via accretion from the intergalactic medium (IGM) and through merging, and lose a fraction of their gas via galactic outflows powered by feedback from starbursts Corresponding author: Dragan Salak dragan@oia.hokudai.ac.jp and/or active galactic nuclei (AGNs).The baryon cycle is believed to regulate how much molecular gas is available for star formation and the growth of the central supermassive black holes (SMBHs) (e.g., Murray et al. 2005;Veilleux et al. 2005Veilleux et al. , 2020;;Hopkins et al. 2012;Zubovas & King 2014;Tumlinson et al. 2017).
Recent observations have revealed the existence of massive (stellar mass M ⋆ ∼ 10 11 M ⊙ ) galaxies with diminished star formation activity already at z ≳ 6, indicating that these objects experienced a vigorous star-burst episode at an earlier epoch followed by quenching (e.g., Straatman et al. 2014;Glazebrook et al. 2017;Girelli et al. 2019;Merlin et al. 2019;Carnall et al. 2020Carnall et al. , 2023;;Forrest et al. 2020;Valentino et al. 2020;Santini et al. 2021;Labbé et al. 2023;Looser et al. 2023;Nanayakkara et al. 2023).When and how this quenching occurred is still debated, but quasar feedback is one of the possible mechanisms considered as a leading internal process to explain the rapid (< 1 Gyr) inside-out quenching of star formation in massive galaxies (Tacchella et al. 2015;Barai et al. 2018;Costa et al. 2018;Spilker et al. 2019;Bischetti et al. 2021;Mercedes-Feliz et al. 2023).Quasars at z ∼ 6, powered by AGN and star formation, are thus believed to be an important evolutionary phase in massive galaxy evolution (e.g., Carilli & Walter 2013;Casey et al. 2014;Lapi et al. 2018).To understand how massive galaxies evolved, it is important to reveal the physical conditions of the interstellar medium (ISM) in quasars at the epoch of reionization (EoR; 6 ≲ z ≲ 20), and this has been at the focus of recent research (e.g., Venemans et al. 2017Venemans et al. , 2018;;Walter et al. 2018;Hashimoto et al. 2019a;Novak et al. 2019;Li et al. 2020a,b;Neeleman et al. 2021;Meyer et al. 2022;Pensabene et al. 2022;Decarli et al. 2023).Many of these quasars exhibit high far-infrared luminosities (L FIR ≳ 10 13 L ⊙ ) that indicate dust heating by extreme star formation and AGN radiation (e.g., Walter et al. 2009;Wang et al. 2013;Leipski et al. 2014;Venemans et al. 2018Venemans et al. , 2020)).Extrapolating from the known properties of local galaxies (Fluetsch et al. 2019;Lutz et al. 2020;Roberts-Borsani 2020), it is expected that energy released in such nuclear activity is sufficient to generate galactic outflows.Since molecular gas is the primary fuel for star formation and SMBH growth, it is important to investigate the molecular phase, which has been largely untraced, at the EoR.
While AGN-driven outflows have been observed extensively at low/moderate redshifts, and the relations between the outflows and the physical properties of host galaxies investigated (e.g., Rupke et al. 2005Rupke et al. , 2021;;Feruglio et al. 2010;Veilleux et al. 2013;Cicone et al. 2014;Fiore et al. 2017;Gowardhan et al. 2018;Förster Schreiber et al. 2019;Lutz et al. 2020), our understanding of quasar feedback in the early Universe has been limited.This is because detecting outflows from the emission lines of standard tracers such as CO relies on identification of broad wings in the spectra that are generally much weaker than the main component of the line profile and requires a high signal-to-noise ratio (e.g., Cicone et al. 2014).Most high-z outflow studies are based on searches of broad wings in the emission spectra of [C II] 158 µm and CO lines (e.g., Decarli et al. 2018;Novak et al. 2020) and extended halos of cold gas (Fujimoto et al. 2019(Fujimoto et al. , 2020)), but unambiguous detections of outflows (distinguished from inflows) using these lines are still rare at the EoR and beyond (Izumi et al. 2021), making it difficult to evaluate the quasar-driven feedback.
The 119 µm doublet can unambiguously reveal the presence of cold molecular outflows and/or inflows through its P-Cygni profile (e.g., Veilleux et al. 2013;Herrera-Camus et al. 2020).Since the energy required for the excitation of OH 2 Π 3/2 from the rotational ground state J = 3/2 to the state J = 5/2 is E/k ≈ 120 K, where k is the Boltzmann constant, cold gas (≲ 100 K) is observed in absorption against a bright continuum source.On the other hand, the gas density required to thermalize the rotational transitions of OH is very high (n H2 ≳ 10 9 cm −3 ), so the transition can be observed in J = 5/2 → 3/2 emission in environments where molecular gas is highly excited (dense warm gas, either through shocks, or because it is exposed to strong far-infrared continuum radiation), such as those in AGNs (Veilleux et al. 2013).
At high redshift, previous works have showed that OH outflows can readily be detected in strongly lensed, dusty star-forming galaxies up to z ≈ 5 (George et al. 2014;Spilker et al. 2018Spilker et al. , 2020a,b),b); there are also reports of two OH detections in quasars at z ≈ 6 (Butler et al. 2023), and one tentative (Herrera-Camus et al. 2020).Interestingly, the results in Spilker et al. (2020a) suggest that OH 119 µm may be a more reliable tracer of line-of-sight outflows at high z than [C II] 158 µm and CO lines.It is therefore of great interest, and motivation of this work, to investigate whether the 119 µm line can provide a good probe of outflows at the EoR.
To search for molecular outflows in EoR quasars, we observed OH 119 µm toward J2054-0005 using the Atacama Large Millimeter/submillimeter Array (ALMA).The quasar was discovered from the Sloan Digital Sky Survey (SDSS) data (Jiang et al. 2008), and later detected by ALMA in continuum as well as [C II] 158 µm, [O III] 88 µm, and CO lines (Wang et al. 2010(Wang et al. , 2013;;Hashimoto et al. 2019a;Venemans et al. 2020).The measurements of these lines have determined its redshift to be z = 6.0391 ± 0.0002.The bolometric luminosity of the source is L bol ≈ 1.2 × 10 47 erg s −1 , corresponding to 3.2 × 10 13 L ⊙ (Farina et al. 2022).The total IR luminosity of L IR ≈ 1.3 × 10 13 L ⊙ suggests a star formation rate (SFR) of ≈ 1900 M ⊙ yr −1 (Hashimoto et al. 2019a), where a Kroupa (2001) initial mass function (IMF) is assumed, although this is an upper limit because of possible AGN contribution to dust heating (Schneider et al. 2015;Di Mascia 2023).As discussed in Section 5.1 below, the contribution of AGN to L IR is estimated to be ≈ 59% in this source, yielding a lower star formation rate of SFR ≈ 800 M ⊙ yr −1 .However, despite the intense star formation and the presence of an AGN, neither [C II] 158 µm, [O III] 88 µm, nor CO lines have revealed outflows in previous observations.Is there no outflow, or is it difficult to detect it with these tracers?Establishing a reliable tracer of molecular outflows is essential for future studies of galaxies at highest redshifts.
The paper is organized as follows.In Section 2, we describe the ALMA observations and data reduction.The resulting continuum image and OH 119 µm spectrum is presented in Section 3.This is followed by an analysis of the OH gas outflow in Section 4, discussion on the outflow's driving mechanism and imprint on galaxy evolution in Section 5, and a summary in Section 6.

OBSERVATIONS AND DATA REDUCTION
The observations were conducted between August 4 and 12 in 2022 during ALMA cycle 8.The antennas of the 12 m array observed toward a single field centered at (α, δ) ICRS = (20 h 54 m 06.s 503, −00 • 05 ′ 14. ′′ 43), which corresponds to the central position of the ALMA 87 µm continuum (Hashimoto et al. 2019a).In most observing runs, 44 antennas were used, but the number ranged from 41 to 46.The array was in configuration C-5 with baselines from 15 m to 1301 m.
The Band 7 receivers were tuned to cover the OH doublet at an observing frequency for the adopted red-shift z = 6.0391.Two spectral windows (upper sideband; USB) were centered at the observing frequencies 356.315GHz and 358.090GHz for the line observations.Since the bandwidth of each of them is 1.875 GHz, this setup makes the two spectral windows next to each other with an overlap of 100 MHz.The 119.23 µm line of the doublet was set to lie in this overlap region, whereas the 119.44 µm line is separated in velocity by ≈ 521 km s −1 .With a frequency resolution of 7.813 MHz, the achieved velocity resolution (average over the bandwidth) is 6.56 km s −1 .To improve the signal-to-noise ratio for the analysis, we smoothed the data cubes to a resolution of 35 km s −1 .There are 44 velocity channels in each window, with one perfectly overlapping channel, and the total effective velocity coverage of the two adjacent windows is 3045 km s −1 with the 119.23 µm line at the center.
The other two spectral windows (lower sideband; LSB) were dedicated to continuum observations.The central observing frequencies were 344.2 GHz and 346.1 GHz, and the bandwidth of each of them was 2 GHz.
The scheduling block was executed 9 times.J2253+1608 (8 data sets) and J1924-2914 (1 data set) were observed for the flux and bandpass calibration, whereas J2101+0341 (all data) was observed for the phase calibration.The total on-source time was 7.3 hours, whereas the total time including calibrator observations and other overheads was 12.7 hours.The uncertainties in the paper are only statistical errors; the absolute flux accuracy in Band 7 is reported to be 10% (Braatz et al. 2021).
The data were reduced using the Common Astronomy Software Applications (CASA) package (The CASA team et al. 2022).Basic calibration was performed with a CASA pipeline resulting in 9 calibrated measurement sets.The calibrated data were then combined and imaged using the CASA task tclean.
The continuum image was reconstructed using the line-free LSB spectral windows in the multi-frequency synthesis mode with the hogbom deconvolver and standard gridding.The weighting was set to briggs with the robust parameter equal to 0.5 (compromise between resolution and sensitivity).To conduct unbiased mask-based image reconstruction, we employed the auto-multitresh tool (Kepley et al. 2020).We also tried tclean in interactive mode, but there was no obvious difference in the result, so we adopted the automatically-created image.The threshold for iterations was set to 2 σ, where σ was calculated using the CASA task imstat on a first-generation clean image masked for the central region where the source is located.The rms sensitivity in the final continuum image is σ = 13 µJy beam −1 .The synthesized beam size (full-width half maximum; FWHM) is (b maj , b min ) = (0.′′ 205, 0. ′′ 176), corresponding to ≈ 1 kpc.For the adopted cosmological parameters, 1 ′′ is equivalent to 5.689 kpc and the luminosity distance to the source is D L = 58.1465Gpc.These are the highest spatial resolution observations of OH toward a quasar at z > 6.
The OH image was reconstructed from USB spectral windows.The visibilities of two spectral windows were processed by tclean separately and continuum was not subtracted.The deconvolver and weighting setups were the same as for the continuum, and the auto-multitresh tool was used.The mean sensitivity over two spectral windows is σ = 0.13 mJy beam −1 in a channel of 35 km s −1 , and the synthesized beam size is (b maj , b min ) = (0.′′ 204, 0. ′′ 174).The two spectral windows were merged using the task imageconcat.
The velocity is expressed in radio definition with respect to the rest frame of the source (z = 6.0391).The frequency that corresponds to zero velocity is ν = ν rest (1 + z) −1 = 357.193± 0.010 GHz, where ν rest = 2514.31640360GHz is the rest frequency of the OH 2 Π 3/2 J = 3/2−5/2, F = 3 − −2 + transition, 1 where F is the quantum number for the total angular momentum of the molecule (including nuclear spin), and " + " and " − " denote the Λ-doubling of energy levels.
The final images were corrected for the primary beam attenuation.The basic parameters of the resulting images are summarized in Table 1.
We begin this section by a presentation of the continuum image.It is followed by a description of the OH spectrum and derivation of its basic properties.
We also found an additional continuum source positioned 2. ′′ 4 west of the quasar and detected at S/N = 8.9 (Figure 2).Since OH is not detected there, it is not clear at this point whether the source is a physical companion or is at different redshift and happens to lie within the solid angle subtended by the primary beam.The projected separation from J2054-0005 corresponds to ≈ 14 kpc if they are at the same redshift.The peak intensity is measured to be 116 ± 13 µJy beam −1 at (α, δ) ICRS = (20 h 54 m 06.s 344, −0 • 05 ′ 14. ′′ 83), and the flux density is S ν (r < 0. ′′ 5) = 377 ± 9 µJy.

OH gas
The OH 119 µm line is robustly detected toward the quasar.Figure 3 shows an integrated OH spectrum (total flux density) extracted from the region where the 123−µm (LSB) continuum is detected at > 3 σ.We selected this broad region because there is a possibility that OH is distributed throughout the galactic disk  2).To guide the eye, the contours are plotted at (5, 50) × σ, where σ = 1.28382 × 10 −5 Jy beam −1 .The image is corrected for the primary beam attenuation, so the noise is increased at the edges.Right.Close-up view of the continuum in J2054-0005.The contours are plotted at (−3, 3, 5, 10, 20, 40, 80, 160, 240) × σ.The beam size is shown at the bottom left as a filled ellipse.
Although there are reasons to assume that OH is not optically thin, e.g., if the gas distribution is clumpy and the absorbing gas does not cover the continuum entirely, the apparent optical depth at the line center, where absorption is maximum, is τ ≈ 0.36.
In the vicinity of the OH doublet, there are CH + (J = 3−2) at the sky frequency of 355.364GHz, and 18 OH (J = 5/2−3/2) at 355.021 GHz, that may be responsible for the decrease in flux that appears at the offset velocity v ≳ 1000 km s −1 .A similar feature attributed to the 18 OH line is found in Mrk 231, indicating an enhanced [ 18 OH]/[OH] abundance due to processing by star formation (Fischer et al. 2010;González-Alfonso et al. 2014).However, the lines are sufficiently separated from OH and unlikely to significantly affect the analysis of the line profile described below.

OH line fitting
To investigate the kinematics of OH gas, we performed a least-squares fitting of the line profile using Python (scipy.optimize.curvefit).Since the line is a doublet, and there may be multiple components (systemic, outflow, or inflow), fitting was conducted using two double-gaussian functions.Although the spectrum is likely to be more complicated than this simple structure, we aimed at limiting the number of free parameters while extracting key quantities related to the analysis of outflows.The fitting constraints were the following: the separation between the doublet lines of a double gaussian is fixed to 521 km s −1 , and their peak values and FWHMs are set to be equal (e.g., Goicoechea & Cernicharo 2002).The continuum intensity within the velocity range covered by the two adjacent USB spectral windows was assumed to be constant.On the other hand, the line intensity, continuum intensity, and the velocity separation of the double gaussians of different components (e.g., systemic and outflow) were set as free parameters.
The results of fitting are shown in Figure 3 and listed in Table 2.We find OH emission, traced by a double gaussian with an FWHM line width of 306 ± 55 km s −1 near the systemic velocity, and a broad absorption feature with FWHM = 1052 ± 234 km s −1 with the peak absorption velocity of v cen = −669 ± 87 km s −1 relative to the systemic velocity.The FWHM of the emission line obtained from fitting is comparable to those of the emission lines of The terminal velocity (maximum extent of the blueshifted wing of the absorption) is at least v max ≈ −1500 km s −1 but may be beyond the spectral coverage.Another indicator of terminal velocity is v 98 , the velocity above which 98% of absorption takes place.This quantity is found to be v 98 = −1574 ± 35 km s −1 .On the other hand, the velocity above which 84% of absorption takes place, is v 84 = −1104 ± 35 km s −1 .The fitted emission components are red-shifted relative to the systemic by 65 ± 15 km s −1 .This is not unusual, as such positive shifts in emission components have been observed in the majority of nearby galaxies that exhibit P-Cygni profiles and likely arise from outflows on the opposite side of the continuum (receding relative to the observer) (Veilleux et al. 2013).The uncertainties above include only those from fitting.The redshift uncertainty expressed in velocity is ≈ 8 km s −1 .
The continuum flux density in the USB spectral windows was found by fitting to be S cont = 6.175 ± 0.092 mJy.This is ≈ 8% higher than the LSB continuum (see Section 3.1), but not unexpected, because the continuum emission at this frequency is in the Rayleigh-Jeans domain.Also, the regions where the fluxes were extracted (LSB continuum within r < 0. ′′ 5, USB continuum where LSB continuum is detected at > 3 σ) are not equal, but are similar in size.Given the fact that there are almost no line-free channels in the USB spectrum and that fitting was done under simple assumptions, we find this to be a reasonable estimate.
Figure 4 shows a continuum-subtracted OH spectrum, including the best fit and fitting residuals.

MOLECULAR GAS OUTFLOW
The absorption line is tracing the OH gas between the continuum source and the observer.Since this is observed in blueshift relative to the host galaxy, the line reveals unambiguous signature of an outflow, as it has been observed toward AGNs and star-forming galaxies at lower z (e.g., Veilleux et al. 2013).For example, the nearby ULIRG Mrk 231 exhibits a similar terminal velocity of ≈ −1500 km s −1 and a P-Cygni line profile, which are attributed to an AGN-driven outflow (Fischer et al. 2010;Spoon et al. 2013).Since we do not detect red-shifted absorption, there is no significant molecular gas inflow in J2054-0005 along the line of sight.In the analysis below, we consider the absorption to be tracing an outflow, and derive the basic properties, such as mass, mass outflow rate, and kinetic energy.

Outflow mass
Assuming that the outflow has the shape of a thin spherical shell of radius r out , the molecular gas mass in the outflow (including helium and heavier elements) can be expressed as (Veilleux et al. 2020) where µ = 1.36 is the mean particle mass per hydrogen nucleus, m H is the hydrogen atom mass (in kilograms), N H = 2N H2 is the column density of hydrogen nuclei (in m −2 ), and Ω is the solid angle ("opening angle") subtended by the shell (Veilleux et al. 2020).Thus, the term Ωr 2 out is the area of the outflowing shell, and the opening angle is given by Ω = 4πf , where f is the dimensionless covering factor, i.e., the fraction of the sphere covered by the outflow as seen from its origin.
The column density N H2 cannot be measured directly, but we can find the column density of OH molecules (N OH ) and then apply a [OH]/[H 2 ] abundance ratio to get N H2 .N OH can be calculated from the measured absorption line under the approximation of local thermodynamic equilibrium (LTE).The OH column density can be expressed as (e.g., Mangum & Shirley 2015) where λ ul = 119.23 µm is the rest-frame wavelength, A ul = 0.1388 s −1 is the Einstein coefficient for the u → l transition J = 5/2 → 3/2 (Schöier et al. 2005), g u = 6 is the statistical weight of the J = 5/2 level, ∆E = E u −E l is the energy difference of the two levels (E l = 0), T ex is the excitation temperature, and Q is the partition function.If τ ≪ 1, the integral is approximately equal to the equivalent width, defined as W = (1−e −τ )dv.Assuming T ex = 50 K, which is equal to the dust temperature (T d = 50 ± 2 K; Hashimoto et al. 2019a), we get Q ≃ 19 (Pickett et al. 1998).The term Q(1 − e −∆E/kTex ) −1 increases by a factor of ≈ 3 from 50 to 150 K.Note that the equivalent width in Equation ( 3) is calculated for a single line of the doublet (W in Table 2) because Q accounts for the Λ-doubling.
The OH abundance has been measured for the Milky Way galaxy (e.g., Goicoechea & Cernicharo 2002;Goicoechea et al. 2006;Nguyen et al. 2018;Rugel et al. 2018) and a number of nearby galaxies where multiple OH lines were detected (Spinoglio et al. 2005;Falstad et al. 2015;Stone et al. 2018), and it has been analyzed in theoretical work employing simulations (Richings & Faucher-Giguère 2018).A relatively high value of [OH]/[H 2 ] = 5×10 −6 is reported for Sgr B2 (Goicoechea & Cernicharo 2002).Although this value is often used to derive H 2 mass from OH, it may be lower at sub-solar metallicities that may be applicable to high-z sources.This is, however, not necessarily the case for evolved systems, as near-solar metallicities have been found in some quasars at z > 6 (Novak et al. 2019;Li et al. 2020b;Onoue et al. 2020).On the other hand, a low abundance of [OH]/[H 2 ] ≈ 1×10 −7 has been reported recently even  4).The quantities SOH and W are given for a single line of the doublet (the total equivalent width is twice the value).The continuum flux density is Scont = 6.175 ± 0.092 mJy.
for various regions in the Galaxy (Nguyen et al. 2018;Rugel et al. 2018).We derive the outflow properties using [OH]/[H 2 ] = 5 × 10 −7 (Table 3) as a reasonable compromise between the two extremes.Choosing this value is also motivated by the fact that it yields an outflow mass and mass outflow rate comparable to those obtained using an empirical relation presented in Section 4.2.2 below.
The radius of the outflow (r out ), as given in Equation 2, is the radius of a thin outflow shell (Veilleux et al. 2020).We adopt a radius of r out = 2 kpc, which is the size of the 123 µm continuum shown in Figure 1 and the region from which the spectrum in Figure 3 is extracted.This radius is consistent with discussion in Section 4.3 below, where it is argued that the OH outflow is extended to at least r ≳ 1 kpc from the center.
For the covering factor, we use f = 0.3 in the analysis below.This value is also adopted in Spilker et al. (2020b) for a sample of star-forming galaxies at redshift z = 4−5, although it can be as large as f ≈ 0.8 (Spilker et al. 2018).Even if we adopt f = 1, the main results of the scaling relations discussed in Section 5 do not change, albeit the mass outflow rates would be larger by a factor of ≈ 3.
The integral in Equation ( 3) is calculated by inserting τ from Equation (1), where S(v) is the profile obtained from gaussian fitting and S cont is a constant (Table 2).The equivalent width of a single line of the doublet is W = 200 km s −1 , yielding M out ≈ 4.9 × 10 9 M ⊙ .The obtained W is larger than that found in almost all nearby ULIRGs (Veilleux et al. 2013) and dusty star-forming galaxies at z = 4−5 (Spilker et al. 2020a), and is largest in a quasar at z > 6 reported to date.The calculated outflow gas mass and other dynamical quantities (derived below) are listed in Table 3.The outflow mass is ∼ 8−16% of the total molecular gas mass (Decarli et al. 2022).

Optically-thin outflow model
Assuming that the outflow is expanding at velocity v out as a thin spherical shell, the mass outflow rate averaged over the outflow lifetime can be expressed as This equation gives a conservative estimate (Maiolino et al. 2012;Lutz et al. 2020;Salak et al. 2020;Veilleux et al. 2020).For the outflow velocity, we adopt v out = 669 km s −1 , based on the center velocity (v cen ) of the absorption feature (see Table 2 and Section 3.3).Although this is the mean value along the line of sight, it is equal to the outflow velocity in the case of a sphericallysymmetric outflow.Taking W = 200 km s −1 (equivalent width of the absorption line), [OH]/[H 2 ] = 5×10 −6 , r out = 2 kpc, and f = 0.3, we obtain a lower limit of Ṁout > 168 M ⊙ yr −1 .Note that if we adopt a low abundance of [OH]/[H 2 ] = 1 × 10 −7 and f = 1, the upper limit of the mass outflow rate becomes Ṁout ≈ 2.8 × 10 4 M ⊙ yr −1 .However, the upper limit yields an outflow mass that exceeds the total molecular gas mass in the host galaxy (Decarli et al. 2022).Clearly, the uncertainty of the mass outflow rate is dominated by the poorly constrained OH abundance.We adopt a moder- The dynamical age of the outflow is defined as the time needed for outflowing gas to travel a distance r out at a constant velocity v out .Using the derived quantities above, the outflow age is t out = r out /v out ∼ 3 × 10 6 yr.The timescale is much shorter than the depletion time due to star formation (∼ 4−8 × 10 7 yr).

Empirical relation
Alternatively, we can circumvent the above assumptions and use an empirical formula for the mass outflow rate discussed in the literature, hoping that it is applicable to the EoR quasar.The "recipe" formula from Herrera-Camus et al.Using the mass outflow rate from Equation ( 6), we calculate the outflow mass M emp out = (r out /v out ) Ṁ emp out and other dynamical quantities.All outflow properties derived from this empirical relation are listed in Table 3.Some caveats of this approach include the fact that Equation (6) only incorporates the equivalent width at v < −200 km s −1 (to exclude the systemic component assuming that it does not exceed this velocity) regardless of the line width of the outflow-tracing absorption line and the mean outflow velocity.None the less, the mass outflow rate and outflow mass obtained using Equation ( 6) are in agreement with those obtained under LTE and moderate OH abundance (Table 3).

Resolved OH absorption and emission
The high angular resolution (≈ 1 kpc) and sensitivity allow us to probe the spatial distribution of OH gas velocity for the first time in a quasar at z > 6.We extracted OH spectra from adjacent rectangular regions of area 0. ′′ 1 × 0. ′′ 1, corresponding to approximately one half of the synthesized beam, and performed double-gaussian fitting in each region using the procedure described in Section 3.3.To obtain successful fits with a minimum number of free parameters, we fit all regions with only two double-gaussians.A moment 1 image of OH that shows the positions of the regions as pixels is shown in Figure 5.The spectra in Figure 6 were extracted from the 12 pixels in Figure 5, labelled (179,178) in the bottom left corner, (181, 181) in the top right corner, etc., and shown as a 3 × 4 profile map.All spectra that could yield reasonable fits are plotted in Figure 6 together with the best fits.The OH doublet absorption was successfully fitted in 11 rectangular regions (Figure 6 and Table 4).The profile could not be well-fitted in the surrounding regions, where the continuum intensity is weaker and OH is not significantly detected.
The absorption line is marginally spatially resolved, which can be inferred from a north-south shift in the peak absorption and emission velocities obtained by fitting.This implies that the distribution of OH gas is not confined to the AGN, which is too compact to be resolved, but extends over a broader (r ≳ 1 kpc) central region.The FWHM line widths (Table 4) are relatively comparable throughout the region (≈ 800−1300 km s −1 ).The peak velocity is not minimum (most negative) at the continuum center, as may be expected from a spherically symmetric outflow emerging from the center, but retains comparable values throughout the map.Although the fitting uncertainties are large, the fitted absorption peak velocities appear to be more negative on the south side compared to the north side, though region (180,181) seems to deviate from this trend.The mean absorption velocities in each row in Figure 6 from north to south are (−648 ± 80, −565 ± 61, −667 ± 47, −806 ± 70) km s −1 , excluding (181,180).These results suggest that the outflow is not uniform and might have the shape of a cone whose axis is inclined with respect to the line of sight.Observations at higher resolution are needed to get a clearer picture of the outflow geometry.
OH is detected in emission at > 3 σ in some regions and exhibits relatively comparable FWHM line widths throughout the region (≈ 200−330 km s −1 ), consistent with reported [C II] 158 µm, µm and [O III] 88 µm, and CO (J = 6 → 5) line widths (Wang et al. 2010(Wang et al. , 2013;;Hashimoto et al. 2019a).This suggests that highly excited (warm or dense, shocked) molecular gas is distributed in the central 2 kpc region, either in the host galaxy or in the outflowing gas (the size of the region where emission could be fitted is ≈ 2 kpc in diameter).The positive peak velocities of the emission lines are generally lower in the north compared to the south.The exception is at (181,178), though the line is only marginally detected there.The mean emission veloc-ities in each row in Figure 6 from north to south are (15 ± 11, 67 ± 12, 101 ± 11, 72 ± 13) km s −1 .
Pixel (181,179) was fitted with an absorption feature and a very narrow emission feature (Figure 6).The latter is likely an artifact, because the line width is too narrow and the emission line is not significantly detected here.Pixel (181,180) was successfully fitted with two different absorption features.Since we did not put constraints on whether the fit should yield absorption or emission, the fitting turned out to be more successful with two absorption features than one absorption and one emission features at this pixel.
In order to investigate the origin of the apparent velocity shift in the OH emission line, we compare the OH data with [C II] 158 µm data.The moment 1 image in Wang et al. (2013) shows that the [C II] 158 µm line exhibits a velocity gradient in the northwest-southeast direction, that is generally consistent with the northsouth velocity increase in the fitted OH emission line spectra, albeit with an offset: the fitted OH emission is systematically redder than [C II] 158 µm.Thus, it is possible that the OH emission follows the velocity field of the bulk gas in the host galaxy traced by [C II] 158 µm.
Figure 7 shows a comparison of the OH 119 µm and [C II] 158 µm spectra (the [C II] data are from #2019.1.00672;S. Fujimoto, in prep.), extracted from the central pixel (maximum 123 µm continuum intensity; pixel size 0. ′′ 034).These are the highest-resolution [C II] data of this source and therefore best for comparison.In addition to the main emission profile, the [C II] line profile exhibits a secondary component on the redshifted side (up to +500 km s −1 ), and there appears to be a blue-shifted component at velocities comparable to the OH outflow velocity (−600 km s −1 ), although it is tentative at this choice of aperture.This is consistent with recent findings that OH 119 µm is a more robust tracer of outflows compared to [C II] 158 µm (Spilker et al. 2020a).
More details of the [C II] observations and results will be presented in Fujimoto et al. (in prep.).

DISCUSSION
In this section, we first estimate the fractional contribution of AGN to the IR luminosity, and then discuss on the driving mechanism of the outflow, the fate of the outflowing gas, and its impact on the host galaxy.

Fractional contribution of AGN to IR luminosity
Estimating the SFR in high-z sources is often done by applying a conversion factor to the far-IR luminos-  ity.However, in quasars, the IR luminosity is produced not just by dust heated by star formation, but also by dust heated by the AGN.To obtain a reliable estimate of SFR, it is therefore important to subtract the fractional contribution of the AGN (e.g., Duras et al. 2017).
Here, we applied a multi-wavelength spectral energy distribution (SED) modeling of the spectrum of J2054-0005 using the latest version of Code Investigating Galaxy Emission (CIGALE; Boquien et al. 2019).CIGALE is a state-of-the-art modeling and fitting code that combines a broad range of components, including a stellar population, AGN, and dust.For each parameter, CIGALE conducts a probability distribution function (PDF) analy-sis, yielding the output value as the likelihood-weighted mean of the PDF.For the SED fitting of J2054-0005, we used the following data: SDSS-z, WISE1, Herschel (PACS, SPIRE), ALMA bands 6 and 7, and the upper limits from Herschel 500 µm and the Very Large Array 1.4 GHz (Bañados et al. 2015;Shao et al. 2019).
We assigned a flexible star formation history composed of a delayed component to account for a recent burst (e.g., Donevski et al. 2020).The assumed values for the main stellar population are set to be between 300 and 750 Myr, with an e-folding time of 90 Myr, and a Chabrier (2003) IMF.The gas-phase metallicity was fixed to be two times lower than the solar value.This is  motivated by the fact that the fitting result was somewhat more successful (in terms of χ 2 ) with this value compared to that based on solar metallicity (e.g., Novak et al. 2019), although both yielded a similar L IR .Dust attenuation was modeled using a modified law from Charlot & Fall (2000), and dust emission was modeled based on Draine & Li (2007).The model assumes a mixture of grains exposed to variable radiation fields.We fixed the dust emission slope to β = 2, and allowed a range of radiation field intensities (10 < U min < 50).This approach can account for a very intense central heating source.
The AGN module used to determine the fractional contribution of the AGN to the total IR luminosity is based on Fritz et al. (2006).The model takes into account three components through radiative transfer: the primary source located in the torus, the scattered emission by dust, and the thermal dust emission.We fixed the optical depth to 2, following some prescriptions from the literature (e.g., Ciesla et al. 2017).For the sake of computational efficiency, we modeled two inclination angles of the torus (30 • and 70 • ; Mountrichas et al. 2019).
The result of the procedure is shown in Figure 8.A reasonably-well fit was obtained, as indicated by a median reduced χ 2 = 0.52 ± 0.11.We found that the total IR luminosity is L IR = (1.34 ± 0.17) × 10 13 L ⊙ , consistent with previous studies (Hashimoto et al. 2019a), and that the fractional contribution of AGN to L IR is f AGN = 0.59 ± 0.08.Thus, the AGN may be responsible for as much as ≈ 59% of L IR in this source.This is comparable to the fractions reported for quasars at high z (e.g., Schneider et al. 2015;Duras et al. 2017;Tripodi et al. 2022).Accounting for this effect, the IR-derived star formation rate becomes SFR = 770 ± 180 M ⊙ yr −1 , calculated using SFR/M ⊙ yr −1 = 1.40 × 10 −10 L ′ IR /L ⊙ , where L ′ IR is the AGN-subtracted IR luminosity (integrated within λ rest = 8 − 1000 µm).The conversion factor corresponds to the Chabrier IMF, and was calculated by dividing the Kroupa IMF factor in Murphy et al. (2011) by 1.06 (e.g., Zahid et al. 2012;Speagle et al. 2014).We adopt this value in the analysis below.
The best fit was obtained for the torus inclination angle of 30 • .Based on low-z studies, this value is consistent with a Type 1 AGN (e.g., Yang et al. 2020).It is also consistent with a relatively broad line width (FWHM ≈ 4890 km s −1 ) of Lyα emission reported in Jiang et al. (2008).

Driving mechanism
Is the outflow driven by star formation or does the AGN feedback from the accretion onto the supermassive black hole (SMBH) play a role?One way to address this problem is to investigate the energy and momentum of the outflow and compare them to the expected input from star formation.
Using the outflow mass and velocity derived in Section 4, we calculate the bulk kinetic energy of the outflowing gas.Adopting the molecular gas mass derived under LTE (Table 3), the kinetic energy is and the power required to drive the molecular outflow (kinetic power) is This is equivalent to ≈ 6.2 × 10 10 L ⊙ , which is ≈ 0.5% of the total infrared luminosity of the source.If other ISM phases (atomic and ionized gas) are present in the outflow, the energy and power are larger.The total momentum of the molecular outflow is By comparison, the momentum injection by a typical core-collapse supernova (SN; mass m 0 ≈ 10 M ⊙ , velocity v 0 ≈ 3000 km s −1 ) is of the order of p 0 ≈ 3 × 10 4 M ⊙ km s −1 , and the total momentum of an outflow driven by SN explosions is p SN ≈ p 0 R SN t SN , where R SN is the SN rate, and t SN is the time interval measured from the onset of SN explosions.Thus, R SN t SN is the total number of SN explosions over the starburst episode.Adopting a relation between the SN rate and star formation rate, R SN /yr −1 ≈ α SN (SFR/M ⊙ yr −1 ), where α SN ≈ 0.01−0.02depends on IMF and assumes continuous star formation (Veilleux et al. 2020), and SFR = 770 M ⊙ yr −1 , we obtain R SN ≈ 15 yr −1 (for α = 0.02).
If the time interval of the SN feedback is equal to the dynamical age of the outflow (t SN = t out ), the total SN momentum becomes p SN ≈ 1.3 × 10 12 M ⊙ km s −1 , produced by R SN t SN ≈ 4.5 × 10 7 SN explosions.Theoretical works suggest that the final momentum input per SN may be even higher (e.g., Kim & Ostriker 2015;Walch & Naab 2015).Moreover, the total momentum could be larger by a factor of two if stellar winds (radiation pressure) from massive stars play a significant role (e.g., Leitherer et al. 1999;Murray et al. 2010).Taking into account the possibility that other gas phases also participate in the outflow, so that the total momentum (molecular, neutral atomic, and ionized gas) is somewhat larger, it seems that star formation activity alone may be approximately sufficient to explain the observed outflow.Although we obtain this result based on a moderate relative OH abundance, a similar conclusion is reached if the empirical relation for the mass outflow rate is used.On the other hand, the mean outflow velocity of 670 km s −1 and the terminal velocity of 1500 km s −1 exceed the outflow velocities typically measured in starforming galaxies, which are found to be 100−500 km s −1 and < 1000 km s −1 , respectively (e.g., Sugahara et al. 2019), although the outflow velocities are found to be higher in more extreme systems (Spilker et al. 2020a).By contrast, terminal velocities of outflows in AGNdominated systems are often found to ≳ 1000 km s −1 and as large as 1500 km s −1 (e.g., Rupke et al. 2005;Sturm et al. 2011;Spoon et al. 2013;Veilleux et al. 2013;Ginolfi et al. 2020).Theoretical studies also reproduce the velocities of ≳ 1000 km s −1 and mass outflow rates of ∼ 10 3 M ⊙ yr −1 in AGN-driven outflows (e.g., Ishibashi et al. 2018;Costa et al. 2018).The results indicate that the AGN feedback (radiation pressure) may play a role in boosting the velocity in J2054-0005.
Generally, there is an increase in v out with L IR , though the relation is not clear for the quasars.One possibility is that this is simply because the AGN contribution in driving the outflow in J2054-0005 is relatively large, making the velocity higher than what it would be if star formation were the only driving mechanism.On the other hand, as discussed in Butler et al. (2023), if the outflows in these quasars are anisotropic (e.g., conical), it is possible that the random orientation results in no correlation because OH in most high-z sources is unresolved.Further observations including emission lines at higher resolution are necessary to clarify the outflow geometry.

Suppression of star formation
The mass loading factor, defined as the ratio of the mass outflow rate to star formation rate, is an indicator of the outflow impact on star formation activity.In J2054-0005, using the AGN-corrected SFR calculated in Section 5.1, we find Using these values, we calculate the depletion time, i.e., the time it takes for the outflow to remove molecular gas from the galactic center region, as The relatively short timescale, as compared to the time required for star formation to consume molecular gas in ordinary star-forming galaxies (∼ 10 9 yr), implies that J2054-0005 is undergoing an episode of rapid quenching of star formation.The depletion time is shorter by an order of magnitude compared to that in nearby ULIRGs where OH outflows have been detected (González-Alfonso et al. 2017).On the other hand, since molecular gas is being evacuated from the galactic center region, the outflow is also quenching the gas reservoir available to feed the SMBH (M BH ≈ 2 × 10 9 M ⊙ ; Farina et al. 2022), thereby limiting its growth.
The result has important implications for the evolution of the central stellar component and its coevolution with the SMBH, which is believed to be the origin of the relation between the SMBH mass and bulge velocity dispersion found in a large sample of galaxies from the local Universe to high redshifts (e.g., Murray et al. 2005;Kormendy & Ho 2013;Izumi et al. 2019).The short depletion timescales inferred from this work, as well as other dynamical properties, agree with recent predictions from models of massive galaxy formation at z > 6 (Lapi et al. 2018;Pantoni et al. 2019).The models predict that the quasar outflow phase is accompanied by significant stellar component increase up to within ∼ 30 Myr, and support the making of quiescent galaxies recently observed at redshifts z ≈ 3−5 (Carnall et al. 2023).Our result and other recent OH and OH + observations (Shao et al. 2022;Butler et al. 2023) indicate that cool gas outflows may be common in quasars at z > 6, but further observations of larger samples are needed to investigate their statistical properties.

Ṁout −L IR relation
Figure 10 shows a comparison of mass outflow rates and the total IR luminosity (L IR ) in high-z quasars and dusty star-forming galaxies (DSFGs) with OH outflow detections.Since the uncertainty of the mass outflow rate is dominated by the OH abundance, we plot the product W v out r out in addition to Ṁout , because it is the directly measured quantity that is proportional to Ṁout in the expanding-shell model under LTE, whereas Ṁout depends on OH abundance, f , and T ex .Here, W is the equivalent width of the outflow component, v out is the center velocity of the absorption line, and r out is the adopted radius of the OH outflow (Section 4.1).For all sources, we assumed that the radius is equal to the size of the continuum-emitting region in the host galaxy from which the OH spectrum was extracted.This is r out ≈ 2.0 kpc for J2054-0005, r out ≈ 2.5 kpc for J2310+1855, and r out ≈ 4.0 kpc for P183+05, based on Figure 2 in Butler et al. (2023).For the DSFGs, the outflow radius is assumed to be the continuum size r out = r cont from Table 2 in Spilker et al. (2020a).W was calculated using Equation ( 4), and we assumed T ex = 50 K and [OH]/[H 2 ] = 5×10 −7 for all sources in calculating Ṁout .
The plot shows that there is a positive relation between the mass outflow rate with L IR , and the power law is log y = k log x + m, where k = 1.61 ± 0.40 and m = −16.4± 5.2, although the scatter is large (R 2 = 0.64).
As discussed in Herrera-Camus et al. ( 2020) and Spilker et al. (2020b), the correlation can also be found if the mass outflow rate is set to be proportional to Ṁout ∝ W √ L IR .This relation circumvents the uncertainty of the quantities such as the outflow radius and OH optical depth.We reproduce a similar relation in Figure 11, although W here is the equivalent width of The ordinate on the left shows the product of measured quantities (equivalent width, outflow velocity, and radius) that is proportional to Ṁout in the expanding-shell model; Ṁout = 6.28 × 10 −3 W voutrout, where W is in km s −1 , vout is in km s −1 , and rout is in kpc.The black circles are dusty starforming galaxies (DSFGs) at z = 4−5 (Spilker et al. 2020a).The data of J2310+1855 and P183+05 are taken from Butler et al. (2023), and W was calculated using Equation ( 4).The ordinate on the right is Ṁout calculated using Equations ( 3) and ( 5).The dotted line is η = 1 and the solid line is the fit.The open square symbol (not included in the fit) shows J2054-0005 corrected for the AGN contribution.The data points for J2310+1855 and P183+05 are not AGN-corrected.
the outflow (absorption) component instead of W v<200 as in their works.The relation is nearly linear (plotted as solid line) with k = 1.13 ± 0.36 and intercept m = −6.3± 4.7 (R 2 = 0.52).
The above analysis yields a positive relation between Ṁout and L IR , though it should be noted that L IR gives an upper limit to the SFR if the AGN fraction is not subtracted, and the fraction may differ among the sources.We have corrected the SFR by subtracting the fractional AGN contribution to dust heating for J2054-0005.Although the sample of quasars is small, Figure 10 suggests that the Ṁout − L IR relation may be less clear when such correction is made.The outflows may be driven by the combined contribution of the starburst and AGN (e.g., Gowardhan et al. 2018), which is responsible for the positive Ṁout −L IR relation and the scatter may be caused by nonuniform outflow geometry, projection effects, and different covering factors.On the other hand, Butler et al. (2023) suggest that if a quasar is unobscured, it has already cleared the material near the nucleus, thereby leaving a less energetic outflow due to inefficient coupling with the surrounding ISM material.Since the AGN in J2054-0005 appears to be Type 1, i.e., unobscured (Section 4.2; Jiang et al. 2008) but the outflow is significantly more powerful compared to those in J2310+1855 and P183+05 (also Type 1), the obscuration effects do not seem to play a major role in outflow energetics.

Gas escaping into the IGM
Previous observations have revealed the presence of atomic and molecular gas in the halos and circumgalactic medium of high-redshift galaxies (e.g., Emonts et al. 2016;Tumlinson et al. 2017;Fujimoto et al. 2020;Cicone et al. 2021;Scholtz et al. 2023).This suggests that the gas was ejected from host galaxies by powerful outflows with terminal velocities that may even exceed the escape velocity.Here, we make a simple analysis to investigate whether the molecular outflow in J2054-0005 is fast enough to transport OH gas into the IGM.
The dynamical mass of J2054-0005, derived from [C II] 158 µm data for a disk inclination angle of 24 • , is estimated to be M dyn ≈ 7.2 × 10 10 M ⊙ (Wang et al. 2013) within the [C II]-emitting region of radius R ≈ 1 kpc.Assuming a spherically symmetric mass distribution, this yields an escape velocity of v esc (R) ≈ 2GM dyn /R ≈ 780 km s −1 at R. The estimate corresponds to a rotational velocity of 560 km s −1 , implying a very massive host galaxy.The escape velocity is similar to the velocities reported for dusty star-forming galaxies at z = 4−5 (Spilker et al. 2020b).Since the obtained value is comparable to the mean velocity of the outflow along the line of sight, a significant fraction (up to ∼ 50%) of the outflowing molecular gas may be able to escape the gravitational potential well and inject metals and dust into the IGM.This is in agreement with re-cent observations of enriched gas in the circumgalactic medium at z ∼ 6 (Wu et al. 2021).

[O III]88/[C II]158 luminosity ratio
Recently, the luminosity ratio of [O III] 88 µm to [C II] 158 µm has drawn attention as it is found to be higher at high-redshift compared to local star-forming galaxies (e.g., Inoue et al. 2016;Laporte et al. 2019;Hashimoto et al. 2019a,b;Pallottini et al. 2019;Tamura et al. 2019;Arata et al. 2020;Bakx et al. 2020;Carniani et al. 2020;Harikane et al. 2020;Lupi et al. 2020;Vallini et al. 2021;Katz et al. 2022;Sugahara et al. 2022;Witstok et al. 2022;Ren et al. 2023) and local dwarf galaxies (Ura et al. 2023).One of the scenarios proposed to explain the observations is the possibility of outflows affecting the covering factor of photodissociation regions (PDRs) (Harikane et al. 2020).In sources with powerful outflows, low-ionization PDRs traced by [C II] 158 µm may be cleared so that their covering factor is decreased relative to that of H II regions traced by [O III] 88 µm.If that is the case, we may expect to see more powerful outflows in sources with a high luminosity ratio.
So far, only two reionization-epoch quasars (J2054-0005 and J2310+1855) with OH outflow detections have been detected also in [C II] 158 µm and [O III] 88 µm (Hashimoto et al. 2019a;Butler et al. 2023).A comparison of their properties shows the following characteristics.(1) The luminosity ratio L [OIII] /L [CII] is ∼ 7 times larger in J2054-0005 (2.1±0.4) compared to J2310+1855 (0.3±0.1).(2) The mass outflow rate in J2054-0005 is ∼ 3 times larger than in J2310+1855 (Figure 10); the mean outflow velocity is also higher (669±87 km s −1 in J2054-0005 compared to 334 ± 14 km s −1 in J2310+1855).On the other hand, J2310+1855 has slightly larger total IR luminosity (1.9 × 10 13 L ⊙ ) compared to J2054-0005 (1.3 × 10 13 L ⊙ ).Although we cannot draw conclusions from only two sources, J2054-0005, with a more powerful outflow, also has a significantly higher luminosity ratio, which supports the scenario of a low PDR covering factor.Previous [C II] 158 µm studies suggest that the inclination angle of a rotating host galaxy in J2054-0005 is relatively low (≈ 24 • ; Wang et al. 2013).If that is the case, and if the outflow is predominantly propagating perpendicular to a rotating disk, it is close to the line of sight, hence the observed velocity is higher compared to that in J2310+1855.

OH emission: highly excited molecular gas
The OH 119 µm line has been detected in emission toward a number of nearby AGNs (Spinoglio et al. 2005;Veilleux et al. 2013).Recently, Butler et al. (2023) reported a detection of OH emission in one quasar at z ≈ 6.However, the line has been detected only in absorption toward a sample of dusty star-forming galaxies at z = 4−5 (Spilker et al. 2020a).
Based on multi-line OH observations, Spinoglio et al. (2005) argue that collisional excitation dominates the population of the 2 Π 3/2 J = 5/2 level that leads to radiative decay and 119 µm emission in the active nucleus of the local Seyfert galaxy NGC 1068.On the other hand, Veilleux et al. (2013) found that the strength of the OH 119 µm absorption relative to emission is correlated with the 9.7 µm silicate strength, an indicator of the obscuration of the nucleus, in their ULIRG sample.Spoon et al. (2013) argue that the OH emission arises from dustobscured central regions and that, except in two outliers, radiative excitation may be dominant.Since the peak of the spectral energy distribution of dust thermal emission in J2054-0005 is close to λ rest = 53 µm, the wavelength that corresponds to the energy difference between the levels 2 Π 1/2 J = 3/2 and 2 Π 3/2 J = 3/2, absorption of the continuum from dust emission could excite the level which would radiatively decay into the ground state.In that case, it is expected that the 2 Π 1/2 J = 3/2 → 1/2 line at λ rest = 163 µm would also be observed in emission.Given that the 120 µm continuum emission is spatially extended, and the fact that OH emission is marginally spatially resolved (Section 4.3), it is likely that the excited OH gas is not confined to the compact AGN, but distributed in a broader (∼ 2 kpc) region.Further multi-line observations are necessary to constrain the excitation mechanism of OH molecules in EoR quasars.
1.The OH 119.23,119.44 µm doublet line and the 120 µm continuum are detected toward the quasar.
The continuum is detected at high signal-to-noise ratio of 260.The OH line exhibits a P-Cygni profile with absorption and emission components.
2. We fitted the OH profile with two double-gaussian functions using a least-squares fitting tool.The fits reveal a blue-shifted absorption component (peak absorption depth τ ≈ 0.36), that unambiguously reveals as outflowing molecular gas, and emission component at near-systemic velocity.The absorption peak velocity is v cen = −669 ± 87 km s −1 , the FWHM line width is 1052 ± 234 km s −1 , and the terminal velocity is v 98 = −1574±35 km s −1 , indicating a fast molecular outflow.This is the first quasar with such high molecular outflow velocity discovered at z > 6.
3. The mass outflow rate, calculated under LTE approximation, OH abundance [OH]/[H 2 ] = 5 × 10 −7 , and assuming the geometry of an expanding thin spherical shell with a covering factor f = 0.3 and radius r out = 2 kpc, is Ṁout ≈ 1700 M ⊙ yr −1 .Using an empirical relation from the literature, we obtain Ṁ emp out ≈ 1500 M ⊙ yr −1 .4. The absorption and emission components of the OH line are marginally spatially resolved in the central 2 kpc, suggesting that the outflow extends over this region.Since the critical density and excitation energy for the upper rotational levels of OH are relatively high (n cr ≳ 10 9 cm −3 , E/k = 120 K), the detection of OH emission implies that molecular gas is highly excited (warm or dense, shocked), possibly by far-IR radiation pumping from dust grains and by collisions with H 2 (e.g., shocks).The OH line is significantly broader compared to [C II] 158 µm in the central 1-kpc region.
5. In order to estimate the fractional contribution of AGN to the total infrared luminosity (L IR ), we performed SED fitting of the spectrum of J2054-0005 using the code CIGALE.The result yields a total infrared luminosity of L IR = (1.34 × 0.17) × 10 13 L ⊙ and suggests that as much as 59% of L IR is produced by dust heated by the AGN and not by star formation.The IR-derived SFR corrected for this effect is estimated to be SFR = 770 ± 180 M ⊙ yr −1 (Chabrier IMF).
6.The mass outflow rate is comparable to the AGNcorrected star formation rate in the host galaxy ( Ṁout /SFR ∼ 2); it is higher compared to other two quasars with OH detections at z > 6 and among the highest at high redshift.At the current mass loss rate, molecular gas is expected to be depleted after only t dep ≈ (2−4) × 10 7 yr, implying rapid quenching of star formation.The dynamical age of the outflow is t out = r out /v out ≈ 3 × 10 6 yr.
7. An analysis of the outflow energetics and terminal velocity indicates that the outflow in J2054-0005 may be powered by the combined effects of the AGN and star formation.This is supported by the fact that we find a positive correlation between Ṁout and total luminosity L IR using a sample of 8 dusty star-forming galaxies at z = 4−5 and 3 quasars at z > 6, that the outflow velocity in J2054-0005 is relatively high compared to the majority of star-forming galaxies at high z, and that as much as 59% of L IR is produced by the AGN.
8. The mean outflow velocity is comparable to the estimated escape velocity.This implies that as much as ∼ 50% of the outflowing molecular gas may be able to escape from the host galaxy and enrich the intergalactic medium with heavy elements.9. We report the discovery of a companion at a projected separation of 2. ′′ 4.This source is detected only in continuum at the significance of 8.9 σ.

Figure 1 .
Figure 1.Left.The 123 µm continuum image.The solid rectangle indicates the position of J2054-0005; a close-up view is shown in the right panel.The dashed rectangle indicates the position of a projected companion (see Figure2).To guide the eye, the contours are plotted at (5, 50) × σ, where σ = 1.28382 × 10 −5 Jy beam −1 .The image is corrected for the primary beam attenuation, so the noise is increased at the edges.Right.Close-up view of the continuum in J2054-0005.The contours are plotted at(−3, 3, 5, 10, 20, 40, 80, 160, 240) × σ.The beam size is shown at the bottom left as a filled ellipse.

Figure 3 .
Figure3.Integrated OH 119 µm spectrum extracted from the region that includes all pixels where the 123−µm continuum is detected at > 3 σ.The dashed blue and dotted red curves are the absorption and emission components, respectively, determined from two double-gaussian line fitting, and the solid magenta curve is their sum (see Section 3.3).The offset velocity is measured with respect to the 119.23 µm line, as indicated by a vertical line.
[C II] 158 µm (243 ± 10 km s −1 ; Wang et al. 2013), [O III] 88 µm (282 ± 17 km s −1 ; Hashimoto et al. 2019a), and CO (J = 6 → 5) (360 ± 110 km s −1 ; Wang et al. 2010).However, [O III] 88 µm is a tracer of ionized gas, and it is unclear whether [C II] 158 µm is traces the same molecular medium as OH does.Only CO is directly comparable to OH as a molecular gas tracer, and the line widths of the two are in agreement with each other.

Figure 4 .
Figure 4. Continuum-subtracted OH 119 µm spectrum with the computed fit (as in Figure 3) and its residuals shown at the bottom.

W
(2020) modified by Spilker et al. (2020b) takes the form Ṁ emp out M ⊙ yr −1 = 1.4 where W v<−200 is the OH 119 µm equivalent width at v < −200 km s −1 .This equivalent width, derived from the observed spectrum, is W v<−200 ≈ 268 km s −1 , yielding a mass outflow rate of Ṁ emp out ≈ 1500 M ⊙ .The value is significantly larger compared to those for the star-forming galaxies at redshift z = 4−5 reported in Spilker et al. (2020b), which have values between 220 and 1290 M ⊙ yr −1 .

Figure 5 .
Figure 5. Moment 1 image of the OH data calculated within the velocity range of [−1505, 420] km s −1 that exhibits OH absorption.The spectra from each pixel are shown in Figure 6.The contours are the continuum as in Figure 1.

Figure 6 .
Figure6.Fitted OH 119 µm spectra (profile map) toward the central ≈ 1.7 × 2.3 kpc 2 region extracted from the pixels in Figure5.Each pixel has the area 0. ′′ 1 × 0. ′′ 1 (approximately one half of the beam size) with offset 0. ′′ 1.The spectra are denoted by the pixel coordinates (x, y) at the top right.The continuum peak is approximately between(180,179) and (180,180).The center velocities of the fitted emission and absorption lines are shown at the bottom right of each spectrum.The narrow emission feature at (181,179) is likely an artifact.

Figure 8 .Figure 9 .
Figure 8. Results of the SED fitting for J2054-0005 using the CIGALE code.The red line labeled "Dust emission" is the dust component heated by star formation, and the orange line labelled "AGN emission" is the contribution from the AGN.The latter makes up 59% ± 8% of the total infrared luminosity.The fit has χ 2 = 0.52 ± 0.11 (reduced).See the text for more details on the parameter setup.1000 2000 3000 4000 5000 SFR [M yr 1 ] of star formation.The total molecular gas mass in this quasar was estimated to be M mol = (3−6) × 10 10 M ⊙ by Decarli et al. (2022) from a variety of tracers (dust, CO, [C I] 609 µm, and [C II] 158 µm emission).

Figure 10 .
Figure10.Mass outflow rate vs. total IR luminosity (LIR).The ordinate on the left shows the product of measured quantities (equivalent width, outflow velocity, and radius) that is proportional to Ṁout in the expanding-shell model; Ṁout = 6.28 × 10 −3 W voutrout, where W is in km s −1 , vout is in km s −1 , and rout is in kpc.The black circles are dusty starforming galaxies (DSFGs) at z = 4−5(Spilker et al. 2020a).The data of J2310+1855 and P183+05 are taken fromButler et al. (2023), and W was calculated using Equation (4).The ordinate on the right is Ṁout calculated using Equations (3) and (5).The dotted line is η = 1 and the solid line is the fit.The open square symbol (not included in the fit) shows J2054-0005 corrected for the AGN contribution.The data points for J2310+1855 and P183+05 are not AGN-corrected.
Figure 11.W √ LIR and total infrared luminosity (LIR).The equivalent width W is the outflow component obtained from fitting the absorption feature.The solid line is the fitted relation.The open square symbol (not included in the fit) shows J2054-0005 corrected for the AGN contribution.The data points for J2310+1855 and P183+05 are not AGNcorrected.

Table 1 .
Image Parameters

Table 2 .
OH Line Fit Parameters

Table 3 .
Molecular Outflow PropertiesNote-The LTE values are calculated using the covering factor f = 0.3, outflow radius rout = 2 kpc, excitation temperature Tex = 50 K, and outflow velocity vout = 669 km s −1 .The empirical relation for the mass outflow rate is given in Equation (6).

Table 4 .
Fitting Results for Resolved Emission and Absorption ComponentsRegion v emi cen [km s −1 ] v abs cen [km s −1 ] FWHM emi [km s −1 ] FWHM abs [km s −1 ]Note-Regions are designated by image pixel numbers (x, y).Each pixel has area 0. ′′ 1 × 0. ′′ 1, which is approximately one half of the synthesized beam.The spectrum at (181,180) could not be fitted with emission.