# FOREVER22: Gas and metal outflow from massive galaxies in protocluster
regions

Naoki Harada1,Hidenobu Yajima1 and Makito Abe1
1Center for Computational Sciences, University of Tsukuba, Ten-nodai, 1-1-1
Tsukuba, Ibaraki 305-8577, Japan
E-mail<EMAIL_ADDRESS>
(Accepted 2023 August 22. Received Received 2023 August 21; in original form
2022 October 09)

###### Abstract

We study gas and metal outflow from massive galaxies in protocluster regions
at $z=3-9$ by using the results of the FOREVER22 simulation project. Our
simulations contain massive haloes with $M_{\rm h}\gtrsim 10^{13}~{}\rm
M_{\odot}$, showing high star formation rates of $>100~{}\rm
M_{\odot}~{}yr^{-1}$ and hosting supermassive black holes with $M_{\rm
BH}\gtrsim 10^{8}~{}\rm M_{\odot}$. We show that the mass loading factor
($\eta_{\rm M}$) sensitively depends on the halo mass and it is $\eta_{\rm
M}=1.2~{}(9.2)$ for $M_{\rm h}=10^{13}~{}(10^{11})~{}{\rm M_{\odot}}$. Once
the halo mass exceeds $\sim 10^{12.5}~{}\rm{\rm M_{\odot}}$, the outflow
velocity of the gas rapidly decreases near a virial radius, and the gas
returns to a galactic centre finally as a fountain flow. Also, the metal
inflow and outflow rates sensitively depend on the halo mass and redshift. At
$z=3$, the inflow rate becomes larger than the outflow one if $M_{\rm
h}\gtrsim 10^{13.0}~{}{\rm M_{\odot}}$. Thus, we suggest that massive haloes
cannot be efficient metal enrichment sources beyond virial radii that will be
probed in future observations, e.g., studies of metal absorption lines with
the Prime Focus Spectrograph on the Subaru telescope.

###### keywords:

galaxies: evolution – galaxies: formation – galaxies: high-redshift

††pagerange: FOREVER22: Gas and metal outflow from massive galaxies in
protocluster regions–FOREVER22: Gas and metal outflow from massive galaxies in
protocluster regions††pubyear: 2023

##  1 Introduction

In the modern theory of cosmological structure formation, galaxies evolve via
mergers and matter accretion from the intergalactic medium (IGM). The growth
rate of galaxies sensitively depends on the environment in the large-scale
structures. In the overdense regions, the matter accretion rate and merger
frequency are expected to be larger than the cosmic mean fields. Thus,
multiple massive/bright galaxies can collectively form in such overdense
regions even in the early Universe. According to the hierarchical structure
formation scenario, such particular regions finally evolve into local galaxy
clusters (Overzier 2016). Therefore, the clustering regions of massive
galaxies in the early Universe are dubbed protoclusters (PCs). Thanks to the
brightness of massive galaxies in PCs, recent galaxy surveys have discovered a
lot of PCs and allowed us to study the galaxy evolution in overdense regions
(e.g., Matsuda et al. 2005; Toshikawa et al. 2018; Higuchi et al. 2019;
Umehata et al. 2019; Miller et al. 2018; Hill et al. 2020). These observations
indicate that a diversity of galaxy populations emerge in PC regions. However,
the origin of the diversity has not been understood well.

It is believed that the baryon cycle, i.e., accretion of matter, outflow, and
recycling of gas, has an inseparable relationship with the evolution of
galaxies. Star formation activities in galaxies are closely related to the
supply of the gas accreted from circumgalactic medium (CGM) and IGM, and from
galaxy mergers. The interstellar gas and metals can be returned into CGM/IGM
via stellar/active galactic nuclei (AGN) feedback. Therefore, the physical
states of CGM/IGM are tightly related to the activities of the galaxies. Some
parts of the ejected components would fall back again into the galaxy due to
its deep gravitational potential. This recycling process can affect the
subsequent star formation histories. One of the central issues in galactic
astronomy is to correctly follow the baryon cycle with the stellar/AGN
feedback. So far, various theoretical works have proposed the numerical
feedback models to reproduce the baryon cycle through the simulation of the
galaxy formation (e.g., Springel & Hernquist 2003; Schaye et al. 2015; Yajima
et al. 2017; Pillepich et al. 2018). However, there is a large room in
parameters for the star formation and the feedback models in cosmological
simulations of galaxy formation. By introducing the models and parameters
appropriately, previous studies have successfully reproduced observational
properties, e.g., the star formation history (e.g., Springel & Hernquist 2003;
Oppenheimer & Davé 2006; Pillepich et al. 2018), the galaxy luminosity and
mass function (e.g., Benson et al. 2003; Oppenheimer et al. 2010; Schaye et
al. 2015; Pillepich et al. 2018), the galaxy size (e.g., Pillepich et al.
2018), the mass-metallicity relation (MZR) (e.g., Finlator & Davé 2008; De
Rossi et al. 2017), and the presence of metals in the diffuse IGM (e.g.,
Aguirre et al. 2001; Oppenheimer & Davé 2006).

Recent observations have been studying the baryon recycling process between
galaxies, circumgalactic medium (CGM) and IGM, that can provide hints to
understand the feedback mechanism related to the enhancement and quenching of
star formation (Tumlinson et al. 2017). Bertone et al. (2007) introduced the
model of recycling in their semi-analytical studies. Also, a lot of
theoretical studies have been devoted to clear the relation between the baryon
recycle and the galaxy evolution, for example, in the viewpoint of gas
dynamics in CGM (Oppenheimer & Davé 2008; Oppenheimer et al. 2010; Ford et al.
2014; Christensen et al. 2016; Hafen et al. 2019; Borrow et al. 2020; Lochhaas
et al. 2020; Mitchell et al. 2020; Hafen 2020; Wright et al. 2021; Mitchell &
Schaye 2022), the impacts of stellar and AGN feedback (Shen et al. 2013;
Tollet et al. 2019; Wright et al. 2020; Donahue & Voit 2022), stellar mass and
star formation rates (SFRs) (Davé et al. 2011), metal enrichment history
(Brook et al. 2014; Muratov et al. 2017), and galaxy mass assembly (Anglés-
Alcázar et al. 2017; Mitchell & Schaye 2022). Among them, the particle
tracking method in smoothed-particle hydrodynamics (SPH) simulations is a
promising way to study and visualize the trace of gas flow in CGM and/or IGM
(Oppenheimer & Davé 2008; Oppenheimer et al. 2010; Ford et al. 2014;
Christensen et al. 2016; Hafen et al. 2019; Borrow et al. 2020). However, the
mass range of galaxies investigated is still limited. Recent discoveries of
high-redshift PC regions and passive galaxies bring the motivation of
understanding the feedback and the outflow of very massive galaxies.

Progenitors of today’s high-mass ”Coma-type” clusters with $M_{\rm
h}>10^{15}~{}{\rm M_{\odot}}$ is thought to have the halo mass of $M_{\rm
h}>10^{13}~{}{\rm M_{\odot}}$ at $z\sim 3$ (Chiang et al. 2013). Simulations
of massive PCs are challenging because of their rareness in space which
requires large simulation volumes. In this paper, we perform cosmological SPH
simulations with a large simulation box of $(714.3~{}\rm cMpc)^{3}$, and study
the baryon cycling in massive dark matter (DM) haloes at redshifts of $3-9$
that are supposed to be progenitors of massive galaxy clusters in the local
Universe.

In Section 2, we explain FOREVER22 project and the star formation, and the
feedback models. The inflow and outflow of gas in massive galaxies are
presented in Section 3. We show the metal outflow rate and spatial
distribution around massive haloes in Section 4. In Section 5, we summarize
our main results.

##  2 Simulation

###  2.1 The FOREVER22 simulation

In this study, we use the result of FOREVER22 project which is a new series of
cosmological hydrodynamic simulations focusing on 10 PC regions (Yajima et al.
2022b). We utilize the modified version of smoothed particle hydrodynamics
(SPH)/$N$-body code GADGET-3 (Springel et al. 2001; Springel 2005). The code
incorporates the star formation, the radiation feedback and supernova (SN)
feedback from the stars, and the AGN feedback from the massive black holes
(BHs). Here we use a standard SPH scheme that can suppress the mixing of
clouds and background gas (see also Agertz et al. 2007; Schaller et al. 2015;
Braspenning et al. 2023). Therefore, the outflow gas can propagate more
easily. On the other hand, we carefully tune the subgrid models to be fit for
use in the cosmological simulation with the standard SPH. Here we summarize
the implementation of the code. More details of the simulations are described
in Yajima et al. (2022b, a).

###  2.2 Star formation

In the FOREVER22 project, the star formation is followed by using the model
developed in Schaye & Dalla Vecchia (2008) that is based on the Kennicutt-
Schmidt law $\dot{\Sigma}_{\ast}=A(\Sigma_{\rm g}/1~{}{\rm M_{\odot}}~{}{\rm
pc}^{-2})^{n}$ of local galaxies. The star formation rate $\dot{m}_{\ast}$ is
evaluated from the total pressure $P$ of the interstellar medium (ISM) as
follows

$\dot{m}_{\ast}=m_{\rm SPH}A(1~{}{\rm M_{\odot}}~{}{\rm
pc}^{-2})^{-n}\left(\frac{\gamma}{G}f_{\rm g}P\right)^{(n-1)/2},$ (1)

where $m_{\rm SPH}$ is the SPH particle mass, $\gamma=5/3$ is the ratio of
specific heats, $f_{\rm g}$ is the gas mass fraction (set to unity), $n$ and
$A$ are the free parameter. During the simulation time step $\Delta t$, the
conversion probability from an SPH particle to a stellar particle is evaluated
by (Schaye & Dalla Vecchia 2008)

$\dot{p}_{\ast}=\min\left(\frac{\dot{m}_{\ast}\Delta t}{m_{\rm
SPH}},1\right).$ (2)

In the simulation, the stellar particle is treated as a star cluster assuming
a Chabrier initial mass function (IMF) with the mass range of $0.1-100~{}{\rm
M_{\odot}}$. We confirmed that the star formation and feedback models of the
FOREVER22 project reproduced observational properties, e.g., the cosmic star
formation rate density and the stellar mass function (Yajima et al. 2022b),
and most parameters referred to the values derived in Schaye et al. (2015).
The slope $n=1.4$ and $A=1.5\times 10^{-4}~{}{\rm{\rm M_{\odot}}~{}{\rm
yr^{-1}}}~{}{\rm kpc}^{-2}$ are the fiducial values, but the slope $n$ is set
to 2.0 if the gas density exceeds the metallicity dependent critical value
$n_{0}(Z/0.002)^{-0.64}~{}{\rm cm}^{-3}$. Also, if the gas density satisfies
$n_{\rm H}>n_{0}$, the effective equation of state with the adiabatic index
$\gamma_{\rm eff}$ is used. To prevent the undesirable gravitational
fragmentation in the unresolved scale, we set $\gamma_{\rm eff}=4/3$. The
critical density $n_{0}$ and the floor temperature at $n_{\rm H}=n_{0}$ are
the parameters, which depend on the resolution of the simulation (see §2.1 in
Yajima et al. 2022b). Here we set $n_{\rm 0}=0.1~{}\rm cm^{-3}$.

The gas temperature is estimated with a net cooling rate derived based on the
collisionally and photoionization equilibrium states of each gas and metal
species using cloudy v07.02 (Ferland 2000). We take into account the line
cooling processes of metals and primordial gas and the cooling by continuum
radiation such as from thermal bremsstrahlung and recombination.

###  2.3 Stellar feedback

Young stars have an impact on the physical states of the interstellar medium
(ISM), the CGM, and the IGM via radiation and SN feedback. These feedback
processes are considered in our simulations as follows.

####  2.3.1 Radiation feedback

First, the simulations take into account the photoionization feedback. We
consider UV radiation only from star clusters with an age younger than
$10^{7}~{}\rm yr$. The photoionized region is evaluated by taking the balance
between the ionizing photon production rate $\dot{N}_{\rm ion}$ of a stellar
particle and the recombination rate of neighbour gas particles as

$\dot{N}_{\rm ion}=\sum_{i}n_{{\rm H},i}^{2}\alpha_{\rm
B}\frac{m_{i}}{\rho_{i}},$ (3)

where $n_{{\rm H},i}$ $m_{i}$ and $\rho_{i}$ are hydrogen number density, mass
and density of $i$-th SPH particle, $\alpha_{\rm B}$ is the case-B
recombination coefficient (Hui & Gnedin 1997). The local recombination rate is
summed up in order of the distance from the stellar particle until the
recombination rate reaches $\dot{N}_{\rm ion}$. The $\dot{N}_{\rm ion}$ of
each stellar particle is evaluated by using a synthesized SED of the Chabrier
IMF considering its age and metallicity. Once the gas is ionized, the
temperature of the photoionized SPH particles is set to be $3\times
10^{4}~{}{\rm K}$ as hot $\rm{H_{II}}$ bubbles, and the star formation is
forbidden in the $\rm{H_{II}}$ region. This feedback can affect the state of
ISM if we have sufficient resolution to resolve the cold gas component.
Therefore, this feedback can play an important role in the suppression of star
formation significantly in the high-resolution series of the FOREVER22
simulations, First0 and First1 runs.

In addition to the photoionization process, the UV photon absorption by dust
gives momentum to gas. The simulation models the feedback of the UV radiation
pressure by giving the momentum kick to the SPH particle. We evaluate the mean
free path of the UV continuum photons against the dust as $l_{\rm
mfp}=1/\kappa_{\rm d}\rho$ where $\kappa_{\rm d}$ is the opacity of dust. If
gas particles are within $l_{\rm mfp}$ from a star cluster, they are
influenced by the radiation pressure as

$\bm{F}_{\rm rad}=\frac{\kappa_{\rm d}\rho L_{\rm UV}}{4\pi
r^{2}c}\frac{\bm{r}}{r},$ (4)

where $r$ is the distance from the stellar particle to the gas particle,
$L_{\rm UV}$ is the luminosity for the UV continuum ($1000-5000~{}$Å).

Besides, we take into account the UV background radiation at $z\lesssim 10$ by
using the model of Haardt & Madau (2001). Once the UV background turns on, the
chemical compositions and the cooling rates can change due to the photo-
ionization process. We consider the modifications of the cooling rates for 11
elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe) under the UV background
following Wiersma et al. (2009a). Also, the self-shielding effect for dens gas
is incorporated. We assume that if the hydrogen number density is larger than
$0.01~{}\rm cm^{-3}$, the flux decreases with a fraction $(n_{\rm
H}/0.01~{}{\rm cm^{-3}})^{-2}$ (Yajima et al. 2012; Rahmati et al. 2013).
Thus, the UV background suppresses the cooling of the low-density gas.

####  2.3.2 Supernova Feedback

The SN feedback is modelled by using the recipe proposed by Dalla Vecchia &
Schaye (2012). In this scheme, the SN energy is injected into neighbour gas
particles as thermal energy. To avoid the over-cooling problem, the heating
gas particles are stochastically chosen and are heated to $10^{7.5}~{}{\rm
K}$. Dalla Vecchia & Schaye (2012) argued that the injected thermal energy is
converted into kinetic one if the cooling time is longer than the local sound
crossing time. Thus, the SN feedback effectively works if the gas density is
less than the following critical value,

$n_{\rm H}\sim 100~{}{\rm cm^{-3}}~{}\left(\frac{T}{10^{7.5}~{}{\rm
K}}\right)^{3/2}\left(\frac{m_{\rm SPH}}{10^{4}~{}{\rm
M_{\odot}}}\right)^{-1/2}.$ (5)

As in Schaye et al. (2015), we introduce a multiplication factor $f_{\rm th}$
to avoid the over-cooling problem even at the high-density star-forming region
($n_{\rm H}>100~{}{\rm cm^{-3}}$),

$f_{\rm th}=f_{\rm th,min}+\frac{f_{\rm th,max}-f_{\rm
th,min}}{1+\left(\frac{Z}{0.1~{}Z_{\odot}}\right)^{n_{\rm
z}}\left(\frac{n_{\rm H,birth}}{n_{\rm H,0}}\right)^{-n_{\rm n}}},$ (6)

where $n_{\rm H,birth}$ is the gas density where a stellar particle forms,
$n_{\rm z}$, $n_{\rm n}$ and $n_{\rm H,0}$ are free parameters. The SN
feedback is likely to be strong in low metallicity and low-density regions
because of inefficient radiative cooling. Therefore, the multiplication factor
is set as the function of metallicity $Z$ and the density $n_{\rm H}$ as
$n_{\rm z}=n_{\rm n}=2/\ln(10)$ and $n_{\rm H,0}=0.67~{}{\rm cm^{-3}}$ (Schaye
et al. 2015). We set $f_{\rm th,max}=2.5$ and $f_{\rm th,min}=0.3$ as the
upper/lower limit of $f_{\rm th}$. The multiplication factor can exceed unity.
This considers the additional feedback processes that are not included in the
simulation, e.g., stellar winds, cosmic rays, and/or more energetic SN yields
than assumed in the simulation. Note that the upper limit of the
multiplication factor is lower than that of previous works ($f_{\rm
th,max}=3.0$ is chosen in Schaye et al. 2015) by considering the additional
radiative feedback channels as described above. With the parameters, our
previous work reproduced properties of observed high-redshift galaxies (Yajima
et al. 2022b). In our simulations, stellar particles continuously release
hydrogen, helium, and metals into the neighbour gas. Based on the prescription
in Wiersma et al. (2009b), we consider the gas and metal productions with
tabulated yields for Types Ia and II SNe, and from asymptotic giant branch
stars.

###  2.4 Black holes

The AGN feedback is expected to suppress the star formation in massive
galaxies via the heating and evacuating ISM (e.g., Dubois et al. 2012; Schaye
et al. 2015). Our simulations include the growth and the feedback processes of
massive black holes (BHs) in galaxies. Since resolving an accretion disc
around a BH is still difficult even in the current galaxy formation
simulations, we introduce subgrid models of the growth and feedback of BHs.
First, a seed BH with the mass of $10^{5}~{}h^{-1}~{}{\rm M_{\odot}}$ is put
at the galactic centre when the halo mass reaches $10^{10}~{}h^{-1}~{}{\rm
M_{\odot}}$. As for the orbital evolution of BH particles, the simple
repositioning technique is adopted to mimic the effect of dynamical friction
on the BH. Recent work suggested that this technique is important to consider
the feedback from supermassive BHs (SMBH) in the large-scale cosmological
hydrodynamic simulations (Bahé et al. 2022). We evaluate the gas accretion
rate onto the BH particle based on the Bondi-Hoyle-Lyttleton model as

$\dot{m}_{\rm Bondi}=\frac{4\pi cGM_{\rm BH}^{2}\rho}{\left(c_{\rm
s}^{2}+v_{\rm rel}^{2}\right)^{3/2}},$ (7)

where $\rho$, $c_{\rm s}$ are the density and the sound speed of gas particle
around a BH, $v_{\rm rel}$ is the relative velocity between the gas and the
BH. We consider the suppression of the gas accretion due to the angular
momentum of the gas as

$\dot{m}_{\rm acc}=\dot{m}_{\rm Bondi}\times{\rm min}(C^{-1}_{\rm visc}(c_{\rm
s}/V_{\rm\phi})^{3},~{}1),$ (8)

where $C_{\rm visc}$ is a numerical parameter regarding the viscosity of the
accretion disc, $C_{\rm visc}=200\pi$ is chosen as the fiducial setup (Schaye
et al. 2015). Using the accretion rate, the growth rate of BHs is estimated by
$\dot{m}_{\rm BH}=(1-f_{\rm r})\dot{m}_{\rm acc}$, where $f_{\rm r}=0.1$ is
the radiative efficiency. The physical properties of neighbour gas are
evaluated by using the nearest 100 gas particles. The fiducial model takes
into account the upper limit of the accretion rate by the Eddington limit,

$\dot{m}_{\rm Edd}=\frac{4\pi GM_{\rm BH}m_{\rm p}}{f_{\rm r}\sigma_{\rm
T}c}.$ (9)

The FOREVER22 project also includes the run with the optional allowing the
super-Eddington accretion, i.e., $\dot{m}_{\rm acc}>\dot{m}_{\rm Edd}$.
Throughout this paper, we investigate only the results from the fiducial
models with the Eddington limit.

In the simulations, during the timestep of $\Delta t$ the released energy is
estimated as $\Delta E=f_{\rm e}f_{\rm r}\dot{m}_{\rm BH}c^{2}\Delta t$, where
$f_{\rm e}=0.15$ is the thermal coupling factor. As described below, the
FOREVER22 project considers two types of feedback mechanisms, thermal (quasar
mode) or kinetic (radio mode) feedback.

####  2.4.1 Quasar mode feedback

In the quasar mode feedback, the released energy $\Delta E$ is injected into
neighbour gas particles by heating the temperature up to $\Delta
T=10^{9}~{}{\rm K}$. Unlike the SN feedback, the gas-particle heated is
selected from the nearest one. If the released energy from a BH is large, the
next particles are heated up in order of the distances. If the BH mass is
smaller than $10^{9}~{}{\rm M_{\odot}}$, the quasar mode feedback alone is
always used.

####  2.4.2 Radio mode feedback

The jet-like kinetic feedback is activated when the BH mass exceeds
$10^{9}~{}{\rm M_{\odot}}$. In radio mode feedback, half of $\Delta E$ is
injected as the kinetic energy, and the other one is used as the thermal
energy. The target gas particle of the radio mode feedback receives the
momentum kick with the velocity of $3000~{}\rm km\;s^{-1}$. The direction axis
of momentum is determined along the total angular momentum vector
$\bm{n}=\bm{L}/|\bm{L}|$ of the neighbour gas particles. Then the direction
$\bm{n}$ or $-\bm{n}$ is randomly chosen.

###  2.5 The Models

The FOREVER22 project focuses on massive 10 PC regions chosen in the (714.3
cMpc)3 simulation volume. The PC regions are explored with the different three
sizes of zoom-in simulation volumes with different resolutions to cover the
wide dynamic ranges; Proto Cluster Region (PCR) run (28.6 cMpc)3, Brightest
proto-Cluster Galaxy (BCG) run ($\sim$ 10 cMpc)3 and First run ($\sim$ 3
cMpc)3. In this study, we refer to the results of PCR runs. The PCR runs allow
us to investigate the large-scale gas structure around galaxies covering a
wide range of their mass. The PCR runs resolve the zoom-in volume with the gas
and DM particle mass resolutions of $m_{\rm SPH}=2.9\times
10^{6}~{}h^{-1}~{}{\rm M_{\odot}}$ and $m_{\rm DM}=1.6\times
10^{7}~{}h^{-1}~{}{\rm M_{\odot}}$. The gravitational softening length is set
to $2.0~{}h^{-1}~{}{\rm kpc}$. 10 PC regions are labeled as PCR0 to PCR9 and
the simulations are performed from $z=100$ to $2$. Most massive haloes in each
PC region reach $\sim 10^{14}~{}{\rm M_{\odot}}$ at $z=2$. We have total 200
snapshots for each run, which is corresponding to the time bin $\Delta t\sim
16~{}{\rm Myr}$. We perform the post-processing analysis for the snapshot data
to explore the gas dynamics and the statistical properties of very massive
galaxies in detail.

Figure 1: Spatial distributions of gas column densities in most massive haloes
at $z=3.0$ in the three protocluster regions: PCR0, PCR1, and PCR2. The halo
masses are $1.2\times 10^{14}~{}\rm{\rm M_{\odot}}$ (PCR0), $8.4\times
10^{13}~{}\rm{\rm M_{\odot}}$ (PCR1), and $5.5\times 10^{13}~{}\rm{\rm
M_{\odot}}/h$(PCR2). The projection depth is 50 $\rm ckpc~{}h^{-1}$ and white
arrows represent streamlines of gas.

##  3 GAS INFLOW AND OUTFLOW IN MASSIVE HALOES

###  3.1 Halo mass dependence of outflow properties

Fig. 1 shows structures of gas flow around most massive haloes at $z=3$ in
PCR0, PCR1, and PCR2 runs. There are large-scale inflow structures along IGM
filaments. Inflow and outflow coexist near star-forming regions, and some
streamlines show fountain structures. Below, we statistically evaluate the
inflow and outflow properties.

We evaluate the inflow and outflow rates as in Nelson et al. (2019). First, we
estimate the radial velocity as $v_{\rm r,i}=\frac{\bm{v}_{\rm
i}(t)\cdot\bm{r}_{\rm i}}{|{\bm{r}_{\rm i}}|}$ , where $\bm{v}_{\rm i}$ is the
velocity vector of $i$-th gas particle, and define the outflow (inflow) if it
is positive (negative). The position vector of $i$-th particle $\bm{r}_{\rm
i}$ is estimated from the centre of mass of stellar component. Then, we
consider all outflowing gas particles of radial distances in the range between
$r$ and $r+\delta r$ and estimate the outflow rate as

$\dot{M}_{\rm out}(r)=\sum_{i=1}^{n}\frac{m_{\rm i}}{\delta r}|v_{\rm r,i}|,$
(10)

where $m_{\rm i}$ is mass of $i$-th gas particle. In the case of inflow, we
consider only gas particles with negative radial velocities. Note that, here
the inflow gas includes accretion from IGM, galaxy mergers, and returning ones
after evacuation from star-forming regions as the outflow.

Figure 2:  Halo mass $M_{\rm h}$ scaling of the mass loading factor $\eta_{\rm
M}$, median of the radial velocity $v_{\rm r-50}$, and the outflow rate
$\dot{M}_{\rm out}$ extracted from top 500 massive haloes in PC0 at $z=3.16$.
Red circles indicate the median and gray shades show the quantiles (25 - 75
percent) in each bin with a bin width of 0.4 dex. Red dashed lines in top
panels are fitting results of medians by the power law, $\eta_{\rm
M}=7.43\times 10^{5}(M_{\rm h}/{\rm M_{\odot}})^{-0.45}$ for the left top
panel and $\eta_{\rm M}=7.24\times 10^{3}(M_{\rm h}/{\rm M_{\odot}})^{-0.27}$
for the right top panel. Blue lines in the middle panels are the escape
velocity estimated assuming the NFW profile. In the left panels, each quantity
is evaluated at radius $r=0.1R_{\rm vir}$ as galaxy-scale outflow whereas at
$r=R_{\rm vir}$ as halo-scale outflow in right panels.

In this work, we measure the inflow and outflow rates at $r=0.1$ or
$1.0\;R_{\rm vir}$. The former ones indicate “galaxy-scale” outflow, whereas
the latter ones do “halo-scale” outflow. Fig. 2 presents the outflow rate
$\dot{M}_{\rm out}$, mass loading factor $\eta_{\rm M}$, and outflow velocity
at $z=3$. To understand the relationship between the feedback and the halo
gravitational potential, we show the halo mass scaling of these values rather
than the stellar mass. Here we consider the top 500 massive haloes identified
by the standard friend-of-friend (FoF) grouping technique (containing central
and satellite galaxies), corresponding to the mass range from $6.3\times
10^{10}~{}{\rm M_{\odot}}$ to $1.0\times 10^{14}~{}{\rm M_{\odot}}$. The mass
loading factor is evaluated as $\eta_{\rm M}=\dot{M}_{\rm out}/{\rm SFR}$. We
find that, at $r=0.1R_{\rm vir}$, $\eta_{\rm M}$ monotonically decreases with
increasing halo mass at the mass range of $M_{\rm h}\lesssim 10^{13}~{}{\rm
M_{\odot}}$, showing $\eta_{\rm M}\sim 8.8$ at $M_{\rm h}=10^{10.8}~{}{\rm
M_{\odot}}$ and $\eta_{\rm M}\sim 1.1$ at $M_{\rm h}=10^{12.8}~{}{\rm
M_{\odot}}$. On the other hand, the trend of $\eta_{\rm M}$ reverses at the
halo mass higher than $10^{13}~{}{\rm M_{\odot}}$ and the most massive halo
shows $\eta_{\rm M}=9.3$. Such a trend is also seen in some previous studies.
Nelson et al. (2019) found an upturn at the stellar mass of $\sim
10^{10.5}~{}{\rm M_{\odot}}$ in Illustris TNG50 simulation at $z=2$, while
Mitchell et al. (2020) showed it occurred at $M_{\rm h}\sim 10^{12}~{}{\rm
M_{\odot}}$ in EAGLE simulation at $z<7.3$. They suggested the upturn was due
to the strong BH feedback in massive haloes. Reasonable AGN feedback models
are still under debate and depend on the resolution and simulation schemes.
The feedback efficiency $f_{\rm e}$ and the energy fraction deposited into the
kinetic or thermal form are quite different among different simulations. Our
model is basically similar to that in EAGLE (Schaye et al. 2015) but includes
the jet-type feedback as in illustrisTNG (Nelson et al. 2019). Note that,
however, the number of such massive haloes is not sufficient for statistically
robust discussion due to random variability within the cosmological
simulations (Genel et al. 2019; Keller et al. 2019).

Using the median values at $M_{\rm h}<10^{13}~{}{\rm M_{\odot}}$, we derive a
fitting with a form, ${\rm log_{10}}\eta_{\rm M}=\alpha{\rm log_{10}}M_{\rm
h}/{{\rm M_{\odot}}}+\beta$. We find the best-fit parameters as $\alpha=-0.45$
and $\beta=5.9$. This slope is intermediate between two characteristic scaling
law in $\eta_{\rm M}-M_{\rm h}$ relation, one is the momentum-conservation
model ($\eta_{\rm M}\propto M_{\rm h}^{-1/3}$) and another is the energy-
conservation one ($\eta_{\rm M}\propto M_{\rm h}^{-2/3}$) (Murray et al. 2005;
Veilleux et al. 2020). The reasonable model sensitively depends on the cooling
efficiency at a post-shock region. The adiabatic case is corresponding to the
energy-conservation model. Our results suggest that the cooling is moderately
effective on the galactic scale.

In the halo scale, $\eta_{\rm M}$ slowly decreases at $M_{\rm h}\lesssim
10^{13}~{}{\rm M_{\odot}}$ and becomes almost flat at $M_{\rm h}\gtrsim
10^{13}~{}{\rm M_{\odot}}$. We perform the fitting at $M_{\rm
h}<10^{13}~{}{\rm M_{\odot}}$ as with the galaxy-scale analysis, and the best-
fit slope is $\alpha=-0.27$. Therefore, we indicate that the outflow
propagates with the momentum-conservation mode in the halo scale. Our results
at $M_{\rm h}\lesssim 10^{13}~{}{\rm M_{\odot}}$ are similar to Mitchell et
al. (2020) in the galaxy scale, but the slope is shallower than their work
($\alpha\sim-1$) in the halo scale. This may indicate that in our simulations
the cooling works efficiently rather than Mitchell et al. (2020) at the halo
scale. Although the number of samples is not sufficient to discuss the
statistical features, our simulation contains the more massive haloes. Our
result shows even massive haloes with $M_{\rm h}\sim 10^{14}~{}{\rm
M_{\odot}}$ can keep $\eta_{\rm M}$ greater than unity that indicates
protoclusters in the early Universe can provide gas and metals into IGM.

Radial velocities of outflow gas are also shown in the figure. In the
galactic-scale, median values of the outflow gas ($v_{\rm r-50}$)
monotonically increase with the halo mass from $38~{}\rm km\;s^{-1}$ at
$M_{\rm h}=10^{10.8}~{}{\rm M_{\odot}}$ to $304~{}\rm km\;s^{-1}$ at
$10^{14}~{}{\rm M_{\odot}}$. One of the characteristic velocities of the
galactic halo is the escape velocity $V_{\rm esc}=(2GM_{\rm h}/R_{\rm
vir})^{1/2}$ which is equal to $\sqrt{2}$ times the circular velocity $V_{\rm
c}$. Assuming $M_{\rm h}\propto R_{\rm vir}^{3}$, we obtain the halo mass
scaling of escape/circular velocities as $\propto M_{\rm h}^{1/3}$. In the
middle panel of Fig. 2, we overplot $V_{\rm esc}$ as the function of $M_{\rm
h}$. Our simulations show that $v_{\rm r-50}$ has a halo mass dependence
similar to $V_{\rm esc}$ whereas the absolute value is smaller. Martin (2005)
have pointed out the tight relationship between the outflow velocity and the
circular velocity. The velocity ratio of $v_{\rm r-50}$ to $V_{\rm c}$ is
$\sim 0.4$ at $10^{10.8}~{}{\rm M_{\odot}}<M_{\rm h}<10^{12.8}~{}{\rm
M_{\odot}}$. For halo-scale outflow, however, the increasing trend disappears
beyond $M_{\rm h}\sim 10^{12}~{}{\rm M_{\odot}}$. The difference between
$V_{\rm esc}$ and $v_{\rm r-50}$ increases with the halo mass, suggesting that
most gas particles cannot escape from massive haloes. This will be discussed
in Section 3.4 and 4.

The bottom panels show outflow rates. As shown in Yajima et al. (2022b),
massive galaxies are on the observed main-sequence (Pearson et al. 2018),
having SFR of $\gtrsim 1000~{}{\rm{\rm M_{\odot}}~{}{\rm yr^{-1}}}$. A similar
amount of gas is evacuated, resulting in $\eta\gtrsim 1$. Given that the gas
accretion rate onto a halo is low, most gas can be lost within the cosmic time
scale. Note that, however, Yajima et al. (2022b) suggested that most massive
haloes in PC regions keep gas-rich states.

Massive stars and BHs are thought to be the two main drivers of the galactic
outflow. Fig. 3 shows the dependence of $\eta_{\rm M}$ on specific star
formation rate (sSFR) and gas accretion rate onto BH for $r=0.1~{}R_{\rm
vir}$. We can find that massive haloes have a lower mass loading factor as
compared with low-mass haloes even if they have similar sSFR or $\dot{M}_{\rm
BH}$. This is consistent with the negative $\eta_{\rm M}-M_{\rm h}$ relation
in Fig. 2. We find that $\eta_{\rm M}$ decreases with sSFR. Low-mass haloes
show frequent fluctuations in star formation histories due to the cycles of
outflows and refreshing of gas (see also, Yajima et al. 2017; Arata et al.
2020; Abe et al. 2021). Then, galaxies can keep high star formation rates as
the halo mass increases. Therefore, the low values of sSFR may indicate the
phase of the suppressed star formation after a starburst. Note that the
fluctuations of SFR can also be shown in the case with a merging resolution
following the structure low-mass galaxies (e.g., Crain et al. 2015).
Considering the time lag between the launching due to the starburst and the
large-scale outflow $0.1R_{\rm vir}/V_{\rm c}\sim
31.5\left(\frac{1+z}{4}\right)^{-3/2}~{}{\rm Myr}$, the anti-correlation can
be explained. The decreasing tendency for $\eta_{\rm M}$ shifts to flat (
orange line) or slightly increasing tendency (shade) at ${\rm sSFR}\gtrsim
1~{}\rm Gyr^{-1}$. A part of galaxies with ${\rm sSFR}\gtrsim 1~{}\rm
Gyr^{-1}$ are thought to be massive galaxies in the starburst phase. As the
result of a large amount of SN energy from stars, the mass loading factor may
keep the constant value slightly less than 10 (as the median value). Contrary
to the result for sSFR, we find no strong dependency of $\eta_{\rm M}$
regarding the BH activities. Here, the BH activity is measured at each time
step. However, it can change rapidly over time scales shorter than the
dynamical time of the galactic outflow. The time-scale difference may weaken
the dependence. Note that the energetic AGN feedback has the potential to
drive the outflowing gas comparable to or greater than that of stars. As we
described in §2.4, the BH particles release the energy of $\Delta E=f_{\rm
e}f_{\rm r}\dot{m}_{\rm BH}c^{2}\Delta t$ during the time interval of $\Delta
t$. Thus, we estimate the releasing energy rate of a BH as

$\displaystyle\dot{E}_{\rm BH}$ $\displaystyle=$ $\displaystyle f_{\rm
e}f_{\rm r}\dot{m}_{\rm BH}c^{2}$ (11) $\displaystyle=$ $\displaystyle
8.5\times 10^{44}~{}{\rm erg~{}s^{-1}}~{}\left(\frac{f_{\rm e}f_{\rm
r}}{0.015}\right)\left(\frac{\dot{m}_{\rm BH}}{{\rm{\rm M_{\odot}}~{}{\rm
yr^{-1}}}}\right)$

On the other hand, the available releasing SN energy per unit mass is
$\epsilon_{\rm SNII}=8.73\times 10^{15}~{}{\rm erg~{}g^{-1}}$ by assuming the
Chabrier IMF and the energy from a single event of SN as $10^{51}~{}{\rm erg}$
(Dalla Vecchia & Schaye 2012). Considering the continuous star formation over
$\sim 10$ Myr (typical lifetime of the massive stars), the releasing energy
rate due to the SN feedback may become as

$\dot{E}_{\rm SN}=\epsilon_{\rm SNII}\times{\rm SFR}=5.5\times 10^{43}~{}{\rm
erg~{}s^{-1}}~{}\left(\frac{{\rm SFR}}{100~{}{\rm{\rm M_{\odot}}~{}{\rm
yr^{-1}}}}\right).$ (12)

In the case of the energy conservation outflow, the releasing energy rate
$\dot{E}$ is equated with the production rate of kinetic energy as
$\frac{1}{2}\dot{M}_{\rm out}V_{\rm out}^{2}=\dot{E}$. On the other hand, in
the case of the momentum conservation mode, the production rate of the
momentum can be estimated as $\dot{M}_{\rm out}V_{\rm
out}=\dot{P}\sim\dot{N}m_{\rm SPH}v_{\rm th}\sim\dot{E}/v_{\rm th}$, where
$v_{\rm th}$ is the thermal velocity of heated gas particles due to the SN/BH
feedback. Given that $V_{\rm out}\sim V_{\rm esc}$, the ratio of the outflow
rates due to the BH and the SN feedback is evaluated as

$\frac{\dot{M}_{\rm out}^{\rm BH}}{\dot{M}_{\rm out}^{\rm
SN}}=\xi\left(\frac{f_{\rm e}f_{\rm r}}{0.015}\right)\left(\frac{\dot{m}_{\rm
BH}/{\rm SFR}}{0.01}\right),$ (13)

where the factor $\xi=15.3$ for the energy conservation case and
$\xi=2.7\left(\frac{10^{7.5}~{}{\rm K}}{\Delta T_{\rm
SN}}\right)^{1/2}\left(\frac{10^{9}~{}{\rm K}}{\Delta T_{\rm
BH}}\right)^{-1/2}$ for the momentum conservation case, respectively. Hence,
the BH feedback would contribute to the outflow rate if $\dot{m}_{\rm BH}/{\rm
SFR}$ is larger than $6.5\times 10^{-4}~{}(3.7\times 10^{-3})$ for the energy
(momentum) conservation case. As indicated in McAlpine et al. (2017), BHs can
grow rapidly once the halo mass exceeds $\sim 10^{12}~{}\rm M_{\odot}$. A
similar trend is shown in FOREVER22 simulation (see Fig. 8 and 11 in Yajima et
al. 2022b). In the numerical simulations, the activity of BH feedback
drastically changes over time and we find no clear dependency between
$\eta_{\rm M}$ and $\dot{M}_{\rm BH}$. The BH growth and activity are quite
different among simulations (e.g., Davé et al. 2019; Dubois et al. 2021).
However, Nelson et al. (2019) showed a similar weak dependency regardless of
the different resolution and BH models. However, the BH feedback has the
potential to boost the $\eta_{\rm M}$ even in the high-mass galaxies with the
mass $M_{\rm h}\gtrsim 10^{12}~{}{\rm M_{\odot}}$ (see Fig. 2).

Figure 3:  left panel: Mass loading factor $\eta_{\rm M}$ as a function of
sSFR and gas accretion rate onto BHs divided by stellar mass at
$r=0.1~{}R_{\rm vir}$. The top 500 massive haloes at $z=3.16$ are analyzed.
Red lines and shades show the median values and quartiles (25 - 75 percent) in
each bin. Each point is classified by color which indicates the halo mass.

###  3.2 Radial profile of inflow and outflow rate

The velocity of the outflow gas gradually decreases as interstellar matter or
CGM is stacked. Once the velocity reaches the escape velocity at a specific
radius, the outflow gas starts to return as the inflow. To study the point, we
examine radial profiles of inflow and outflow rates. Fig. 4 presents the
specific inflow and outflow rates defined as $f_{\rm
s-in,s-out}\equiv\dot{M}_{\rm in,out}(r)/M_{\rm h}$. We divide haloes into
four mass bins and measure the median values. $f_{\rm s-out}$ of the most
massive group (top-left-hand panel) declines at $r\gtrsim R_{\rm vir}$. This
indicates a part of gas returns into a halo at $r\sim R_{\rm vir}$ as a halo-
scale fountain flow. As the radial distance increases, the difference between
$f_{\rm s-in}$ and $f_{\rm s-out}$ increases and reaches around ten times. In
cases of lower mass haloes, the outflow rates keep constant or increase even
beyond $r\sim R_{\rm vir}$, suggesting that metal-enriched gas from star-
forming regions is provided into IGM. The specific outflow rates exceed $\sim
3\times 10^{-2}~{}\rm Gyr^{-1}$. Therefore, the low-mass haloes can lose most
gas within the cosmic age at $z\sim 3$ considering the gas fraction $M_{\rm
gas}/M_{\rm h}\sim\Omega_{\rm b}/\Omega_{\rm M}$ if no gas inflow exists.
However, we confirm that $f_{\rm s-in}$ is higher than $f_{\rm s-out}$ overall
radius within $2R_{\rm vir}$. Therefore, most high-redshift galaxies grow in a
gas-rich state. Mitchell et al. (2020) also showed the radial profiles of the
outflow gas within a viral radius at $z=0$. Our simulation results show
similar structures that the outflow rate increases with the radial distance
within a virial radius. Note that, the outflow rate is higher than the figure
of the radial profile in Mitchell et al. (2020) by a factor of few because of
the different redshifts. Here, by expanding the range of the radial distance
and the halo mass, we show that the rates of the outflow gas rapidly decrease
just beyond the virial radii in the case of massive haloes.

The structures of the inflow and the outflow gas are likely to be reflected in
properties of absorption lines in SEDs of background objects like QSOs or
bright galaxies, which will be one of the main topics carried out by future
observational missions, e.g., the Subaru Prime Focus Spectrograph. Wright et
al. (2020) suggested that the baryon accretion rate sensitively depended on
the sub-grid physics like stellar or AGN feedback. In particular, the
accretion rate for low-mass and low-redshift galaxies is sensitive to the sub-
grid models because of the efficient galaxy outflows. We will investigate the
reasonable sub-grid models via comparisons with upcoming observational data in
future work.

Figure 4:  Radial profiles of inflow and outflow rates divided by halo masses
at $z=3.0$. The inflow/outflow rates are evaluated at $r=R_{\rm vir}$. Red
solid and blue dashed lines are the median values of outflow and inflow rates
in each halo mass range. Shades represent the quartiles (25-75 percent).

###  3.3 Redshift evolution of inflow and outflow

Fig. 5 shows the redshift evolution of the specific gas inflow/outflow rate at
$r=R_{\rm vir}$. We pick up haloes with the stellar mass larger than
$10^{8}~{}{\rm M_{\odot}}$ and divided them into three groups according to
their halo mass $M_{\rm h}$ at each redshift. Overall, we find that there are
no clear differences between the three groups. All cases indicate an
increasing trend in $f_{\rm s-in}$ between $3<z<9$. Dekel et al. (2013); Dekel
& Mandelker (2014) developed a simple toy model on the galaxy evolution and
showed the mass accretion rate as a function of the scale factor,
$a=(1+z)^{-1}$. Given the Einstein-de Sitter universe, the specific mass
accretion rate in the galaxy by the dark matter shows a monotonical increase
with $z$ and is proportional to $(1+z)^{n}$, $n=5/2$. Dekel et al. (2013)
assumed that the specific gas inflow rate is also expressed with the same
relation and confirmed it by their simulations (Dekel & Mandelker 2014). In
the cases with the halo masses of $<10^{12}~{}{\rm M_{\odot}}$, the inflow
rates increase with redshift as $(1+z)^{2.5}$, which indicates the haloes grow
via mergers with gas-rich haloes simply. In the case of massive haloes with
$M_{\rm h}>10^{12}~{}{\rm M_{\odot}}$, the slope is somewhat shallower. This
can be due to matter accretion or halo mergers with a low baryon fraction. The
outflow rates, on the other hand, is remarkably different from the inflow and
it does not show clear evolution. It keeps fairly constant with $z$ between
$3<z<9$, whereas it shows some small variations, the increasing trend at
$z\sim 3$, and the decreasing feature at high redshifts. Since the difference
between the inflow and the outflow produce the change of baryon mass $M_{\rm
b}$ in haloes, the large difference at high redshift suggests active growth in
haloes. At a lower redshift of $z\sim 3$, however, the difference becomes
small and the baryon mass growth may be suppressed. Note that, $f_{\rm s-out}$
of the massive halo group decreases at $z\gtrsim 6$. This is likely due to the
deep gravitational potential that hosts most gas against the SN feedback.

Figure 5:  Redshift evolution of outflow and inflow rates normalized by halo
masses for three halo groups: $10^{10}<M_{\rm h}/{\rm M_{\odot}}<10^{11}$
(left panel), $10^{11}<M_{\rm h}/{\rm M_{\odot}}<10^{12}$ (center panel), and
$10^{12}<M_{\rm h}/{\rm M_{\odot}}$ (right panel). Red and blue lines are the
median values of outflow and inflow rates in each halo group, and shades
represent the quartiles (25-75 percent).

###  3.4 Trajectories of outflowing gas particles

Thanks to the Lagrangean format of our code, we can track outflow gas
particles (Oppenheimer & Davé 2008; Ford et al. 2014; Christensen et al. 2016;
Anglés-Alcázar et al. 2017; Borrow et al. 2020; Hafen 2020). Fig. 6 shows
radial distances of the outflow gas particles of two sample haloes with the
halo masses of $1.3\times 10^{14}$ (Halo0) and $1.1\times 10^{12}~{}{\rm
M_{\odot}}$ (Halo37). We follow the motions of gas particles from $z=4.8$ to
$z=2.9$, corresponding to the time period of $1.0$ Gyr. We classify target
particles into two groups to study galaxy scale and halo scale outflow. For
the galaxy-scale samples (group A), all gas particles with initial radial
positions between $0.05$ and $0.1R_{\rm vir}$ are selected. For the halo scale
(group B), we consider gas particles between $0.5$ and $0.6~{}R_{\rm vir}$.
The black and blue thick lines respectively represent median values for groups
A and B at each time, while thin lines show trajectories of randomly selected
10 particles. As shown in the top panel of Fig. 6, in Halo0, both median
values reach $r\sim 0.8~{}R_{\rm vir}$ at $t\sim 1.8~{}\rm Gyr$ and keep
staying a similar distance. This indicates that most outflow gas contributes
to hot halo gas. We find that a part of gas particles return to a galactic
centre after reaching the radius of $\gtrsim 0.5~{}R_{\rm vir}$, i.e., the
large-scale fountain flows (see the black thin line in the top panel of Fig.
6). This feature indicates that most outflow gas have outflow velocities
similar to halo escape velocities but it does not exceed them. Thus, in the
case of massive haloes, the outflow rate decreases near the virial radius as
shown in Fig. 4. In Halo37 (bottom panel in Fig. 6), on the other hand, most
gas particles exceed the virial radius within $\sim 0.3$ Gyr and reach $\sim
2.5~{}R_{\rm vir}$ at $t\sim 1$ Gyr finally. Thus, we suggest that the baryon
cycle between the galaxy and CGM sensitively depends on the host halo mass,
which will be investigated with future observations with the Prime Focus
Spectrograph on the Subaru telescope.

Figure 6:  Radial motion of outflowing gas particles in the most massive halo
(panel (a), $M_{\rm h}\sim 1.3\times 10^{14}~{}{\rm M_{\odot}}$ at $z=2.89$)
and Halo37 (panel (b), $M_{\rm h}\sim 1.1\times 10^{12}~{}{\rm M_{\odot}}$ at
$z=2.89$) from $z=4.79$ to $2.89$, corresponding to $0.96$ Gyr as the time
interval. At the initial step ($z=4.79$), some outflow particles are
categorized into two groups depending on their initial distances:
$0.05<r/R_{\rm vir}<0.1$ (black lines) and $0.5<r/R_{\rm vir}=0.6$ (blue
lines). Thick lines show median values of all outflow particles in each
category. Thin lines represent the trajectories of 10 particles in each group.

We measure the maximum distances of outflow gas particles ($R_{\rm max}$).
Fig. 7 shows the mass dependence of $R_{\rm max}$ by using 44 haloes with the
mass range from $3.66\times 10^{11}~{}{\rm M_{\odot}}$ to $2.83\times
10^{13}~{}{\rm M_{\odot}}$ at $z=3.0$. Here we pick the gas particles
distributed within $0.1~{}R_{\rm vir}$ initially and consider their
trajectories with the redshift range from $z=4.0$ to $3.0$. We confirm that
these 44 haloes are always major progenitors of the haloes at $z=3$, i.e.,
they do not experience mergers with heavier haloes, because it is difficult to
distinguish fountain flows and accretion onto merger companions. We find that
most gas in haloes with a mass larger than $10^{12}~{}{\rm M_{\odot}}$ do not
go beyond virial radii and $R_{\rm max}$ distributes at $\sim R_{\rm vir}$ or
$\lesssim 0.5~{}R_{\rm vir}$. On the other hand, some haloes of $M_{\rm
h}<10^{12}~{}{\rm M_{\odot}}$ show the large $R_{\rm max}$. Note that, in the
case of haloes of $M_{\rm h}<10^{12}~{}{\rm M_{\odot}}$, there is a large
dispersion with the range between $R_{\rm max}\sim 0.02$ and $2.3$ that can be
induced gas structures around star-forming regions. Borrow et al. (2020)
showed that a large fraction of gas can move beyond $1~{}\rm Mpc$ even in the
case of massive galaxies by $z=0$. They suggested that the jet feedback could
induce the outflow. Therefore, the migration distance of the gas can depend on
the redshift and the feedback models. If young stars are embedded in dense
clouds even at the end of the lifetime of massive stars $\sim 10~{}\rm Myr$,
thermal energy released by SN feedback can be lost via radiative cooling
before the galactic wind is launched. The haloes with the low values of
$R_{\rm max}$ are likely to keep star formation smoothly, while the star
formation of the ones with the large $R_{\rm max}$ can be quenched until the
gas is recovered. The quenched time scale can be $\tau\sim R_{\rm vir}/V_{\rm
c}\sim 315~{}\rm Myr$ at $z\sim 3$.

Figure 7:  Maximum distances of outflow gas particles at $3.0<z<4.0$ for 44
massive haloes. The halo mass is measured at $z=4.0$. Only outflow particles
with the initial positions of $0<r/R_{\rm vir}<0.1$ are chosen at $z=4.0$.
Black-filled circles and error bars represent the median values of the outflow
particles and the quartiles (25-75 percent).

##  4 Metal outflow

In this section, we investigate metal outflow from galaxies. The metal-
enriched CGM can be a powerful tool to reveal the feedback and the star
formation processes via comparisons between observations and theoretical
models (e.g., Tumlinson et al. 2017; Fujimoto et al. 2019). Also, the metal
distribution in IGM can reflect the star formation history over the cosmic age
(e.g., Cowie & Songaila 1998; Aguirre et al. 2001). The complex nature of the
metal recycling and the distribution would be related to not only the star
formation but also the AGN activities (Moll et al. 2007). Previous numerical
simulations investigated the outflow of metals from galaxies while considering
the variation of halo mass. For instance, Christensen et al. (2018) performed
high-resolution cosmological hydrodynamic simulations to investigate the
history of metal recycling from dwarf to Milky Way-sized galaxies. They found
that the gas outflow effectively removed metals from galactic discs and the
metal was distributed dispersedly even beyond virial radii, while more massive
galaxies can retain more metals within the virial radii. Recently, Pandya et
al. (2021) performed the state-of-art cosmological hydrodynamic simulation and
argued the properties of the gas/metal outflow. They found that the metal
loading factor for low-mass galaxies could be an order of unity, which implies
that most of the metals were ejected with the outflow. Also, for the low-mass
galaxies, the ratio of metal loading factors for ISM and halo regions
exhibited similar values, suggesting the efficient supply of the metals to CGM
and IGM. However, the ratio drops for massive galaxies and the smaller
fraction of the metals leaves the haloes. In this study, we extend the mass
range of the galaxies up to $\gtrsim 10^{13}~{}{\rm M_{\odot}}$ and
investigate the properties of metals between the galaxies and CGMs.

Fig. 8 represents the metal inflow and outflow rates at $r=R_{\rm vir}$ for
500 massive haloes at $z=3,4,5$, and $6$. The inflow and outflow rates are
calculated as

$\dot{M}^{\rm metal}_{\rm in,out}=\sum_{i=1}^{n}\frac{m_{\rm i}Z_{\rm
i}}{\delta r}|v_{\rm r,i}(r_{\rm i},t)|.$ (14)

Similar to the case of the gas flow, the metal outflow rate increases with the
halo mass. At $z=3$, massive haloes induce the metal outflow larger than
$1.0~{}{\rm{\rm M_{\odot}}~{}{\rm yr^{-1}}}$, corresponding to the gas outflow
of $\dot{M}_{\rm out}\gtrsim 1000~{}{\rm{\rm M_{\odot}}~{}{\rm yr^{-1}}}$ with
the metallicity $Z\sim 0.1~{}Z_{\odot}$. The mass loading factors for the
massive galaxies show an order of unity at the virial radii (see Fig. 2). The
metal loading factor is evaluated as $\sim 10^{-3}$. Therefore, a small
fraction of the metals are ejected into IGM in the case of massive galaxies.

As expected from the results in the previous section, the metal outflow rate
does not increase significantly at $M_{\rm h}>10^{12}~{}{\rm M_{\odot}}$. On
the other hand, the metal inflow rate simply increases with the halo mass via
merger processes and accretion of metal-enriched gas around galaxies. In
particular, the IGM around the massive haloes is metal-enriched earlier, which
is likely due to the metal outflow from neighbour low-mass haloes. Therefore,
the metal inflow rate becomes high at $M_{\rm h}\gtrsim 10^{12}~{}{\rm
M_{\odot}}$. The metal outflow rate overcomes the inflow rate in the case of
low-mass haloes and the trend is reversed at massive ones. Therefore, the
process of metal enrichment within galaxies would be distinguishable depending
on the halo mass. For low-mass haloes, the metal contents within the virial
radius can be dominated by the metals produced in-situ star formation via
supernovae. As the halo mass increases, the fraction of the metals supplied
from outside the halo via accretion and mergers can be dominant. Christensen
et al. (2018) investigated the galaxies with the mass range of $M_{\rm h}\sim
10^{9.5-12}~{}{\rm M_{\odot}}$ at $z=0$ and pointed out the similar mass
trend.

Fig. 9 shows net metal outflow rates, i.e., $\dot{M}^{\rm metal}_{\rm
out}-\dot{M}^{\rm metal}_{\rm in}$. We find that the net metal outflow changes
from positive to negative as the halo mass increases. The critical halo mass
increases as the redshift decreases and it is $\sim 10^{13.0}~{}{\rm
M_{\odot}}$ ($\sim 10^{11.2}~{}{\rm M_{\odot}}$) at $z=3$ ($z=6$). Because of
deeper gravitation potentials at higher redshifts, the gas cannot escape from
the haloes efficiently compared to the low-redshift haloes that have the same
masses. The moderately massive galaxies with $M_{\rm h}\sim 5\times
10^{12}~{}{\rm M_{\odot}}$ have the net metal outflow rates of $\gtrsim
0.1~{}{\rm{\rm M_{\odot}}~{}{\rm yr^{-1}}}$. Also, the net-metal outflow rate
of higher redshift is lower than that of lower redshift. This is likely due to
lower metallicities of IGM and ISM at higher redshifts. Galaxies with the halo
masses of $\gtrsim 10^{12}~{}{\rm M_{\odot}}$ show the metallicities with
$\sim 0.1-0.2~{}Z_{\odot}$ at $z=6$ and it increases by a factor of $\sim
1.5-2.0$ at $z=3$. Thus, low-mass haloes at lower redshifts can be efficient
metal-enrichment sources of IGM.

Figure 8:  Metal outflow and inflow rates evaluated at virial radii. Red and
blue lines are the median values of outflow and inflow rates in each halo mass
range. Shades represent the quartiles (25-75 percent).  Figure 9:  Net outflow
rate, $|\dot{M}^{\rm metal}_{\rm out}-\dot{M}^{\rm metal}_{\rm in}|$, where
$\dot{M}^{\rm metal}_{\rm out}$ and $\dot{M}^{\rm metal}_{\rm in}$ are taken
from Fig. 8. Filled and open symbols show the cases of positive $\dot{M}^{\rm
metal}_{\rm out}-\dot{M}^{\rm metal}_{\rm in}$ (outflow dominant) and negative
ones (inflow dominant) respectively.

Fig. 10 presents the radial profiles of metallicities for haloes in PCR0
region from $z=5.0$ to $3.0$. The metallicity monotonically increases with the
cosmic time, indicating the accumulation of metals in the CGM region, and also
decreases as the radial distance increases. First, metals are produced near
the centre of the galaxy and then carried to CGM and IGM via the stellar/AGN
feedback (Boissier & Prantzos 1999). We make fitting as
$\frac{Z}{Z_{\odot}}\propto\left(\frac{r}{R_{\rm vir}}\right)^{\alpha}$ for
$0.5<r/R_{\rm vir}<1$. The slope becomes shallower with increasing halo mass,
and it is $-1.06$ for $10^{10}<M_{\rm h}/{\rm M_{\odot}}<10^{11}$, $-0.85$ for
$10^{11}<M_{\rm h}/{\rm M_{\odot}}<10^{12}$, $-0.77$ for $10^{12}<M_{\rm
h}/{\rm M_{\odot}}<10^{13}$ and $-0.76$ for $10^{13}<M_{\rm h}/{\rm
M_{\odot}}$. The shallower slopes of massive haloes may reflect that metal
outflows stall near virial radii and accumulate into halo gases as suggested
in the previous section. The distributions of massive haloes with $M_{\rm
h}>10^{12}~{}{\rm M_{\odot}}$ get shallower as the redshift decreases, while
that for lower halo masses does not change significantly with time. This
reflects the different efficiencies of the metal accumulation in the CGM
depending on the nature of the outflow gas, i.e., the fountain flow or
escaping from host haloes. We note that the AGN-driven outflow would have the
potential to affect the metal distribution. For instance, Taylor & Kobayashi
(2015) argued that the powerful AGN-driven outflow can remove the gas and
metals from galaxies and enhance the metallicity of CGM compared to the result
without AGN feedback. As we mentioned in §3.1, BHs grow rapidly and give
feedback to surrounding gas once the halo mass exceeds $\sim 10^{12}~{}{\rm
M_{\odot}}$ in our simulations. Therefore, the radial profile would also be
shallower if the AGN-driven outflow activates in the low-redshift high-mass
galaxies and reduces the metals at the central region, though the impact might
strongly depend on the modeling of AGN feedback. The metallicity profiles are
also related to the mass-metallicity relation (MZR) suggested in observations
(e.g., Maiolino et al. 2008; Mannucci et al. 2009) and simulations (e.g.,
Brooks et al. 2007; Ma et al. 2016; Collacchioni et al. 2020). At $r=R_{\rm
vir}$ and $z=3.0$, for example in our results, $Z/Z_{\odot}$ increases from
$0.023$ at $10^{10}<M/{\rm M_{\odot}}<10^{11}$ (upper-left panel) to $0.061$
at $10^{12}<M_{\rm h}/{\rm M_{\odot}}<10^{13}$ (lower-left panel), but
decreases to $0.055$ in massive haloes with $10^{13}<M_{\rm h}/{\rm
M_{\odot}}$ (lower-right panel). Such a peaking behavior in the MZR is also
suggested for galaxies at $z\lesssim 1$(Collacchioni et al. 2020).

The radial metal distribution can be investigated via observations of
absorption lines in SEDs of background objects. However, observational sight-
lines are too sparse to probe the radial distributions because the background
sources are frequently limited to bright QSOs due to sensitivities of current
observational facilities (e.g., Steidel et al. 2010; Lehner et al. 2014).
Recently, Méndez-Hernández et al. (2022) showed the metal distributions from
stacking spectra using background galaxies. Future observational missions like
the Subaru PFS will be able to detect metal absorption lines in SEDs of many
background galaxies and allow us to study the radial metal distributions. With
the upcoming data, we will do comparison studies and investigate reasonable
star formation and feedback models in future work. Note that, there are
dispersions of the metal distributions in the same mass range as in the
figure. Therefore, statistical studies would be required.

Figure 10:  Radial profiles of gas metallicity. Blue dot-dashed, Black solid
and red dashed lines show median values of all haloes in each mass range at
$z=3,4$ and $5$. The shades represent the ranges of 25 - 75 percent.

##  5 Summary

In this paper, we investigate gas and metal outflows from massive galaxies in
protocluster regions by using the results of the FOREVER22 simulation project
that follows the formation of protocluster regions at $z\geq 2$. Our
simulations contain very massive haloes of $\gtrsim 10^{13}~{}{\rm M_{\odot}}$
even at high redshifts $z\gtrsim 2$. Also, there are starburst galaxies with
${\rm SFR}\gtrsim 1000~{}{\rm{\rm M_{\odot}}~{}{\rm yr^{-1}}}$ and
supermassive black holes with ${\rm M_{\rm BH}}\gtrsim 10^{8}~{}{\rm
M_{\odot}}$ in the protocluster regions at $z\sim 3$. Here we study the
relation between the outflows and the halo mass or SFR, or the gas accretion
rate onto a BH. The main conclusions are summarized as follows.

  * •

The gas outflow rate increases with the halo mass and it reaches $\sim
1000~{}{\rm{\rm M_{\odot}}~{}{\rm yr^{-1}}}$ for $M_{\rm h}\gtrsim
10^{13}~{}{\rm M_{\odot}}$. The mass loading factor decreases with the halo
mass, it is $\eta\sim 1(10)$ for $M_{\rm h}\sim 10^{13}(10^{11})~{}{\rm
M_{\odot}}$.

  * •

Massive haloes with $M_{\rm h}>10^{12.5}~{}{\rm M_{\odot}}$ show the radial
profiles of the gas outflow rates that decrease significantly at virial radii.
Therefore the outflow from star-forming regions contributes to halo gas or
large-scale fountain flow in the case of massive haloes.

  * •

The metal outflow depends on the halo mass and redshift. Considering the metal
inflow and outflow, total metal masses in massive galaxies of $M_{\rm
h}\gtrsim 1\times 10^{13}~{}{\rm M_{\odot}}$ at $z=3$ increase with time,
while lower-mass haloes can lose metal via galactic winds.

Thus, we suggest the origin of CGM with metals sensitively depends on the halo
mass. Future observational missions with James Webb Telescope or the Prime
Focus Spectrograph on the Subaru telescope will be able to investigate spatial
distributions of gas and metal around high-redshift galaxies.

## Acknowledgments

The numerical simulations were performed on the computer cluster, XC50 in
NAOJ, and Trinity at Center for Computational Science in University of
Tsukuba. This work is supported in part by MEXT/JSPS KAKENHI Grant Number
17H04827, 20H04724, 21H04489, 22H00149, JST FOREST Program, Grant Number
JPMJFR202Z, ABC project research, Grant Number AB041008 (HY).

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abe et al. (2021) Abe M., Yajima H., Khochfar S., Dalla Vecchia C., Omukai K., 2021, MNRAS, 508, 3226
  * Agertz et al. (2007) Agertz O., et al., 2007, MNRAS, 380, 963
  * Aguirre et al. (2001) Aguirre A., Hernquist L., Schaye J., Katz N., Weinberg D. H., Gardner J., 2001, ApJ, 561, 521
  * Anglés-Alcázar et al. (2017) Anglés-Alcázar D., Faucher-Giguère C.-A., Kereš D., Hopkins P. F., Quataert E., Murray N., 2017, MNRAS, 470, 4698
  * Arata et al. (2020) Arata S., Yajima H., Nagamine K., Abe M., Khochfar S., 2020, MNRAS, 498, 5541
  * Bahé et al. (2022) Bahé Y. M., et al., 2022, MNRAS, 516, 167
  * Benson et al. (2003) Benson A. J., Bower R. G., Frenk C. S., Lacey C. G., Baugh C. M., Cole S., 2003, ApJ, 599, 38
  * Bertone et al. (2007) Bertone S., De Lucia G., Thomas P. A., 2007, MNRAS, 379, 1143
  * Boissier & Prantzos (1999) Boissier S., Prantzos N., 1999, MNRAS, 307, 857
  * Borrow et al. (2020) Borrow J., Anglés-Alcázar D., Davé R., 2020, MNRAS, 491, 6102
  * Braspenning et al. (2023) Braspenning J., Schaye J., Borrow J., Schaller M., 2023, MNRAS, 523, 1280
  * Brook et al. (2014) Brook C. B., Stinson G., Gibson B. K., Shen S., Macciò A. V., Obreja A., Wadsley J., Quinn T., 2014, MNRAS, 443, 3809
  * Brooks et al. (2007) Brooks A. M., Governato F., Booth C. M., Willman B., Gardner J. P., Wadsley J., Stinson G., Quinn T., 2007, ApJ, 655, L17
  * Chiang et al. (2013) Chiang Y.-K., Overzier R., Gebhardt K., 2013, ApJ, 779, 127
  * Christensen et al. (2016) Christensen C. R., Davé R., Governato F., Pontzen A., Brooks A., Munshi F., Quinn T., Wadsley J., 2016, ApJ, 824, 57
  * Christensen et al. (2018) Christensen C. R., Davé R., Brooks A., Quinn T., Shen S., 2018, ApJ, 867, 142
  * Collacchioni et al. (2020) Collacchioni F., Lagos C. D. P., Mitchell P. D., Schaye J., Wisnioski E., Cora S. A., Correa C. A., 2020, MNRAS, 495, 2827
  * Cowie & Songaila (1998) Cowie L. L., Songaila A., 1998, Nature, 394, 44
  * Crain et al. (2015) Crain R. A., et al., 2015, MNRAS, 450, 1937
  * Dalla Vecchia & Schaye (2012) Dalla Vecchia C., Schaye J., 2012, MNRAS, 426, 140
  * Davé et al. (2011) Davé R., Oppenheimer B. D., Finlator K., 2011, MNRAS, 415, 11
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827
  * De Rossi et al. (2017) De Rossi M. E., Bower R. G., Font A. S., Schaye J., Theuns T., 2017, MNRAS, 472, 3354
  * Dekel & Mandelker (2014) Dekel A., Mandelker N., 2014, MNRAS, 444, 2071
  * Dekel et al. (2013) Dekel A., Zolotov A., Tweed D., Cacciato M., Ceverino D., Primack J. R., 2013, MNRAS, 435, 999
  * Donahue & Voit (2022) Donahue M., Voit G. M., 2022, Phys. Rep., 973, 1
  * Dubois et al. (2012) Dubois Y., Devriendt J., Slyz A., Teyssier R., 2012, MNRAS, 420, 2662
  * Dubois et al. (2021) Dubois Y., et al., 2021, A&A, 651, A109
  * Ferland (2000) Ferland G. J., 2000, in Arthur S. J., Brickhouse N. S., Franco J., eds, Revista Mexicana de Astronomia y Astrofisica Conference Series Vol. 9, Revista Mexicana de Astronomia y Astrofisica Conference Series. pp 153–157
  * Finlator & Davé (2008) Finlator K., Davé R., 2008, MNRAS, 385, 2181
  * Ford et al. (2014) Ford A. B., Davé R., Oppenheimer B. D., Katz N., Kollmeier J. A., Thompson R., Weinberg D. H., 2014, MNRAS, 444, 1260
  * Fujimoto et al. (2019) Fujimoto S., et al., 2019, ApJ, 887, 107
  * Genel et al. (2019) Genel S., et al., 2019, ApJ, 871, 21
  * Haardt & Madau (2001) Haardt F., Madau P., 2001, in Neumann D. M., Tran J. T. V., eds, Clusters of Galaxies and the High Redshift Universe Observed in X-rays. (arXiv:astro-ph/0106018)
  * Hafen (2020) Hafen Z. Z., 2020, PhD thesis, Northwestern University
  * Hafen et al. (2019) Hafen Z., et al., 2019, MNRAS, 488, 1248
  * Higuchi et al. (2019) Higuchi R., et al., 2019, ApJ, 879, 28
  * Hill et al. (2020) Hill R., et al., 2020, MNRAS, 495, 3124
  * Hui & Gnedin (1997) Hui L., Gnedin N. Y., 1997, MNRAS, 292, 27
  * Keller et al. (2019) Keller B. W., Wadsley J. W., Wang L., Kruijssen J. M. D., 2019, MNRAS, 482, 2244
  * Lehner et al. (2014) Lehner N., O’Meara J. M., Fox A. J., Howk J. C., Prochaska J. X., Burns V., Armstrong A. A., 2014, ApJ, 788, 119
  * Lochhaas et al. (2020) Lochhaas C., Bryan G. L., Li Y., Li M., Fielding D., 2020, MNRAS, 493, 1461
  * Ma et al. (2016) Ma X., Hopkins P. F., Faucher-Giguère C.-A., Zolman N., Muratov A. L., Kereš D., Quataert E., 2016, MNRAS, 456, 2140
  * Maiolino et al. (2008) Maiolino R., et al., 2008, A&A, 488, 463
  * Mannucci et al. (2009) Mannucci F., et al., 2009, MNRAS, 398, 1915
  * Martin (2005) Martin C. L., 2005, ApJ, 621, 227
  * Matsuda et al. (2005) Matsuda Y., et al., 2005, ApJ, 634, L125
  * McAlpine et al. (2017) McAlpine S., Bower R. G., Harrison C. M., Crain R. A., Schaller M., Schaye J., Theuns T., 2017, MNRAS, 468, 3395
  * Méndez-Hernández et al. (2022) Méndez-Hernández H., et al., 2022, A&A, 666, A56
  * Miller et al. (2018) Miller T. B., et al., 2018, Nature, 556, 469
  * Mitchell & Schaye (2022) Mitchell P. D., Schaye J., 2022, MNRAS, 511, 2948
  * Mitchell et al. (2020) Mitchell P. D., Schaye J., Bower R. G., Crain R. A., 2020, MNRAS, 494, 3971
  * Moll et al. (2007) Moll R., et al., 2007, A&A, 463, 513
  * Muratov et al. (2017) Muratov A. L., et al., 2017, MNRAS, 468, 4170
  * Murray et al. (2005) Murray N., Quataert E., Thompson T. A., 2005, ApJ, 618, 569
  * Nelson et al. (2019) Nelson D., et al., 2019, MNRAS, 490, 3234
  * Oppenheimer & Davé (2006) Oppenheimer B. D., Davé R., 2006, MNRAS, 373, 1265
  * Oppenheimer & Davé (2008) Oppenheimer B. D., Davé R., 2008, MNRAS, 387, 577
  * Oppenheimer et al. (2010) Oppenheimer B. D., Davé R., Kereš D., Fardal M., Katz N., Kollmeier J. A., Weinberg D. H., 2010, MNRAS, 406, 2325
  * Overzier (2016) Overzier R. A., 2016, A&A Rev., 24, 14
  * Pandya et al. (2021) Pandya V., et al., 2021, MNRAS, 508, 2979
  * Pearson et al. (2018) Pearson W. J., et al., 2018, A&A, 615, A146
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Rahmati et al. (2013) Rahmati A., Pawlik A. H., Raicevic M., Schaye J., 2013, MNRAS, 430, 2427
  * Schaller et al. (2015) Schaller M., et al., 2015, MNRAS, 451, 1247
  * Schaye & Dalla Vecchia (2008) Schaye J., Dalla Vecchia C., 2008, MNRAS, 383, 1210
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Shen et al. (2013) Shen S., Madau P., Guedes J., Mayer L., Prochaska J. X., Wadsley J., 2013, ApJ, 765, 89
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 312
  * Springel et al. (2001) Springel V., Yoshida N., White S. D. M., 2001, New Astronomy, 6, 79
  * Steidel et al. (2010) Steidel C. C., Erb D. K., Shapley A. E., Pettini M., Reddy N., Bogosavljević M., Rudie G. C., Rakic O., 2010, ApJ, 717, 289
  * Taylor & Kobayashi (2015) Taylor P., Kobayashi C., 2015, MNRAS, 452, L59
  * Tollet et al. (2019) Tollet É., Cattaneo A., Macciò A. V., Dutton A. A., Kang X., 2019, MNRAS, 485, 2511
  * Toshikawa et al. (2018) Toshikawa J., et al., 2018, PASJ, 70, S12
  * Tumlinson et al. (2017) Tumlinson J., Peeples M. S., Werk J. K., 2017, ARA&A, 55, 389
  * Umehata et al. (2019) Umehata H., et al., 2019, Science, 366, 97
  * Veilleux et al. (2020) Veilleux S., Maiolino R., Bolatto A. D., Aalto S., 2020, A&A Rev., 28, 2
  * Wiersma et al. (2009a) Wiersma R. P. C., Schaye J., Smith B. D., 2009a, MNRAS, 393, 99
  * Wiersma et al. (2009b) Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., Tornatore L., 2009b, MNRAS, 399, 574
  * Wright et al. (2020) Wright R. J., Lagos C. d. P., Power C., Mitchell P. D., 2020, MNRAS, 498, 1668
  * Wright et al. (2021) Wright R. J., Lagos C. d. P., Power C., Correa C. A., 2021, MNRAS, 504, 5702
  * Yajima et al. (2012) Yajima H., Choi J.-H., Nagamine K., 2012, MNRAS, 427, 2889
  * Yajima et al. (2017) Yajima H., Nagamine K., Zhu Q., Khochfar S., Dalla Vecchia C., 2017, ApJ, 846, 30
  * Yajima et al. (2022a) Yajima H., Abe M., Fukushima H., Ono Y., Harikane Y., Ouchi M., Hashimoto T., Khochfar S., 2022a, arXiv e-prints, p. arXiv:2211.12970
  * Yajima et al. (2022b) Yajima H., et al., 2022b, MNRAS, 509, 4037