Apparent Effect of Dust Extinction on the Observed Outflow Velocity of Ionized Gas in Galaxy Mergers

In this study, we examine photoionization outflows during the late stages of galaxy mergers, with a specific focus on the relation between the observed velocity of outflowing gas and the apparent effects of dust extinction. We used the N-body/smoothed particle hydrodynamics code ASURA for galaxy merger simulations. These simulations concentrated on identical galaxy mergers featuring supermassive black holes of 108 M ⊙ and gas fractions of 30% and 10%. From the simulation data, we derived velocity and velocity dispersion diagrams for the active galactic nuclei (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, 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 [O iii] λ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.


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 (∼62%) yutaninm@gmail.com 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 ox ≤ -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]) AB > 7.0 and infrared bright F 22µ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 velocitydispersions of [OIII] λ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(Saitoh et al. , 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.

GALAXY MERGER SIMULATIONS
We performed galaxy-merger simulations using the Nbody/SPH code ASURA (Saitoh et al. 2008(Saitoh et al. , 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 * ρ gas /t f f , where C * is a dimensionless star-formation efficiency parameter set to 0.0099, ρ gas is the local gas density, and t f f 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 H > 1000 cm −3 .(2) The gas temperature T gas < 100 K. (3) ∇ • v SP H < 0, where v SP H denotes the velocity of an SPH particle.(4) ∆Q < 0, where ∆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 K -10 8 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 GM BH /(r 2 acc + ϵ 2 ) 1/2 , where ϵ is the gravitational softening length.We assume ϵ = 3.0 pc and r acc = 6.0 pc.
The isotropic AGN-feedback was implemented by providing thermal energy (∆E = η AGN ṀBH c 2 ) weighted by the spline function to the gas particles surrounding the sink particle.ṀBH is the gas mass accretion rate of a BH particle per step at r acc and η AGN is a free  parameter representing the energy-loading efficiency.In this study, η AGN = 0.2% was assumed based on the discussions reported in Kawaguchi et al. (2020).

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 (∼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 highvelocity outflow observed by Toba et al. (2017).The SMBH masses in their samples were distributed at approximately 10 8 M ⊙ ; therefore, we set the SMBH mass in our models to 10 8 M ⊙ .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 ∼ 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 µ gas = M gas /(M gas + M bul (r xy < 1 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 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 θ disk denotes the angle between the gas disk of the secondary system and xy-plane.We considered three cases with θ disk = 0 deg, 30 deg, and 60 deg.In the following discussion, G3T0 denotes the case of θ 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.

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 ⊙ to 10 13 L ⊙ , whereas for a gas fraction 10%, the AGN luminosity ranges from 10 10 L ⊙ to 10 12 L ⊙ .

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] λ5007 ionization parameters (-3 < logU < -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 ν = L ν /4πr 2 ) was obtained from Cloudy's AGN table (Ferland et al. 2017).The AGN SED was assumed to be  The third row shows the total bolometric luminosity of the two AGNs (LAGNs = Σ 0.1 * Ṁ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.
The AGN bolometric luminosity depends on the mass accretion rate on the sink particle.For the optical depth (τ ν ) expressed in Equation 2, we consider scattering and absorption by dust, photoionization of hydrogen, and Thomson scattering: where σ dust is the extinction cross section of MRN dust (Laor & Draine 1993), σ H is the photoionization crosssection of hydrogen (Osterbrock & Ferland 2006), and σ T is 6.65 × 10 25 cm 2 , which is the Thomson scattering cross-section.N 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 v 2 r + C 2 s , where v r is radial velocity from SMBH and C s is sound speed.If the velocity of a gas particle is exceeds 2GM BH /r, then its velocity contributes to Equations 5 and 6.We assume that the [OIII] λ5007 line intensity is proportional to n H 2 because the [OIII] λ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 (5) where v losi is the line-of-sight velocity of each SPH particle in a rest frame centered on sink particle and τ 5007 equals to σ dust (λ = 5007 Å)N H,los (θ 0 , ϕ 0 ).N H,los (θ 0 , ϕ 0 ) is the column density from each SPH particle to r = 6 kpc in the (θ 0 , ϕ 0 ) direction.The θ 0 is the angle from the rotation axis of the primary system (i.e., z-axis in Figure 1).We use N H,los (θ 0 , ϕ 0 ) as the approximate column density only for the line of sight within the solid angle Ω = 2π(1 − cos(π/12)) sr from (θ 0 ,ϕ 0 ).We have fixed ϕ 0 = 0 deg (i.e., in the x-axis direction in Figure 1) and calculated for four cases where θ 0 is 0 deg, 30 deg, 60 deg, and 90 deg (i.e., covering the range from θ = -15 deg to 105 deg and ϕ -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 ϕ direction did not change significantly on that scale (see Figure 2).We choose 360 lines of sight within Ω = 2π(1 − cos(π/12)) for each θ 0 with equal sin θdθdϕ in each SMBH.As there are two SMBHs in a model, we plotted 2 × 4×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 7is not satisfied.
where h i is the kernel radius of SPH particles.As the accuracy of SPH particles depends on the local parti-cle 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 × 4 × 360 lines of sight were selected for each snapshot in all the models.Thus, for the G3T0 model, at the maximum 181×2880 points were plotted on the VVD diagram.

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

Apparent effect of dust extinction
Figure 5 plots the VVD diagram with τ 5007 = 0 expressed in Equations 5 and 6. Figure 5 confirms that a symmetric distribution with respect to vlos = 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 τ 5007 = 0 for v z > 0 particles and τ 5007 = ∞ for v z ≤ 0 with rest frame for SMBH.In Figure 6, for simplicity, VVD diagrams are plotted from a faceon view (0 ≤ θ(deg) ≤ 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 (θ ≤ 45 deg).where dust extinction was considered (i.e., τ 5007 ̸ = 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 ≃ 0, and for G3T0 with a gas fraction of 30%, it is drastically reduced for z ≲ 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 σ 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 (σ 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 τ 5007 = ∞ for v z ≤ 0 particles with rest frame for SMBH.In contrast, for v z > 0 particles, dust extinction was considered (i.e, τ 5007 ̸ = 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.  2017) reported that some IR bright DOGs with fast [OIII] λ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 σ v los = 1000 km s −1 and vlos = -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 pseudoobservations in this study.For comparison with more detailed observations, both infrared and ionizing emission line pseudo-observations should be performed.

model limitation
We performed identical galactic-core merger simulations.In our models, we fixed the initial mass of SMBH at 10 8 M ⊙ .This is because Toba et al. (2017) suggested that the SMBH mass of IR bright DOGs is approximately 10 8 M ⊙ .DOGs that contain even heavier SMBHs than 10 8 M ⊙ 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 ⊙ are considered to be approximately 100 times rarer than SMBH at 10 8 M ⊙ (e.g.Rosas-Guevara et al. 2016).Therefore, we consider that an SMBH mass of 10 8 M ⊙ 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 ⊙ .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.

SUMMARY AND CONCLUSIONS
We simulated of late-stage galaxy mergers using the N-body/SPH code ASURA (Saitoh et al. 2008(Saitoh et al. , 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. Late-stage galaxy mergers exhibit strong ionized outflows with velocity dispersions surpassing 500 km s −1 .(Figure 4).
2. The mean observed velocity (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. The velocity dispersion (σ v los ) of the ionized outflows are also affected by dust extinction.The [OIII] λ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  The colored VVD diagrams for gas fraction of 30% models and gray VVD diagrams for gas fraction of 10% models.The mean value for vlos of each σv los bin is indicated by square for gas fraction 30% and triangle symbols for gas fraction 10%.
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(Saitoh et al. , 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.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×10 4 M ⊙ , and the gas particle mass resolution was 2×10 3 M ⊙ 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.

Figure 1 .
Figure 1.The schematic of the merger model

Figure 2 .
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.

Figure 3 .
Figure3.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 (LAGNs = Σ 0.1 * Ṁ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.

Figure 4 .
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.Nsnap is the total number of snapshots in the VVD diagram; all points in the VVD diagram are 360×Nsnap.LAGN is the mean luminosity of the snapshots used in each panel.

Figure 5 .
Figure5.VVD diagrams for galaxy merger models with τ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 Figure4.The mean value for vlos of each σv los bin is indicated by a circle (τ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 τ5007 = 0.

Figure 6 .
Figure 6.VVD diagram with τ5007 = 0 for vz > 0 particles and τ5007 = for ≤ 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 vlos of each σ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 .
Figure 7. VVD diagram with τ5007 = ∞ vz ≤ 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 vz > 0 particles.The gray plot denotes the VVD considering the dust extinction for all ionized outflow gas particles.The mean value for vlos of each σ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.

Figure 9 Figure 8 .
Figure 9 shows the VVD diagrams depicted in Figure 4, classified with respect to the AGN luminosity.The upper row in Figure9shows 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 Figure5.Thus, the AGN lu-

Figure 10 .
Figure 10.Comparison of the theoretical 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 vlos of each σv los bin is indicated by square for gas fraction 30% and triangle symbols for gas fraction 10%.

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