Galactic Winds across the Gas-Rich Merger Sequence II. Lyman Alpha Emission and Highly Ionized O VI and N V Outflows in Ultraluminous Infrared Galaxies

This paper is the second in a series aimed at examining the gaseous environments of z$\le$0.3 quasars and ultraluminous infrared galaxies (ULIRGs) as a function of AGN/host galaxy properties across the merger sequence. This second paper focuses on the Ly$\alpha$ emission and O VI and N V absorption features, tracers of highly ionized gas outflows, in ULIRGs observed with HST/COS. Ly$\alpha$ emission is detected in 15 out of 19 ULIRGs, and 12 of the 14 clear Ly$\alpha$ detections show emission with blueshifted velocity centroids and/or wings. The equivalent widths of the Ly$\alpha$ emission increase with increasing AGN luminosities and AGN bolometric fractions. The blueshifts of the Ly$\alpha$ emission correlate positively with those of the [O III] emission, where the latter traces the ionized gas outflows. The Ly$\alpha$ escape fractions tend to be slightly larger in objects with stronger AGN and larger outflow velocities, but they do not correlate with nebular line reddening. Among the 12 ULIRGs with good continuum signal-to-noise ratios, O VI and/or N V absorption features are robustly detected in 6 of them, all of which are blueshifted, indicative of outflows. In the combined ULIRG $+$ quasar sample, the outflows are more frequently detected in the X-ray weak or absorbed sources. The absorption equivalent widths, velocities and velocity dispersions of the outflows are also higher in the X-ray weak sources. No other strong correlations are visible between the properties of the outflows and those of the AGN or host galaxies.


INTRODUCTION
Major mergers of gas-rich galaxies, both near and far, are the paradise for magnificent starbursts and rapid growth of supermassive black holes. In the local universe, the majority of the ultraluminous infrared galaxies (ULIRGs) are mergers of gas-rich galaxies. The merger process drives the gas and dust to the central region of the system, fueling the (circum)nuclear starbursts and the rapid accretion of the supermassive black tic medium (e.g., Martin 2005;Rupke et al. 2005;Martin 2006;Tremonti et al. 2007;Martin & Bouché 2009;Sturm et al. 2011;Rupke & Veilleux 2013a;Veilleux et al. 2013a,b;Cicone et al. 2014;Veilleux et al. 2014;Tombesi et al. 2015;Rupke et al. 2017;Liu et al. 2019;Fluetsch et al. 2019Fluetsch et al. , 2020Lutz et al. 2020;Veilleux et al. 2020, for a review).
While the cooler, neutral and/or molecular phases on larger scale ( kpc) often dominate the outflow energetics, it is the hotter, ionized phase of the wind that serves as the best probe for the driving mechanism of these winds. ULIRG F11119+3257, arguably the best example so far, possesses a massive, galactic scale (1-10 kpc) molecular and neutral-gas outflow apparently driven by the fast (>0.1 c), highly ionized (Fe XXV and Fe XXVI at ∼7 keV) nuclear wind (Tombesi et al. 2015Veilleux et al. 2017). While this result is intriguing, the faintness of the majority of ULIRGs at ∼ 7 keV, unlike many quasars, has impeded a statistically meaningful study of this phenomenon in most ULIRGs with current X-ray facilities.
The superb far-ultraviolet (FUV) sensitivity of the Cosmic Origins Spectrograph (COS) on the Hubble Space Telescope (HST) provides a powerful alternative tool for such study in the low-z universe. Rest-frame FUV spectroscopy has enabled a comprehensive study of the multi-phase nature of outflows, built upon the abundant spectral features arising from the high-ionization, low-ionization, and neutral phases of the outflowing gas (e.g., Chisholm et al. 2015;Tripp et al. 2011;Heckman et al. 2015;Hamann et al. 2019;Arav et al. 2020). Up to now, only about a dozen ULIRGs have been studied with HST/COS data, but the results are fascinating. In Mrk 231, highly blueshifted Lyα emission (with respect to systemic velocity) is observed to coincide in velocity with the highly blueshifted absorption features tracing the fast outflow in this galaxy, suggesting an outflow-related origin for the Lyα emission (Veilleux et al. 2013a(Veilleux et al. , 2016. With a larger sample of 11 ULIRGs, Martin et al. (2015, hereafter M15) have shown that prominent, blueshifted Lyα emission down to −1000 km s −1 exists in about half of the objects, and they argued that the blueshifted Lyα emission originates from the clumps of gas condensing out of hot winds driven by the central starbursts (Thompson et al. 2016). In addition, blueshifted absorption features from high-ionization species like O 5+ and/or N 4+ (114 and 77 eV are needed to produce these ions, respectively) and low-ionization species like Si + and Fe + are also detected in a few objects, providing unambiguous evidence of outflowing gas.
Despite the tantalizing evidence of FUV-detected outflows in the ULIRGs described above, the sample ex-amined so far is small and incomplete, where AGNdominated ULIRGs and matched quasars are largely missing. To address this issue, we have selected a more complete sample of ULIRGs and quasars to systematically study the gaseous environments along the merger sequence, from ULIRGs to quasars. In Veilleux et al. (2022) (hereafter Paper I), we presented the results from the first part of our study, focused on the highly-ionized gas outflows, traced by O VI 1032, 1038 and N V 1238, 1243 absorption features, in a sample of 33 local quasars. We found that the O VI and N V outflows are present in ∼61% of the sample, and the incidence rate and equivalent widths (EWs) of these highly ionized outflows are higher among X-ray weak or absorbed sources. Similarly, the flux-weighted outflow velocity dispersions are also the highest among the X-ray weakest sources. However, no significant correlation is visible between the flux-weighted outflow velocities/velocity dispersions and the other properties of the quasars and host galaxies.
In this paper, we report the results from an analysis of the Lyα emission and O VI and N V absorption features of the 21 ULIRGs in the sample 1 , expanding on the results from Paper I by considering the combined ULIRG + quasar sample. In Sec. 2, we describe the HST/COS observations of the ULIRG sample, the reduction of the data sets, and the ancillary data from the literature. In Sec. 3, we present the analysis of the FUV spectra of the ULIRGs, focusing on Lyα emission in the first part and the O VI and N V absorption features in the second. In Sec. 4, we discuss the potential key factors that control the observed Lyα properties, and in Sec. 5, we examine the incidence rates and properties of the O VI and N V outflows. In Sec. 6, we search for trends between the O VI and/or N V outflow properties and the AGN/galaxy properties in the ULIRG+quasar sample. In Sec. 7, we summarize the main results of this paper. Throughout the paper, we adopt a ΛCDM cosmology with H 0 = 75 km s −1 Mpc −1 , Ω m = 0.3, and Ω Λ = 0.7.

HST/COS G130M Observations of ULIRGs
Our sample is selected based on three criteria: (1) They are part of the 1-Jy sample of 118 local ULIRGs with z < 0.3 and |b| > 30 • (hence modest Galactic extinctions; Kim & Sanders 1998b); (2) In order to address the role of AGN feedback in these systems, the ULIRGs

Note-Column
(1): Name of the object; Column (2): Short names of the objects adopted in this paper; Column (3): Redshift. The label on the upper-right corner of each entry indicates the reference of the redshift. a: Based on the narrow optical emission lines from Martin et al. (2015). b: The best redshift adopted by NASA/IPAC Extragalactic Database a . c: Based on the narrow optical emission lines from Rupke et al. (2017). d: Based on the narrow optical emission lines From Rupke & Veilleux (2013b). e: Based on the optical emission lines from Kim & Sanders (1998a); Column (4): Log of bolometric luminosity in solar units based on. For ULIRGs, we assumed L bol = 1.15 L IR where L IR is the 8-1000 µm infrared luminosity from Kim & Sanders (1998a). For the QUEST quasars, we assumed L bol = 7 L(5100Å) + L IR based on Netzer et al. (2007); Column (5): Optical spectral type from Veilleux et al. (1999): S1 means Seyfert 1, S2 means Seyfert 2, L means LINER, HII means star-forming. Column (6): Interaction class from Veilleux+02: I-First approach, II-First contact, III(a/b)-Pre-merger (Wide binary/Close binary), IV(a/b)-Merger (Diffuse/Compact), V-Old Merger, Iso-Isolated, Tpl-triplet. Column (7): Fraction (in percentage) of the bolometric luminosity produced by the AGN, based on the mean values derived in Veilleux et al. (2009b). The typical uncertainties are 10-15%. Column (8): FUV AB magnitudes from GALEX. Column (9): Soft (0.5-2 keV) X-ray luminosity. For X-ray related quantities (Column (9), (10), (11) and (13)), different rows = different observations dates. The references for these quantities are described in Sec. 2.3.2. The errors are omitted whenever they are not available from the literature or public data; Column (10) Hard (2-10 keV) X-ray luminosity; Column (11): The X-ray absorbing column density in units of 10 22 cm −2 ; Column (12): Monochromatic continuum luminosity at the rest-frame 1125Å; Column (13): X-ray-to-UV index, α U V,X , as defined in Sec. 2.3.2.
(3) They are the FUV-brightest ULIRGs of the 1-Jy sample with FUV magnitudes (AB) m F U V 21. These criteria result in a sample of 21 objects (Table 1): 15 of which were observed through the HST cycle 26 program (PID:15662; PI: Sylvain Veilleux), and the remaining 6 objects have archival COS G130M spectra of sufficient quality from three programs (PID: 12533, PI: C. Martin; PID: 12569, PI: S. Veilleux; PID: 12038: PI: J. Green). Among these 6 objects, QSO-B0157+001, 3C 273, and Mrk 231 were also studied in Paper I as they meet the criteria for QUEST quasars, and F01004−2237, Z11598−0112, and F12072−0444 were also studied in M15. In the following sections, we adopt the short names listed in Table 1 when referring to the objects in our sample.
The Cycle 26 HST/COS spectra presented in this paper were obtained in TIME-TAG mode through the PSA using the medium resolution FUV grating, G130M. Four focal plane offset positions were adopted to reduce the impact of fixed-pattern noise associated with the microchannel plate. We got all four FP-POS settings for all targets, except for targets F04103, F14070, F21219, and F23233 with a central wavelength of 1291Å. For these objects, we followed the COS2025 recommendations and used FP-POS = 3 and 4 to get equal exposures for segments A and B. The wavelength setting was adjusted according to the redshift of the target and was selected to optimize the number of strong lines that can be observed with G130M. At all but the highest redshifts (z < 0.2), Lyα, the high-ionization lines from N V, and low-ionization lines from Si II 1260 and Si III 1206 fit within the wavelength coverage of the data. At the highest redshifts (z 0.20-0.27), we lose N V and Lyα but gain Lyβ and O VI 1032, 1038. The wavelength coverage for individual observations is summarized in Table  2.

HST/COS Data Reduction
Among the 6 objects with archival data, we retrieved the fully reduced spectra for 5 of them from the Hubble Legacy Spectrum Archive (Peeples et al. 2017), and obtained the fully reduced spectrum of F12072−0444 published in M15 from C. L. Martin. For the other 15 newly-observed objects presented in this paper, the raw data were processed and combined by the CALCOS pipeline v3.3.10. CALCOS corrects the data for instrumental effects, assigns a vacuum wavelength scale, and extracts flux-calibrated spectra. It applies a heliocentric correction to the final x1d files for each exposure, and combines the individual exposures to a single spectra when possible.
The COS aperture is filled with emission from geocoronal airglow, so the observed wavelengths of these lines are independent of target position in the PSA. By inspecting the velocity offsets of theoretical and observed wavelength of strong geocoronal lines, we can examine the potential systematic errors in the wavelength calibration. The measured velocity offsets are randomly distributed with absolute values < 30 km s −1 , and we take these measurements as typical errors of the wavelength calibration from the pipeline.
Finally, all spectra are corrected for foreground Galactic extinctions from Schlafly & Finkbeiner (2011) and the reddening curve with R V =3.1 of Fitzpatrick (1999).

General Physical Properties
For the bolometric luminosities of our sources, we adopt L bol = 1.15 L IR where L IR is the 8-1000 µm infrared luminosity retrieved from Kim & Sanders (1998b), except for the three sources that are also QUEST quasars, which are quoted from Paper I where we assume L bol = 7 L(5100Å) + L IR based on Netzer et al. (2007). The L(5100Å) is the continuum luminosity λL λ at 5100Å rest wavelength and L IR is the 1 -1000 µm infrared luminosity (the details can be found in the notes of Table 1 in Paper I).
The optical spectral classifications are quoted from Veilleux et al. (1999): S1 means Seyfert 1, S2 means Seyfert 2, L means LINER, HII means star-forming galaxies. The interaction classes (or merger classes) are from Veilleux et al. (2002): I-First approach, II-First contact, III(a/b)-Pre-merger (Wide binary/Close binary), IV(a/b)-Merger (Diffuse/Compact), V-Old Merger, Iso-Isolated, Tpl-triplet. The fraction of the bolometric luminosity produced by the AGN, or simply AGN fractions, are the average values derived in Veilleux et al. (2009a). The AGN luminosities are defined as the bolometric luminosities multiplied by the AGN fractions.

X-ray Data
Published X-ray data and measurements exist for 14 out of the 21 objects: Those for the 3 quasars (QSO-B0157+001, 3C 273, Mrk 231) also studied in Paper I are from Teng & Veilleux (2010); Teng et al. (2014); Veilleux et al. (2014); Ricci et al. (2017) Those for the remaining 11 sources are from a series of X-ray studies of ULIRGs and quasars (Teng et al. 2005;Teng & Veilleux 2010). Following Paper I, the X-ray weakness of AGN/quasars can be described with the X-ray to optical spectral index, α OX ≡ 0.372log[(F (2 keV)/F (3000Å)] (e.g., Brandt et al. 2000). While α OX is measured for most nearby quasars, there are virtually no published measurements for the ULIRGs in our sample. Instead, we have defined an alternative X-ray to FUV spectral index, α U V,X , based on the ratio of the soft X-ray (0.5-2 keV) flux to the FUV flux from GALEX These results are listed in Table 1.
For the QUEST quasars studied in Paper I, there is a clear positive correlation between α U V,X and α OX (see Fig. 1), with a p-value 2×10 −5 from the Kendall tau test (the null hypothesis is no correlation), which demonstrates that α U V,X is indeed a good surrogate for α OX . For the quasars without published 0.5-2 keV flux from Chandra or XMM-Newton and/or FUV flux from GALEX, we convert their α OX listed in Paper I to α U V,X by adopting the relation α U V,X = 2.17α OX +2.26 , which is obtained from a linear fit to the quasars with both α U V,X and α OX measurements.

Optical Spectra
The Gemini/GMOS IFU spectra from Rupke & Veilleux (2013b), Rupke et al. (2017), or the SDSS spec-2.0 1.5 OX 2 1 0 UV, X Figure 1. The X-ray to optical spectral index αOX versus the X-ray to FUV spectral index αUV,X for the quasars in Paper I. The two indices are defined as αOX ≡ 0.372log[(F (2 keV)/F (3000Å)] and αUV,X ≡ log[F (0.5 − 2 keV)/F (F U V )], respectively. The solid line is a linear fit to the data points. The errors on αOX and αUV,X are uncertain, and largely associated with the uncertainties in the analyses of the X-ray spectra, as described in, e.g., Teng & Veilleux (2010). The cross in the upper left corner of the figure indicates ± 0.1 dex errors adopted in the Kendall tau test and the fit.
tra (Eisenstein et al. 2011) are adopted as the optical spectra for our objects by default unless explicitly stated otherwise (the GMOS data are adopted by default whenever both GMOS and SDSS data are available for the same object). The long-slit, optical spectra of all objects but 3C 273 are retrieved from Veilleux et al. (1999), and the optical spectrum of 3C 273 is retrieved from Buttiglione et al. (2009). They are adopted as the default optical spectra whenever the Gemini/GMOS and SDSS spectra are not available.
For these spectra, the continua are modeled with either stellar population synthesis (SPS) models (González Delgado et al. 2005) adopting pPXF (Cappellari 2017) or 4th-order polynomials and/or power-law function with customized Python software utilizing LM-FIT ), on a case-by-case basis. The properties of the [O iii] λ5007 and Hα emission lines are then measured from these continuum-subtracted optical spectra.

AGN Fractions of Starburst-dominated ULIRGs
Most of the key AGN and host galaxy properties of the QUEST quasars from Paper I and the starburstdominated ULIRGs from M15 are tabulated in these papers. One exception is the AGN fractions of the ULIRGs in M15, which are estimated based on the IRAS-based 25 µm to 60 µm flux ratios (F25/F60) listed in Table  1 of M15. Specifically, the AGN fraction is calculated adopting the best-fit to the trend presented in Fig. 36(c) in Veilleux et al. (2009a): (1) where x = log(F 25/F 60).

RESULTS FROM THE HST DATA ANALYSIS
In this section, we present the major results from our analysis of the HST/COS spectra. First, the properties regarding Lyα emission are examined in Sec. 3.1; Second, the measurement of FUV continuum luminosity is briefly described in Sec. 3.2; Finally, the properties of O VI 1032, 1038 and N V 1238, 1243 absorbers are discussed in Sec. 3.3.

Detection Rates
The Lyα transition falls within the wavelength range of the observations for 19 of the 21 objects. Among the 19 objects, the Lyα emission is detected (S/N > 3) in 15 of them, including F07599 where the Lyα emission is heavily affected by deep, broad and narrow N V 1238, 1243 absorption features (and perhaps also broad Lyα absorption). After Lyα, N V 1238, 1243 and O VI 1032, 1038 are the most frequently detected emission lines in our sample (notice that our G130M spectra do not cover Si IV and C IV). Descriptions about the presence of emission/absorption features other than Lyα in our objects may be found in Appendix A.

Line Profiles
The Lyα profiles of the 15 Lyα detections are presented in Fig. 2. For these objects, the Lyα emission is the most prominent feature in the observed spectral range, except for F07599, where the Lyα line is severely suppressed by a highly blueshifted broad absorption line (BAL) and several less blueshifted narrow N V 1238, 1243 absorption features. As a result, a relatively robust measurement of the Lyα profile of F07599 is impossible, and this source is left out from the analyses to characterize the Lyα profiles in the following sections.
The flux and EWs of the entire Lyα profiles are measured in wavelength windows customized for each object based on their line widths. The local continuum of each object is determined by fitting the line-free windows adjacent to the Lyα features with power-law or loworder polynomials (order ≤ 2). The contamination from nearby N V emission is subtracted for sources Mrk 1014, F05189, 3C 273, and F21219, where the N V doublet are modeled as Gaussian profiles. The foreground absorption features and absorption lines from other species at the systemic velocity are interpolated over with cubic splines. For the non-detections, the 3-σ upper limits on the flux of Lyα are estimated in a velocity window of −1000 to +1000 km s −1 .
We adopt a non-parametric approach to characterize the Lyα profiles quantitatively: the velocities v 50 and v 80 (velocities at 50 and 80 percentiles of the total flux calculated from the red side of the line), the line width W 80 (which encloses the central 80 percent of the total flux), and the line asymmetry A 91 ([v 90 +v 10 ]/W 80 ), where v 90 and v 10 are the velocities at 90 and 10 percentiles of the total flux calculated from the red side of the line) are the primary measurements adopted for our analyses in the following sections. All velocities are calculated with respect to the systemic redshifts listed in Table 1. Table 3 summarizes the results. A visualization of these non-parametric measurements is shown in Fig. 3.
Blueshifted wings of the Lyα emission are often seen across our sample (12 out of the 14 objects with robustly measured Lyα profiles show v 80 < 0 and line asymmetry A 91 < 0), while the P-Cygni-like profile (blueshifted absorption accompanied by redshifted emission) is only seen in F15250 (in F16156, a weak blueshifted wing is seen in the Lyα emission blueward of the strong blueshifted absorption feature, and the overall Lyα profile is thus not P-Cygni-like). This is in clear contrast with the prevalence of P-Cygni-like profiles seen in lowredshift star-forming galaxies (e.g., Wofford et al. 2013). The blueshift of the velocity centroid of the Lyα line is also prevalent in our sample, where 11 out of the 14 objects with Lyα detections show v 50 < 0. Moreover, in two Type 1 sources, 3C 273 and Mrk 1014, the broad Lyα emission (accompanied by broad N V 1238, 1243 emission, when the spectrum covers this region) shows line width typical for the broad emission line region (BELR) of AGN, which is often the case for low-redshift Type 1 AGN and quasars (e.g., Shull et al. 2012). Finally, Lyα absorption features near systemic velocities are clearly seen in 9 objects.

Profiles
The observed Lyα profiles in these dusty ULIRGs are likely affected by complex radiative transfer effects due to the resonant nature of the transition. While it is impossible to derive the intrinsic profile of the resonant Lyα line with our data, the non-resonant optical emission lines (e.g. [O iii] λ5007 forbidden line, Hα recombination line) provide a point of reference since they are less affected by absorption and scattering. Qualitatively, we may therefore infer the extent to which Lyα photons are absorbed and/or scattered, by comparing the Lyα and optical line profiles.
In Fig. 4 and 5, we plot the Lyα profiles in comparison with the [O iii] λ5007 profiles and the Hα profiles. All of these line profiles are continuum-subtracted and rescaled for better visualization. In total, 11 objects with [O iii] λ5007 observations and 15 objects with Hα observations show Lyα detections. The measurements of the [O iii] λ5007 and Hα profiles are summarized in Table 4. At first glance, one similarity between the Lyα profiles and the non-resonant line profiles is the occurrence of blueshifted emission line wings in many objects.  Blueshifted [O iii] λ5007 (and, to a lesser extent, Hα) emission is a tell-tale signature of ionized gas outflows, so the blueshifted Lyα emission may also arise from the outflowing gas. In addition, Lyα emission is generally broader than or comparable in width to the [O iii] λ5007 and Hα emission. Otherwise, the Lyα profiles are diverse among the objects in our sample and no apparent trend is seen between the profiles of Lyα and the non-resonant optical lines. We will discuss these results further in Sec. 4.2.

Lyα Escape Fraction
At first, it may be surprising that significant Lyα emission is observed from these dusty ULIRGs given the implied huge optical depth to Lyα photons. In this section, we quantify the escape of Lyα photons by calculating the Lyα escape fraction.
Under Case B recombination, the ionized region is optically thick to the Lyman series, and the intrinsic F(Lyα)/F(Hα) flux ratio depends only on the electron density and temperature. Adopting the low density limit where n e << n e,crit = 1.55 × 10 4 cm −3 , the col-lisions can be safely neglected. The intrinsic Lyα flux is then predicted to be 8.1 times the intrinsic Hα flux (Hummer & Storey 1987;Draine 2011).
Galaxies usually show Lyα-to-Hα flux ratios much less than the predicted Case B ratio. We can therefore describe the suppression of Lyα photons by comparing the observed Lyα flux to the intrinsic values indicated by the Hα emission. The escape fraction of Lyα photons can thus be defined as The intrinsic Hα flux is calculated from the observed Hα flux, corrected for the nebular reddening using the Cardelli et al. (1989) reddening curve, namely and the E(B − V ) is calculated from the Balmer decrement:  We set the intrinsic Hα/Hβ ratio to 3.1 for all but one objects as they show optical spectral features consistent with AGN activity. The only exception is F01004: it is located in the star-forming region in the BPT and VO87 diagnostic line ratio diagrams (Baldwin et al. 1981;Veilleux & Osterbrock 1987;Osterbrock & Ferland 2006), and we set the intrinsic Hα/Hβ ratio of this object to 2.87.
The Hα fluxes are measured from spectra gathered from literature, as described in Sec. 2.3.3, and can be divided into three groups: (1) Gemini/GMOS IFU observations; (2) SDSS spectra; (3) long-slit spectra from Veilleux et al. (1999). Some of our objects are point sources or show very compact morphology in the narrowband Hα images (based on the GMOS data) or R-band images (from SDSS or Veilleux et al. 1999), so no aperture corrections are needed for their Hα flux in the calculation of Lyα escape fraction. However, for the more extended objects, the aperture difference between the COS FUV spectroscopy and the optical observations need to be taken into account, as the throughput of the 2.5 COS aperture drops sharply beyond the central 0.4 . In practice, the aperture correction is negligible if at least one of the following three criteria is met: (1) it is a Type 1 AGN; (2) the PSF contribution to its overall flux in the R-band image is more than 50% based on the measurements in Veilleux et al. (2002); (3) the R-band effective radius is less than 1 based on the measurements in Veilleux et al. (2002). For the other objects, aperture corrections are needed to account for the rapid decrease of COS throughput at large radius as mentioned above. Specifically, we calculate the aperture correction factors for each group of optical observations separately: (i) for the GMOS IFU observations, we generate the Hα flux maps based on the data cube, and use the COS throughput function to vignette the Hα flux maps within the region with radius r < 1.25 . The aperture correction factor is then the vignetted Hα flux within r < 1.25 (corresponding to the COS aperture) divided by the original Hα flux within the same aperture; (ii) for the SDSS spectra, we adopt the r-band images as surrogates for the Hα flux maps. We then vignette the r-band images within the same COS throughput func-tion, and the aperture correction factor is the vignetted r-band flux within r < 1.25 divided by the original rband flux within SDSS aperture (D=3 or D=2 ); (iii) for the long-slit spectra from Veilleux et al. (1999), we follow the same logic as adopted for the SDSS spectra but use the R-band images in Kim et al. (2002). The aperture correction factor is thus the vignetted R-band flux within r < 1.25 divided by the original R-band flux within the extraction region of the long-slit spectra (2 ×5 kpc).
The aperture correction factors adopted in the calculations and resulting Lyα escape fractions are listed in Tables 4 and 3, respectively.

Continuum Luminosity at 1125Å
As a surrogate for the FUV continuum luminosity adopted in Paper I, the monochromatic luminosities at rest-frame 1125Å, log(λL 1125 ), are measured whenever the continuum is detected, with a bandpass of 20Å. These results are recorded in Table 1.

FUV Absorption Features
The focus of this section is the strongest metal absorption features detected in our objects, O VI 1032, 1038 and N V 1238, 1243, tracers of the highly ionized gas in these systems. Only 12 out of the 21 objects have continuum S/N in the vicinity of O VI 1032, 1038 and/or N V 1238, 1243 that are high enough (S/N 10 in a 500 km s −1 window) to allow for the detection of corresponding absorption lines. Out of these 12 objects, 6 objects show O VI and/or N V absorption features associated with the galaxy (velocity centroid < 13000 km s −1 and not from intervening systems), and the velocity centroids of these absorption features are all blueshifted. One more object, F15250, may display a N V absorption feature, but the doublet is so close to a group of geo-coronal emission lines that the N V 1239 transition is heavily contaminated and no robust measurements of the N V feature can be made. Our estimates for the EW and centroid velocity of the N V 1242 absorber alone are ∼0.3Å and −500 km s −1 , respectively, which has not taken into account the infilling from the N V 1238, 1243 emission. This source is excluded from the discussions related to the absorption features in the following sections.
The properties of these detected O VI and/or N V absorption features vary wildly: F07599 shows a >25000 km s −1 wide N V 1238, 1243 BAL accompanied by narrower absorbers at smaller velocities, whereas F23060 shows relatively narrow and shallow N V 1238, 1243 absorption features on top of the N V 1238, 1243 emission. Overall, the absorption features in F07599 and F01004 fall in the BAL category (velocity width > 2000 km s −1 ), while all other absorption features are classified as narrow absorption lines (NAL; velocity width < 500 km s −1 ). The object-by-object description of these absorption features is given in Appendix A.
To quantify the strength of these absorbers, we follow the same procedure as in Paper I. First, we fit the continuum and/or broad emission lines (Lyα, O VI, N V) with low-order polynomials or Gaussian profiles. After the spectra are normalized by the best-fits from the continuum and/or broad emission line fits, these absorbers are quantified using a non-parametric approach, where we measure the total velocity-integrated EWs of the outflowing absorbers in the object's rest frame, the weighted average outflow velocity and the weighted outflow velocity dispersion, The results are summarized in Table 5.

Evidence of Partial Covering for the O VI and N V Doublets
The profiles of the resolved O VI and N V doublets may be used to derive the basic characteristics of the absorbing cloud -background source system. For the absorption features in Z11598 and F13218, there is evidence of partial covering: the optical depth ratio of the spectrally resolved doublet deviates from the theoretical expectation from a simple model where the foreground cloud is illuminated by the background point source with a 100% covering fraction.
For the cases where the optical depths of the doublets (proportional to λf osc , the product of wavelength and oscillator strength) differ by a factor of ∼2 (like O VI 1032O VI , 1038O VI and N V 1238O VI , 1243, and the continuum intensity is normalized to unity, the coverage fraction (or covering factor, C f ) as a function of velocity may be obtained. Following Hamann et al. (1997a), in the simple situation where the two transitions of the doublet do not overlap with each other, we have I 1 and I 2 are the normalized intensities of the weaker and stronger absorption lines, respectively, and C f is the covering factor. Then the optical depth as a function of the velocity can be written as The column density of the ion can then be obtained by integrating the optical depth over the velocity adopting The resulting values of N ion are listed in Table 5.

Voigt Profile Fitting of O VI and N V Absorbers
A popular approach to quantify the absorption features is to fit them with the product of individual components assuming Voigt profiles for the optical depth distribution in frequency (or velocity, wavelength) space. As discussed in Sec. 3.3.1, there is evidence for partial covering in a few objects. To account for this, we also  (1): Object and absorption feature name; Column (2)-(4): Velocity-integrated EW, Weq (eq. 5), average depth-weighted velocity vwtavg (eq. 6), and average depth-weighted velocity dispersion σrms (eq. 7) of the absorption lines; Column (5): Ion column densities obtained from the analysis of the absorption doublet with partial covering model as described in Sec. 3.3.1; Column (6): Velocity-weighted covering fractions calculated from the same analysis for ion column densities in Column (5); Column (7): Number of components in the best-fits from the Voigt profile fits as described in Sec. 3.3.2. Notice that the O VI absorption in F01004 and the N V absorption of F07599 are BAL, and the 1-component fit is only tentative/experimental, with the aim to capture the overall absorption profile. Column (8): Ion column densities from the Voigt profile fits. The values for individual components from the best-fit model are separated by comma. This is also True for Column (9); Column (9): Covering fractions from the Voigt profile fits. The flag "f" in parenthesis indicates that the covering fraction is fixed to the corresponding value in the fits.
include a constant covering fraction parameter to each component of the model. The final model of the normalized intensity can then be written as C f , τ , ν, N , b, f osc , and σ 0 are the covering fraction parameter, optical depth, frequency, ion column density, Doppler parameter, oscillator strength, and cross section, respectively. Φ(ν|b, z) is the normalized Voigt profile. The adopted atomic parameters are taken from Morton (2003).
In the fitting procedures, the model described above is further convolved with the line-spread function of HST/COS tabulated on the HST/COS website 2 . We adopt a customized software built on the non-linear least-squares fit implemented in LMFIT to search for the best-fit model. In our software, a velocity component is added to the model if the Bayesian Information Criterion (BIC) (Ivezić et al. 2019) decreases, and this pro-cess stops when the minimum BIC value is found. The model with the lowest BIC value is then chosen as the best-fit to the data, which is also confirmed by a visual inspection. In addition, we have tested this by manually fitting the absorption line profile with n+1 components when n components are required by the best fit. The change in total column density is in general 0.1 (in logarithm), within the uncertainty of total column density derived from the best-fit model. The uncertainties of the best-fit parameters are calculated from the 1-σ (68.3%) confidence interval adopting the conf interval function of LMFIT, which takes into account the covariances between blended absorption components.
The best fits for all absorbers clearly detected in our sample are shown in Fig. 6. The covering factor, C f , is fixed to unity in two objects: for F21219, the fit is unable to break the degeneracy between the covering factor and ion column density; for F23060 no evidence of partial covering is suggested by the data. For Z11598 and F13218, the ion column densities are reported as lower limits due to the saturation of the absorption features. Additionally, the fits for the N V BAL in F07599 and O VI BAL in F01004 are highly uncertain given the model parameter degeneracy caused by the saturation  and smoothness of the absorption feature, and the large uncertainties in the continuum determination. Following Paper I, our fitting scheme assumes that the velocity dependence of the optical depth can be parameterized as the sum of discrete independent components with Voigt profiles and constant covering fractions, and that the individual absorbers simply overlap with each other along the line-of-sight. In reality, the C f is probably a more complex function of velocities (Arav et al. 2005(Arav et al. , 2008(Arav et al. , 2013, and the absorbing material may completely overlap, partially overlap or not overlap each other. While our approach cannot account for these details, it is sufficient to meet our primary goal, which is to characterize the overall strength and kinematics of these absorbers, and put constraints on the ion column density.
The column densities and covering fractions from these fits are summarized in Table 5. For sources Z11598, F13218 and F21219, the results are consistent with those obtained from the analysis of the absorption doublets with partial covering model as described in Sec. 3.3.1. Further discussion based on these fits is postponed until Sec. 5.

ORIGIN OF THE Lyα EMISSION
Intuitively, strong Lyα emission is not expected in dusty ULIRGs due to the huge optical depth. The origin of Lyα emission in our objects is therefore worth investigating. In this section, we focus on three factors that may help with the production and/or escape of Lyα photons in our ULIRG sample: (i) AGN: The AGN activity may intrinsically produce more Lyα photons than what we have assumed. For example, the gas density may be so high in the broad line region (and perhaps also narrow line region) of AGN that collisional excitation becomes important in promoting Lyα emission (e.g., Dijkstra 2017). Radiation from the AGN may also ionize the gas and destroy the dust in the Lyα-emitting regions and material along the lineof-sight, and therefore reduce the overall opacity to the Lyα emission. (ii) Outflow: The blueshifted Lyα emission may come from the outflowing gas. The velocity offset between the outflow and the interstellar medium decreases the optical depth to the Lyα emission radiated from the fast-moving gas. The outflow can also create low opacity pathways for Lyα photons by clearing out the gas and dust. In addition, the outflow may have broken out of the dusty ISM so that Lyα photons can escape freely. (iii) Dust: The Lyα photons are heavily affected by complicated, dust-related radiative transfer effects, which directly affect the observed properties of the Lyα emission in ULIRGs.
In our analyses, we adopt Kendall tau correlation tests to examine potential correlations between the properties of Lyα emission and those of the AGN, outflows, and dust reddening. Specifically, we adopt the method from Isobe et al. (1986) to compute the Kendall tau correlation coefficient. The p-value of null hypothesis (no correlation) is calculated to show the statistical significance of the correlation. This method can handle censored data, which is the case for EW Lyα and f esc (Lyα). We use the implementation of pymccorrelation (Privon et al. 2020), which perturbs the data with Monte Carlo method to compute the error in the correlation coefficient (Curran 2014). To expand the dynamic ranges of the variables in the analyses, by default, we also include the results of the starburst-dominated ULIRGs from M15, whenever possible. The results from these correlation tests are summarized in Table 6. 4.1. Effect of the AGN First, we examine the possible link between the strength of the AGN (more specifically, the AGN luminosities, L AGN , and AGN fractions, f AGN ; see Table 1) and the properties of Lyα emission. As shown in Fig. 7, the EWs of Lyα (EW Lyα ) increase with both L AGN and f AGN . The p-values are ∼0.002 and ∼0.023 based on the Kendall tau tests, respectively. In addition, L AGN may also correlate with the kinematic properties of Lyα (v 80 , v 50 A 91 and W 80 of Lyα; p-values 0.06-0.10). As an example, v 80 of Lyα are plotted against L AGN and f AGN in Fig. 8. These weak trends are mainly driven by the Type 1 sources with small v 80 and large W 80 : indeed, the correlations are not statistically significant (p-values > 0.1) when the Type 1 sources are excluded.
Next, we explore the behavior of the Lyα escape fraction, f esc (Lyα), with the strength of the AGN, as shown in Fig. 9. There are possible positive correlations between L AGN and f esc (Lyα) (p-value 0.062), and between f AGN and f esc (Lyα) (p-value 0.068). However, we note that these correlations are no longer significant (p-values 0.2 and 0.4, respectively) if we exclude the data points of the starburst-dominated ULIRGs from M15.

Effect of the Outflow
The prevalence of blueshifted Lyα emission line profiles in our sample, as stated in Sec. 3.1, hints at a potential link between outflows and the Lyα emission in our objects. Therefore, we start by simply checking whether there is a correlation between the EW Lyα and v 80 of Lyα, and confirm a positive result (p-value 0.02). However, we note that this trend disappears (p-values > 0.1) if we do not include the starburst-dominated objects from M15. Similarly, we have also checked the potential correlation between f esc (Lyα) and v 80 of Lyα, but no statistically significant trend is present (p-value 0.13).

Connection with the Ionized Outflow in Emission
The blueshift of the non-resonant, forbidden emission [O iii] λ5007 line in galaxies is strong evidence for ionized gas outflows . To investigate the connection between the blueshift of the Lyα emission and [O III] outflowing gas, it is thus natural to determine whether the kinematic properties derived from the  Table 6).
Additionally, as shown in Fig. 11 For the sake of completeness, we also briefly discuss the strong Hα line emission in these objects. We note that v 80 , W 80 , and A 91 of Hα are positively correlated with those of Lyα (p-value 0.001-0.039; see Table 6), whereas no correlation is seen between v 50 of Hα and Lyα (p-value 0.437). While blueshifted Hα emission may indicate outflowing gas, Hα is likely dominated by the emission from the broad emission line region (BELR) in Type 1 sources (e.g. 3C 273, Mrk 1014; see Fig. 5). In addition, the nearby [N ii] λλ6548, 6583 emission lines add uncertainty to the kinematic measurements based on Hα. Therefore, we do not use the Hα-based kinematics to examine the link between the blueshifts of Lyα and the ionized outflows in the remainder of the paper.

Connection with the O VI and/or N V Outflows
As discussed in Sec. 3.3, highly ionized O VI and/or N V outflows are detected in the HST/COS spectra of our ULIRGs. Among the 9 objects with both Lyα measurements and continuum S/N high enough to allow for There is no clear difference in the mean values and overall ranges of the Lyα escape fractions between the objects with and without O VI and/or N V outflows. The mean value of the Lyα EWs in the objects without outflow detection in O VI and/or NV is higher than those with outflows by ∼60%. However, these results are based on a very small sample, so no statistically robust conclusion may be drawn.

Connection with the Neutral Phase Outflows
Neutral gas outflows, traced by the blueshifted Na I D λλ5890, 5896 absorption features, are often detected in ULIRGs (e.g., Rupke et al. 2005;Rupke & Veilleux 2011, 2013bRupke et al. 2017). As shown in Fig. 12, three of our objects have both blueshifted Lyα emission and blueshifted interstellar Na I D absorption features with similar kinematics. Specifically, for F05189, v 50 of Lyα emission is similar to that of the blueshifted component of the Na I D absorption (∼−400 km s −1 ). For F11119, the blueshifted peak of the Lyα emission and v 50 of Na I D absorption have a similar velocity of ∼−800 km s −1 . For Mrk 231, v 50 of Na I D absorption (∼−5000 km s −1 ) is close to the velocity of the peak of the blueshifted wing of Lyα emission. Overall, if the blueshifted Lyα emission is indeed tracing the outflowing gas, as argued in Section 4.2.1, the results above hint at a possible connection between the outflowing gas traced by blueshifted Lyα and the outflowing neutral gas in these three objects. For instance, Lyα could be scattered off of the outflowing neutral gas traced by Na I D. However, given the very small number of objects where the Lyα − Na I D comparison was possible, these results may not apply to all objects in the sample.

Effect of Dust
Complex dust-related radiative transfer processes may shape the observed Lyα emission. In the following, we explore how dust and its distribution within the galaxy may affect the escape of Lyα photons qualitatively, by examining the relation between the nebular line color excess, E(B − V ) (reflecting dust reddening of the lineemitting gas) and f esc (Lyα). A more quantitative analysis on this issue requires a careful modelling of the radiative transfer processes, which is beyond the scope of this paper.
In Fig. 13, we plot f esc (Lyα) as a function of E(B − V ) for both our objects and the starburst-dominated ULIRGs from M15. Also shown in the figure is the expected f esc (Lyα) given the continuum attenuation at the Lyα transition derived from the values of E(B − V ) based on the Cardelli et al. (1989) reddening curve. In the figure, there is a lack of correlation between the E(B − V ) and f esc (Lyα) (p-value 0.371 from the Kendall tau test, see Table 6), which is inconsistent with the naive expectation that f esc (Lyα) should decrease with increasing E(B − V ). In addition, the Lyα emission in 6 out of the 14 ULIRGs with clear Lyα detections show f esc (Lyα) much higher than the values expected from the reddening curve. Similar phenomena have been seen in nearby Lyα-emitting galaxies with E(B − V ) 0.3, where Lyα emission is on average several times stronger than expected (Scarlata et al. 2009;Atek et al. 2014;Hayes et al. 2014).
The two phenomena described above may be caused by the fact that our observations reflect the integrated properties of regions with various dust extinction within one galaxy. The observed emission lines come from both the dusty regions (more likely located in the central parts of the systems) and dust-free/less dusty regions such as broad line region of the AGN, off-nuclear H II region, diffused ionized gas, tidally stripped gas during the merging process or even outflowing gas at large distance from the nucleus. The observed Lyα flux and E(B − V ) are the integrated values of all these regions with different weights, which may therefore lead to large scatters in the relation between E(B−V ) and f esc (Lyα). The higher f esc (Lyα) at a given E(B − V ) may be explained if dusty regions bias the overall E(B − V ) to a value higher than those of less dusty regions where most of the observed Lyα emission comes from (this scenario has been pointed out in previous work; e.g., Atek et al. 2014  from two ionized gas clouds: one with zero dust reddening, and one with increasing dust reddening. The Hα luminosity ratios of the two regions are also varied. The overall E(B − V ) and f esc (Lyα) integrated over the two gas clouds are shown in Fig. 14. At large E(B − V ), the f esc (Lyα) can be easily larger than the expected values based on the dust extinction curve, adopting case B. In addition, the overall distribution of the data points can also mimic the scatter/non-correlation seen in Fig. 13. Similarly, the enhancement of f esc (Lyα) may also be caused by the internal geometry of the ISM and the dust distribution within it that, together, change the behavior of Lyα photons with respect to dust attenuation. At least two possible solutions have been proposed by previous works: (i) Atek et al. (2009) andFinkelstein et al. (2009) have invoked the Neufeld (1991) geometry where dust is embedded within the H I clumps of a multiphase ISM, and the scattering of Lyα photons prevents them from encountering dust. However, radiative transport simulations show that this effective "boost" of Lyα is very difficult unless parameters are carefully fine-tuned (Laursen et al. 2013;Duval et al. 2014). The predicted increase of Lyα EW with measured attenuation is not observed in our sample. (ii) Instead, Scarlata et al. (2009) argue for a scenario that is also built upon a clumpy dust distribution. It neither requires preservational scattering as in scenario (i) nor predicts that Lyα EW increase with E(B − V ). Alternatively, the higher-than-expected f esc (Lyα) in the 6 sources may be caused by higher intrinsic Lyα emission. Such deviation from case B may be caused by the high density gas within the BELR of the AGN, where the collisional excitation enhances the intrinsic Lyα emission. Nevertheless, following this logic, it may be odd that the higher-than-expected f esc (Lyα) are mostly seen in Type 2 sources rather than Type 1 sources where the Lyα enhancement due to the dense BELR should be more prominent. For example, in the Type 1 source 3C 273, while the broad Lyα emission line resembles the broad Hα emission line (see Fig. 5) and is thus likely originated from the BELR, the f esc (Lyα) in 3C 273 is close to the expected value under case B condition.

Overall Trends
In short, the analyses above indicate that the EWs of Lyα are more closely related to the strength of AGN activity (e.g., L AGN , f AGN ), while the blueshifts of Lyα emission are more closely linked to those of non-resonant optical emission lines tracing ionized gas outflows. It is likely that the AGN activity governs the overall production of Lyα emission and the outflowing gas generates the blueshifted Lyα emission.
Additionally, the Lyα escape fractions, f esc (Lyα), tend to be slightly higher in sources with stronger AGN and faster outflows. Nevertheless, the f esc (Lyα) does not correlate with the dust reddening, E(B − V ), and 6 out of the 14 objects show f esc (Lyα) higher than the expectation from attenuation adopting case B conditions and the extinction curve from Cardelli et al. (1989).

Origin of the O VI and N V Absorbers
We now turn our attention to the O VI and N V absorption features. Given the general blueshifts of these features, they are most likely tracing gas driven out of these galaxies by the starburst and/or AGN. Other possible origins include tidal debris from the galaxy merger, intervening circumgalactic medium (CGM), and stellar absorption.
Following Paper I, the characteristics of O VI and N V absorption features that indicate a quasar-driven wind origin include (1) line profiles that are broad and smooth compared to the thermal line widths (10-20 km s −1 for N 4+ and O 5+ ) ions at temperature T 10 4.5 -10 5.5 K, (2) line ratios of the doublets O VI λ1032/λ1038 and N V λ1238/λ1243 that imply partial covering of the quasar emission source, and perhaps (3) large column densities in these high-ionization ions and/or high O VI/H I column density ratio (Hamann et al. 1997b,a;Tripp et al. 2008;Hamann et al. 2019).
All blueshifted O VI and N V absorbers detected in the 6 ULIRGs meet criterion (1) except for the 2 narrow (σ ∼ 30 km s −1 ) components out of the 3 components in the best-fit O VI profile of F13218 and the 2 narrow (σ ∼ 20 − 30 km s −1 ) components out of the 4 components in the best-fit N V profile of F21219 4 . In addition, the O VI and N V absorption features in Z11598 and F13218 meet criterion (2), and the deep BAL features in F01004 and F07599 suggest potentially high column densities in agreement with criterion (3). So, in the end, all 6 ULIRGs show certain O VI and N V absorption features consistent with quasar-driven winds. This translates into an apparent outflow incidence rate of ∼50%, close to that of the quasar sample in Paper I. We will expand on this topic in Sec. 6.
In contrast, the redshifted component of the O VI 1032, 1038 absorption feature in Z11598 may arise from tidal debris or infalling gas. It is not likely associated with the CGM since such relatively strong N V absorption line is rarely found in CGM studies at the low redshift (Werk et al. 2016), but with several notable exceptions (Ding et al. 2003;Lehner et al. 2009;Savage et al. 2010;Tripp et al. 2011;Muzahid et al. 2015;Rosenwasser et al. 2018;Gatkine et al. 2019;Zahedy et al. 2020). In addition, we cannot rule out the possibility that the two narrow components (|v| 500 km s −1 , σ ∼ 30 km s −1 ) in the O VI profile of F13218 come from the turbulent ISM and/or CGM of the system.

Location and Energetics of the Outflows
4 The narrow line widths of these components do not rule out their outflow origin. AGN outflows can also show combinations of smooth broad components with narrow comps superimposed (e.g., Yuan et al. 2002)  Note-Column (1): Object and absorption feature names; Column (2)-(4): Estimated lower limits of mass, momentum, and kinetic energy outflow rates as described in Sec. 5.2. A radial distance of 0.1 pc, an ionization correction factor of 0.2 (e.g., Tripp & Savage 2000), and a solar metallicity are adopted in the calculation. The absorption features in F01004 and F07599 are broad absorption lines (BAL), and the lower limits listed are just a rough estimation from a tentative/experimental single component fit to the BAL.
Here we start with the constraints on the location of the O VI and N V absorbers detected in our ULIRG sample. In Z11598 and F13218 (and perhaps also all other 4 ULIRG with detected absorbers), the O VI and/or N V absorption features are deeper than the FUV continuum level, implying that part of the O VI and N V emission produced in the BELR is also absorbed. These absorbers are thus located outside the BELR which has a scale of r BELR 0.1( λL λ (1350Å) 2 × 10 46 erg s −1 ) 0.55 pc (16) (e.g., Kaspi et al. 2005Kaspi et al. , 2007Bentz et al. 2013). The equation (16) above is based on the C IV-traced BELR size luminosity relation from equation (2) in Kaspi et al. (2007), and the denominator, 2 × 10 46 erg s −1 , corresponds to the λL λ (1350Å) of 3C 273. This typical scale is also consistent with the VLTI/GRAVITY result on the size of the BELR of 3C 273 (r 0.12± 0.03 pc; Gravity Collaboration et al. 2018).
As for the upper limit on the radial distance of the outflows, qualitative constraints exist for the two objects with evidence of partial covering, Z11598 and F21219. The distance of the outflows in these two sources cannot be significantly larger than the size of the region where the continuum radiation comes from. The evidence of partial covering in these two objects also sets interesting limits on the size of the absorbing cloud: If the absorbing material is a single uniform cloud with a 100% filling factor, the cloud size should thus be 0.1 pc. However, if the absorbing material is made of multiple clouds, the partial covering may instead reflect the small filling factor of the clouds, and the sizes of individual clouds may be much smaller than 0.1 pc.
The mass, momentum, and kinetic energy outflow rates of the highly ionized outflow detected in our sample are estimated using the following equations: In the equations above, Q is an approximate global outflow covering factor quoted from Paper I, based on the incidence of mini-BALs in the SDSS quasars (Trump et al. 2006;Knigge et al. 2008;Gibson et al. 2009;Allen et al. 2011); R is the radial distance of the outflow and the value of 0.1 pc is a place-holder adopted for illustrative purposes; N H is the column density of hydrogen, and v 1000 is the outflow velocity in units of 1000 km s −1 obtained from the Voigt profile fits described in Sec.

3.3.2.
Note that the Voigt profile fits are highly uncertain for the O VI BAL in F01004, and the N V BAL in F07599, as stated in Sec. 3.3.2. Therefore, only rough estimations of the ion column densities may be obtained from the fits. For Z11598 and F13218, the O VI and N V absorption features are also saturated despite their much narrower line widths, so the corresponding ion column densities from the fits are also uncertain. Nevertheless, the results derived from the Voigt profile fits are consistent with those derived from analyzing the absorption doublet with partial covering model, as described in Sec. 3.3.1 (see Table 5). For the N V absorption in F21219, the ion column density acquired from the analysis with partial covering model is a bit higher than that from the Voigt profile fits, which is expected as the partial covering factor, while fixed to unity in the Voigt profile fit, should be less than unity.
Next, the ion column densities are converted into hydrogen column densities. The metal abundance and ionization correction factor needed for the conversion can, in principle, be determined from elaborate photoionization modeling when multiple absorbers from both highand low-ionization species are present (e.g., Arav et al. 2013;Haislmaier et al. 2021), but this information is not available for our objects. For simplicity, we adopt a solar abundance (while super-solar metallicities are also reported in the literature; e.g., Moe et al. 2009) and a ionization correction factor of 0.2 (which is a conservative upper limit reported in the literature; e.g., Tripp & Savage 2000) in the calculations. As discussed at the beginning of this section, the radial distances of these absorbers are largely unconstrained, other than the fact that the absorbers in Z11598 and F13218 are located outside of the BELR (r 0.1 pc). Adopting these aforementioned values, the obtained mass, momentum, and energy outflow rates are likely lower limits as reported in Table 7. In general, they are modest ( 1%) compared with the star formation rates, AGN radiation momenta, and AGN luminosities of the systems. However, for the highly saturated O VI BAL in F01004 and N V BAL in F07599, the column density of the outflowing gas and thus the outflow energetics are probably severely underestimated.

O VI AND N V ABSORBERS IN THE COMBINED ULIRG + QUASAR SAMPLE
In this section, we explore the properties of the highly ionized O VI and N V absorbers along the ULIRG-QSO merger sequence (e.g., see Sanders et al. 1988;Hopkins et al. 2009;Veilleux et al. 2009a, and references therein) by combining the 11 ULIRGs with high enough continuum S/N to allow for O VI and N V absorption detections (ULIRG F15250 is excluded from the analysis given that the N V absorption feature is highly uncertain due to contamination from geo-coronal emission) and the 30 quasars presented in Paper I (3 of the 33 quasars in Paper I overlap with the ULIRG sample and are thus already included). In total, 6 ULIRGs (see Sec. 3.3) and 17 quasars (see Section 7.1 in Paper I) show O VI and N V absorbers indicative of quasar-driven outflows. In addition, to be consistent with the analyses in Paper I and to maximize the sample size for better statistics, we also include the narrow O VI and/or N V absorption features in the 3 quasars from Paper I that do not meet our criteria for quasar-driven outflows. The absorption features in these three objects are relatively narrow (σ rms 10-30 km s −1 ) and are redshifted in two of these objects.

Overall Incidence Rates and Regressions
Based on the measurements listed in Table 5, while quite uncertain due to the small sample size, we can estimate the incidence rate of absorption features in the ULIRG-only sample. Adopting the β distribution (Cameron 2011) used in Paper I, we obtain an incidence rate of ∼55% (1-σ range: ∼40%-68%) for the detection of O VI or N V or both absorption features. This is similar to the rate of ∼61% (1-σ range: ∼52%-68%) obtained in Paper I, which is based on the quasar sample alone. In the combined ULIRG + quasar sample, the overall incidence rate of O VI or N V or both absorption features is ∼63% (1-σ range: ∼55%-70%).
Next, we explore how the incidence rates and properties (velocity-integrated EWs, W eq , depth-weighted velocities, v wtavg , and velocity dispersions, σ rms ) of the O VI and N V absorption features depend on the AGN and host galaxy properties of our objects, adopting the β distribution above and regressions described below. Following Paper I, we apply linear regressions adopting the Bayesian model in LINMIX ERR . We use the Metropolis-Hastings sampler and a single Gaussian to represent the distribution of the parameters. LINMIX ERR allows censored values for dependent variables (y), which is the case for W eq . The only exceptions involve N H (X-ray), where both x and y values are censored, in which case we adopt the Kendall tau correlation test described at the beginning of Sec. 4. For both methods, we calculated the correlation coefficients r and their 1-σ errors. A perfect correlation gives r = 1 and a perfect anti-correlation gives an r = −1. A sample with no correlation at all gives r = 0. In addition, we have computed the significance of a correlation, P(r), as the fraction of correlation coefficients r ≤ 0 (r ≥ 0) for a positive (negative) correlation 5 . For LINMIX ERR, the distribution of r are acquired from the posterior distribution, while for pymccorrelation, they are obtained from the Monte Carlo perturbations.
For the regressions, we do not take O VI and N V data as independent measurements. Therefore, when measurements are available for both doublets in a given source, we take the average of the measurements (either detection or limit) from the two lines. If only one line is detected, we use the measurement for the detection. In addition, when multiple X-ray measurements exist for a source, we take the average of them. Errors in L bol , L IR , L IR /L BOL , L FIR /L BOL , and α U V,X are unknown or largely uncertain, so for the regressions we fix their errors to ±0.1 dex 6 . For λL 1125 , we ignore the negligible statistical measurement errors. As examples, the final 5 From LINMIX ERR, it is technically difficult to calculate the classic p-value in null hypothesis significance testing. Nevertheless, our definition of P(r) can describe the significance of the correlation similarly. Like p-value < 0.05, P(r) < 0.05 also indicates a statistically significant correlation: it suggests that the possibility is 95% for the correlation coefficient r to be larger (smaller) than 0, in the case of a positive (negative) correlation. Note that this P(r) was called p-value in Paper I, which is abandoned in this Paper II to avoid ambiguity. 6 Note that the results of the regression analysis do not change even if we adopt larger errors of up to ±0.5 dex. data points adopted for the regressions are shown in the inset panels of Fig. 15 and 16.

Dependence on the X-ray Properties
The dependence of the incidence rates on several primary AGN/host galaxy properties are listed in Table  8. The regression results between the AGN/host galaxy properties and the properties of the absorption features (W eq , v wtavg , and σ rms ) are listed in Table 9. In brief, we find that the incidence rates of the absorption features do not depend on the AGN/host galaxy properties, such as the bolometric luminosities, AGN luminosities, AGN fractions, IR luminosities, FIR luminosities, FIRto-bolometric luminosity ratios, and FUV luminosities (note that only a few key quantities are listed in Table  8). Similarly, no statistically significant (P(r) < 0.05, |r| >> 0) trends are seen between these AGN/galaxy properties and the properties of the absorption features. These negative results are largely consistent with those in Paper I based on the quasar-only sample.
Nevertheless, as discussed in Paper I, the incidence rate and properties of the absorption features do correlate with several X-ray properties of the sources. We examine these trends with the ULIRG + quasar sample below. In general, we find that the incidence rate of the absorption features is higher in the X-ray weak (relative to their UV luminosities, as quantified by α U V,X described in Sec. 2.3.2) or absorbed sources (Table 8). We also find that the W eq , v wtavg and σ rms of the absorption features correlate with α U V,X (see Table 9).
The incidence rates of the absorption features (either O VI or N V or both) are 88% (1-σ range: 76%-92%) for the objects with α U V,X < −1.3 7 and 45% (1-σ range: 35%-56%) for those with α U V,X ≥ −1.3. Such difference is statistically significant with a p-value of ∼0.014, adopting the scipy.stats implementation of the Fisher exact test with a null hypothesis that galaxies with α U V,X < − 1.3 and α U V,X ≥ −1.3 are equally likely to show O VI or N V absorption features. Additionally, as shown in Fig. 15, the objects with lower α U V,X tend to have larger W eq , smaller v wtavg , and larger σ rms , where the regressions give correlation coefficients r of −0.59 +0.15 −0.11 , 0.48 +0.15 −0.19 and −0.50 +0.18 −0.16 , respectively. These results in general confirm or strengthen those from Paper I based on the quasar-only sample and α OX : the incidence rates in Paper I were found to be 75% (1-σ range: 59%-83%) versus 55% (1-σ range: 44%-65%) for sources with α OX < −1.6 and α OX ≥ −1.6 (This difference was not considered significant since the p-value was ∼0.45, adopting the same Fisher test on the difference between the incidence rates), and the correlation coefficients r for the trends between the α OX and the W eq , v wtavg , and σ rms , were −0.62 +0.17 −0.13 , 0.31 +0.21 −0.24 and −0.55 +0.20 −0.15 , respectively. Similarly, the incidence rate of absorption features for sources with N H (X-ray) > 10 22 cm −2 is 76% (1-σ range: 65%-83%), whereas the rate for those with lower N H (Xray) is 25% (1-σ range: 17%-41%). This result is almost identical to that obtained in Paper I. Moreover, as shown in Fig. 16, the W eq of the absorption features may be higher in objects with higher N H (X-ray), which is consistent with the result from Paper I: The correlation coefficients r are 0.24 +0.03 −0.03 for the ULIRG + quasar sample and 0.19 +0.03 −0.03 for the quasar-only sample, respectively. As for the v wtavg and σ rms of the absorption features, their lack of dependence on the N H (X-ray) found in Paper I remains.

Radiation Pressure as the Most Plausible
Wind-Driving Mechanism Overall, the results from the analyses of the combined ULIRG + quasar sample reinforce the main conclusions of Paper I: (i) The incidence rate and properties of the O VI and N V absorption features (i.e., EWs, outflow velocities and outflow velocity dispersions) are positively correlated with the X-ray weakness of the sources. (ii) The incidence rate of these absorption features is higher in sources with larger X-ray absorbing column densities. The EWs of absorption features may also be higher in such sources.
This dependence of the incidence rate, EWs and kinematic properties of these outflows on the X-ray weakness and/or absorbing columns of the sources can best be explained if these outflows are radiatively driven. As discussed in detail in Section 7.3 of Paper I, the combined radiative force ("force multiplier"; Arav & Li 1994) is greatly suppressed when the gas is over-ionized by the extreme-ultraviolet (EUV)/X-ray photons, becoming too transparent to be radiatively accelerated effec-    Figure 16. Same as Fig. 15 but as function of the logarithm of the X-ray absorbing column densities NH (X-ray).
tively. The successful launching of a radiatively-driven wind thus depends on whether this ionizing EUV/X-ray radiation is shielded and/or intrinsically weak.
In the first case, the over-ionized material may serve as a radiative shield to soften the ionizing spectrum enough so that the outflow material downstream can (2) indicate relatively strong correlations with P(r) < 0.05 and |r| 0.5. Column (3): Number of data points. Column (4): Probabilities for the correlation coefficients r to be ≤0 (for positive correlation) and ≥0 (for negative correlation), as defined in the 3rd paragraph in Sec. 6.1. Column (5): Correlation coefficients r and their 1-σ errors.
be effectively accelerated (Murray et al. 1995;Proga & Kallman 2004;Proga 2007;Sim et al. 2010). However, as mentioned in Paper I, the predicted strong near-UV absorption features near systemic velocity produced by the shielding material (e.g., Hamann et al. 2013) are in general not observed in our sample. In the second case, it is proposed that the X-ray emission in weak-lined "wind-dominated" quasars are intrinsically faint and unabsorbed Wu et al. 2011;Luo et al. 2015;Veilleux et al. 2016). In our sample, several sources with fast O VI/N V outflows show evidence of intrinsically weak X-ray emission, including F07599, PG1001 and PG1004 (Luo et al. 2013(Luo et al. , 2014. More sensitive hard X-ray (>10 keV) observations of our sample would shed light on the exact origin of this X-ray weakness.
Apart from the aforementioned trends with X-ray weakness, the new data on the ULIRGs do not add significantly more support to the radiatively driven wind scenario. The lack of positive trends between the maximum velocities of the outflows and the optical, UV, bolometric luminosities or the Eddington ratios is likely due to the limited dynamic range of properties of the combined ULIRG+quasar sample, and noise in the predicted correlations associated with projection effects and variance in the launching radius and efficiency of the radiative acceleration associated with the complex microphysics of the photon interaction with the clouds (Paper I). Furthermore, we find no other case of line-locking among the outflows of ULIRGs (line-locking was observed in the outflows of two quasars in Paper I). Lastly, as in Paper I, no evidence is present that radiation pressure on dust grains is an important contributor to the radiative acceleration in our sample (the outflow properties do not correlate with the mid-, far-, and total (1 -1000 µm) infrared excesses). While the alternative thermal wind and "blast wave" models cannot be formally ruled out by our data, these models cannot readily explain the observed connection between the outflow properties and X-ray weakness and absorbing column densities (interested readers are referred to Section 7.3 of Paper I for a more detailed discussion of these models).

The Effects of Stochasticity of AGN-Outflow Activity
Intuitively, the lack of correlations between the properties of the O VI and/or N V outflows and those of the AGN is unexpected given that these winds are driven by the AGN. Together with the result that the outflow incidence rate in the ULIRG sample is virtually identical to that in the quasar sample, it may imply that the launching of these quasar-driven outflows is a stochastic phenomenon throughout the late merger stages. This is consistent with the picture that the triggering of AGN activity has a significant chaotic/random component in local gas-rich mergers and AGN (e.g., Davies et al. 2007;Veilleux et al. 2009a). Given this stochasticity of AGN activity, Veilleux et al. (2009a) warns that a sample size of 50-100 may be needed to detect any trends with merger phase. Time delays between bursts in AGN activity, the ejection of the material driven this AGN activity, and the detection of the ejected material on pc and kpc scales also likely complicate this picture .  Table 3, and Section 3.1.
• The equivalent widths of Lyα increase with increasing AGN fractions and AGN luminosities. The strength of the Lyα emission is therefore correlated with that of the AGN. See Fig. 7, Fig. 9, and Section 4.1.
• The blueshifted line centroids and/or wings of the Lyα emission correlate with those of the nonresonant optical emission lines. The 80-percentile velocities v 80 of Lyα are positively correlated with those of [O iii] λ5007 (or Hα), with the highest statistical significance among all kinematic properties measured from the data. This suggests that the blueshifted wings of Lyα emission are physically linked to the ionized outflowing gas. There is also a possible connection between the blueshifted Lyα emission lines and the blueshifted Na I D λλ5890, 5896 absorption lines tracing the cool neutral-atomic gas outflows, although the sample size (3) in this case is very limited. See Fig. 10, Fig. 12, and Section 4.2.
• For 6 of the 14 objects with clear Lyα detections, the Lyα escape fractions, calculated as the observed Lyα flux divided by the intrinsic Lyα flux expected from the extinction-corrected Hα flux, are higher than the values expected under Case B recombination adopting Cardelli et al. (1989) reddening law. Weak, positive correlations exist between the Lyα escape fractions and the AGN strength (e.g. L AGN , f AGN ) or outflow velocities (e.g. −v 80, [O III] ). See Fig. 9, Fig. 11, Section 4.1, and Section 4.2.
• Among the 12 objects with good continuum S/N, at least 6 objects show clear O VI and/or N V absorbers. The velocity centroids of these absorbers are all blueshifted and show large ranges of depth-weighted velocities (from ∼ −12690 to −170 km s −1 ) and depth-weighted velocity dispersions (from ∼100 to 4600 km s −1 ). They are likely tracing quasar-driven outflows based on their broad and smooth profiles, as well as the evidence for partial covering in several objects. The implied incidence rate of highly ionized gas outflows in our ULIRG sample (∼50%) is similar to that of the QUEST quasars in Paper I (∼60%). See Table 5, Table 8, Section 3.3, and Section 5.1.
• The locations of these O VI and N V outflows are not well constrained, although they are probably located outside of the broad emission line regions since the absorption features are deeper than the underlying continuum level in at least two (and perhaps all six) ULIRGs. The lower limits on the power and momenta of these outflows, based on conservative values of the metal abundances, ionization corrections, and radial distances of the outflowing material, are generally modest compared with the radiative luminosities and momenta of the central energy source (AGN+starburst). See Table 7 and Section 5.2.
• When combining the results on the ULIRGs presented in this paper with those on the QUEST quasar sample from Paper I, we find that the incidence rates of O VI and/or N V absorption features are higher in the X-ray weak sources with smaller X-ray-to-UV indices, α U V,X . Specifically, the incidence rate of either O VI or N V or both absorption features is 88% (1-σ range: 76%-92%) in objects with α U V,X < −1.3, and 45% (1-σ range: 35%-56%) in objects with α U V,X ≥ −1.3. Similarly, the equivalent widths, weighted outflow velocities, and weighted velocity dispersions of these features are higher in the X-ray weak sources.
These results reinforce the main conclusions of Paper I and favor radiative acceleration as the dominant driving mechanism. See Fig. 15, Table 8,  Table 9, and Section 6.2.
• As found in Paper I, the incidence rate of O VI or N V or both absorption features for sources with X-ray absorbing column densities N H (X-ray) > 10 22 cm −2 is larger (76%; 1-σ range: 65%-83%) than the rate among those with lower N H (Xray) (25%; 1-σ range: 17%-41%). The equivalent widths of the O VI and N V absorption features may also be higher in sources with larger X-ray absorbing column densities. See Fig. 16, Table 8, Table 9, and Section 6.2.
• Apart from the aforementioned correlations with the X-ray properties of the sources, the properties of the outflows do not correlate with those of the AGN/host galaxies along the late-stage merger sequence (i.e. from AGN-dominated ULIRGs to quasars). Since the incidence rate of outflows found in our AGN-dominated ULIRGs is also virtually the same as that in the quasars, these results suggest that the launching of these quasar-driven outflows is stochastic throughout the late merger stages. A rigorous exploration of the outflow properties along the merger sequence would require a larger ( 50-100) sample that covers equally well the pre-merger and late-merger stages of ULIRGs. See Section 6.4. F15250+3608: The Lyα line shows a P-Cygni-like profile. N V emission is also detected. Emission and blueshifted absorption from the N V 1242 line is visible, whereas the N V 1239 transition overlaps with the strong geocoronal emission nearby so that no measurements of it can be made. As a result, the overall properties of the N V doublet is highly uncertain. Our estimates for the EW and centroid velocity of the N V 1242 absorber alone are ∼0.3Å and −500 km s −1 , respectively, which has not taken into account the infilling from the N V 1238, 1242 emission. There