Discovery of a strong ionized-gas outflow in an AKARI-selected Ultra-luminous Infrared Galaxy at z = 0.5

TL;DR

This paper reports the discovery of a powerful ionized-gas outflow in an AKARI-selected ULIRG at z=0.5, showing significant outflow velocities and extending over 4 kpc, with implications for galaxy feedback processes.

Contribution

It presents the first detailed observation of a strong outflow in an AKARI-selected ULIRG at intermediate redshift, highlighting its properties and potential impact on galaxy evolution.

Findings

The outflow velocity reaches up to 1830 km/s.

The outflow extends to approximately 4 kpc from the galaxy center.

The galaxy exhibits a high star formation rate of 990 M_sun/yr.

Abstract

In order to construct a sample of ultra-luminous infrared galaxies (ULIRGs, with infrared luminosity, $L_{\rm IR} > 10^{12}$ L$_{\odot}$) at 0.5 < z < 1, we are conducting an optical follow-up program for bright 90-$\mu$m FIR sources with a faint optical (i < 20) counterpart selected in the AKARI Far-Infrared Surveyor (FIS) Bright Source catalog (Ver.2). AKARI-FIS-V2 J0916248+073034, identified as a ULIRG at z = 0.49 in the spectroscopic follow-up observation, indicates signatures of an extremely strong outflow in its emission line profiles. Its [OIII] 5007\AA\ emission line shows FWHM of 1830 km s$^{-1}$ and velocity shift of -770 km s$^{-1}$ in relative to the stellar absorption lines. Furthermore, low-ionization [OII] 3726\AA\ 3729\AA\ doublet also shows large FWHM of 910 km s$^{-1}$ and velocity shift of -380 km s$^{-1}$. After the removal of an unresolved nuclear component, the long-slit spectroscopy 2D image possibly shows that the outflow extends to radius of 4 kpc. The mass outflow and energy ejection rates are estimated to be 500 M$_{\odot}$ yr$^{-1}$ and $4\times10^{44}$ erg s$^{-1}$, respectively, which imply that the outflow is among the most powerful ones observed in ULIRGs and QSOs at 0.3 < z < 1.6. The co-existence of the strong outflow and intense star formation (star formation rate of 990 M$_{\odot}$ yr$^{-1}$) indicates that the feedback of the strong outflow has not severely affect the star-forming region of the galaxy.

Figures

Tables

Table 1. Table 1: Parameters of the CIGALE SED fitting.

Star formation history
OSP e-folding time $\tau$ (Myr)	500
Age of OSP $t_{\mathrm{main}}$ (Myr)	100, 200, 400, 600
Age of YSP $t_{\mathrm{burst}}$ (Myr)	10, 20, 40, 80
Mass fraction of YSP $f_{\mathrm{burst}}$	0.01, 0.03, 0.05, 0.1, 0.3, 0.5
Single stellar population (Bruzual & Charlot, 2003)
Initial mass function	Salpeter (1955)
Metallicity	0.02
Dust attenuation (Calzetti et al., 2000)
E(B-V) of YSP	0.2–2.0, per 0.2
Reduction factor $f_{\mathrm{att}}$	0.5
Dust re-emission (Dale et al., 2014)
Dust heating slope ${\alpha }_{\mathrm{SF}}$	1.125, 1.625, 2.125
AGN emission (Fritz et al., 2006)
$R_{\max }/R_{\min }$ of dusty torus	60
Optical depth at 9.7 $\mu$m ${\tau }_{9.7}$	1.0, 3.0, 6.0
Radial dust distribution $\beta$	$-$0.5
Angular dust distribution $\gamma$	0.0
Opening angle of torus $\theta$	100.0∘
Angle of viewing axis $\psi$	0.001∘, 89.990∘
$f_{\mathrm{IR}}^{\mathrm{AGN}}=L_{\mathrm{IR}}^{\mathrm{AGN}}/L_{\mathrm{IR}}^{\mathrm{tot}}$	0.01–0.09, per 0.01
	0.1–1.0, per 0.05
Radio non-thermal emission (Colina & Perez-Olea, 1992)
FIR / radio coefficient $q_{\mathrm{FIR}}$	2.58
Radio power-law slope ${\alpha }_{\mathrm{rad}}$	0.7

Table 2. Table 2: Estimation of galaxy properties from the best fit SED ⋆ .

SFR estimated from $L_{\mathrm{IR}}^{\mathrm{dust}}$	$990\pm 44$ M⊙ yr-1
Stellar mass $M_{\star }$	$9.46\pm 1.69\times {10}^{10}$ M⊙
Stellar un_att luminosity $L_{\mathrm{un}\_\mathrm{att}}^{\mathrm{star}}$	$6.86\pm 0.31\times {10}^{12}$ L⊙
$E\left(B-V\right)$ of YSP	$1.49\pm 0.10$
SF IR luminosity $L_{\mathrm{IR}}^{\mathrm{dust}}$	$5.75\pm 0.26\times {10}^{12}$ L⊙
Total IR luminosity $L_{\mathrm{IR}}^{\mathrm{tot}}$	$6.13\pm 0.25\times {10}^{12}$ L⊙
AGN bolometric luminosity $L_{\mathrm{bol}}^{\mathrm{AGN}}$	$4.65\pm 0.39\times {10}^{11}$ L⊙
$f_{\mathrm{IR}}^{\mathrm{AGN}}=L_{\mathrm{IR}}^{\mathrm{AGN}}/\left(L_{\mathrm{IR}}^{\mathrm{AGN}}+L_{\mathrm{IR}}^{\mathrm{dust}}\right)$	$6.27\pm 0.66\%$
$f_{\mathrm{bol}}^{\mathrm{AGN}}=L_{\mathrm{bol}}^{\mathrm{AGN}}/\left(L_{\mathrm{bol}}^{\mathrm{AGN}}+L_{\mathrm{un}\_\mathrm{att}}^{\mathrm{star}}\right)$	$6.34\pm 0.65\%$
⋆in the case of ${\alpha }_{\mathrm{SF}}=1.125$ and $q_{\mathrm{FIR}}=2.60$

Equations

$$\sigma \equiv \frac{1}{2\sqrt{2\ln 2}}\sqrt{{\left(\frac{\Delta v}{c}{\lambda }_{\mathrm{line}}\right)}^2+\Delta {\lambda }_{\mathrm{ins}}^2},$$

$$F={\int }_{-\infty }^{v_F}f(v^{\prime })\mathrm{d}{v}^{\prime }/{\int }_{-\infty }^{\infty }f(v^{\prime })\mathrm{d}{v}^{\prime },$$

$$L_{\mathrm{bol}}^{\mathrm{AGN}}=L_{\mathrm{central}}^{\mathrm{AGN}}+L_{\mathrm{scattering}}^{\mathrm{AGN}}+L_{\mathrm{thermal}}^{\mathrm{AGN}},$$

$$L_{\mathrm{bol}}^{\mathrm{host}}=L_{\mathrm{attenuated}}^{\mathrm{star}}+L_{\mathrm{IR}}^{\mathrm{dust}}=L_{un\_attenuated}^{\mathrm{star}}.$$

$$M_{\mathrm{ion}}$$

$$L_{\mathrm{H}\alpha }$$

$$\frac{M_{\mathrm{ion}}}{10^8{\mathrm{M}}_{\odot }}=3.330\times \frac{L_{\mathrm{H}\alpha }}{10^{43}\mathrm{erg}{\mathrm{s}}^{-1}}{\left(\frac{n_{\mathrm{e}}}{100\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1}.$$

$$L_{[\mathrm{O}\mathrm{III}]}$$

$$\frac{M_{\mathrm{ion}}}{10^8{\mathrm{M}}_{\odot }}=4.613\times \frac{L_{[\mathrm{O}\mathrm{III}]}}{10^{44}\mathrm{erg}{\mathrm{s}}^{-1}}{\left(\frac{n_{\mathrm{e}}}{100\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1}.$$

$$E(B-V)=-\frac{2.5}{k_{\mathrm{H}\alpha }-k_{\mathrm{H}\beta }}{\log }_{10}\left[\frac{L_{\mathrm{H}\alpha ,\mathrm{obs}}/L_{\mathrm{H}\beta ,\mathrm{obs}}}{L_{\mathrm{H}\alpha ,\mathrm{int}}/L_{\mathrm{H}\beta ,\mathrm{int}}}\right],$$

$$L_{\mathrm{cor}}=L_{\mathrm{obs}}\times 10^{0.4\times k_{\lambda }\times E(B-V)}.$$

$${\dot{M}}_{\mathrm{out}}=M_{\mathrm{out}}{v}_{\mathrm{out}}\times \frac{3}{R_{\mathrm{out}}},{\dot{E}}_{\mathrm{out}}=\frac{1}{2}{\dot{M}}_{\mathrm{out}}{v}_{\mathrm{out}}^2,$$

$$\frac{L_{\rm wind}}{10^{46}\ \rm erg\ s^{-1}}=3.0\times\Big{(}\frac{v_{\rm out}}{1000\ \rm km\ s^{-1}}\Big{)}^{3}\Big{(}\frac{R_{\rm out}}{10\ \rm kpc}\Big{)}^{2}\frac{n_{0}}{1\ \rm cm^{-3}},$$

$$\frac{L_{\mathrm{Edd}}}{10^{46}\mathrm{erg}{\mathrm{s}}^{-1}}\equiv 1.26\times \frac{M_{\mathrm{BH}}}{10^8{\mathrm{M}}_{\odot }}.$$

$${\log }_{10}(\frac{M_{\mathrm{BH}}}{10^8{\mathrm{M}}_{\odot }})=1.05\times {\log }_{10}(\frac{M_{\star }}{10^{11}{\mathrm{M}}_{\odot }})-\mathrm{0.55.}$$

Full text

1 Introduction

Ultra-luminous infrared galaxies (ULIRGs, with infrared luminosity, $L_{\rm IR}>10^{12}$ L*⊙) are a population of the most intensely star-forming galaxies in the local universe, with star formation rates (SFR) of $10^{2}$–$10^{3}$ M⊙* yr*-1*(Sanders & Mirabel, 1996; Rowan-Robinson, 2000). They are thought to represent rapidly growing phase of massive galaxies before quenching of their star formation. In a widely accepted scenario of massive galaxy formation, a ULIRG is thought to be formed after merging of two gas-rich disk galaxies. An intense star formation in the nucleus is triggered by a rapid feeding of gas into the central region. This star formation results in a formation of dust particles, which, in turn, absorb rest-frame optical and ultraviolet (UV) emission from young blue stars, and re-radiate at far-infrared (FIR) / sub-millimeter wavelengths. Following the feeding of gas into the nuclear region, an active galactic nucleus (AGN) emerges (Sanders et al., 1988; Hopkins et al., 2008). As AGN heats the dusty materials around the supermassive blackhole (SMBH), the spectral energy distribution (SED) of the ULIRG shows a clear mid-infrared (MIR) bump since the temperature of dust heated by AGN is higher than that heated by the star formation (Imanishi et al., 2007; Alonso-Herrero et al., 2012; Lopez-Rodriguez et al., 2013; Ichikawa et al., 2014). The vigorous starburst and / or AGN in a ULIRG can induce a strong outflowing wind, which would blow out the gas and dust and terminate the activity of the galaxy. It is commonly accepted that a powerful AGN is required to drive outflows to high velocity, e.g., $v_{\rm max}\geq 500$ km s*-1* ($v_{\rm max}=|v_{\rm shift}|+v_{\rm dispersion}$, Rupke et al. (2005); Westmoquette et al. (2012); Arribas et al. (2014); Harrison et al. (2014)).

Recently, signatures of such outflowing gas in multi-phase have been found in various ULIRGs (neutral, Rupke & Veilleux (2011), Perna et al. (2015); ionized, Soto et al. (2012), Rodríguez Zaurín et al. (2013); and molecular, Veilleux et al. (2013), Saito et al. (2017)). Those observations support the idea that outflow plays an important role in the transition of a ULIRG from an extreme starburst to a quiescent galaxy. One widely employed diagnostic to identify the ionized outflowing gas is the broad [Oiii] 5007Å emission line profiles (Christopoulou et al., 1997; Tadhunter et al., 2001; Zamanov et al., 2002). As a forbidden (collisionally excited) transition, the [Oiii] line is a good tracer of the kinematics of gas in Hii regions formed by hot, young stars in a massive starburst, or gas in the narrow line region (NLR) formed by an AGN (Osterbrock & Ferland, 2006), which is found to be extended over parsecs to several kiloparsecs by the spatially resolved long-slit and integral-field spectroscopy observations (e.g., Westmoquette et al. (2012); Liu et al. (2013); Harrison et al. (2012)). Compared to the [Oiii] line with relatively high ionization potential (IP, 35.12 eV), the [Oii] 3726Å 3729Å doublet with lower IP (13.62 eV) can trace ionized gas with lower ionization degree, which shows more extended structure than gas traced by [Oiii] line (Collins et al., 2009). Several works reported that [Oii] emission line also shows outflow signatures, although it usually does not display the extremely broad feature seen in [Oiii] emission line (Zakamska & Greene, 2014; Perna et al., 2015; Toba et al., 2017). Broad [Oii] line profiles possibly imply that the ionized outflow could couple to the interstellar medium (ISM) at larger radius, i.e., over galaxy-wide scales ($\sim 10$ kpc).

In this paper we report a discovery of one ULIRG selected from AKARI 90-$\mu$m FIR survey, AKARI-FIS-V2 J0916248+073034 (hereafter J0916a), which indicates signatures of an extremely strong outflow in both [Oiii] and [Oii] emission line profiles but possibly with low AGN contribution to its bolometric luminosity. The paper is organized as follows. We present the photometric and spectroscopic observations and data reduction in Section 2. The methods and results of optical spectral analyses and multi-band SED fitting are shown in Section 3 and Section 4, respectively. We discuss the properties of the extreme outflow in Section 5. Section 6 summarizes the conclusion. Throughout the paper, the rest frame wavelengths are given in the air. We adopt the cosmological parameters $H_{0}=$ 70 km s*-1* Mpc*-1*, $\Omega_{m}=0.3$ and $\Omega_{\Lambda}=0.7$.

2 Observation and Data Reduction

In order to construct a sample of ULIRGs at intermediate redshifts ($0.5<z<1$), we are conducting an optical follow-up program for 90-$\mu$m FIR sources in the AKARI Far-Infrared Surveyor (FIS) Bright Source catalog (Ver.2, hereafter FISBSCv2, Yamamura et al., in prep.) by utilizing the Sloan Digital Sky Survey (SDSS) optical imaging data. We focus on optically-faint sources, which are expected to be at relatively high redshifts ($z\sim 1$) and luminous in the IR wavelength range. The detection limit of FISBSCv2 catalog reaches 0.44 Jy at 90 $\mu$m, which is deeper than the previous FIR all-sky survey by the Infrared Astronomical Satellite (IRAS). Furthermore, thanks to the smaller PSF of AKARI (spatial resolution of FWHM $=1$–$1.5^{\prime}$, Doi et al. (2015)) than IRAS, we can achieve reliable identifications even with faint ($i\sim 20$) optical counterparts. Additionally, recently released Wide-field Infrared Survey Explorer (WISE) all-sky survey catalog (Wright et al., 2010) in the MIR wavelength range (3.4 (W1), 4.6 (W2), 12 (W3), and 22 (W4) $\mu$m) with a spatial resolution of FWHM $=6.1$–$12.0^{\prime\prime}$, is useful to narrow down the positional uncertainty of the AKARI FIR sources. The detection limit of the WISE survey is significantly deeper (1 mJy at 12 $\mu$m) than the AKARI FIR catalog for sources with typical IR SED (e.g., Polletta et al. (2007)). Therefore, the AKARI FIR survey with the WISE pinpointing and SDSS optical photometric data is unique in identifying rare luminous FIR galaxies with faint optical counterpart. The cross-matched sample covers $\sim 7000$ deg2, which is 10 times wider than the Herschel Astrophysical Terahertz Large Area Survey (H-ATLAS) FIR survey at the similar depth (Smith et al., 2017).

We started spectroscopic follow-up of the optically-faint FIR sources using FOCAS on the Subaru telescope. Eight objects are observed in a service program in S17A (S17A0216S, PI:Masayuki Akiyama). For each object, the long slit spectra were obtained from two exposures, with integration time of 420 seconds per exposure. The 300B grism with SY47 filter were used. The configuration provides spectral resolution of $\rm R=500$ at 5500 Å with slit width of 0.5*′′* and covers 4600–9200 Å in the first-order of SY47. The data is reduced through bias-subtraction, flat-fielding, cosmic rays removal and wavelength calibration, before extraction of an individual spectrum using the IRAF package focasred. The wavelength calibration accuracy and spectral resolution measured using the positions and widths of the night-sky emission lines are 1 Å and 10.6 Å (FWHM), respectively, in the entire spectrum. The flux calibration is done with the spectra of standard stars, Feige34 and GD153, taken after the observations of the main targets. The atmospheric reddening was corrected using the Maunakea mean atmospheric extinction curve from Buton et al. (2013). We also corrected for the atmospheric oxygen and water vapor absorptions by comparing the observed spectra of standard stars with those archived in the CALSPEC database111http://www.stsci.edu/hst/observatory/crds/calspec.html. The slit-loss effect was estimated to be approximately 50 % by comparing the integrated flux to the SDSS photometry in the g- and i-bands, and corrected in the reduced spectrum. The seeing during the observation was 0.5*′′*.

All of the eight objects are identified as emission line galaxies at $0.3<z<0.6$. One of the objects, J0916a with i-band magnitude of 20.29 mag shows emission lines at $z=0.49$ with features of a strong outflow. Since the H$\alpha$-[Nii] complex of J0916a lies around 9800 Å in the observed frame, we also examine the observed spectrum in the 9200–10000 Å range after correcting for the contamination of the second-order spectrum, using the throughput of the first- and second-order of the 300B+SY47 setting222https://www.subarutelescope.org/Observing/Instruments/FOCAS/spec/

efficiency.html and the observed first-order spectrum below 5000 Å. The optical image of the object is shown in Figure 1. The object has a neighboring galaxy J0916b with i-band magnitude of 21.43 mag. During the observation of this object the slit was aligned to the direction along the two galaxies.

We collected the photometric data for J0916a from SDSS (u, g, r, i, z, 3500–9100 Å, DR7, Abazajian et al. (2009)), WISE (W1–W4, 3.4–22 $\mu$m, Wright et al. (2010)), AKARI (Wide-S, 90 $\mu$m, Yamamura et al., in prep.) and the Very Large Array Faint Images of the Radio Sky at Twenty-cm (VLA FIRST Survey, Becker et al. (1995)). WISE images of J0916a show one foreground star at the edge of the AKARI beam in the short wavelength bands (W1 and W2 panels in Figure 2), but we expect that it does not contaminate the photometry of AKARI  because the star has blue infrared color and the photometry in W4 is dominated by J0916a.

3 Spectroscopic Properties

3.1 Spectral Fitting and Line profiles

The reduced spectrum of J0916a is shown in Figure 4, which covers the wavelength range from [Oii] doublet to H$\alpha$-[Nii] complex. In order to examine the emission line properties of J0916a, firstly the stellar continuum with absorption lines are subtracted from the galaxy spectrum using the Penalized Pixel-Fitting stellar kinematics extraction code (pPXF), which is a public code for extracting the stellar kinematics and stellar population from absorption line spectra of galaxies (Cappellari & Emsellem, 2004; Cappellari, 2012). The code uses the Medium resolution INT Library of Empirical Spectra (MILES) library (Vazdekis et al., 2010), which contains single stellar population synthesis models and covers the full range of the optical spectrum with a resolution of 2.3 Å (FWHM). The best luminosity-weighted fitting result is shown as red curve in Figure 4 , and a zoomed-in spectrum covering strong Balmer absorption lines (H$\epsilon$, H8, H9, etc.) is shown in Figure 5. The systemic redshift is estimated to be $0.4907\pm 0.0001$. In order to determine the uncertainty in the systemic redshift, we generate artificial spectra by adding random noise with scatter equal to the standard deviation of the fitting residuals and then perform the same fitting procedure for the artificial spectra using pPXF. The procedures are repeated for 100 times and the scatter is adopted as the uncertainty of the measurement. 333 For the uncertainty of the systemic velocity determination, we also consider the difference from the velocity of the narrow H$\alpha$ line. We take the velocity of narrow H$\alpha$ emission line in relative to stellar absorption lines as the upperlimit of the uncertainty of the systemic velocity, which is measured as approximately 90 km s*-1*, or 0.0003 in redshift.

The measured velocity dispersion of the stellar component is about 370 km s*-1* after correction for the instrumental profile. The luminosity-weighted age distribution indicates a main stellar population of 2 Gyr and a starburst population of 70 Myr.

The emission line profile fitting is carried out using the Python Spectroscopic Toolkit (PySpecKit) developed by Ginsburg & Mirocha (2011). At first, we examine the profiles of the [Oiii] 4959Å 5007Å doublet since they are bright and isolated from other emission lines. We tie the central wavelengths and dispersions of the Gaussian profiles for the [Oiii] doublet, and fix the intensity ratio of the two lines to the theoretical value of 1:2.92 (Osterbrock & Ferland, 2006). The velocity shift and FWHM of [Oiii] line is $-770$ km s*-1* and $1830$ km s*-1*, respectively, if a pair of tied Gaussian profiles is used to fit [Oiii] doublet. The residual shows a blue excess with shift of over $2000$ km s*-1* but the excess is within $3\sigma$ level, thus an additional Gaussian profile is not adopted in the following fitting procedure. No narrow Gaussian component is required to explain the profile of the [Oiii] doublet.

We then fit the [Oii] doublet and H$\alpha$-[Nii] complex regions. The intensity ratio of [Nii] 6548Å 6583Å doublet is fixed to the theoretical value of 1:3.06 (Osterbrock & Ferland, 2006). We also fix the intensity ratio of [Oii] 3726Å 3729Å doublet to 1:1.31. The line ratio is calculated with the empirical function of Sanders et al. (2016) by assuming electron temperature of $10^{4}$ K and electron density of 100 cm*-3*, which are the typical photoionization equilibrium temperature and electron density of Hii regions (Sanders et al., 2016) and AGN NLR (Netzer, 1990). We simultaneously fit the forbidden lines, by tying the centers and dispersions of the multiple emission lines for each component444Note that tying the dispersion here means the assumption that for a certain component every emission line shows the same velocity dispersion after correcting for the instrumental broadening. Thus for each component the width $\sigma$ of the Gaussian profiles is defined as:

$\displaystyle\ \ \ \sigma\equiv\frac{1}{2\sqrt{2\ln 2}}\sqrt{\left(\frac{\Delta v}{c}\lambda_{\rm line}\right)^{2}+\Delta\lambda_{\rm ins}^{2}},$

  with the same velocity FWHM $\Delta v$ for each line, where the $\lambda_{\rm line}$ and $\Delta\lambda_{\rm ins}$ are the observed central position of the emission lines and instrumental broadening FWHM, respectively. . In the case where only a narrow Gaussian profile ($<500$ km s*-1*) is adopted for [Oii] line, a residual excess of 6–7$\sigma$ arises on the blue edge of the [Oii] line, which indicates that another broad component is required to fit this blue wing. Owing to the blending of H$\alpha$ and [Nii] lines, it is hard to determine whether the broad [Nii] lines exist or not. Since IP of [Nii] (14.53 eV) is close to IP of [Oii] (13.62 eV), we assume [Nii] line has the same profile as the [Oii] line and consider the case with broad components of [Nii] doublet (Figure 6) as a nominal fitting result, hereafter. In order to reduce the degeneracy, the shift and width of broad [Oii] and [Nii] lines are tied to the values measured from the [Oiii] line.

{longtable}

c—ccccc Line properties determined in the best fitting profile (bottom panel of Figure 6). Name $v_{50}$ (km s*-1*) FWHM (km s*-1*) $w_{80}$ (km s*-1*) Flux ($10^{-16}$erg s*-1* cm*-2*) $L_{\rm obs}$ ($10^{41}$erg s*-1*)

[Oii] 3726Å 3729Å  $-380\pm 90$ $910\pm 50$ $1690\pm 50$ $7.6\pm 0.2$ $7.0\pm 0.2$

[Neiii] 3869Å  $-670\pm 110$ $1800\pm 110$ $1910\pm 70$ $3.3\pm 0.3$ $3.0\pm 0.2$

H$\beta$  $-50\pm 90$ $310\pm 50$ $380\pm 50$ $0.9\pm 0.2$ $0.9\pm 0.2$

[Oiii] 5007Å  $-770\pm 100$ $1830\pm 70$ $1940\pm 70$ $11.6\pm 0.3$ $10.7\pm 0.3$

H$\alpha$  $-240\pm 90$ $360\pm 40$ $2010\pm 140$ $13.8\pm 1.0$ $12.8\pm 0.9$

[Nii] 6583Å  $-380\pm 90$ $720\pm 50$ $1680\pm 50$ $17.2\pm 0.8$ $16.0\pm 0.7$

\endhead

In order to describe the multi-component profiles, following Zakamska et al. (2014, 2016) and Perna et al. (2015), we use velocity defined from the normalized cumulative distribution:

[TABLE]

where $f(v)$ is the flux per velocity unit from the best fit model, $F$ equals to the fraction of flux with velocity $v\leq v_{F}$. The velocity shift is defined to be $v_{50}$. The width comprising 80 percent of the total flux, $w_{80}\equiv v_{90}-v_{10}$, is equivalent to $2.563\ \sigma$ for a single Gaussian profile. We also define the total FWHM as the width at half maximum of the total profile consisting of both of the narrow and broad components. The $v_{50}$, FWHM, and $w_{80}$ of each emission line are listed in Table 3.1. The uncertainty for the velocity shifts and widths are estimated using a Monte Carlo method. We generate artificial spectra by adding random noise with scatter equal to the standard deviation of the fitting residuals and then perform the same fitting procedure for the artificial spectra. The procedures are repeated for 100 times and the $1\sigma$ dispersion is used as the measurement error for each parameter.

Signatures of an extremely strong outflow are indicated in the emission line profiles. The [Oiii] line shows FWHM of $1830\pm 70$ km s*-1* with velocity shift of $-770\pm 100$ km s*-1* in relative to the stellar absorption lines 555 The uncertanty of the velocity shift is estimated to be $\sqrt{\sigma_{\rm shift,line}^{2}+\sigma_{\rm system}^{2}}$, where $\sigma_{\rm shift,line}$ and $\sigma_{\rm system}$ are the scatters for emission lines and the systemic velocity, respectively. We take $\sigma_{\rm system}=90$ km s*-1*, which is the deviation between the velocities of stellar absorption lines and narrow H$\alpha$ emission line, and can be considered as the upperlimit of the uncertainty of the systemic velocity. . More interestingly, the low-ionization [Oii] doublet also shows large FWHM of $910\pm 50$ km s*-1* with velocity shift of $-380\pm 90$ km s*-1*. The velocities of [Oiii] and [Oii] indicate the object shows one of the most powerful outflows observed among ULIRGs / DOGs (Toba et al. 2017; see discussion of Section 5.1) at $z<1$. The extremely large width of the low-ionization [Oii] line probably suggests that the strong outflow can extend beyond the nuclear / bulge region of the galaxy.

3.2 Long-slit spectroscopy 2D image of the emission lines

In order to examine the spatial distribution of the outflow along the slit direction, we plot long-slit spectroscopy 2D image of the emission lines after subtracting the stellar continuum components. As shown in Figure 1, one companion galaxy (hereafter J0916b) exists at approximately $2^{\prime\prime}$ ($\sim 11$ kpc) from J0916a. Therefore the stellar contribution from both of the two galaxies are considered. We extract the integrated stellar spectrum of J0916b using the same method as for the J0916a. It is hard to fit the stellar continuum for each line of the data due to the low S / N ratio at the outer regions. Thus, we assume the inner and outer regions of each galaxy have the same stellar continuum, and apply a Gaussian profile to normalize the spatial distribution. The Gaussian profile is determined using the integrated flux distribution in line-free regions. The emission line components can be obtained after subtracting the stellar components of these two galaxies. Three zoom-in views around [Oii], [Oiii], and H$\alpha$-[Nii] emission lines are shown in Figure 7. Only regions with S/N $>3\sigma$ are shown in each panel. The noise ($\sigma$) is estimated to be the scatter of the off-source region in the spectroscopy 2D image.

J0916b lies at the same redshift ($z=0.49$) as J0916a. Among the three zoom-in views in Figure 7, J0916b only has relatively luminous [Oii] line, which is connected to the extended [Oii] component of J0916a. The [Oii] ‘bridge’ component possibly indicates a direct interaction between the two galaxies, which is consistent with the tidal optical morphology (Figure 1), and also consistent with the evolutionary scenario that ULIRGs are thought to be the descendants of galaxy mergers. On the other hand, the outflowing gas of J0916a shows extended structure ([Oii] and [Oiii] panels in Figure 7), which supports a possible evidence of the galaxy-wide outflow. We discuss the velocity structure and ionization source further in Section 5.1.

4 SED Analyses with multi-band photometric data

4.1 Building templates

The offset between J0916a and J0916b ($\sim$ 2*′′) is smaller than the spatial resolution of WISE (6.1–12.0′′) and AKARI (1.0–1.5′*). Since J0916b is fainter in the optical and has only weak emission lines, which suggests no strong AGN or star formation activity, we assume that all of the observed MIR and FIR fluxes are associated with J0916a.

We perform a SED fitting for J0916a using the Code Investigating GALaxy Emission (CIGALE, Noll et al. (2009); Serra et al. (2011)) of version 0.12. Stellar continuum in the UV-optical as well as dust emission in the IR band are modeled in CIGALE, assuming the energy balance between absorption and re-emission by interstellar dust. AGN component is also modeled with the energy balance between absorption and re-emission by the dusty torus. The total IR luminosity between 1 $\mu$m and 1000 $\mu$m, $L_{\rm IR}$, is defined as the sum of the emissions from dust heated by stars and AGN. Several parameters need to be set to generate templates. The observed photometric data points are then fitted with the templates. The best fit template is selected with $\chi^{2}$ minimization in the CIGALE code.

In order to model the stellar continuum, the star formation history (SFH) needs to be assumed. We adopt two stellar components, i.e., an old stellar population (OSP) with an exponentially declining SFR and a young stellar population (YSP) with a constant SFR (Buat et al., 2011, 2014). In order to avoid too many parameters in the SED fitting, only three parameters, i.e., the ages of OSP and YSP as well as the mass fraction of YSP ($f_{\rm burst}$), are used as free parameters to reproduce SFH (see Table 1). The values of the age parameters are selected based on the estimation from luminosity-weighted fitting result of stellar continuum (Figure 4). The e-folding time ($\tau$) of OSP is fixed to be 500 Myr, assuming a typical quenching timescale for local galaxies (Lian et al., 2016). The single stellar population model of Bruzual & Charlot (2003) with the initial mass function of Salpeter (1955) is adopted. Metallicity is fixed to be the solar metallicity of 0.02.

Calzetti et al. (2000) dust attenuation law and Dale et al. (2014) dust emission model are adopted in the SED fitting. In the dust attenuation model, a reduction factor is applied to the color excess of the old stellar population. The reduction factor is defined as $f_{\rm att}\equiv E(B-V)_{\rm OSP}/E(B-V)_{\rm YSP}=A_{\rm V}^{\rm ISM}/A_{\rm V}^{\rm BC}$, under the assumption that the young stars are embedded in dense molecular birth clouds (BC) and old stars are in diffuse ISM (Calzetti et al., 2000; Charlot & Fall, 2000; da Cunha et al., 2010). In order to reduce the number of free parameters in the SED fitting, we fix $f_{\rm att}=0.5$, following Buat et al. (2012) and Lo Faro et al. (2017). In the dust emission model of Dale et al. (2014), the dust is assumed to be exposed to a distribution of starlight intensities described by a power-law function: $dM_{\rm dust}/dU=U^{-\alpha_{\rm SF}}$. The exponent $\alpha_{\rm SF}$ determines the relative contributions of different SED templates, which are comprised of calculated emission from small and large dust grains, as well as Polycyclic Aromatic Hydrocarbon (PAH) whose spectrum is derived from the average MIR spectrum of star-forming galaxies observed by the Infrared Space Observatory (ISO) satellite. A lower $\alpha_{\rm SF}$ indicates stronger heating intensities with the FIR emission peak at shorter wavelengths, and corresponds to more active star-forming activities in a galaxy (Dale et al., 2014). Due to the lack of detections in the sub-millimeter bands, it is hard to constrain the FIR emission peak for J0916a. Therefore we only test three values of $\alpha_{\rm SF}$ based on the results in the literature, i.e., $\alpha_{\rm SF}\sim 1.1$, which corresponds to a SED peaks at around 60 $\mu$m in the rest frame with a translated grey-body temperature of 40 K, from an IRAS 60-$\mu$m selected ‘warm’ ULIRG sample at $z=0.1$ (Symeonidis et al., 2011); $\alpha_{\rm SF}\sim 1.6$ for a SED with a peak at 90–100 $\mu$m and temperature of about 30 K, from the $Spitzer$ / MIPS 70-$\mu$m selected ‘cold’ ULIRG sample at $0.1<z<1.2$ (Symeonidis et al., 2011) and the 90-$\mu$m selected ULIRGs in the AKARI Deep Field-South (ADF-S) survey at $0.4<z<1.2$ (Małek et al., 2017); and $\alpha_{\rm SF}\sim 2.1$ for SEDs of normal galaxies in ADF-S field (Małek et al., 2017).

The Fritz et al. (2006) model is adopted to generate AGN emission templates. The model consists of three different components: the emission directly from the central source, the scattered component, and the thermal radiation of a dusty torus heated by the central source. Six parameters were used by Fritz et al. (2006) to describe the structure and geometry of the dusty torus and to calculate the radiation transfer. However, in order to reduce the number of free parameters in the SED fitting with limited photometric data points, we fix the radial ($\beta=-0.5$) and angular ($\gamma=0.0$) dust distribution, the opening angle of the torus ($\theta=100.0^{\circ}$) as well as the ratio of the radii ($R_{\rm max}/R_{\rm min}=60$) following Małek et al. (2017). Only two extreme values of the viewing angle of the torus $\psi$, 0.001*∘* and 89.990*∘*, are considered, which corresponds to type-2 and type-1 AGNs, respectively. Silicate optical depth $\tau_{9.7}$ and AGN fraction $f_{\rm IR}^{\rm\scriptscriptstyle AGN}$, which is defined as the ratio of AGN IR luminosity to the total IR luminosity of the galaxy, are chosen as free parameters in the SED fitting procedure with ranges of 1.0–6.0 and 0.01–1.0, respectively. Although AGN IR emission is also considered in the Dale et al. (2014) model, we exclude this component because the templates are mainly from type-1 AGNs and do not have the absorption features in the MIR band which are observed in type-2 AGNs.

Since J0916a was also detected at 1.4 GHz in the VLA FIRST survey, we consider the radio non-thermal component in the SED fitting. Two parameters, FIR / radio correlation coefficient ($q_{\rm FIR}$) and slope of the power-law synchrotron emission ($\alpha_{\rm rad}$) are fixed to be 2.58 and 0.70, respectively, based on the model of Colina & Perez-Olea (1992), in which the non-thermal radio emission from radio supernovae and supernova remnants are the main contributor to the radio luminosity in the ULIRGs. A fixed $q_{\rm FIR}$ results in the constraint on the ratio between the FIR and radio luminosities of J0916a. We also test other values of $q_{\rm FIR}$ from the literature, e.g., 2.1 for sub-millimeter galaxies (SMGs) and ULIRGs at $1<z<3$ (Kovács et al., 2006), 2.3 for radio-identified, flux-limited ($S_{\rm 60\mu m}\geq 2$ Jy) IRAS sources (Yun et al., 2001), and 2.8 for the local spirals and starburst galaxies (Condon & Broderick, 1991). If we fit the SED with changing $q_{\rm FIR}$ of 2.1, 2.3, 2.58, and 2.8, the reduced $\chi^{2}$ for the best fit cases are 11.4, 6.3, 0.8, and 3.7, respectively, with the estimated $L_{\rm IR}$ from interstellar dust varying between $2.6\times 10^{12}$ and $6.4\times 10^{12}$ erg s*-1*. The value of $q_{\rm FIR}$ is hard to be constrained with the current limited data points, especially in the FIR wavelength range. Therefore we select $q_{\rm FIR}=2.58$, since it results in the smallest $\chi^{2}$. The values and ranges of the free parameters are summarized in Table 1.

In order to estimate the uncertainties of the output parameters, we perform the same fitting procedure to the simulated mock SED. The mock SEDs are generated by adding random values with the same scatters as the measurement errors to the observed photometric data points in each band. We repeat the fitting procedure for 100 mock SEDs and use resulting standard deviations as the $1\sigma$ uncertainties of the output parameters.

4.2 The best fit properties from the SED fitting

The best fit parameters and SED model of J0916a are shown in Table 2 and Figure 8, respectively. The resulting reduced $\chi^{2}=0.79$ suggests that the model with relatively large number of free parameters (8) does not overfit the observation with limited number of data points (11).

The total IR luminosity between 1 $\mu$m and 1000 $\mu$m in the rest frame reaches $L_{\rm IR}$ $=6.13\pm 0.25\times 10^{12}$ L*⊙. 94 % of the IR luminosity comes from the re-emission of the interstellar dust heated by star formation, and the component dominates the FIR peak. Using the empirical relationship of SFR (M⊙* yr*-1*) $=4.5\times 10^{-44}\ L_{\rm IR}^{\rm dust}$ (erg s*-1*) (Kennicutt, 1998), the SFR is estimated to be 990 M*⊙* yr*-1*. With the estimation of stellar mass of $9.46\times 10^{10}$ M*⊙*, J0916a lies above the star-forming main sequence at $0.0<z<0.7$ (Tasca et al., 2015), and the deviation indicates strong starburst activity in this object. The reduced $\chi^{2}$ for the cases where $\alpha_{\rm SF}$ equals to 1.125, 1.625, and 2.125 are 0.79, 4.09, and 22.87, respectively, which suggests that J0916a is similar with the ‘warm’ ULIRG samples with relatively large grey-body temperature ($\sim 40$ K, Symeonidis et al. (2011)).

In CIGALE the contribution of AGN component (green solid curve in Figure 8) is defined as $f^{\rm\scriptscriptstyle AGN}_{\rm IR}\equiv L_{\rm IR}^{\rm\scriptscriptstyle AGN}/L_{\rm IR}^{\rm tot}$, where $L_{\rm IR}^{\rm\scriptscriptstyle AGN}=L_{\rm thermal}^{\rm\scriptscriptstyle AGN}$ is the thermal emission of dusty torus, and the fraction is estimated to be $6.27\pm 0.66\,\%$. We can also estimate the AGN contribution using the bolometric luminosity of the AGN, which is calculated as the sum of the luminosities of the three components in the model of Fritz et al. (2006):

[TABLE]

where $L_{\rm central}^{\rm\scriptscriptstyle AGN}$ is the isotropic luminosity of the central energy source integrated in the range of 0.001–20 $\mu$m, $L_{\rm scattering}^{\rm\scriptscriptstyle AGN}$ and $L_{\rm thermal}^{\rm\scriptscriptstyle AGN}$ are luminosities of scattering and thermal emission from dusty torus, respectively, integrated in the range of 1–1000 $\mu$m. In the best fit SED of J0916a, the AGN emission only dominates the MIR band and does not extend to the range shorter than 1 $\mu$m, indicating that the radiation directly from the central accretion disc is obscured in the UV-optical range. The obscuration of AGN is consistent with the extended optical morphology (see Figure 1). The bolometric luminosity of the host galaxy can be estimated as:

[TABLE]

Therefore, we obtain $f^{\rm\scriptscriptstyle AGN}_{\rm bol}\equiv L_{\rm bol}^{\rm\scriptscriptstyle AGN}/(L_{\rm bol}^{\rm\scriptscriptstyle AGN}+L_{\rm bol}^{\rm host})=$ $6.34\pm 0.65\,\%$, which is consistent with $f^{\rm\scriptscriptstyle AGN}_{\rm IR}$. The bolometric luminosity of obscured AGN is assumed to be the sum of total IR luminosity and X-ray luminosity in several works (Vasudevan et al., 2010; Lusso et al., 2011). We do not have the X-ray constraint on the AGN activity of J0916a. However, the intrinsic 2–10 keV to AGN bolometric luminosity is 0.03 %–0.81 % in the AGN hosted ULIRGs (Teng et al., 2015) or 4 %–8 % in obscured AGNs (Lusso et al., 2011). Thus the lack of X-ray observation will not affect the estimation of $L_{\rm bol}^{\rm\scriptscriptstyle AGN}$ of J0916a.

The resulting contribution of AGN is significantly smaller than the average AGN fraction observed among ULIRGs (35 %–40 %, Veilleux et al. (2009); 40 %–50 %, Goto et al. (2010); 20 %, Ichikawa et al. (2014); 19 %, Małek et al. (2017)) and even smaller than that of LIRGs (12 %, Alonso-Herrero et al. (2012); 10 %, Buat et al. (2015); 12 %, Małek et al. (2017)). In order to test the necessity of the AGN component in the SED fitting, we perform the same fitting procedure with AGN contribution fixed to be 0 %, and the result is shown in Figure 9. The reduced $\chi^{2}$ significantly increase to 7.50, which suggests that the AGN component is necessary to reproduce the MIR photometry.

5 Discussion

5.1 Outflow structure and ionization source

It is commonly found that highly ionized gas, e.g., [Oiii] and [Neiii], tends to show strong outflows. J0916a shows extreme outflow not only in the [Oiii] emission line, but also in the [Oii] emission line. In order to compare the broad line widths of both [Oiii] and [Oii] emission lines with galaxies in the local universe, we collect the samples in the literature with both [Oiii] and [Oii] measurements. We examine the line widths of 1.8 million SDSS DR7 galaxies at $z\sim 0.1$, whose absorption and emission profiles are measured by the MPA-JHU group666https://wwwmpa.mpa-garching.mpg.de/SDSS/DR7/. Only the objects with $S/N>7$ for both [Oiii] and [Oii] emission lines are selected, and their velocity distributions ($w80$) are shown in Figure 10 with contours. We separate the sample into star-forming galaxies (SFGs) and Seyfert 2 galaxies (Sy2s) following the classification of Brinchmann et al. (2004) with the “Baldwin, Phillips & Terlevich” (BPT) diagram (Baldwin et al., 1981; Osterbrock & Ferland, 2006). The distributions of SFGs and Sy2s show peak at 100 km s*-1* and 200 km s*-1*, respectively. J0916a lies far beyond the velocity range of both SFG and Sy2 samples in the local universe.

Recently several works have been presented on strong outflows detected in dust obscured quasars / galaxies at high redshifts. Perna et al. (2015) and Zakamska et al. (2016) reported the discovery of very broad (FWHM $=1500\textrm{--}3000$ km s*-1*) and strongly blueshifted [Oiii] emission line in the spectra of obscured quasars at $z\sim 2$. Toba et al. (2017) analyzed the ionized gas properties for 36 IR-bright DOGs at $0.3<z<1.0$, of which 25 objects are classified as ULIRGs/HyLIRGs, and found that the IR-bright DOGs show relatively strong outflows compared to Sy2s at $z<0.3$. The dispersions of [Oiii] and [Oii] emission lines of these objects are also shown in Figure 10. In order to compare with the velocity dispersion measured in the same method, we obtain SDSS spectra of the DOG sample in Toba et al. (2017) and re-preform spectral fitting for the DOGs with significant detections of both [Oiii] and [Oii] lines to measure the $w80$ dispersions. The powerful outflow of J0916a characterized by the broad ionized oxygen lines is among the fastest ones compared with the ULIRGs/DOGs at intermediate redshifts, and even comparable to the luminous high-redshift QSOs.

In order to examine the spatial extent of the ionized outflow, we perform spectral fitting to every line of the stellar continuum subtracted long-slit spectroscopy image shown in Figure 7. The pixel sampling of $0.2^{\prime\prime}$ corresponds to a physical scale of 1.3 kpc at the object. The radial distributions of flux (left), velocity shifts (middle) and FWHM (right) for each emission line are shown in Figure 11. The [Oii] line shows a more extended structure than stellar continuum, while [Oiii] and H$\alpha$ emission lines show relatively compact flux distribution. The velocity profiles show an extent to about 5 kpc with velocity shift larger than 500 km s*-1* and FWHM larger than 1000 km s*-1*, especially in the North-East side, indicating the outflow is extended up to 5 kpc scale. However, there also exists the possibility that the observed [Oiii] line with high-velocity at 5 kpc comes from the smeared-out flux of an unresoved core component, e.g., compact star-forming region or AGN NLR, by the PSF, which could lead to an overestimation of the size of outflow.

In order to evaluate a possible effect of an unresolved component on the estimation of the size of the outflow region, firstly we plot the flux radial distribution per velocity window with a width of 300 km s*-1* (Figure 12) to check if high-velocity [Oiii] line can be detected in the outskirt region. The unresolved core component is assumed to follow the profile of PSF and contribute to 100% of the observed [Oiii] flux in the central position. The residuals of observed [Oiii] flux minus the PSF-convolved core component (green solid curves in Figure 12) can be used to trace the extended component. We fit the residuals using one Gaussian profile on each of the South-West and North-East side, and find that the extended component with a velocity shift of 600–900 km s*-1* can be detected to about 3.5 kpc with S/N $\sim 7$. The noise ($\sigma$) is estimated to be the scatter of the off-source region in the 2D spectroscopy image. If we take $5\sigma$ as the detection limit, then the maximum radius of the extended component is 4–5 kpc with a velocity shift of 600–900 km s*-1*.

Assuming the original observed spectrum of the central pixel as the representative spectrum of the unresolved core component, we further subtract the PSF-convolved core component from the spectroscopy 2D image. The flux of the central pixel is entirely reduced. The spectroscopy 2D image of extended [Oii], [Oiii], and H$\alpha$-[Nii] emission lines after removal of the PSF-convolved core component are shown in Figure 13. We re-perform spectral fitting per pixel along the long-slit direction, and the measured velocity shifts and FWHM of [Oiii] and [Oii] lines are shown in Figure 14. The results indicate that, even after the removal of the PSF-convolved core component, the [Oiii] line with high velocity shift ($\sim$ 800 km s*-1*) and FWHM ($\sim$ 1200 km s*-1*) can still be detected to 4 kpc scale on the North-East side. As for the South-West side, the velocity shift of [Oiii] decreases from $\sim$ 600 km s*-1* to $\sim$ 300 km s*-1* at 2–4 kpc. The difference between the velocity shift on North-East side and South-West side would imply that the outflow is within an inclined structure. Although the velocity shift of [Oiii] is lower on the South-West side, the FWHM of [Oiii] is still over 1000 km s*-1*up to 4 kpc. Figure 15 shows the [Oiii] spectra of the unresolved core component (left), and the extended spectra after removal of the core component at 4 kpc on both sides (middle and right), respectively.

The more extended flux distribution as well as the lower ionization potential of [Oii] suggest that the outflow shown by [Oii] emission line can be extended farther than [Oiii] line. However, Figure 14 shows that the profiles of velocity shifts and FWHM of [Oii] line is more compact than those of [Oiii]. Furthermore, the velocity FWHM of extended [Oii] gas at 3–4 kpc decreases from 1000 km s*-1to 600 km s-1*after the removal of unresolved core component, indicating that the highest-velocity outflowing [Oii] gas is dominated by the core spectrum. Therefore, although the [Oii] gas is more spatially extended, the compact velocity profiles and core-dominated high-velocity outflowing gas imply that the extended [Oii] gas in the outer region of the galaxy has been less affected by the outflow than [Oiii] emitting gas.

The [Nii] / H$\alpha$ and [Oiii] / H$\beta$ emission line ratios within the outflow region ($r<4$  kpc) are higher than those of outside region. If we plot the line ratios of J0916a on the BPT diagram, then all the data points within outflow region are located on AGN (Seyfert 2) sequence (Figure 16). It is widely accepted that the high-velocity outflows shown by emission lines with high-IP, e.g., [Oiii] and [Neiii], are driven by AGN radiation or wind. On the other hand, the low-IP (13.62 eV) [Oii] line can also be significantly contributed from the star-forming regions, thus it is interesting to examine whether the fast [Oii] outflows are associated with the star formation activity. We use the MPA-JHU SDSS galaxy catalog to find the possible ionization source of fast [Oii] outflows observed among them. The distributions of galaxies with low and high FWHM of [Oii] emission line on the BPT diagram are shown in Figure 16. Galaxies with FWHM larger than 600 km s*-1* tend to show higher [Nii] / H$\alpha$ and [Oiii] / H$\beta$ line ratios, and mostly lie on the AGN sequence, suggesting that the fast [Oii] outflowing gas is ionized by AGN radiation.

5.2 Mass outflow and energy ejection rates

In the recent works, the mass of the ionized outflowing gas has been estimated by either Balmer lines (H$\alpha$, e.g., Soto et al. (2012), Arribas et al. (2014); H$\beta$, e.g., Harrison et al. (2014), Carniani et al. (2015)) or [Oiii] 5007Å line (e.g., Cano-Díaz et al. (2012); Kakkad et al. (2016)). Since the H$\beta$ emission of J0916a is weak, we estimate the mass of the ionized outflowing gas of J0916a using H$\alpha$ and [Oiii] lines. Following Carniani et al. (2015), the mass of the ionized gas can be expressed as:

[TABLE]

where $m_{\rm p}$ is the proton mass; $n_{\rm He{\scriptscriptstyle III}}$, $n_{\rm H{\scriptscriptstyle II}}$, and $n_{\rm e}$ are the number density of helium ions, hydrogen ions, and electrons, respectively; $V$ is the volume of the considered region and $f_{\rm H{\scriptscriptstyle II}}$ the volume filling factor of the ionized hydrogen gas. With Equation 4 and the luminosity of H$\alpha$ emission line:

[TABLE]

the mass of the ionized gas can be estimated as:

[TABLE]

In the derivation we assume (a) helium gas is fully ionized and $n_{\rm He{\scriptscriptstyle III}}=0.1\times n_{\rm H{\scriptscriptstyle II}}$, thus $n_{\rm e}=n_{\rm H{\scriptscriptstyle II}}+2n_{\rm He{\scriptscriptstyle III}}=1.2\times n_{\rm H{\scriptscriptstyle II}}$; (b) $\langle n_{\rm e}^{2}\rangle=\langle n_{\rm e}\rangle^{2}$, i.e., all of the ionized gas clouds have the same density (Cano-Díaz et al., 2012). The emissivity $j_{\rm H\alpha}=3.536\times 10^{-25}$ erg s*-1* cm*-3* is calculated using PyNeb (Luridiana et al., 2015) under the typical temperature ($T_{\rm e}=10^{4}$ K) and electron density ($n_{\rm e}=100$ cm*-3*) of Hii regions and AGN NLRs. The mass of the ionized gas can be similarly derived from [Oiii] emission line, whose luminosity can be given by:

[TABLE]

where $n_{\rm O\scriptscriptstyle III}$ is the density of $\mathrm{O^{2+}}$ ions and $n_{\rm Oion}$ the total density of oxygen ions. In addition to the assumptions (a)-(b) adopted in the calculation with H$\alpha$, following Cano-Díaz et al. (2012) and Carniani et al. (2015) we also assume (c) $n_{\rm Oion}=n_{\rm O\scriptscriptstyle III}$, i.e., most of the oxygen ions are in its doubly ionized form777If we consider the observed flux ratio of [Oiii] and [Oii]  the actual $n_{\rm O\scriptscriptstyle III}/n_{\rm Oion}$ is approximately 0.7, which will result in a larger estimation of ionized gas mass. In order to compare the results with those in the literature (e.g., Harrison et al. (2012); Kakkad et al. (2016); Zakamska et al. (2016)), we keep $n_{\rm Oion}=n_{\rm O\scriptscriptstyle III}$ in the proceeding calculations.; (d) $n_{\rm Oion}/n_{\rm H{\scriptscriptstyle II}}=n_{\rm O}/n_{\rm H}$, i.e., oxygen and hydrogen gas have the same ionization degree, and $n_{\rm O}/n_{\rm H}=[n_{\rm O}/n_{\rm H}]_{\odot}=730\pm 100\ \mathrm{parts\ per\ million}$, where $[n_{\rm O}/n_{\rm H}]_{\odot}$ is the solar oxygen abundance (Centeno & Socas-Navarro, 2008); (e) $f_{\rm O{\scriptscriptstyle III}}V=f_{\rm H{\scriptscriptstyle II}}V$, i.e., the regions from which oxygen and hydrogen lines are emitted have the same size. The emissivity $j_{\rm[O{\scriptscriptstyle III}]}=3.497\times 10^{-21}$ erg s*-1* cm*-3* is calculated with PyNeb under the same condition as that assumed for H$\alpha$ ($n_{\rm e}=100$ cm*-3*, $T_{\rm e}=10^{4}$ K). Combining Equation 4 and Equation 7, we obtain the mass of the ionized gas as:

[TABLE]

The mass of the ionized outflowing gas is proportional to the intrinsic luminosity of outflowing components of H$\alpha$ and [Oiii], and inversely proportional to the electron density. In order to obtain the intrinsic luminosity, the amount of dust extinction need to be determined. The color excess $E(B-V)$ can be estimated from the H$\alpha$ / H$\beta$ emission line ratio through the formula:

[TABLE]

where the attenuation curve, $k_{\lambda}=A_{\lambda}/E(B-V)$, which is calculated empirically for star-forming galaxies, resulting in $k_{\rm H\alpha}=3.33$ and $k_{\rm H\beta}=4.60$ (Calzetti et al., 2000). We assume the intrinsic value of the Balmer decrement, $L_{\rm H\alpha,int}/L_{\rm H\beta,int}$, of 2.86, which corresponds to the electron temperature of $10^{4}$ K and electron density of 100 cm*-3* assuming Case B recombination (Osterbrock & Ferland, 2006). The ratio is commonly assumed for star-forming galaxies in the literature (Groves et al., 2012; Domínguez et al., 2013). Assuming the intrinsic value, we derive the color excess as $E(B-V)=1.01\pm 0.30$ from the narrow components of H$\alpha$ and H$\beta$, which is close to the typical amount of dust attenuation in the local U/LIRGs (Veilleux et al., 1995; Alonso-Herrero et al., 2006; García-Marín et al., 2009) and DOGs (Hwang & Geller, 2013). The SED fitting yields higher color excess $E(B-V)_{\rm YSP}=1.49\pm 0.10$ for the young stellar population. The structure of the dust extinction in ULIRGs was found to be patchy and clumpy (García-Marín et al., 2009; Piqueras López et al., 2013), and $E(B-V)_{\rm YSP}=1.49$ possibly reflects the dust attenuation in dense molecular clouds where the young stars are formed. Since the relationship between the dust extinction of the stellar continuum and emission lines has not been concluded yet (Puglisi et al., 2016; Lo Faro et al., 2017), we adopt $E(B-V)=1.01$ from the Balmer decrement to correct for the dust extinction of the outflowing gas, with the formula:

[TABLE]

Considering the projection effect in the line of sight, and for the comparison to the results in the literature, the total luminosity of [Oiii] including both of the unresolved core component and extended component is used to estimate the mass of the outflowing gas 888The integrated [Oiii] flux at 2.5–6.5 kpc contributes to approximately 35% of the total [Oiii] flux of the galaxy. The estimated gas mass will decrease by 0.45 dex if only the gas in the spatially extended component is considered in the calculation.. The luminosity of broad H$\alpha$ is obtained from the spectral fitting decomposition with both of the broad H$\alpha$ and broad [Nii] components (see discussion in Section 3). The [Sii] 6716Å 6731Å doublet of J0916a lies in the observed wavelength range where the contamination from night sky lines is severe, thus the electron density is assumed to be $n_{\rm e}=100$ cm*-3*, which is a conventionally employed value for the estimation of ionized outflowing gas mass (Liu et al., 2013; Kakkad et al., 2016; Toba et al., 2017). With the dust extinction corrected luminosity and the assumption of $n_{\rm e}=100$ cm*-3*, the mass of the ionized outflowing gas in J0916a is estimated to be $M_{\rm out,H\alpha}=4.7\times 10^{8}$ M*⊙* using the broad component of H$\alpha$, or $M_{\rm out,[O{\scriptscriptstyle III}]}=3.1\times 10^{8}$ M*⊙* using the [Oiii] emission line. The $M_{\rm out,H\alpha}/M_{\rm out,[O{\scriptscriptstyle III}]}\sim 1.5$ is consistent with the average ratio ($\sim 3$) of AGN outflows (Fiore et al., 2017). The difference between $M_{\rm out,H\alpha}$ and $M_{\rm out,[O{\scriptscriptstyle III}]}$ can be due to the assumptions employed in Equation 8, e.g., $n_{\rm Oion}=n_{\rm O\scriptscriptstyle III}$ possibly results in the underestimation of $M_{\rm out,[O{\scriptscriptstyle III}]}$.

In gas clouds with constant density, the mass outflow rate can be derived as $\dot{M}_{\rm out}=M_{\rm out}\times\dot{V}/V=M_{\rm out}v_{\rm out}\times A/V$, where $v_{\rm out}=w_{80}/1.3$ is the bulk outflow velocity (Liu et al., 2013), and $A/V$ is the surface-area-to-volume ratio. Assuming the ionized outflowing gas in a spherically symmetric sector (Fiore et al., 2017), we have $A/V=\Omega R_{\rm out}^{2}/(\Omega R_{\rm out}^{3}/3)=3/R_{\rm out}$, where $\Omega$ is the opening angle and the $R_{\rm out}$ is the radius of outflow region. We adopt $R_{\rm out}\sim 4$ kpc, which is estimated from the core-component-removed spectroscopy 2D image at $5\sigma$ level. Then the mass outflow rate $\dot{M}_{\rm out}$ and energy ejection rate $\dot{E}_{\rm out}$ can be derived as:

[TABLE]

with the results of $560$ M*⊙* yr*-1* and $4.1\times 10^{44}$ erg s*-1* using broad H$\alpha$ line, or $370$ M*⊙* yr*-1* and $2.7\times 10^{44}$ erg s*-1* using [Oiii] line, respectively. Those values reflect the instantaneous feedback effect from AGN and / or star formation activity in the radius $R_{\rm out}$ at the observed time.

The estimation of $\dot{M}_{\rm out}$ and $\dot{E}_{\rm out}$ is based on a few assumptions and the results are dependent on the definition of outflow velocity and kinetic power, the size of outflow region, the correction for the PSF smearing effect, the correction for dust extinction, as well as the electron density. The derived $\dot{E}_{\rm out}$ are consistent with those estimated with another conventional definition of kinetic power $\dot{E}_{\rm out}=\dot{M}_{\rm out}\left(\Delta v^{2}+3\sigma^{2}\right)/2$ (Holt et al., 2006; Harrison et al., 2014), where $\Delta v=(v_{05}+v_{95})/2$ ($v_{05}$ and $v_{95}$ are defined from Equation 1) and $\sigma=w_{80}/2.355$, within a factor of 0.6. Bae et al. (2017) proposed an outflow velocity estimate as $v_{\rm out}^{\star}\equiv 2\times\sqrt{v^{2}+\sigma^{2}}$, where $v$ and $\sigma$ are the flux-weighted averaged shift and dispersion of the emission lines. If we calculate $\dot{M}_{\rm out}$ and $\dot{E}_{\rm out}$ with the definition of $v_{\rm out}^{\star}$, the results increase by a factor of 1.5. $\dot{M}_{\rm out}$ and $\dot{E}_{\rm out}$ will decrease by 0.45 dex if only the gas in the spatially extended component is considered in the calculation. The uncertainty in the extinction correction using $E(B-V)=1.01\pm 0.30$ can lead to the variations of $\dot{M}_{\rm out}$ and $\dot{E}_{\rm out}$ by 0.4 and 0.5 dex for H$\alpha$ and [Oiii] lines, respectively. The assumption of the electron density, which is thought as the main source of uncertainties in the calculation of ionized gas mass and hence the kinetic power of the outflow, was widely discussed in the literature (Rodríguez Zaurín et al., 2013; Fiore et al., 2017; Harrison et al., 2018). The typical value of $n_{\rm e}$ for the AGN NLR is 50–1500 cm*-3* (Peterson, 1997). 100–500 cm*-3* was assumed in several papers (Cano-Díaz et al., 2012; Harrison et al., 2014; Carniani et al., 2015; Kakkad et al., 2016; Fiore et al., 2017; Toba et al., 2017) and roughly consistent with the result of [Sii] ratio of a luminous obscured quasar (Perna et al., 2015). Rupke & Veilleux (2013) and Liu et al. (2013) employed a lower density of 10 cm*-3* for ionized outflowing gas clouds in a thin shell, which was found in the superbubble in NGC 3079 (Cecil et al., 2001). On the contrary, Holt et al. (2011), Rose et al. (2018), and Kawaguchi et al. (2018) reported that the density of outflowing gas can reach as high as $10^{3}$ to $10^{5}$ cm*-3*. We have no direct constraint on $n_{\rm e}$ for J0916a. However, the critical density of [Oii] ($\sim 10^{4}$ cm*-3*, Osterbrock & Ferland (2006)) can be considered as a proxy for the upper limit of $n_{\rm e}$. The large uncertainty of $n_{\rm e}$ can result in the uncertainty of the estimated gas mass and kinetic power of the outflow by one to two orders of magnitude.

A method assuming an energy-conserving bubble in a uniform medium was also widely used to estimate the maximum of energy ejection rates (Heckman et al., 1990; Veilleux et al., 2005; Nesvadba et al., 2006; Harrison et al., 2012, 2014). In the scenario (Castor et al., 1975; Weaver et al., 1977), the fast inside wind induced by the central engine release energy into a hot ($T\geq 10^{6}$ K), low-density, and coronal shocked gas bubble, where the radiative cooling loss is insignificant, by a constant rate, $L_{\rm wind}$. The adiabatically expanding shocked bubble sweeps up a warm ($T\sim 10^{4}$ K), dense, and thin shell into the ISM in the host galaxy. The mechanical luminosity of the inside AGN or stellar wind, $L_{\rm wind}$, can be estimated with the formula:

[TABLE]

where $n_{0}$ is the gas density of the un-disturbed ambient ISM, and we adopt $n_{0}=0.5$ cm*-3* following Nesvadba et al. (2006) and Harrison et al. (2014); $R_{\rm out}$ and $v_{\rm out}$ are the same as defined in Equation 11. According to the calculation of Weaver et al. (1977), 55 % of the total feedback energy is released into the swept-up gas region, and 40 %–70 % of the shell energy is converted to the kinetic energy of the swept-up gas as it collapse into a thin shell through radiative cooling. Therefore the kinetic power of the swept-up ISM corresponds to 22 %–38 % of the feedback energy, i.e., $\dot{E}_{\rm out}\simeq 0.30\times L_{\rm wind}$. For the sake of clarity, the kinetic energy ejection rates derived from the uniform spherical sector model and the energy-conserving bubble model, are named $\dot{E}_{\rm k,out}^{\rm sph}$ and $\dot{E}_{\rm k,out}^{\rm bub}$, respectively. The $\dot{E}_{\rm k,out}^{\rm bub}$ from [Oiii] lines is $2.6\times 10^{45}$ erg s*-1*, which should be considered as the upper limit of the kinetic energy ejection rate because in the scenario it is assumed that all the ambient ISM within $R_{\rm out}$ is entrained in the outflow and no filling factor is employed to reflect the clumpiness of the ambient gas. If we adopt filling factor of 0.2, which is a typical value for warm neutral medium (Hi) and warm ionized medium (Hii) of the Galaxy (Brinks, 1990), the estimation of $\dot{E}_{\rm k,out}^{\rm bub}$ will decrease by approximately one order of magnitude.

The derived mass outflow rate (400–500 M*⊙* yr*-1*, from the uniform spherical sector model) is comparable to the star formation rate (990 M*⊙* yr*-1*), leading to a mass loading factor of $\eta=\dot{M}_{\rm out}/\textrm{SFR}\sim 0.5$, which informs the efficiency of the outflow to remove gas from star-forming region. The co-existence of the strong outflow and intense star formation indicates that the feedback of the outflow has not severely affect the star-forming region of the galaxy, although the outflow would be extend to the radius of 4 kpc. The probable explanation can be that the axis of the outflow is perpendicular to the disc of the host galaxy. It is also possible that the galaxy stays in the intermediate stage in which the feedback just becomes effective and begins to sweep out the ISM in the entire galaxy, i.e., on-set of a quenching of star formation.

In order to discuss the dominant power source of the outflow further, we can compare the energy input power from the AGN and star formation to the energy ejection rate (Cicone et al., 2014; Kakkad et al., 2016; Harrison et al., 2018). The coupling efficiency $\epsilon_{\rm f}$ is defined as the fraction of energy from AGN or star formation that couples to the surrounding gas. In the simulations of galaxy evolution with AGN feedback, $\epsilon_{\rm f,AGN}^{\rm m}$ is assumed to be from 0.5 % (assuming hot gas fraction of 10 %, Hopkins & Elvis (2010)) to 20 % (employed for AGN in the low-accretion mode, Weinberger et al. (2017)), with a typical value of 5 %, which is estimated in DiMatteo et al. (2005) to reproduce the normalization of the local $M_{\rm BH}$–$\sigma$ relation. As for feedback from star formation activity, the power provided by a starburst can be estimated with the formula, $P_{\rm SF}$ (erg s*-1*) $=7\times 10^{41}$ SFR (M*⊙* yr*-1*), by assuming that gas has a solar-metallicity and the mass-loss rate and mechanical luminosity are constant beyond 40 Myr (Veilleux et al., 2005). Combining with the empirical relationship of SFR (M*⊙* yr*-1*) $=4.5\times 10^{-44}\ L_{\rm IR}^{\rm dust}$ (erg s*-1*) (Kennicutt, 1998), we obtain $\epsilon_{\rm f,SF}^{\rm m}=P_{\rm SF}/L_{\rm IR}^{\rm dust}=3.15\,\%$. From the SED fitting, we obtain the AGN bolometric luminosity, $L_{\rm AGN}=1.78\times 10^{45}$ erg s*-1*, as well as the star formation contributed IR luminosity, $L_{\rm SF}=L_{\rm IR}^{\rm dust}=2.20\times 10^{46}$ erg s*-1*. Therefore the observed kinetic coupling efficiency is $\dot{E}_{\rm k,out}^{\rm sph}/L_{\rm AGN}\sim 20\,\%$ and $\dot{E}_{\rm k,out}^{\rm sph}/L_{\rm SF}\sim 2\,\%$, respectively. Hence, both AGN and star formation are sufficient to drive the strong outflow, but larger coupling efficiency would be required for the AGN feedback case.

In order to compare the energy ejection rate of J0916a to that reported in the literature, we collect the outflow samples from the local U/LIRGs and QSOs (Soto et al. (2012); Arribas et al. (2014); Harrison et al. (2014)), the DOGs and QSOs at intermediated redshifts (Liu et al., 2013; Brusa et al., 2015; Kakkad et al., 2016; Toba et al., 2017), to the ULIRGs and QSOs at high redshifts (Harrison et al. (2012); Carniani et al. (2015); Zakamska et al. (2016); Bischetti et al. (2017)). In majority of the samples, the $L_{\rm AGN}$ is estimated using multi-band SED fitting, except for Arribas et al. (2014)999The AGN contributions to the total luminosity of the galaxies in the sample of Arribas et al. (2014) were not presented, we apply an average AGN fraction of 20 % as mentioned in Section 4.2 of Arribas et al. (2014)., Liu et al. (2013)101010Liu et al. (2013) estimated $L_{\rm AGN}$ using the rest 12 $\mu$m luminosity, $\nu L_{\rm 12\mu m}$, and the bolometric correction of 9 from Richards et al. (2009)., and Carniani et al. (2015)111111Carniani et al. (2015) estimated $L_{\rm AGN}$ using the empirical relation $L_{\rm AGN}\simeq 6\times\lambda L_{5100}$.. For the galaxies in all of the reference samples, the masses of ionized outflowing gas are estimated using Equation 6 for Balmer lines or Equation 8 for [Oiii] line, and then the energy ejection rates, $\dot{E}_{\rm k,out}^{\rm sph}$ and $\dot{E}_{\rm k,out}^{\rm bub}$, are re-calculated using Equation 11 and Equation 12, respectively. For illustrative purposes, we assume that in the uniform spherical sector model, all of the galaxies have the same electron density ($n_{\rm e}=100$ cm*-3*) and electron temperature ($T_{\rm e}=10^{4}$ K). In the energy-conserving bubble model, we adopt the kinetic fraction of feedback energy, i.e., $\dot{E}_{\rm out}\simeq 0.30\times L_{\rm wind}$, and employ a filling factor of 0.2 to reflect the clumpiness of the ambient gas. The results are also shown in Figure 17.

J0916a shows one of the strongest outflow among the galaxies at $z<1.6$. However, if we consider the bolometric luminosity of AGN ($L_{\rm AGN}$) from the SED fitting, the $L_{\rm AGN}$ of J0916a corresponds to only 1 %–10 % of $L_{\rm AGN}$ of the galaxies with similar energy ejection rate, which would suggest the largest coupling ratio ($\dot{E}_{\rm out}/L_{\rm AGN}$) among the plotted objects in Figure 17. In order to investigate the status of the AGN activity, we also estimate the Eddington ratio $\lambda_{\rm Edd}\equiv L_{\rm AGN}/L_{\rm Edd}$ of J0916a. The Eddington luminosity is defined with (e.g., Rybicki & Lightman (1979)):

[TABLE]

As for J0916a, the black hole mass $M_{\rm BH}$ can be estimated from the total stellar mass $M_{\star}$ with an empirical relationship between $M_{\rm BH}$ and $M_{\star}$ in the local universe (Reines & Volonteri, 2015):

[TABLE]

From the SED fitting we obtain the estimation of $M_{\star}=9.46\pm 1.69\times 10^{10}$ M*⊙. Hence the black hole mass and Eddington luminosity are derived as $M_{\rm BH}=2.66\times 10^{7}$ M⊙* and $L_{\rm Edd}=3.35\times 10^{45}$ erg s*-1*, respectively. The Eddington ratio $\lambda_{\rm Edd}=0.6$ of J0916a lies between those of ULIRGs / DOGs, which show a super-Eddington accretion (3.35, Kawakatu et al. (2007)121212Note that Kawakatu et al. (2007) calculated the ratio of $L_{\rm IR}^{\rm tot}/L_{\rm Edd}$ instead of $L_{\rm AGN}/L_{\rm Edd}$. The nuclear starbursts ($<1$ kpc) can contribute a fraction ($<30\,\%$) of the infrared luminosity for their ULIRG sample. ; 2–3, Assef et al. (2015); $\sim$ 3, Toba et al. (2017)), and nearby QSOs, which show a sub-Eddington value (0.09, Kawakatu et al. (2007)).

The extremely strong outflow associated with less luminous AGN in J0916a is puzzling. It is possible that the AGN bolometric luminosity ($L_{\rm AGN}$) of J0916a is underestimated in the SED fitting. We compare the SED-based $L_{\rm AGN}$ and [Oiii]-based $L_{\rm AGN}$ for this ULIRG and the reference sources in Figure 17, and the comparisons are shown in Figure 18. The bolometric correction factor of [Oiii] is 600 for extinction-corrected sources and 3500 for extinction-uncorrected sources (Kauffmann & Heckman, 2009; Lamastra et al., 2009). The luminosity ratio, $\log{(L_{\rm AGN,[O{\scriptscriptstyle III}]}/L_{\rm AGN,SED})}$, is $0.2\pm 0.6$ for the reference sources, and 1.3 for J0916a, indicating the $L_{\rm AGN}$ of J0916a estimated from the extinction-corrected [Oiii] luminosity is about 20 times higher than that estimated from the AGN-contributed IR luminosity. If we assume the [Oiii]-based $L_{\rm AGN}$, then J0916a (shown as open stars in Figure 17) becomes close to the most luminous high-$z$ QSOs with the most powerful outflows in the $\dot{E}_{\rm out}$–$L_{\rm AGN}$ diagram. The Eddington ratio $\lambda_{\rm Edd}\equiv L_{\rm AGN,[O{\scriptscriptstyle III}]}/L_{\rm Edd}\sim 10$ is similar with that of the DOGs/ULIRGs with super-Eddington accretion (Kawakatu et al., 2007; Assef et al., 2015; Toba et al., 2017). The [Oiii]-based $L_{\rm AGN}$ requires that the entire IR luminosity is owing to AGN, which is hard to be reproduced by the current SED fitting templates, where the FIR peak is mainly contributed by star formation heated dust. The hard X-ray observation is necessary to determine the intrinsic power of the possible hidden AGN in J0916a.

Another possible explanation for the large outflow coupling ratio is that J0916a is in the phase where the central AGN is less active than its peak epoch. The low IR luminosity of AGN, which originates from dusty torus in the vicinity (pc-scale) of SMBH, would imply the possibility that AGN lies in a fading status; while the observed extreme [Oiii] and [Oii] outflows would reflect a historical effect of the central engine during its preceding active phase, due to the time-lag between AGN activity in a nuclear region and outflow in an ionization cone (kpc-scale; Harrison (2017)). The traveling timescale of the outflow in J0916a can be estimated with $R_{\rm out}/v_{\rm out}$, which results in $\sim 3$ Myr. Theoretical simulations (Novak et al., 2011; Gabor & Bournaud, 2014) and observational investigations (Hickox et al., 2014) suggest that the accretion rates of SMBHs and AGN luminosities can vary by several orders of magnitude in timescales of $\sim 1$ Myr or less. Additionally, recent works reported a population of AGNs called ‘dying AGNs’, which show AGN signatures in large spatial scale (e.g., radio jets and / or bright [Oiii] line in the kpc-scale NLRs), but lack the features in small scale (e.g., weak or lack of X-ray and / or MIR emission), and imply the transient stage that the central engine was active in the past, but now seems quenched (Schawinski et al. (2010); (Ichikawa et al., 2016, 2017, 2019); Schirmer et al. (2016)). Both AGN variability and AGN quenching can be one of possible mechanisms of the fading AGN scenario of J0916a.

It is also possible that star formation activity at least partially contributes to the strong outflow. While several studies presented that AGNs are required to accelerate high-velocity winds (e.g., $v_{\rm max}=|v_{\rm shift}|+\textrm{FWHM}/2\geq 500$ km s*-1*, Rupke et al. (2005); Westmoquette et al. (2012); Arribas et al. (2014); Harrison et al. (2014)), recently high-velocity outflows traced using Mgii 2796Å 2803Å doublet absorption lines in galaxies with intense star formation activity but weak or no AGN activity have been reported (Diamond-Stanic et al., 2012; Bradshaw et al., 2013; Sell et al., 2014; Heckman & Borthakur, 2016). In these galaxies, the outflow originate from a compact star-forming region with effective radius of several hundred parsecs. Heckman & Borthakur (2016) found that the maximum outflow velocity ($v_{02}$, 10–2500 km s*-1*) correlates with SFR (0.1–600 M*⊙* yr*-1*) and star formation rate density (SFRD, 0.1–5000 M*⊙* yr*-1* kpc*-2*). J0916a also lies on the velocity–SFR relationship, with a fast outflow ($v_{02}\sim 1600$ km s*-1* for [Oiii]) and an intense star formation activity (SFR $=990$ M*⊙* yr*-1*) similar with the most extreme starbursts in the sample of Sell et al. (2014) and Heckman & Borthakur (2016). Thus we can not rule out the possibility that star formation activity contributes to the fast outflow in J0916a, since the size of the star-forming region has not yet been constrained based on the current observations.

6 Conclusion

In the spectroscopic follow-up observations of AKARI-selected FIR-bright optically-faint objects, one ULIRG, J0916a, indicates signatures of an extremely strong outflow in its emission line profiles. Both of high- and low-ionization potential lines, e.g., [Oiii] and [Oii], show large velocity dispersions and shifts in relative to the stellar absorption lines. The velocity dispersions and shifts correspond to one of the fastest outflow among ULIRGs / DOGs at $0.3<z<1.0$, and are comparable to the obscured quasars at $z\sim 2$. After the correction for a maximally possible contribution from an unresolved core component, the long-slit spectroscopy 2D image suggests that the outflow could extend to radius of 4 kpc. However, the co-existence of the strong outflow and vigorous starburst (SFR = 990 M*⊙* yr*-1*) suggests that the star formation has not yet been quenched by the outflow. The $L_{\rm bol}^{\rm\scriptscriptstyle AGN}$ of J0916a estimated from the SED fitting is only 5 %–10 % of $L_{\rm bol}^{\rm\scriptscriptstyle AGN}$ estimated from extinction-corrected [Oiii] luminosity, which leads to a large uncertainty in determinating the status of the AGN hidden in this ULIRG. Further observations are required in order to reveal the properties and origins of the outflow, e.g., hard X-ray observation to directly detect the intrinsic AGN radiation, spectroscopy observation with higher spatial and spectral resolution to determine the electron density in the outflowing gas, integral-field spectroscopy observation to investigate the structure of the outflow, and sub-millimeter observation to determine whether the warm ionized outflow really affect the cold molecular gas reservoirs.

{ack}

This research is based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. This research is based on observations with AKARI, a JAXA project with the participation of ESA. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory / California Institute of Technology, funded by the National Aeronautics and Space Administration. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State / Notre Dame / JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University. MA is supported by JSPS KAKENHI (17H06129). YT is supported by Grant-in-Aid for JSPS Research Fellow (Grant No.18J01050). KI is supported by Program for Establishing a Consortium for the Development of Human Resources in Science and Technology, Japan Science and Technology Agency (JST) and JSPS KAKENHI (18K13584).