Apparent effect of dust extinction on the observed outflow velocity of ionized gas in galaxy mergers

We investigated the effect of a galaxy merger on the ionized gas outflow strength at the AGN origin. We performed a galaxy–merger simulation focusing on the galactic nucleus at a kpc scale in the late stages of galaxy mergers. We assume that the bulge-core ($\sim$kpc) scale structure is not considerably distorted in the early stages of galaxy mergers. The parameters of pre-merger system are reported in Table 1. The pre-merger system comprises a supermassive BH, stellar bulge, and gas disk with a kpc-scale radius. We are interested in the high-velocity outflow observed by Toba et al. (2017). The SMBH masses in their samples were distributed at approximately 10${}^{8}M_{\odot}$; therefore, we set the SMBH mass in our models to 10${}^{8}M_{\odot}.$ The stellar bulge mass was determined using the Magorian relation. Because we focus on gas-rich systems, the mass profile of the stellar bulge is based on star-forming galaxies with a redshift of $\sim$ 2. According to observational studies of star-forming galaxies, we have set the effective radius at 2.2 kpc and the Sersic index at 2.0 for the stellar bulge (Barro et al., 2017; Paulino-Afonso et al., 2017). The distribution function for the stellar bulge was obtained from AGAMA, which provides a self–consistent N-body system (Vasiliev, 2019). The total gas mass is set such that the gas fractions $\mu_{gas}=M_{gas}/(M_{gas}+M_{bul}(r_{xy}<1\,{\rm kpc}))$ are 10% and 30%, respectively, where $r_{xy}$ is the radial length in cylindrical coordinates with its origin at the SMBH. $M_{bul}(r_{xy}<1\,{\rm kpc})$ is the bulge mass within the bulge of a 1–kpc radius. The gas fraction is based on observations of star forming galaxies Erb et al. (2006).

We simulated collisions between identical systems. As the schematic in Figure 1 shows, the systems were separated from each other by 5 kpc and given relative velocities of 100 pc/Myr in the y-direction. We are interested in the late stages of galaxy mergers, where the two galactic nuclei share a common envelope (Stierwalt et al., 2013). Given that the effective radius of the bulge was 2.2 kpc and the initial radius of the disk was 1 kpc, we separated it by 5 kpc to prevent contact and initial mass accretion to the nucleus. In a collisional system, the initial spin angular momentum vector of the gas disk satisfies

where $\theta_{disk}$ denotes the angle between the gas disk of the secondary system and xy-plane. We considered three cases with $\theta_{disk}$ = 0 deg, 30 deg, and 60 deg. In the following discussion, G3T0 denotes the case of $\theta_{disk}=0$ deg for G3 with a gas fraction of 30% and coalescence between identical systems. We simulated six models: G1T0, G1T30, G1T60, G3T0, G3T30, and G3T60 as listed in Table 2.

Figure 2 shows the evolution of the gas density and temperature distribution for G3T0 and G3T60. As the system approaches, the gas is compressed and sink particles are buried in the dense gas region. The mass accretion rate to the sink particle is enhanced in this stage. In addition, the two right columns in Figure 2 show the temperature distribution of the gas, and the rightmost column shows the temperature distribution of the gas particles that received AGN-feedback energy within Myr. The temperature distributions indicate that the hot gas particles are blown vertically up from the disk by the AGN-feedback.

Figure 3 shows the evolution of SMBH binary distance, total bolometric AGN-luminosity, and star formation rate for the three models. The top row of Figure 3 shows that the SMBH binary orbits are similar in all three models with the same gas fraction. The middle row in Figure 3 shows that there was no significant difference in AGN luminosity among the three models with the same gas fraction (models G1 and G3). In all three models, the AGN activity and star formation activity were enhanced around the first and second pericenters of the BH-binary orbit. AGN luminosity and star formation rate tended to be higher with a gas fraction of 30% than with 10%. For a gas fraction of 30%, the AGN bolometric luminosity ranges from $10^{11}L_{\odot}$ to $10^{13}L_{\odot}$, whereas for a gas fraction 10%, the AGN luminosity ranges from $10^{10}L_{\odot}$ to $10^{12}L_{\odot}$.

We used the velocity and velocity-dispersion (VVD) diagrams to investigate the relation between the galaxy merger and ionized outflow velocity. VVD diagrams have been used to investigate the statistical properties of the ionized outflow velocity in AGNs (Woo et al., 2017; Rakshit & Woo, 2018). We constructed the VVD diagrams using the intrinsic velocities of SPH particles derived from the simulated galaxy mergers data to avoid the high computational cost of 3D ionizing emission line pseudo-observations.

In a 3D polar coordinate system centered on each sink particle, we selected photoionized gas particles from those above 8000 K by limiting them to typical [OIII] $\lambda$5007 ionization parameters (-3 $<$ $\log$U $<$ -1).

where $N_{e}$ is the number density of SPH particles at temperatures above 8000 K. The ionization parameter U is defined for an SPH particle as the ratio of the number density of ionizing photons with energy $>$ 13.6 eV to the electron number density $N_{e}$ (e.g., Wada et al., 2018). We consider the radiation emitted from the AGN as the central point source. The AGN spectral energy distribution (SED) ($f_{\nu}=L_{\nu}/4\pi r^{2}$) was obtained from Cloudy’s AGN table (Ferland et al., 2017). The AGN SED was assumed to be

Here, $\alpha_{uv}=-0.5$, $T_{BB}=1.5\times 10^{5}$ K, $kT_{IR}=0.01$ Ryd, $\alpha_{x}$ = -1.0, and b is a constant yielding $\alpha_{ox}$ = -1.4. The AGN bolometric luminosity depends on the mass accretion rate on the sink particle. For the optical depth ($\tau_{\nu}$) expressed in Equation 2, we consider scattering and absorption by dust, photoionization of hydrogen, and Thomson scattering:

where $\sigma_{\rm dust}$ is the extinction cross section of MRN dust (Laor & Draine, 1993), $\sigma_{\rm H}$ is the photoionization cross–section of hydrogen (Osterbrock & Ferland, 2006), and $\sigma_{\rm T}$ is $6.65\times 10^{25}$cm², which is the Thomson scattering cross–section. $N_{{\rm H},agn}$ is the column density from the SMBH to the SPH particle. In addition, we extracted gas particles not bound to the SMBH potential. Because we were interested in the outflow component, we used the radial velocity in a 3D polar coordinate system centered on each sink particle.

From the gas particles selected using the U parameter (-3 $<\log{U}<$ -1), we selected gas particles that were not trapped by an SMBH potential and calculated the line-of-sight velocity and velocity dispersion. To assess whether gas particles are trapped by SMBH potential, the gas particles velocities are calculated by $\sqrt{v^{2}_{r}+C_{s}^{2}}$, where $v_{r}$ is radial velocity from SMBH and $C_{s}$ is sound speed. If the velocity of a gas particle is exceeds $\sqrt{2GM_{BH}/r}$, then its velocity contributes to Equations 5 and 6. We assume that the [OIII] $\lambda$5007 line intensity is proportional to ${n_{H}}^{2}$ because the [OIII] $\lambda$5007 line is the collision–excited line. In the cartesian coordinate system centered on each sink particle, we calculate mean velocity and velocity dispersion as

where ${v_{los}}_{i}$ is the line-of-sight velocity of each SPH particle in a rest frame centered on sink particle and $\tau_{5007}$ equals to $\sigma_{dust}(\lambda=5007{\rm\AA}){N_{{\rm H},los}}(\theta_{0},\phi_{0})$. ${N_{{\rm H},los}}(\theta_{0},\phi_{0})$ is the column density from each SPH particle to r = 6 kpc in the ($\theta_{0},\phi_{0}$) direction. The $\theta_{0}$ is the angle from the rotation axis of the primary system (i.e., z-axis in Figure 1). We use ${N_{{\rm H},los}}(\theta_{0},\phi_{0})$ as the approximate column density only for the line of sight within the solid angle $\Omega=2\pi(1-\cos(\pi/12))$ sr from ($\theta_{0}$,$\phi_{0}$).

We have fixed $\phi_{0}$ = 0 deg (i.e., in the x-axis direction in Figure 1) and calculated for four cases where $\theta_{0}$ is 0 deg, 30 deg, 60 deg, and 90 deg (i.e., covering the range from $\theta$ = -15 deg to 105 deg and $\phi$ -15 deg to 15 deg). In our simulations, the OIII ionization outflow of the AGN origin was on the scale of a few hundred pc or less (see Figure 8), and the gas density distribution in the $\phi$ direction did not change significantly on that scale (see Figure 2). We choose 360 lines of sight within $\Omega=2\pi(1-\cos(\pi/12))$ for each $\theta_{0}$ with equal $\sin{\theta}d\theta d\phi$ in each SMBH. As there are two SMBHs in a model, we plotted $2\times 4\times$360 points per snapshot on the VVD diagram for all models.

To assess the relation between galaxy mergers and ionized outflows accurately, we define the merger phase based on the distance between binary BHs. We are interested in the stage in which the gas disks merge with each other. Because the initial gas disk radius is 1 kpc and BH accretion radius is 6 pc (see Table 1), the merger phase is defined as the distance between the binary BHs between 2 kpc and 12 pc. The reason for setting 2 kpc is that if the distance between SMBHs is less than 2 kpc, two bulge systems with an effective radius of 2.2 kpc strongly gravitationally interact and the gas disk with a radius of 1 kpc is incredibly distorted and the mass accretion rate to the SMBHs increases.

We plotted our models on a VVD diagram for snapshots during the merger phase. However, the following two cases are not plotted in the VVD diagram: 1) when $L_{agn}$ = 0, and 2) when Equation 7 is not satisfied.

where $h_{i}$ is the kernel radius of SPH particles. As the accuracy of SPH particles depends on the local particle density, this condition eliminates snapshots that are heavily influenced by noisy SPH particles. In our calculations, the gravity softening was 3 pc; therefore, we set the threshold at four times the gravity softening (12 pc). For G3T0 model, the merger phase ranged from 13.0 to 30.1 Myr, and snapshots were obtained every 0.1 Myr. Therefore, at the maximum 181 snapshots were plotted on the VVD diagram. $2\times 4\times 360$ lines of sight were selected for each snapshot in all the models. Thus, for the G3T0 model, at the maximum 181$\times$2880 points were plotted on the VVD diagram.

Figure 4 shows VVD diagrams for merger models with gas fractions of 30% and 10%. These VVD diagrams show an ionization outflow component with a velocity dispersion exceeding 500 km s^-1 regardless of the gas fraction, confirming the high-velocity outflow from the galaxy merger simulations. In addition, the model with a gas fraction of 30% is more tilted toward the blue-shifted side than that with a gas fraction of 10%. Indeed, the percentage of samples that were blue-shifted more than -100 km s^-1 was 40.7% for the 30% gas fraction and 30.0% for the 10% gas fraction model. According to Bae & Woo (2016), the tilt in the VVD diagram is produced by dust extinction against the redshift component of the dipole outflow. In this study, we investigated the effect of dust extinction on the line-of-sight average velocity and velocity dispersion (equations 5 and 6) during the late stages of galaxy mergers.

Figure 9 shows the VVD diagrams depicted in Figure 4, classified with respect to the AGN luminosity. The upper row in Figure 9 shows the VVD diagram for a gas fraction of 30%, whereas the bottom row shows the diagram for a gas fraction of 10%. The tilt toward the blueshift of the VVD diagram is not significantly changed against AGN luminosity, compared to the effect of dust extinction in Figure 5. Thus, the AGN luminosity does not significantly affect the tilt of the VVD diagram.

Toba et al. (2017) reported that some IR bright DOGs with fast [OIII] $\lambda$5007 outflows exhibit velocity dispersion and velocity offsets higher than 500 km s^-1. We compare the VVD diagrams of our model with Toba et al. (2017) in Figure 10. Figure 10 shows that our model covers most sources observed in Toba et al. (2017). However, in Toba et al. (2017), a few objects with $\sigma_{v_{los}}\,$= 1000 km s^-1 and $\bar{v}_{los}\,$= -1000 km s^-1 are observed, and the sources was tilted more strongly than in our model. The tilt of the VVD diagram was determined by the dust extinction based on the results of this study. This implies that the high-velocity outflow sources observed in Toba et al. (2017) may be even more obscured than those in our results. This can be reproduced using a galaxy merger model with a higher gas fraction. We were unable to select DOGs from our simulation data because we did not conduct infrared pseudo-observations in this study. For comparison with more detailed observations, both infrared and ionizing emission line pseudo-observations should be performed.

We performed identical galactic-core merger simulations. In our models, we fixed the initial mass of SMBH at 10${}^{8}M_{\odot}$. This is because Toba et al. (2017) suggested that the SMBH mass of IR bright DOGs is approximately 10${}^{8}M_{\odot}$. DOGs that contain even heavier SMBHs than 10${}^{8}M_{\odot}$ are often observed as Hot DOGs (i.e., Wu et al., 2018). In addition to IR bright DOGs, observational studies have suggested that Hot DOGs are also accompanied by very fast ionization flows (e.g., Jun et al., 2020). However, at redshift 2, SMBHs of 10${}^{9}M_{\odot}$ are considered to be approximately 100 times rarer than SMBH at 10${}^{8}M_{\odot}$ (e.g. Rosas-Guevara et al., 2016). Therefore, we consider that an SMBH mass of 10${}^{8}M_{\odot}$ is the most appropriate initial condition for a galaxy merger simulation in considering the high-velocity outflows associated with buried AGNs.

In addition, it is believed that buried AGNs are formed during the late stages of major merger (e.g. Dey & Ndwfs/MIPS Collaboration, 2009; Narayanan et al., 2010; Yamada et al., 2021; Yutani et al., 2022; Dougherty et al., 2023). Our identical major merger simulation in this study is based on a major merger scenario. However, IR bright DOGs could be also formed by a minor merger with a galaxy, accompanied by an SMBH of about 10${}^{7}M_{\odot}$. We suspect that as long as the central nuclei are buried in dust, the apparent effect of dust extinction suggested in this study on the observed outflow velocity would also be the case. The cases of BH masses are not idendical and minor mergers will be examined elsewhere.