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

The continuum emission (λ_rest = 123 μm), extracted from the LSB, is detected toward the quasar with a high signal-to-noise ratio of S/N = 260 (Figure 1). The emission is spatially resolved, although strongly concentrated in the center. The flux density in the region within a radius of 05 of the brightest pixel is S_ν (r < 05) = 5.723 ± 0.009 mJy. We estimated the peak coordinates by two-dimensional Gaussian fitting of a region of radius 05 centered at the brightest pixel. Using CASA's imfit, we obtained ${(\alpha ,\delta )}_{\mathrm{ICRS}}=({20}^{{\rm{h}}}{54}^{{\rm{m}}}06\buildrel{\rm{s}}\over{.} 501,-0^\circ 05^{\prime} 14\buildrel{\prime\prime}\over{.} 44)$. Since S/N is very high, the positional uncertainty is determined by the absolute astrometric accuracy of ALMA observations in Band 7, which is 5% of the synthesized beam size (≈10 mas). By comparison, the peak of the SDSS optical z-band image is at ${(\alpha ,\delta )}_{\mathrm{ICRS}}=({20}^{{\rm{h}}}{54}^{{\rm{m}}}06\buildrel{\rm{s}}\over{.} 486,-0^\circ 05^{\prime} 14\buildrel{\prime\prime}\over{.} 50)$ with an uncertainty of ≈470 mas, and the peak position of the [O iii] 88 μm is at ${(\alpha ,\delta )}_{\mathrm{ICRS}}=({20}^{{\rm{h}}}{54}^{{\rm{m}}}06\buildrel{\rm{s}}\over{.} 503,-0^\circ 05^{\prime} 14\buildrel{\prime\prime}\over{.} 48)$ with an uncertainty of ≈48 mas (Hashimoto et al. 2019b). The 123 μm dust continuum, ionized gas traced by [O iii] 88 μm, and z-band optical peak positions are thus in agreement within their respective uncertainties. The peak intensity obtained from the Gaussian fitting is 3.308 ± 0.014 mJy beam⁻¹, and the size (FWHM) of the central region where the emission is concentrated, deconvolved from the beam, is estimated to be $({d}_{\mathrm{maj}},{d}_{\min })=(0\buildrel{\prime\prime}\over{.} 1567\pm 0\buildrel{\prime\prime}\over{.} 0017,0\buildrel{\prime\prime}\over{.} 1321\pm 0\buildrel{\prime\prime}\over{.} 0022)$ at a position angle of 171°, corresponding to 890 ± 10 pc for the major axis.

Zoom In Zoom Out Reset image size

We also found an additional continuum source positioned 24 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 it is at the same redshift. The peak intensity is measured to be 116 ± 13 μJy beam⁻¹ at ${(\alpha ,\delta )}_{\mathrm{ICRS}}=({20}^{{\rm{h}}}{54}^{{\rm{m}}}06\buildrel{\rm{s}}\over{.} 344,-0^\circ 05^{\prime} 14\buildrel{\prime\prime}\over{.} 83)$, and the flux density is S_ν (r < 05) = 377 ± 9μJy.

Zoom In Zoom Out Reset image size

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 traced by dust. The OH profile is dominated by a broad absorption feature at negative velocities, and emission at near-systemic velocities, exhibiting a typical P-Cygni profile, such as the one observed toward the local ULIRG Mrk 231 (Fischer et al. 2010). The line is very broad: what appears to be either OH absorption or emission that extends over a continuous velocity from −1500 to +1000 km s⁻¹.

Zoom In Zoom Out Reset image size

We can estimate the optical depth τ as long as the line is not completely opaque. If emission at velocities of the absorption line is negligible, the observed flux density is S(v) = S_cont e^−τ , where S_cont is the continuum flux density; hence,

$$\begin{eqnarray}&&\tau (v)=-\mathrm{ln}\left[\displaystyle \frac{S(v)}{{S}_{\mathrm{cont}}}\right].\end{eqnarray} \tag{ 1 }$$

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.364 GHz, and ¹⁸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⁻¹. A similar feature attributed to the ¹⁸OH line is found in Mrk 231, indicating an enhanced [¹⁸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.

To investigate the kinematics of OH gas, we performed a least-squares fitting of the line profile using Python (scipy.optimize.curve_fit). 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⁻¹, 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⁻¹ near the systemic velocity, and a broad absorption feature with FWHM = 1052 ± 234 km s⁻¹, with the peak absorption velocity of v_cen = −669 ± 87 km s⁻¹ relative to the systemic velocity. The FWHM of the emission line obtained from fitting is comparable to those of the emission lines of [C ii] 158 μm (243 ± 10 km s⁻¹; Wang et al. 2013), [O iii] 88 μm (282 ± 17 km s⁻¹; Hashimoto et al. 2019b), and CO (J = 6 → 5; 360 ± 110 km s⁻¹; Wang et al. 2010). However, [O iii] 88 μm is a tracer of ionized gas, and it is unclear whether [C ii] 158 μm 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.

Table 2. OH Line Fit Parameters

	Absorption (Outflow)	Emission (Systemic)
Peak value Smax [mJy]	−1.039 ± 0.078	1.35 ± 0.27
Integrated flux ${{ \mathcal S }}_{\mathrm{OH}}$ [Jy km s−1]	−1.16 ± 0.27	0.44 ± 0.12
Center velocity vcen [km s−1]	−669 ± 87	65 ± 15
v84 [km s−1]	−1104 ± 35	...
v98 [km s−1]	−1574 ± 35	...
Standard deviation σv [km s−1]	446 ± 99	129 ± 23
FWHM line width [km s−1]	1052 ± 234	306 ± 55
Equivalent width W [km s−1]	200 ± 44	−66 ± 19

Note. All quantities are calculated from Gaussian fits. Here, $\mathrm{FWHM}=\sqrt{8\mathrm{ln}2}{\sigma }_{v}$, ${{ \mathcal S }}_{\mathrm{OH}}=\sqrt{2\pi }{\sigma }_{v}{S}_{\max }$, and W is calculated according to Equation (4). The quantities ${{ \mathcal S }}_{\mathrm{OH}}$ and W are given for a single line of the doublet (the total equivalent width is twice the value). The continuum flux density is S_cont = 6.175 ± 0.092 mJy.

Download table as:  ASCIITypeset image

The terminal velocity (maximum extent of the blueshifted wing of the absorption) is at least ${v}_{\max }\approx -1500\,\mathrm{km}\,{{\rm{s}}}^{-1}$ but may be beyond the spectral coverage. Another indicator of terminal velocity is v₉₈, the velocity above which 98% of absorption takes place. This quantity is found to be v₉₈ = −1574 ± 35 km s⁻¹. On the other hand, the velocity above which 84% of absorption takes place is v₈₄ = −1104 ± 35 km s⁻¹. The fitted emission components are redshifted relative to the systemic by 65 ± 15 km s⁻¹. 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⁻¹.

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 < 05, 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 the 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.

Zoom In Zoom Out Reset image size

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)

$$\begin{eqnarray}&&{M}_{\mathrm{out}}=\mu {m}_{{\rm{H}}}{N}_{{\rm{H}}}{\rm{\Omega }}{r}_{\mathrm{out}}^{2},\end{eqnarray} \tag{ 2 }$$

where μ = 1.36 is the mean particle mass per hydrogen nucleus, m_H is the hydrogen atom mass (in kilograms), ${N}_{{\rm{H}}}=2{N}_{{{\rm{H}}}_{2}}$ is the column density of hydrogen nuclei (in per squared meter), and Ω is the solid angle ("opening angle") subtended by the shell (Veilleux et al. 2020). Thus, the term ${\rm{\Omega }}{r}_{\mathrm{out}}^{2}$ 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}_{{{\rm{H}}}_{2}}$ cannot be measured directly, but we can find the column density of OH molecules (N_OH) and then apply a [OH]/[H₂] abundance ratio to get ${N}_{{{\rm{H}}}_{2}}$. 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)

$$\begin{eqnarray}&&{N}_{\mathrm{OH}}=\displaystyle \frac{8\pi }{{\lambda }_{{ul}}^{3}{A}_{{ul}}}\displaystyle \frac{Q}{{g}_{u}}\displaystyle \frac{{e}^{{E}_{l}/{{kT}}_{\mathrm{ex}}}}{1-{e}^{-{\rm{\Delta }}E/{{kT}}_{\mathrm{ex}}}}\int \tau (v){dv},\end{eqnarray} \tag{ 3 }$$

where λ_ul = 119.23 μm is the rest-frame wavelength, A_ul = 0.1388 s⁻¹ 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. 2019b), we get Q ≃ 19 (Pickett et al. 1998). The term $Q{(1-{e}^{-{\rm{\Delta }}E/{{kT}}_{\mathrm{ex}}})}^{-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 (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₂] = 5 × 10⁻⁶ is reported for Sgr B2 (Goicoechea & Cernicharo 2002). Although this value is often used to derive H₂ mass from OH, it may be lower at subsolar 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. 2020a; Onoue et al. 2020). On the other hand, a low abundance of [OH]/[H₂] ≈ 1 × 10⁻⁷ has been reported recently even for various regions in the Galaxy (Nguyen et al. 2018; Rugel et al. 2018). We derive the outflow properties using [OH]/[H₂] = 5 × 10⁻⁷ (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.

Table 3. Molecular Outflow Properties

Quantity	LTE	Empirical Relation
Column density NOH [cm−2]	7.5 × 1015	...
OH abundance [OH]/[H2]	5 × 10−7	...
Column density ${N}_{{{\rm{H}}}_{2}}\,[{\mathrm{cm}}^{-2}]$	1.5 × 1022	...
Mass Mout [M⊙]	4.9 × 109	4.5 × 109
Mass outflow rate ${\dot{M}}_{\mathrm{out}}$ [M⊙ yr−1]	1700	1500
Mass loading factor η	2.2	1.9
Depletion time tdep [yr]	(2–4) × 107	(2–4) × 107
Kinetic energy Eout [erg]	2.2 × 1058	2.0 × 1058
Power ${\dot{E}}_{\mathrm{out}}\,[{L}_{\odot }]$	6.2 × 1010	5.7 × 1010

Note. The LTE values are calculated using the covering factor f = 0.3, outflow radius r_out = 2 kpc, excitation temperature T_ex = 50 K, and outflow velocity v_out = 669 km s⁻¹. The empirical relation for the mass outflow rate is given in Equation (6).

Download table as:  ASCIITypeset image

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. (2020a) 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),

$$\begin{eqnarray}&&W=-\int \mathrm{ln}\left[\displaystyle \frac{S(v)}{{S}_{\mathrm{cont}}}\right]{dv},\end{eqnarray} \tag{ 4 }$$

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⁻¹, yielding M_out ≈ 4.9 × 10⁹ M_⊙. The obtained W is larger than that found in almost all nearby ULIRGs (Veilleux et al. 2013) and DSFGs at z = 4–5 (Spilker et al. 2020b), 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).

4.2.1. 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

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{out}}={M}_{\mathrm{out}}\displaystyle \frac{{v}_{\mathrm{out}}}{{r}_{\mathrm{out}}}.\end{eqnarray} \tag{ 5 }$$

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⁻¹, 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 spherically symmetric outflow. Taking W = 200 km s⁻¹ (equivalent width of the absorption line), [OH]/[H₂] = 5 × 10⁻⁶, r_out = 2 kpc, and f = 0.3, we obtain a lower limit of ${\dot{M}}_{\mathrm{out}}\gt 168\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. Note that, if we adopt a low abundance of [OH]/[H₂] = 1 × 10⁻⁷, and f = 1, the upper limit of the mass outflow rate becomes ${\dot{M}}_{\mathrm{out}}\approx 2.8\times {10}^{4}\,{M}_{\odot }\,{\mathrm{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 moderate abundance of [OH]/[H₂] = 5 × 10⁻⁷, and f = 0.3 in the analysis below.

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⁶ yr. The timescale is much shorter than the depletion time due to star formation (∼4–8 × 10⁷ yr).

4.2.2. 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. (2020) modified by Spilker et al. (2020a) takes the form

$$\begin{eqnarray}&&\left(\displaystyle \frac{{\dot{M}}_{\mathrm{out}}^{\mathrm{emp}}}{{M}_{\odot }\,{\mathrm{yr}}^{-1}}\right)=1.4\left(\displaystyle \frac{{W}_{v\lt -200}}{\mathrm{km}\,{{\rm{s}}}^{-1}}\right){\left(\displaystyle \frac{{L}_{\mathrm{IR}}}{{10}^{12}{L}_{\odot }}\right)}^{1/2}+180,\end{eqnarray} \tag{ 6 }$$

where W_v<−200 is the OH 119 μm equivalent width at v < −200 km s⁻¹. This equivalent width, derived from the observed spectrum, is W_v<−200 ≈ 268 km s⁻¹, yielding a mass outflow rate of ${\dot{M}}_{\mathrm{out}}^{\mathrm{emp}}\approx 1500\,{M}_{\odot }$. The value is significantly larger compared to those for the star-forming galaxies at redshift z = 4–5 reported in Spilker et al. (2020a), which have values between 220 and 1290 M_⊙ yr⁻¹.

Using the mass outflow rate from Equation (6), we calculate the outflow mass ${M}_{\mathrm{out}}^{\mathrm{emp}}=({r}_{\mathrm{out}}/{v}_{\mathrm{out}}){\dot{M}}_{\mathrm{out}}^{\mathrm{emp}}$ 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⁻¹ (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. Nonetheless, 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).

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 01 × 01, 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, labeled (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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Fitting Results for Resolved Emission and Absorption Components

Region	${v}_{\mathrm{cen}}^{\mathrm{emi}}$	${v}_{\mathrm{cen}}^{\mathrm{abs}}$	FWHMemi	FWHMabs
	(km s−1)	(km s−1)	(km s−1)	(km s−1)
179, 178	98 ± 19	−827 ± 147	202 ± 56	1000 ± 557
179, 179	101 ± 15	−665 ± 108	295 ± 62	1368 ± 338
179, 180	61 ± 16	−579 ± 83	278 ± 67	1057 ± 212
179, 181	41 ± 19	−512 ± 148	282 ± 80	849 ± 350
180, 178	75 ± 22	−814 ± 75	230 ± 64	842 ± 269
180, 179	102 ± 15	−588 ± 85	327 ± 60	1224 ± 223
180, 180	74 ± 19	−551 ± 89	330 ± 82	1095 ± 206
180, 181	−11 ± 12	−785 ± 61	96 ± 31	1084 ± 243
181, 178	44 ± 28	−778 ± 132	240 ± 93	963 ± 436
181, 179	...	−747 ± 36	...	933 ± 155
181, 180	...	−559 ± 50, −1296 ± 75	...	621 ± 125, 468 ± 134

Note. Regions are designated by image pixel numbers (x, y). Each pixel has area 01 × 01, which is approximately one-half of the synthesized beam. The spectrum at (181, 180) could not be fitted with emission.

Download table as:  ASCIITypeset image

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⁻¹). 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) kms⁻¹, 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⁻¹), consistent with reported [C ii] 158 μm, and [O iii] 88 μm, and CO (J = 6 → 5) line widths (Wang et al. 2010, 2013; Hashimoto et al. 2019b). 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 those in the south. The exception is at (181, 178); though, the line is only marginally detected there. The mean emission velocities in each row in Figure 6 from north to south are (15 ± 11, 67 ± 12, 101 ± 11, 72 ± 13) kms⁻¹.

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, which is generally consistent with the north–south 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 preparation), extracted from the central pixel (maximum 123 μm continuum intensity; pixel size 0034). 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⁻¹), and there appears to be a blueshifted component at velocities comparable to the OH outflow velocity (−600 km s⁻¹), 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. 2020b).

Zoom In Zoom Out Reset image size

More details of the [C ii] observations and results will be presented in S. Fujimoto et al. (2023, in preparation).

Estimating the SFR in high-z sources is often done by applying a conversion factor to the far-IR luminosity. 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 multiwavelength 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) analysis, 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 χ²) 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\lt {U}_{\min }\lt 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 χ² = 0.52 ± 0.11. We found that the total IR luminosity is L_IR = (1.34 ± 0.17) × 10¹³ L_☉, consistent with previous studies (Hashimoto et al. 2019b), 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⁻¹, calculated using $\mathrm{SFR}/{M}_{\odot }\,{\mathrm{yr}}^{-1}=1.40\times {10}^{-10}{L}_{\mathrm{IR}}^{{\prime} }/{L}_{\odot }$, where ${L}_{\mathrm{IR}}^{{\prime} }$ 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.

Zoom In Zoom Out Reset image size

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⁻¹) of Lyα emission reported in Jiang et al. (2008).

Is the outflow driven by star formation, or does the AGN feedback from the accretion onto the 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

$$\begin{eqnarray}&&{E}_{\mathrm{out}}=\displaystyle \frac{1}{2}{M}_{\mathrm{out}}{v}_{\mathrm{out}}^{2}\approx 2.2\times {10}^{58}\mathrm{erg},\end{eqnarray} \tag{ 7 }$$

and the power required to drive the molecular outflow (kinetic power) is

$$\begin{eqnarray}&&{\dot{E}}_{\mathrm{out}}=\displaystyle \frac{1}{2}{\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}^{2}\approx 2.4\times {10}^{44}\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 8 }$$

This is equivalent to ≈6.2 × 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

$$\begin{eqnarray}&&{p}_{\mathrm{out}}={M}_{\mathrm{out}}{v}_{\mathrm{out}}\approx 3.3\times {10}^{12}\,{M}_{\odot }\,\mathrm{km}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 9 }$$

By comparison, the momentum injection by a typical core-collapse supernova (SN; mass m₀ ≈ 10 M_⊙, velocity v₀ ≈ 3000 km s⁻¹) is of the order of p₀ ≈ 3 × 10⁴ M_⊙ km s⁻¹, and the total momentum of an outflow driven by SN explosions is ${p}_{\mathrm{SN}}\approx {p}_{0}{R}_{\mathrm{SN}}{t}_{\mathrm{SN}}$, where R_SN is the SN rate, and t_SN is the time interval measured from the onset of SN explosions. Thus, ${R}_{\mathrm{SN}}{t}_{\mathrm{SN}}$ is the total number of SN explosions over the starburst episode. Adopting a relation between the SN rate and star formation rate, ${R}_{\mathrm{SN}}/{\mathrm{yr}}^{-1}\approx {\alpha }_{\mathrm{SN}}(\mathrm{SFR}/{M}_{\odot }\,{\mathrm{yr}}^{-1})$, where ${\alpha }_{\mathrm{SN}}\approx 0.01\mbox{-}0.02$ depends on IMF and assumes continuous star formation (Veilleux et al. 2020), and SFR = 770 M_⊙ yr⁻¹; we obtain ${R}_{\mathrm{SN}}\approx 15\,{\mathrm{yr}}^{-1}$ (for α = 0.02). If the time interval of the SN feedback is equal to the dynamical age of the outflow (${t}_{\mathrm{SN}}={t}_{\mathrm{out}}$), the total SN momentum becomes ${p}_{\mathrm{SN}}\approx 1.3\times {10}^{12}\,{M}_{\odot }\,\mathrm{km}\,{{\rm{s}}}^{-1}$, produced by ${R}_{\mathrm{SN}}{t}_{\mathrm{SN}}\approx 4.5\times {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 2 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⁻¹ and the terminal velocity of 1500 km s⁻¹ exceed the outflow velocities typically measured in star-forming galaxies, which are found to be 100–500 and <1000 km s⁻¹, respectively (e.g., Sugahara et al. 2019), although the outflow velocities are found to be higher in more extreme systems (Spilker et al. 2020b). By contrast, the terminal velocities of outflows in AGN-dominated systems are often found to ≳1000 km s⁻¹ and as large as 1500 km s⁻¹ (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⁻¹ and mass outflow rates of ∼10³ M_☉ yr⁻¹ in AGN-driven outflows (e.g., Costa et al. 2018; Ishibashi et al. 2018). The results indicate that the AGN feedback (radiation pressure) may play a role in boosting the velocity in J2054-0005.

In Figure 9, we show the mean line-of-sight outflow velocity plotted against the total IR luminosity (L_IR) for 8 DSFGs at z = 4–5 and 3 quasars at z > 6. Here, L_IR is obtained by integrating the flux over λ_rest = 8–1000 μm (Section 5.1; Shao et al. 2019; Spilker et al. 2020b), except P183+05, for which L_IR = 1.41L_FIR (Venemans et al. 2020). For J2054-0005, we used L_IR derived from the CIGALE fitting. The plot also shows the star formation rate, calculated using SFR/M_☉ yr⁻¹ = 1.40 × 10⁻¹⁰ L_IR/L_☉, where a Chabrier (2003) IMF is assumed. The outflow in J2310+1855 was also detected in OH⁺ (1₁ ← 0₁) absorption (Shao et al. 2022), and the absorption velocity is comparable to that of OH plotted here.

Zoom In Zoom Out Reset image size

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.

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

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{\dot{M}}_{\mathrm{out}}}{\mathrm{SFR}}\sim 2,\end{eqnarray} \tag{ 10 }$$

implying efficient suppression of star formation. The total molecular gas mass in this quasar was estimated to be M_mol = (3–6) × 10¹⁰ M_☉ by Decarli et al. (2022) from a variety of tracers (dust, CO, [C i] 609 μm, and [C ii] 158 μm emission). 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

$$\begin{eqnarray}&&{t}_{\mathrm{dep}}=\displaystyle \frac{{M}_{\mathrm{mol}}}{{\dot{M}}_{\mathrm{out}}}\approx (2\mbox{-}4)\times {10}^{7}\,\mathrm{yr}.\end{eqnarray} \tag{ 11 }$$

The relatively short timescale, as compared to the time required for star formation to consume molecular gas in ordinary star-forming galaxies (∼10⁹ 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⁹ 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 a 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.

Figure 10 shows a comparison of mass outflow rates and the total IR luminosity (L_IR) in high-z quasars and DSFGs with OH outflow detections. Since the uncertainty of the mass outflow rate is dominated by the OH abundance, we plot the product Wv_out r_out in addition to ${\dot{M}}_{\mathrm{out}}$, because it is the directly measured quantity that is proportional to ${\dot{M}}_{\mathrm{out}}$ in the expanding-shell model under LTE, whereas ${\dot{M}}_{\mathrm{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. (2020b). W was calculated using Equation (4), and we assumed T_ex = 50 K, and [OH]/[H₂] = 5 × 10⁻⁷ for all sources in calculating ${\dot{M}}_{\mathrm{out}}$. The plot shows that there is a positive relation between the mass outflow rate with L_IR, and the power law is $\mathrm{log}y=k\mathrm{log}x+m$, where k = 1.61 ± 0.40, and m = −16.4 ± 5.2, although the scatter is large (R² = 0.64).

Zoom In Zoom Out Reset image size

As discussed in Herrera-Camus et al. (2020), Spilker et al. (2020a), the correlation can also be found if the mass outflow rate is set to be proportional to ${\dot{M}}_{\mathrm{out}}\propto W\sqrt{{L}_{\mathrm{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 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² = 0.52).

Zoom In Zoom Out Reset image size

The above analysis yields a positive relation between ${\dot{M}}_{\mathrm{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 ${\dot{M}}_{\mathrm{out}}\mbox{-}{L}_{\mathrm{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 ${\dot{M}}_{\mathrm{out}}\mbox{-}{L}_{\mathrm{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.

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¹⁰ 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}_{\mathrm{esc}}(R)\approx \sqrt{2{{GM}}_{\mathrm{dyn}}/R}\approx 780\,\mathrm{km}\,{{\rm{s}}}^{-1}$ at R. The estimate corresponds to a rotational velocity of 560 km s⁻¹, implying a very massive host galaxy. The escape velocity is similar to the velocities reported for DSFGs at z = 4–5 (Spilker et al. 2020a). 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 recent observations of enriched gas in the circumgalactic medium at z ∼ 6 (Wu et al. 2021).

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; Hashimoto et al. 2019a, 2019b; Laporte et al. 2019; 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, the 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. 2019b; 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⁻¹ in J2054-0005 compared to 334 ± 14 km s⁻¹ in J2310+1855). On the other hand, J2310+1855 has slightly larger total IR luminosity (1.9 × 10¹³ L_⊙) compared to J2054-0005 (1.3 × 10¹³ 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.

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 DSFGs at z = 4–5 (Spilker et al. 2020b).

Based on multiline OH observations, Spinoglio et al. (2005) argue that collisional excitation dominates the population of the ²Π_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 dust-obscured central regions and that, except in two outliers, radiative excitation may be dominant. Since the peak of the SED of dust thermal emission in J2054-0005 is close to λ_rest = 53 μm, the wavelength that corresponds to the energy difference between the levels ²Π_1/2 J = 3/2, and ²Π_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 ²Π_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 multiline observations are necessary to constrain the excitation mechanism of OH molecules in EoR quasars.