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

Naomichi Yutani  Kagoshima University, Graduate School of Science and
Engineering, Kagoshima 890-0065, Japan<EMAIL_ADDRESS>Yoshiki Toba
National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan  Department of Astronomy, Kyoto University, Kitashirakawa-
Oiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan  Academia Sinica Institute of
Astronomy and Astrophysics, 11F of Astronomy-Mathematics Building, AS/NTU,
No.1,
Section 4,Roosevelt Road, Taipei 10617, Taiwan  Research Center for Space and
Cosmic Evolution, Ehime University, 2-5 Bunkyo-cho, Matsuyama, Ehime 790-8577,
Japan  Keiichi Wada Kagoshima University, Graduate School of Science and
Engineering, Kagoshima 890-0065, Japan  Research Center for Space and Cosmic
Evolution, Ehime University, 2-5 Bunkyo-cho, Matsuyama, Ehime 790-8577, Japan
Hokkaido University, Faculty of Science, Sapporo 060-0810, Japan

###### Abstract

In this study, we examine photoionization outflows during the late stages of
galaxy mergers, with a specific focus on the relation between observed
velocity of outflowing gas and the apparent effects of dust extinction. We
used the N-body/smoothed particle hydrodynamics (SPH) code ASURA for galaxy
merger simulations. These simulations concentrated on identical galaxy mergers
featuring supermassive black holes (SMBHs) of 108 M⊙ and gas fractions of 30%
and 10%. From the simulation data, we derived velocity and velocity dispersion
diagrams for the AGN-driven ionized outflowing gas. Our findings show that
high-velocity outflows with velocity dispersions of 500 km s-1 or greater can
be observed in the late stages of galactic mergers. Particularly, in buried
AGNs, both the luminosity-weighted outflow velocity and velocity dispersion
increase owing to the apparent effects of dust extinction. Owing to these
effects, the velocity–velocity dispersion diagrams display a noticeable blue-
shifted tilt in models with higher gas fractions. Crucially, this tilt is not
influenced by the AGN luminosity but emerges from the observational impacts of
dust extinction. Our results imply that the observed high-velocity [OIII]
$\lambda$5007 outflow exceeding 1000 km s-1 in buried AGNs may be linked to
the dust extinction that occurs during the late stages of gas-rich galaxy
mergers.

Supermassive black holes, galaxy mergers, active galaxies, active galactic
nuclei, N-body simulations, hydrodynamical simulations

††software: ASURA (Saitoh et al., 2008, 2009; Saitoh & Makino, 2013), AGAMA
(Vasiliev, 2019)

##  1 Introduction

Galaxy mergers are key phenomena in understanding the evolution of galaxies
and supermassive black holes (SMBHs). Observational studies have suggested
that galaxy mergers are related to the activity of galaxy centers (Sanders &
Mirabel, 1996; Farrah et al., 2002; Fan et al., 2016; Gao et al., 2020).
Theoretically, galaxy mergers can enhance BH accretion and trigger active
galactic nuclei (AGNs), which are often buried in a large amount of gas and
dust during the late stages of galaxy mergers (e.g., Narayanan et al., 2010;
Blecha et al., 2018; Kawaguchi et al., 2020; Yutani et al., 2022). Indeed, Fan
et al. (2016) reported a high merger fraction ($\sim$62%) in heavily obscured
quasars selected using a wide-field infrared survey explorer (WISE). In the
late stages of galaxy mergers, as AGN becomes more active, AGN feedback
becomes stronger. AGN feedback blows away the dust and gas around the AGNs,
causing the buried AGNs to evolve into quasars in the late stages of galaxy
mergers (Hopkins et al., 2008; Dey & Ndwfs/MIPS Collaboration, 2009). This
feedback process is characterized by an ionized gas outflow. Thus, the
properties of ionized outflows during the late stages of galaxy mergers are
crucial for understanding the formation of SMBHs and quasars.

Weak X-ray quasars ($a_{\rm ox}\leq$ -1.7) have been suggested to be
associated with strong highly ionized outflows with high equivalent widths
(Laor & Brandt, 2002; Veilleux et al., 2022). Toba et al. (2017) studied the
ionized gas properties of 36 objects with extreme optical/infrared color (i –
[22])${}_{\rm AB}>$ 7.0 and infrared bright $F_{22\mu{\rm m}}$ $>$ 3.8 mJy
sources (i.e., IR-bright dust-obscured galaxies), selected by the Sloan
digital sky survey and WISE. They showed that most IR-bright dust-obscured
galaxies exhibited larger velocity offsets and larger velocity-dispersions of
[OIII] $\lambda$5007 than those observed in the Type-2 Seyfert galaxies. The
ionized outflow velocities of some IR-bright DOGs exceeded 1000 km s-1.
Theoretical models and observations suggest that IR-bright galaxies are often
in the late stages of galaxy merging (Narayanan et al., 2010; Yamada et al.,
2021; Yutani et al., 2022).

Bae & Woo (2016) studied the velocity and velocity dispersion diagram
(hereafter VVD diagram) of typical AGN using the three-dimensional biconical
outflow model with uniform density. In their models, dust extinction by
razor–thin disks was considered. They showed that dust extinction strongly
affected the ionized outflow velocity integrated into the line of sight. They
found that the extinction by the dust plane increased the line-of-
sight–integrated velocity offset because the redshift component of biconical
outflow was affected by dust extinction. This is an apparent effect, and does
not indicate a change in intrinsic outflow velocity. The apparent effect of
dust extinction is more important for buried AGNs during the late stages of
galaxy mergers.

In addition, Woo et al. (2017) suggested that the bolometric luminosity of
AGNs and velocity of the ionized outflow was correlated. The relation between
AGN luminosity and outflow velocity is controversial owing to the relatively
small number of statistical samples; however, if there is a positive
correlation between the two quantities in the late stages of the galaxy
merger, stronger AGN-feedback may increase the intrinsic velocity of ionized
outflows.

Considering the relation between ionized outflow velocity and the apparent
effects of dust attenuation or AGN activity, strong ionized outflows may
accompanying AGNs in the late stages of galaxy mergers is plausible. However,
AGNs buried in the later stages of galaxy mergers are expected to have higher
scale heights of the dust torus and more complex outflow density distributions
than the typical AGNs considered in Bae & Woo (2016). To understand the
ionized outflow associated with buried AGNs, conducting galaxy merger
simulations with sufficiently fine spatial resolution based on 3D numerical
fluid dynamics simulations is necessary. In this study, we use the
N-body/smoothed particle hydrodynamic (N-body/SPH) simulation code ASURA
(Saitoh et al., 2008, 2009; Saitoh & Makino, 2013) implementing isotropic
thermal AGN-feedback.

The remainder of this paper is organized as follows: In Section 2, we describe
the galaxy merger simulation models and the corresponding results. In Section
3, we summarize the VVD diagram analysis method and corresponding results. In
Section 4, we discuss the high-velocity ionized outflow associated with the
buried AGNs identified in the observations based on our simulation results.
Finally, in Section 5 summarize the results of this study.

Table 1: Initial parameters of the pre-merger system

Namea | ${M_{BH}}^{\rm b}$ | ${M_{bul}}^{\rm c}$ | ${M_{\rm gas}}^{\rm d}$ | $\epsilon^{\rm e}$ | ${r_{acc}}^{\rm f}$ | ${r_{h}}^{\rm g}$ | ${R_{\rm gas}}^{\rm h}$ | ${\Delta M_{bul}}^{\rm i}$ | ${\Delta M_{\rm gas}}^{\rm j}$ | ${\mu_{gas}}^{\rm k}$
---|---|---|---|---|---|---|---|---|---|---
| [$10^{8}M_{\odot}$] | [$10^{11}M_{\odot}$] | [$10^{10}M_{\odot}$] | [pc] | [pc] | [kpc] | [kpc] | [$10^{4}M_{\odot}$] | [$10^{3}M_{\odot}$] |
G1 | 1.0 | 1.0 | 0.26 | 3.0 | 6.0 | 2.2 | 1.0 | 2.5 | 0.65 | 10%
G3 | 1.0 | 1.0 | 1.0 | 3.0 | 6.0 | 2.2 | 1.0 | 2.5 | 2.5 | 30%

Note. — (a) Name of the pre-merger system. (b) Mass of the sink, (c) star, and
(d) gas particles, respectively. (e) Gravitational softening. (f) Accretion
radius of sink particle. (g) Effective radius of stellar bulge Sersic profile.
(h) Outer edge radius of uniform-density gas disk. (i) Mass resolution of star
particles, (j) and gas particles. (k) Gas fraction below 1 kpc.

Figure 1: The schematic of the merger model

##  2 galaxy merger simulations

We performed galaxy–merger simulations using the N-body/SPH code ASURA (Saitoh
et al., 2008, 2009). The dynamics of the interstellar medium in a merger
system were calculated using density-independent fomulation of smoothed
particle hydrodynamics (Saitoh & Makino, 2013). The gravity fields were
calculated using the tree algorithm. The tolerance parameter was set to 0.5.

In this study, star formation and supernova explosions were implemented
according to Saitoh et al. (2008). The star formation rate (SFR) is determined
as $C_{*}\rho_{gas}/t_{ff}$, where $C_{*}$ is a dimensionless star–formation
efficiency parameter set to 0.0099, $\rho_{gas}$ is the local gas density, and
$t_{ff}$ is the free–fall time. SPH particles that satisfy all four of the
following conditions are converted to star particles in a stochastic manner.
(1) The hydrogen number density $n_{\rm H}\,>\,1000\rm{\,cm^{-3}}$. (2) The
gas temperature $T_{gas}\,<\,100\rm{\,K}$. (3)
$\nabla\cdot{\mathbf{v}}_{SPH}\,<\,0$, where ${\mathbf{v}}_{SPH}$ denotes the
velocity of an SPH particle. (4) $\Delta Q\,<\,0$, where $\Delta Q$ denotes
the thermal energy received by an SPH particle in a time step. The initial
mass function of the stellar particles is set to Salpeter (1955). A
probabilistic Type II supernova explosion was introduced based on Okamoto et
al. (2008). Optically thin radiative cooling with solar metallicity is assumed
for a gas at $10\rm{\,K}\,$–$\,10^{8}\rm{\,K}$ (Wada et al., 2009). Far-
ultraviolet radiation and photoelectric heating have also been considered.

SMBHs can acquire the mass of gas particles within an accretion radius
$r_{acc}$ satisfying the following two conditions: 1) The kinetic energy of
SPH particles is smaller than the gravitational energy 2) The angular momentum
of an SPH particle is less than
$J_{acc}\,=\,r_{acc}\sqrt{GM_{BH}/(r_{acc}^{2}\,+\,\epsilon^{2})^{1/2}}$,
where $\epsilon$ is the gravitational softening length. We assume $\epsilon$ =
3.0 pc and $r_{acc}$ = 6.0 pc.

The isotropic AGN-feedback was implemented by providing thermal energy
($\Delta E$ = $\eta_{AGN}\dot{M}_{BH}c^{2}$) weighted by the spline function
to the gas particles surrounding the sink particle. $\dot{M}_{BH}$ is the gas
mass accretion rate of a BH particle per step at $r_{acc}$ and $\eta_{AGN}$ is
a free parameter representing the energy-loading efficiency. In this study,
$\eta_{AGN}$ = 0.2% was assumed based on the discussions reported in Kawaguchi
et al. (2020).

###  2.1 Merger models

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

$L_{z}^{prime}=L_{z}^{second}\cos{\theta_{disk}}-L_{x}^{second}\sin{\theta_{disk}},$
(1)

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.

Table 2: Free parameters of the six models

Model name |  |  | ${\mu_{gas}}^{\rm a}$ |  |  | ${\theta_{disk}}^{\rm b}$ |
---|---|---|---|---|---|---|---
G1T0 |  |  | 10% |  |  | 0 deg |
G1T30 |  |  | 10% |  |  | 30 deg |
G1T60 |  |  | 10% |  |  | 60 deg |
G3T0 |  |  | 30% |  |  | 0 deg |
G3T30 |  |  | 30% |  |  | 30 deg |
G3T60 |  |  | 30% |  |  | 60 deg |

Note. — (a) The gas fraction within 1 kpc, i.e.,
$\mu_{gas}=M_{gas}/(M_{gas}+M_{bul}(r_{xy}<1\,{\rm kpc}))$. (b) Angle between
gas disk of secondary system and xy-plane.

Figure 2: Map of gas density (1st and 2nd columns from the left) and gas
temperature (3rd and 4th columns from the left). The rightmost column shows
only gas particles that received AGN-feedback energy within a Myr. The top
three rows show G3T0, while the bottom three rows show G3T60.

###  2.2 Galaxy merger simulation results

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}$.

Figure 3: The first row shows the distance between binary BHs. The second row
denotes the star formation rage (SFR). The third row shows the total
bolometric luminosity of the two AGNs
($L_{AGNs}=\Sigma\,0.1*\dot{M}_{sink}c^{2}$). In the first to third rows, the
left column is for G1T0, G1T30, and G1T60; the right column is for G3T0,
G3T30, and G3T60.

##  3 Velocity and velocity dispersion diagram

###  3.1 Integral of photoionized outflow velocity

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

$U=\frac{\int^{\infty}_{\nu_{0}}L_{\nu}e^{-\tau_{\nu}}/h\nu d\nu}{4\pi
r^{2}cN_{e}},$ (2)

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

$f_{\nu}=\nu^{\alpha_{uv}}\exp(h\nu/kT_{BB})\exp(-kT_{IR}/h\nu)+b\nu^{\alpha_{x}},$
(3)

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:

$\displaystyle\tau_{\nu}=\sigma_{\nu,dust}N_{{\rm H},agn}+$
$\displaystyle\sigma_{\nu,H}N_{{\rm H},agn}(<8000{\rm K})$
$\displaystyle+\sigma_{T}N_{{\rm H},agn}(\geq 8000{\rm K})$ (4)

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}$cm2, 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

$\displaystyle\overline{v}_{los}(\theta,\phi)$
$\displaystyle=\frac{\Sigma_{i}{n_{H}}_{i}^{2}e^{-{\tau_{5007}}_{i}}{v_{los}}_{i}(\theta,\phi)}{\Sigma_{i}{n_{H}}_{i}^{2}e^{-{\tau_{5007}}_{i}}}$
(5) $\displaystyle\overline{\sigma}^{2}_{v_{los}}(\theta,\phi)$
$\displaystyle=\frac{\Sigma_{i}{n_{H}}_{i}^{2}e^{-{\tau_{5007}}_{i}}{v_{los}}_{i}^{2}(\theta,\phi)}{\Sigma_{i}{n_{H}}_{i}^{2}e^{-{\tau_{5007}}_{i}}}-\overline{v}^{2}_{los}(\theta,\phi),$
(6)

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.

Figure 4: Comparison of VVD diagrams for gas fraction of 30% (upper row), and
gas fraction of 10% (bottom row). The upper panel includes three models,
namely G3T0, G3T30, and G3T60. The lower panel includes three models, namely
G1T0, G1T30, and G1T60. The dotted line in the histogram for each axis
indicates the mean value, and the gray background color indicates the standard
deviation. $N_{snap}$ is the total number of snapshots in the VVD diagram; all
points in the VVD diagram are 360$\times N_{snap}$. $\overline{L}_{AGN}$ is
the mean luminosity of the snapshots used in each panel.  Figure 5: VVD
diagrams for galaxy merger models with $\tau_{5007}$ = 0 in equations 5 and 6
for a gas fraction of 30% (upper panel) and a gas fraction of 10% (bottom
panel), respectively. The gray plot denotes the VVD considering the dust
extinction in Figure 4. The mean value for $\bar{v}_{los}\,$of each
$\sigma_{v_{los}}\,$bin is indicated by a circle ($\tau_{5007}$=0) and square
(gray plot). The dotted line in the histogram for each axis indicates the mean
value, while the gray background color indicates the standard deviation for
galaxy merger models with $\tau_{5007}$ = 0.

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.

$\displaystyle\tilde{h}=\frac{\Sigma_{i}{n_{H}}_{i}^{2}e^{-{\tau_{5007}}_{i}}h_{i}}{\Sigma_{i}{n_{H}}_{i}^{2}e^{-{\tau_{5007}}_{i}}}\leq
12\,{\rm pc},$ (7)

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.

###  3.2 VVD diagrams of gas-rich galaxy merger models

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 6:  VVD diagram with $\tau_{5007}$ = 0 for $v_{z}$ $>$ 0 particles and
$\tau_{5007}$ = $\infty$ for $v_{z}$ $\leq$ 0 with rest frame for SMBH for a
G3T0 with a gas fraction of 30% (upper panel) and a G1T0 with a gas fraction
of 10% (bottom panel), respectively. The gray plot denotes the VVD considering
the dust extinction for all ionized outflow gas particles. The mean value for
$\bar{v}_{los}\,$of each $\sigma_{v_{los}}\,$bin is indicated by a circle for
the colored VVD diagram and square for the gray VVD diagram. The dotted line
in the histogram for each axis indicates the mean value while the gray
background color indicates the standard deviation of the colored VVD diagram.
Figure 7:  VVD diagram with $\tau_{5007}$ = $\infty$ for $v_{z}$ $\leq$ 0 with
rest frame for SMBH for a G3T0 with a gas fraction of 30% (upper panel) and a
G1T0 with a gas fraction of 10% (bottom panel), respectively. Here dust
extinction is considered for $v_{z}$ $>$ 0 particles. The gray plot denotes
the VVD considering the dust extinction for all ionized outflow gas particles.
The mean value for $\bar{v}_{los}\,$of each $\sigma_{v_{los}}\,$bin is
indicated by a circle for the colored VVD diagram and square for the gray VVD
diagram. The dotted line in the histogram for each axis indicates the mean
value while the gray background color indicates the standard deviation for the
colored VVD diagram.

####  3.2.1 Apparent effect of dust extinction

Figure 5 plots the VVD diagram with $\tau_{5007}$ = 0 expressed in Equations 5
and 6. Figure 5 confirms that a symmetric distribution with respect to
$\bar{v}_{los}\,$= 0. Since the radiation emitted from all ionized gases is
optically thin in Figure 5, the redshift component will cancel out the
blueshift component. Figure 6 shows the VVD diagram with $\tau_{5007}$ = 0 for
$v_{z}$ $>$ 0 particles and $\tau_{5007}$ = $\infty$ for $v_{z}$ $\leq$ 0 with
rest frame for SMBH. In Figure 6, for simplicity, VVD diagrams are plotted
from a face-on view (0 $\leq$ $\theta$(deg) $\leq$ 45) for coplanar merger
models (G3T0 and G1T0). Figure 6 shows the effect of dust extinction for the
redshift component on the VVD diagram. We confirmed that the tilt in the VVD
diagram was produced by dust extinction for the redshift component. This
result is consistent with those Bae & Woo (2016).

Figure 8 plots the probability density function (hereafter, PDF) of equation 8
for the six models (i.e., G3T0, G3T30, G3T60, G1T0, G1T30, and G1T60) from
face-on view ($\theta\leq$ 45 deg).

$\displaystyle\tilde{z}=\frac{\Sigma_{i}{n_{H}}_{i}^{2}e^{-{\tau_{5007}}_{i}}z_{i}}{\Sigma_{i}{n_{H}}_{i}^{2}e^{-{\tau_{5007}}_{i}}},$
(8)

where dust extinction was considered (i.e., $\tau_{5007}\neq 0$). Figure 8
shows that the PDF with a gas fraction of 30% exhibits a systematic offset of
z $>$ 0 (observer side). Although for G1T0 with a gas fraction of 10%, the PDF
is symmetric for z $\simeq$ 0, and for G3T0 with a gas fraction of 30%, it is
drastically reduced for z $\lesssim$ 0\. The offset was caused by dust
extinction of the redshift component. As shown in Figure 6, when the redshift
component was obscured by dust, the VVD was tilted. Therefore, when the gas
fraction is 10%, the VVD diagram is symmetrical because it is optically thin
with respect to the redshift component (i.e., beyond the disk).

In addition, dust extinction is important for $\sigma_{v_{los}}\,$with a gas
fraction of 30% (see Figure 5). For the models with a gas fraction of 30%,
Figure 5 shows that $P$($\sigma_{v_{los}}\,$$>$ 500 km s-1) is 0.37 times
lower without dust extinction than that with dust extinction. This was caused
by the density gradient in the line of sight (i.e., the vertical density
gradient). In Figure 5, the dense gas in the optically thick region of the gas
disk affects Equations 5 and 6 by the square of the density. Such high-density
gas regions have lower temperatures and velocity dispersions; therefore, the
velocity dispersion in Figure 5 is systematically lower.

Figure 7 shows that the VVD diagram with $\tau_{5007}$ = $\infty$ for $v_{z}$
$\leq$ 0 particles with rest frame for SMBH. In contrast, for $v_{z}$ $>$ 0
particles, dust extinction was considered (i.e, $\tau_{5007}\neq 0$, which is
an important difference from Figure 6). Figure 7 shows a systematically higher
velocity dispersion than that shown in Figures 5 and 6. As high-density
regions are strongly affected by dust extinction, high-temperature, high-
velocity, and optically thin regions will be more significant. Consequently,
in the optically thin region, a relatively dense region was observed to be the
most violent collision–excited. Therefore, in buried AGNs, the velocity
dispersion of the observed emission lines is higher, owing to dust extinction.
Note that in a typical AGN, the smaller the effect of dust extinction, the
more the velocity dispersion increases, owing to the contribution of the
redshift component (Bae & Woo, 2016). The increase in velocity dispersion due
to dust extinction is not important in a typical AGN, as shown in the VVD
diagram for a gas fraction of 10% in Figure 5.

Figure 8: Probability density function of equation 8 for the merger phase of
models with gas fraction of 30% which are G3T0, G3T30, and G3T60 (blue) and
10% which are G1T0, G1T30, and G1T60 (orange) from face-on view ($\theta\leq$
45 deg). The vertical solid and dashed lines correspond to the medians of the
histgrams with gas fraction of 30% and 10%, respectively. Figure 9: VVD
diagrams classified with respect to AGN-luminosity for a gas fraction of 30%,
namely, G3T0, G3T30, and G3T60 (upper row), and gas fraction of 10%, namely,
G1T0, G1T30, and G1T60 (bottom row). The AGN luminosity of individual SMBHs
was used, not the combined luminosity of each pair as reported in Figure 3.
For a gas fraction of 30%, the threshold was set to 1012 L⊙, while for a gas
fraction of 10%, the threshold was set to 1011 L⊙. The mean value for
$\bar{v}_{los}\,$of each $\sigma_{v_{los}}\,$bin is indicated by square (all
AGN-luminosity, i.e., Figure 4) and triangle symbols (each AGN-luminosity
bin). The dotted line in the histogram for each axis indicates the mean value
while the gray background color indicates the standard deviation.

##  4 Discussion

###  4.1 AGN-luminosity and VVD diagrams

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.

###  4.2 Comparison with Toba et al. (2017)

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.

Figure 10:  Comparison of the theoretical model of this study with the VVD
diagram of Toba et al. (2017) (i.e., red markers). The colored VVD diagrams
for gas fraction of 30% models and gray VVD diagrams for gas fraction of 10%
models. The mean value for $\bar{v}_{los}\,$of each $\sigma_{v_{los}}\,$bin is
indicated by square for gas fraction 30% and triangle symbols for gas fraction
10%.

###  4.3 model limitation

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.

##  5 summary and conclusions

We simulated of late-stage galaxy mergers using the N-body/SPH code ASURA
(Saitoh et al., 2008, 2009; Saitoh & Makino, 2013). Our goal is to explore the
velocity characteristics of ionized outflows during the late stages of galaxy
mergers. Our study yields the following key findings:

  1. 1.

Late-stage galaxy mergers exhibit strong ionized outflows with velocity
dispersions surpassing 500 km s-1. (Figure 4).

  2. 2.

The mean observed velocity ($\bar{v}_{los}\,$) of the ionized outflows is
significantly influenced by dust extinction. In our models, higher gas
fractions result in a greater tilt of the VVD diagram toward a blueshift
(Figures 4 and 5).

  3. 3.

The velocity dispersion ($\sigma_{v_{los}}\,$) of the ionized outflows are
also affected by dust extinction. The [OIII] $\lambda$5007 emission lines from
dense gas are attenuated by dust grains, leading to higher velocity
dispersion, enabling the observation of more diffused and hotter gas
particles. This trend holds particular importance for buried AGNs. (Figure 5)

Although the luminosity and equivalent width of emission lines were not
discussed in this study, we intend to perform pseudo-observations of ionized
emission lines part of a future research.

We are grateful to Takayuki Saito for providing us with the ASURA code (Saitoh
et al., 2008, 2009; Saitoh & Makino, 2013). Numerical computations were
performed on a Cray XC50 at the Center for Computational Astrophysics,
National Astronomical Observatory of Japan. We also thank the anonymous
referee for constructive comments on the text.

##  Appendix A Convergence of resolution

We performed galaxy merger simulations with a higher resolution model (here
after HG3T0) than with the G3T0. Table 3 reports the model parameters for
HG3T0. The spatial resolution was 2 pc, the stellar particle mass resolution
was 1.25$\times 10^{4}M_{\odot}$, and the gas particle mass resolution was
2$\times$10${}^{3}M_{\odot}$ for the HG3T0 model. Figure 11 shows the distance
between the binary BHs, total AGN bolometric luminosity, and star formation
rate for the HG3T0 and G3T0 models. There is no significant difference between
the two models. For star formation rates, the results obtained with the HG3T0
model are systematically higher due to its higher spatial resolution; the
difference in AGN luminosity before 15 Myr is due to the randomness of the gas
clumps formed by the “fluctuations” in the initial gas disk.

Table 3: Initial parameters of the pre-merger system

Namea | ${M_{BH}}^{\rm b}$ | ${M_{bul}}^{\rm c}$ | ${M_{\rm gas}}^{\rm d}$ | ${\epsilon}^{\rm e}$ | ${r_{acc}}^{\rm f}$ | ${r_{h}}^{\rm g}$ | ${R_{\rm gas}}^{\rm h}$ | ${\Delta M_{bul}}^{\rm i}$ | ${\Delta M_{\rm gas}}^{\rm j}$ | ${\mu_{gas}}^{\rm k}$
---|---|---|---|---|---|---|---|---|---|---
| [$10^{8}M_{\odot}$] | [$10^{11}M_{\odot}$] | [$10^{10}M_{\odot}$] | [pc] | [pc] | [kpc] | [kpc] | [$10^{4}M_{\odot}$] | [$10^{3}M_{\odot}$] |
HG3 | 1.0 | 1.0 | 1.0 | 2.0 | 4.0 | 2.2 | 1.0 | 1.25 | 2.0 | 30%

Note. — (a) Name of the pre-merger system. (b) Mass of sink, (c) star, and (d)
gas particles, respectively. (e) Gravitational softening (f) Accretion radius
of sink particle. (g) Effective radius of stellar bulge Sersic profile. (h)
Outer edge radius of uniform-density gas disk. (i) Mass resolution of star
particles, (j) and gas particles. (k) Gas fraction within 1kpc.

Figure 11: Comparison between G3T0 and HG3T0 with respect to the distance BH-
binary (left panel), total AGN luminosity (center panel), and star formation
rate (right panel).

## References

  * Bae & Woo (2016) Bae, H.-J., & Woo, J.-H. 2016, ApJ, 828, 97, doi: 10.3847/0004-637X/828/2/97
  * Barro et al. (2017) Barro, G., Faber, S. M., Koo, D. C., et al. 2017, ApJ, 840, 47, doi: 10.3847/1538-4357/aa6b05
  * Blecha et al. (2018) Blecha, L., Snyder, G. F., Satyapal, S., & Ellison, S. L. 2018, MNRAS, 478, 3056, doi: 10.1093/mnras/sty1274
  * Dey & Ndwfs/MIPS Collaboration (2009) Dey, A., & Ndwfs/MIPS Collaboration. 2009, in Astronomical Society of the Pacific Conference Series, Vol. 408, The Starburst-AGN Connection, ed. W. Wang, Z. Yang, Z. Luo, & Z. Chen, 411.  https://arxiv.org/abs/0905.4531
  * Dougherty et al. (2023) Dougherty, S. L., Harrison, C. M., Kocevski, D. D., & Rosario, D. J. 2023, MNRAS, doi: 10.1093/mnras/stad1300
  * Erb et al. (2006) Erb, D. K., Steidel, C. C., Shapley, A. E., et al. 2006, ApJ, 646, 107, doi: 10.1086/504891
  * Fan et al. (2016) Fan, L., Han, Y., Fang, G., et al. 2016, ApJ, 822, L32, doi: 10.3847/2041-8205/822/2/L32
  * Farrah et al. (2002) Farrah, D., Serjeant, S., Efstathiou, A., Rowan-Robinson, M., & Verma, A. 2002, MNRAS, 335, 1163, doi: 10.1046/j.1365-8711.2002.05698.x
  * Ferland et al. (2017) Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mexicana Astron. Astrofis., 53, 385, doi: 10.48550/arXiv.1705.10877
  * Gao et al. (2020) Gao, F., Wang, L., Pearson, W. J., et al. 2020, A&A, 637, A94, doi: 10.1051/0004-6361/201937178
  * Hopkins et al. (2008) Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356, doi: 10.1086/524362
  * Jun et al. (2020) Jun, H. D., Assef, R. J., Bauer, F. E., et al. 2020, ApJ, 888, 110, doi: 10.3847/1538-4357/ab5e7b
  * Kawaguchi et al. (2020) Kawaguchi, T., Yutani, N., & Wada, K. 2020, ApJ, 890, 125, doi: 10.3847/1538-4357/ab655a
  * Laor & Brandt (2002) Laor, A., & Brandt, W. N. 2002, ApJ, 569, 641, doi: 10.1086/339476
  * Laor & Draine (1993) Laor, A., & Draine, B. T. 1993, ApJ, 402, 441, doi: 10.1086/172149
  * Narayanan et al. (2010) Narayanan, D., Dey, A., Hayward, C. C., et al. 2010, MNRAS, 407, 1701, doi: 10.1111/j.1365-2966.2010.16997.x
  * Okamoto et al. (2008) Okamoto, T., Nemmen, R. S., & Bower, R. G. 2008, MNRAS, 385, 161, doi: 10.1111/j.1365-2966.2008.12883.x
  * Osterbrock & Ferland (2006) Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Paulino-Afonso et al. (2017) Paulino-Afonso, A., Sobral, D., Buitrago, F., & Afonso, J. 2017, MNRAS, 465, 2717, doi: 10.1093/mnras/stw2933
  * Rakshit & Woo (2018) Rakshit, S., & Woo, J.-H. 2018, ApJ, 865, 5, doi: 10.3847/1538-4357/aad9f8
  * Rosas-Guevara et al. (2016) Rosas-Guevara, Y., Bower, R. G., Schaye, J., et al. 2016, MNRAS, 462, 190, doi: 10.1093/mnras/stw1679
  * Saitoh et al. (2008) Saitoh, T. R., Daisaka, H., Kokubo, E., et al. 2008, PASJ, 60, 667, doi: 10.1093/pasj/60.4.667
  * Saitoh et al. (2009) —. 2009, PASJ, 61, 481, doi: 10.1093/pasj/61.3.481
  * Saitoh & Makino (2013) Saitoh, T. R., & Makino, J. 2013, ApJ, 768, 44, doi: 10.1088/0004-637X/768/1/44
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161, doi: 10.1086/145971
  * Sanders & Mirabel (1996) Sanders, D. B., & Mirabel, I. F. 1996, ARA&A, 34, 749, doi: 10.1146/annurev.astro.34.1.749
  * Stierwalt et al. (2013) Stierwalt, S., Armus, L., Surace, J. A., et al. 2013, ApJS, 206, 1, doi: 10.1088/0067-0049/206/1/1
  * Toba et al. (2017) Toba, Y., Bae, H.-J., Nagao, T., et al. 2017, ApJ, 850, 140, doi: 10.3847/1538-4357/aa918a
  * Vasiliev (2019) Vasiliev, E. 2019, MNRAS, 482, 1525, doi: 10.1093/mnras/sty2672
  * Veilleux et al. (2022) Veilleux, S., Rupke, D. S. N., Liu, W., et al. 2022, ApJ, 926, 60, doi: 10.3847/1538-4357/ac3cbb
  * Wada et al. (2009) Wada, K., Papadopoulos, P. P., & Spaans, M. 2009, ApJ, 702, 63, doi: 10.1088/0004-637X/702/1/63
  * Wada et al. (2018) Wada, K., Yonekura, K., & Nagao, T. 2018, ApJ, 867, 49, doi: 10.3847/1538-4357/aae204
  * Woo et al. (2017) Woo, J.-H., Son, D., & Bae, H.-J. 2017, ApJ, 839, 120, doi: 10.3847/1538-4357/aa6894
  * Wu et al. (2018) Wu, J., Jun, H. D., Assef, R. J., et al. 2018, ApJ, 852, 96, doi: 10.3847/1538-4357/aa9ff3
  * Yamada et al. (2021) Yamada, S., Ueda, Y., Tanimoto, A., et al. 2021, ApJS, 257, 61, doi: 10.3847/1538-4365/ac17f5
  * Yutani et al. (2022) Yutani, N., Toba, Y., Baba, S., & Wada, K. 2022, ApJ, 936, 118, doi: 10.3847/1538-4357/ac87a2