How Metals Are Transported In And Out Of A Galactic Disk: Dependence On The Hydrodynamic Schemes In Numerical Simulations

Metallicity is a fundamental probe for understanding the baryon physics in a galaxy. Since metals are intricately associated with radiative cooling, star formation, and feedback, reproducing the observed metal distribution through numerical experiments will provide a prominent way to examine our understandings of galactic baryon physics. In this study, we analyze the dependence of the galactic metal distribution on the numerical schemes and quantify the differences in the metal mixing among modern galaxy simulation codes (the mesh-based code Enzo and the particle-based codes Gadget-2 and Gizmo-PSPH). In particular, we examine different stellar feedback strengths and an explicit metal diffusion scheme in particle-based codes, as a way to alleviate the well-known discrepancy in metal transport between mesh-based and particle-based simulations. We demonstrate that a sufficient number of gas particles are needed in the gas halo to properly investigate the metal distribution therein. Including an explicit metal diffusion scheme does not significantly affect the metal distribution in the galactic disk but does change the amount of low-metallicity gas in the hot-diffuse halo. We also find that the spatial distribution of metals depends strongly on how the stellar feedback is modeled. We demonstrate that the previously reported discrepancy in metals between mesh-based and particle-based simulations can be mitigated with our proposed prescription, enabling these simulations to be reliably utilized in the study of metals in galactic halos and the circumgalactic medium.


INTRODUCTION
Metals in galaxies are fundamental probes to understand the baryon physics in the galactic ecosystem. Since metals are mainly produced in stars and released to the interstellar medium (ISM) via supernova (SN) explosion, the metal distribution provides important information on the stellar lifecycle and stellar feedback. Furthermore, metals are not only the passive tracers of the stellar feedback but also one of the main coolants -e.g., C and O lines in the warm and neutral component of the ISM (Dalgarno & McCray 1972) -thus playing a crucial role in star formation. Indeed, the chemical enrichment of a galaxy is the consequence of an intricate interplay of baryonic processes in galaxies, including star for-mation, gas inflow and outflow, and turbulence. Therefore, the metal distribution offers crucial constraints on the processes of galactic evolution.
Observationally, the metal abundance of a galaxy is determined by the flux ratio of line emissions induced by young OB starlight (Tremonti et al. 2004;Nagao et al. 2006;Liu et al. 2008) or by the diffuse ionized gas (DIG; Haffner et al. 1999;Sanders et al. 2017;Kumari et al. 2019). Observations have also revealed a tight global correlation between the stellar mass and gas-phase metallicity (Z; defined as a mass fraction of metals over total gas) such that more massive galaxies are more metal-enriched throughout a wide range of galaxy masses and redshifts (MZR; Lequeux et al. 1979;Tremonti et al. 2004;Mannucci et al. 2010;Sánchez et al. 2019). Moreover, it has been found that the scatter of this relation decreases when star formation rates (SFRs) are also considered , giving rise to the so-called fundamental metallicity relation (FMR) -a plane Z = Z (SFR, M ) in the three-dimensional space of metallicity (Z), stellar mass (M ) -covering a wide mass range down to dwarf galaxies (Lara-López et al. 2010;Yates et al. 2012;Cresci et al. 2019). Meanwhile, spatially resolved ISM observations have revealed a negative radial metallicity gradient; in other words, metals are more abundant in the central region than the outer region (Zaritsky et al. 1994;Swinbank et al. 2012;Jones et al. 2013). This gradient can be explained by the inside-out disk growth that produces a negative stellar population gradient (Matteucci & Francois 1989;Boissier & Prantzos 1999;Prantzos & Boissier 2000). However, a positive or little gradient has also been reported Troncoso et al. 2014;Carton et al. 2018). In particular, interacting galaxies may exhibit lower metallicity in the central region due to the tidally-induced inflow of primordial gas Perez et al. 2011;Troncoso et al. 2014).
In order to test modern theories of galactic baryon physics and metal transport, numerical experiments have been widely used. Simulations have demonstrated that the observed correlation between stellar mass and metallicity can arise naturally in a hierarchical structure formation scenario when the stellar mass growth is regulated by its own negative feedback or mergers (e.g., de Rossi et al. 2007;Davé et al. 2011). The SN-driven outflows are shown to play a key role, especially in low-mass galaxies, as metal-enriched material can escape the galactic potential wells more efficiently, keeping their metallicity low (e.g., Larson 1974;Tremonti et al. 2004;Brooks et al. 2007;De Rossi et al. 2017). Simulations also found that the galactic environment can alter metal distribution. For example, the inflow of the pristine gas dilutes the metal distribution in disks (e.g., Finlator & Davé 2008;Davé et al. 2011); ram pressure can be exerted by the surrounding gas on to the metal-enriched outflow resulting in the ejecta contained within the galactic halo for an extended amount of time (e.g., Ferrara et al. 2005).
Hydrodynamic simulations have also been utilized to investigate more detailed properties related to metals, such as the metal diffusion coupled with the ISM turbulence and the metallicity distribution function (MDF), which are essential to understand the contamination of pristine gas (Pan et al. 2013). As for reproducing the metal mixing in simulations, however, both particle-based codes and mesh-based codes have their own issues. In Lagrangian particle-based schemes such as the smoothed particle hydrodynamics (SPH) approach, metals do not mix between particles unless an explicit diffusion physics is included as a subgrid model. On the other hand, Eulerian mesh-based codes inherently allow metals to diffuse. However, the artificial diffusion required for stable hydrodynamical solutions may over-mix fluids in simulations with insufficient resolution (e.g., Pan et al. 2013;Springel 2016), especially in systems moving relative to the grid (e.g., Pontzen et al. 2020). Naturally, several authors have reported different methods to mitigate the weaknesses in these numerical approaches. For particle-based codes, the turbulent diffusion scheme -calculated with velocity dispersion (e.g., Greif et al. 2009;Revaz et al. 2016) or velocity shear (e.g., Shen et al. 2010;Brook et al. 2012;Su et al. 2017;Escala et al. 2018) following Smagorinsky (1963) are widely used to estimate the effects of subgrid diffusion and to investigate the MDF (see also Hu & Chiang (2020) that took a different approach). Some authors smooth the metallicity over the SPH kernel when they compute, e.g., the cooling rates for gas particles, not actually performing but mimicking the diffusion of metals (Wiersma et al. 2009). For mesh-based codes, the over-mixing in unresolved eddies has been addressed with a turbulent diffusion model based on the probability distribution function (PDF) method (e.g., Pan et al. 2013;Sarmento et al. 2017), a different subgrid model based on a partial differential equation for energy density (e.g., Schmidt et al. 2014;Schmidt 2015), or a "velocityzeroed" initial conditions (Pontzen et al. 2020).
Since metal distribution is highly sensitive to the hydrodynamic scheme and the diffusion model, careful attention on these numerical methods is imperative for studying chemodynamical processes in a galaxy using numerical simulations. The numerical galaxy formation community has collectively responded to this need over the years, and one such effort was the code comparison project AGORA (Kim et al. 2014(Kim et al. , 2016). The AGORA Collaboration tested the reproducibility of numerical experiments using common initial conditions for a dark matter-only cosmological simulation (Kim et al. 2014) and an idealized, isolated galaxy (Kim et al. 2016), providing insights into both the similarities and differences between contemporary simulation codes. While reporting solid convergence in many galactic properties, Kim et al. (2016) also pointed out a discrepancy in the metal distribution between mesh-based and particle-based codes in a test with an idealized disk galaxy. For example, the metal content of the hot-diffuse halo gas is captured in mesh-based codes, whereas, by design, gas particles are scarce in the halo in particle-based simulations. With neither the halo gas nor an explicit metal mixing scheme included (a design choice in Kim et al. 2016), metal-enriched gas particles tend to stay only near the dense star-forming regions in particle-based simulations (see Figures 32 and 33 in Kim et al. 2016). Therefore, in this paper, following up on the AGORA results, we aim to investigate what causes these inter-code discrepancies in metal distribution and how they can be alleviated. We quantitatively compare the metal distribution using various hydrodynamic schemes and analyze several factors that could lead to the discrepancy between mesh-based and particle-based codes -such as the absence of gas particles in the halo region or the lack of diffusion schemes in the particle-based codes. We will examine the metal distribu- NOTE-The parameters for the gas halo are applicable only for the GADGET-2 and GIZMO-PSPH runs in which the gas halo is included (i.e., those with "GasHalo" in their run names in Table 2; see Section 2.2.2). In mesh-based code ENZO, the gas halo is included as a uniform medium around the disk (see Section 2.2.2). All other parameters follow the default disk galaxy initial condition provided by the AGORA Project (Kim et al. 2016). For more information on the parameters listed above, see Section 2.2. tion in GADGET-2 simulations with different hydrodynamic methods and feedback strengths and compare its metal distribution with that of an ENZO simulation. We will also test our proposed prescription to alleviate the inter-code discrepancy in both GADGET-2 and GIZMO-PSPH and discuss its general applicability to particle-based codes. This paper is structured as follows. In Section 2, we describe the simulation codes, the initial conditions, and the physics models adopted. Section 3 compares the metal distributions resulting from different simulation setups (with GADGET-2 and ENZO), focusing on the effect of the initial conditions and the stellar feedback model. In Section 4, we discuss the effect of the metal diffusion scheme in particlebased simulations. Then, in Section 5, we test our proposed prescription to alleviate the inter-code discrepancy in another particle-based code GIZMO-PSPH and discuss its general applicability. Finally, we summarize our findings and discuss future work in Section 6.

Simulation Codes
In this study, we adopt three gravito-hydrodynamics codes widely used in numerical galaxy formation. We briefly explain the physics and key runtime parameters in each code.
2.1.1. Mesh-based Code: ENZO ENZO is an Eulerian 3-dimensional structured mesh code with the adaptive mesh refinement (AMR) capability. The particle-mesh method is used to compute the gravitational interaction (Hockney & Eastwood 1988) while gas dynamics is solved using the 3rd-order accurate piecewise parabolic method (PPM; Colella & Woodward 1984). In this paper, the ENZO simulation uses a 64 3 initial root grid to cover a (1.311 Mpc) 3 simulation box, achieving an 80 pc spatial resolution with eight additional levels of AMR. A cell is refined by a factor of 2 when the mass of the cell is above m gas,IC = 8.593 × 10 4 M in gas mass, or 8 × m ,IC = 8 × 3.437 × 10 5 M in collisionless particle mass. Other adopted schemes are largely in line with the recent numerical studies using ENZO (e.g., Kim et al. 2019;Shin et al. 2020).
2.1.2. Particle-based Code: GADGET-2 and GIZMO-PSPH GADGET-2 is a tree-particle-mesh SPH code developed by Springel (2005), utilizing a standard density-entropy SPH that manifestly conserves energy, entropy, momentum, and angular momentum. Yet, the purely density-based SPH scheme may give rise to fictitious pressure on the interface between two media with extreme density contrast (Agertz et al. 2007) or damped subsonic turbulence (Bauer & Springel 2012). Later variants of GADGET-2 have improved the modelings of complex flows, shocks, or instabilities. In this study, we test the original GADGET-2 rather than any specific variants to focus on the fundamental properties of SPH so that our proposed prescription to alleviate the inter-code discrepancy could be applied to later variants. We also test another particle-based code, GIZMO (Hopkins 2015), which includes various hydrodynamic solvers treating the volume components of simulation differently: densityentropy formalism, pressure-energy formalism, meshless finite mass (MFM), meshless finite volume (MFV), or Eulerian fixed grid schemes. For the present study, we experiment with the pressure-energy SPH (hereafter PSPH; Hopkins 2013), which better captures the instability on the surface between different phases than the traditional SPH codes. This means that, in the treatment of fluids, the performance of our chosen GIZMO(-PSPH) matches that of a contemporary SPH code such as GADGET-3 widely utilized in the community.
In both GADGET-2 and GIZMO-PSPH, we use the cubic spline kernel (Hernquist & Katz 1989) for the softening of the gravitational force with the desired number of neighboring particles N ngb = 32. We adopt the Plummer equivalent gravitational softening length ε grav of 80 pc and allow the hydrodynamic smoothing length to reach a minimum of 0.2ε grav . In both codes, we include the radiative cooling, star formation and feedback following Kim et al. (2016). We also implement an explicit metal diffusion scheme following Hopkins

Initial Condition
We have adopted the disk initial condition (IC) provided by the AGORA Project (Kim et al. 2016) that models an idealized Milky Way-mass disk galaxy of M 200, crit = 1.074 × 10 12 M in isolation. The IC includes a dark matter halo that follows the Navarro-Frenk-White profile (NFW; Navarro et al. 1997), an exponential disk of stars and gas, and a stellar bulge following the Hernquist profile (Hernquist 1990). 1 The detailed structural parameters of the IC are listed in Table 1.
In what follows, we discuss, in particular, the gas distribution in the original AGORA IC and our modified IC.

Gas Distribution In The Original AGORA IC
For the mesh-based code ENZO, the gas density field of the disk is initialized with an exact analytic formula where r is the cylindrical radius, z is the vertical distance from the disk plane, r d = 3.432 kpc, z d = 0.1 r d , and ρ 0 = M d, gas /(4πr 2 d z d ) with M d, gas = 8.593 × 10 9 M . On the other hand, for the particle-based codes GADGET-2 and GIZMO-PSPH, the AGORA IC provides a text file of the initial positions of gas particles in the disk that follow Eq.(1). Notably, the gas particles are absent in the halo region. For both mesh-based and particle-based codes, the initial temperature in the disk is set to 10 4 K and the initial metallicity to Z disk = 0.02041. Finally, for the mesh-based code only, the gas density in the halo is set to a constant n H = 10 −6 cm −3 to avoid a zero value in the cells, with an initial metallicity of Z halo = 10 −6 Z disk and temperature of 10 6 K.

Gas Halo In The Modified AGORA IC
In the original AGORA IC employed in Kim et al. (2016) to model an idealized Milky Way-mass disk galaxy, a gas halo is included only in mesh-based codes, albeit with a negligible density n H = 10 −6 cm −3 as mentioned in the previous section. In this paper, we demonstrate that this small difference can cause substantial discrepancies in the baryonic properties of a galaxy during the 500 Myr evolution, especially in the halo. In this section, we will first argue that a sufficientlyresolved gas halo is necessary for particle-based simulations to properly model the metal transport and mixing in the circumgalactic medium (CGM), and explain the modified IC adopted in some of our simulations.
Consider a situation in which the disk gas is being pushed by strong SN winds and is on the verge of leaving the disk's ISM into the halo. Without any gas particles in the halo (e.g., in the particle-based simulations described in Kim et al. 2016), the supersonic gas outflows do not experience any pressure that impedes its motion. Hence the outflow continues to move into the vacuum while losing little or no momentum. On the contrary, in simulations with a gas halo, the outflow will be slowed down or sometimes even severely suppressed. As we will demonstrate in later sections, simulators or those who analyze simulations should be cautioned that this discrepancy may cause substantial deviations in the baryonic properties in the ISM and CGM. In addition, without the gas halo, the metals in disk gas particles have no particle to diffuse into in the halo region. Therefore, until the halo is populated with (a few) gas particles expelled from the disk by SN winds, the halo rarely becomes metal-enriched.
For these reasons, in this study, we test a different IC to model an isolated galaxy with the particle-based codes GADGET-2 and GIZMO-PSPH, that allocates additional gas particles in the halo region. Gas particles are placed in the halo following the NFW profile in a way that they approximately match the initial halo gas density in the mesh-based code's IC, n H = 10 −6 cm −3 . We match the mass of an individual gas particle in the halo to that in the disk, m gas,IC = 8.593 × 10 4 M , resulting in a total of 4,000 gas particles in the halo (see Table 1). The initial metallicity is set to Z halo = 10 −6 Z disk and the temperature to 10 6 K. The particlebased simulations that utilize this revised IC are denoted with "GasHalo" in their run names in Table 2. The initial metal distributions of the original AGORA IC and the revised IC are displayed in Figure 1, along with that for the mesh-based codes. In the Gad2 run's IC (identical to the particle-based codes' ICs in Kim et al. 2016) the halo is free of gas (vertical height z > 5 kpc), while the Gad2-GasHalo and Enzo runs feature nonzero gas density and metals in the halo region.

Baryon Physics
We consider all the baryon physics that are relevant in the process of galaxy formation by closely following the previous AGORA disk comparison (Kim et al. 2016), along with an optional scheme for explicit metal diffusion.

Cooling, UV Background, Jeans Pressure Support
The radiative cooling and heating rates for the gas are calculated with the GRACKLE library (Smith et al. 2017). We adopt GRACKLE's ionization equilibrium mode with the Haardt & Madau (2012) UV background radiation at z = 0 -i.e., the gas cooling rate is determined by its density, temperature and metallicity in the ionization levels satisfying the equilibrium state using CLOUDY (Ferland et al. 2013). In addition, we include the Jeans pressure floor to avoid any artificial collapse and numerical fragmentation (Truelove et al. 1997). The Jeans pressure is determined as where γ = 5/3 is the adiabatic index, G is the gravitational constant, and ρ gas is the gas density. Here, ∆x is equivalent to the spatial resolution (or its proxy) carried by each simulation code -that is, the finest cell size in ENZO, the smoothing length h sml in GADGET-2, and the radius of the "effective volume" of a cell in GIZMO-PSPH, (4π/(3N ngb )) 1/3 h sml . Correspondingly, we set the controlling parameter N Jeans to 4.0, 4.2 and 6.3 for ENZO, GADGET-2 and GIZMO-PSPH, respectively, to provide a similar amount of pressure across the codes. These choices of N Jeans are in line with Kim et al. (2016, see their Section 3.1), and lead to model galaxies producing similar stellar masses of ∼ 10 9 M in the first 500 Myr (see Table 2; see also    Table 2). The Gad2 run does not contain any gas in the halo region while the initial gas distribution in Gad2-GasHalo is approximately comparable with that in Enzo. See Section 2.2.2 for more information.

Star Formation and Feedback
Gas parcels that are denser than a threshold, n H = 10 cm −3 , spawn stars at a rate following the local Schmidt law where ρ is the stellar density, t ff = (3π/(32 Gρ gas )) 1/2 is the local free-fall time, and ε = 0.01 is the star formation efficiency per free-fall time. 5 Myr after their formation, star particles inject thermal energy, mass, and metals into their surrounding ISM in an attempt to describe Type II SN explosions. Following Chabrier (2003) initial mass function (IMF), we assume that for stars with a mass range of 8 -40 M , a single SN event occurs per erevery 91 M of stellar mass formed, releasing 2.63 M of metals and 14.8 M of gas (including metals). To probe the difference in the efficiency of stellar feedback between mesh-based and particlebased codes -especially in the context of metal transport -we test various thermal energy values of the stellar feedback: 10 51 , 2×10 51 , and 3×10 51 ergs per SN, labelled as the Gad2-GasHalo/+TFB2/+TFB3 run (see Table 2).

Explicit Metal Diffusion In Particle-based Simulations
In particle-based simulations, once a metal field is assigned to a particle in the IC, its value never changes unless the particle is directly affected by a SN bubble. This means that a naive particle-based approach does not capture the interparticle diffusion of metals. Evidently, many particle-based code groups studying metal transport in a galaxy-scale simulation have devised a way to model how the metal is mixed. The diffusion scheme has also been shown indispensable to  Table 2 for the list of our simulations, and Section 3.1 for more information on this figure. match the observed scatters of metal element abundances, such as alpha and r-process elements (see, e.g., Revaz et al. 2016;Escala et al. 2018;Dvorkin et al. 2020). In this study, we consider an explicit turbulent metal diffusion scheme in GADGET-2 and GIZMO-PSPH, and compare the metal distribution in the galactic disk and halo, with and without the scheme.
We adopt the metal diffusion scheme used in Hopkins et al. (2018) and Escala et al. (2018), which itself is based the Smagorinsky-Lilly model (Smagorinsky 1963;Shen et al. 2010). In brief, the model estimates the subgrid diffusion effect driven by the velocity shear between particles, assuming that the local diffusivity is dependent on the velocity shear and the resolution scale. That is, the metal diffusion between the particles in a shear motion with respect to each other is where M i is the scalar field (metallicity) of the i-th particle, h is the effective measurement scale (which we choose to set to the SPH kernel size h sml ), || · || is the Frobenius norm, and C d is the diffusion coefficient that is proportional to the Smagorinsky-Lilly constant, calibrated by numerical simulations based on the Kolmogorov theory. And the symmetric trace-free tensor S i j is given by where i, j = {x, y, z}, v i is the velocity vector for each gas particle, and, x i is the spatial coordinate. Thus, the turbulent diffusion becomes efficient in large eddies where the shear drives local fluid instabilities. We caution the readers that this simplistic diffusion model is always dissipative, and the backscattering from small to large scales is not possible. The single controlling parameter in this model, the diffusion coefficient C d , has difficulty in properly describing various types of turbulent flows, either. For example, the model may overestimate the diffusion in the laminar flows (Rennehan et al. 2019). Additionally, Colbrook et al. (2017) showed that the turbulent diffusivity is dependent on the scale of eddies. Despite these limitations, we adopt the Smagorinsky-Lilly model to show that including such a simple diffusion model can mitigate the known discrepancy between mesh-based and particle-based codes. Different authors have used different values for C d , ranging from 0.003 in Escala et al. (2018) to 0.05 in Shen et al. (2010). 2 Given such a wide range of values found in the literature, here we test C d = 0.006, 0.02, and 0.06 which are labelled as the Gad2-diff0.3, Gad2-diff1, and Gad2-diff3 run, respectively (see Table 2).

METAL DISTRIBUTION IN HALO: DEPENDENCE ON INITIAL GAS DISTRIBUTION AND FEEDBACK STRENGTH
Using a suite of simulations listed in Table 2, we compare how the spatial distribution of metals differs between different types of simulations -first focusing on how the existence of a halo gas and the different feedback strengths change the extent of transported metals in particle-based simulations. Figure 3. 500 Myr snapshots of our isolated disk simulations using different ICs, stellar feedback, and diffusion schemes. Face-on (1st row), edge-on (2nd row) and wider edge-on (3rd row) projections of metal density, and density-weighted edge-on projection of the vertical velocity (4th row). The GADGET-2 simulation without the halo gas -the Gad2 run (see Table 2); the same runtime condition as the particle-code runs in Kim et al. (2016) -substantially differs from the mesh-based Enzo run in the metal distribution in the halo. In contrast, another GADGET-2 simulation, but this time with a gas halo, more feedback energy, and explicit metal diffusion scheme -the Gad2-GasHalo+TFB2+diff1 run -is comparable with the Enzo run. See Table 2 for the list of our simulations, and Section 3.2 for more information on this figure.

Comparison of Star Formation Rates
Metals in galaxies are produced by stars. Therefore, in order to compare the spatial distribution of metals between different simulations, it is necessary to establish a baseline in which all simulations exhibit similar star formation histories in the timespan considered. In Figure 2, we show the SFRs in simulations of different hydrodynamic solvers, different ICs, different feedback strengths, and different diffusion coefficients. As noted in Section 2.3.1, the Jeans pressure support for each code (ENZO, GADGET-2, GIZMO-PSPH) are set in such a way that the runs show similar SFRs and the value N Jeans is in line with the previous AGORA comparison (Kim et al. 2016). As a result, despite the differences in the simulation setup, most of the runs analyzed in this article exhibit similar star formation histories within a few tens of percents at all times, totaling a stellar mass of ∼ 10 9 M in the first 500 Myrs (see the rightmost column of Table 2). 3 Note that even the seemingly important differences such as in the stellar feedback strength or in the diffusion coefficient introduce only marginal changes in SFRs in Figure 2. For example, the run with higher thermal energy suppresses star formation only slightly when compared with the one with lower energy (compare the Gad2-GasHalo/+TFB2/+TFB3 runs in Table 2 and Figure 2). Meanwhile, the coefficient of metal diffusion does not substantially affect the produced stellar mass (compare the Gad2-(GasHalo+TFB2)+diff0.3/diff1/diff3 runs in Table 2). Since all the simulation produce a similar amount of stars -and thus metals -we can now conjecture that any difference in the spatial distribution of metals is due to the Here, the cylindrical profile is made for the disk gas only, defined as the region where the vertical distance from the disk plane z is less than 5 kpc. Most of the runs exhibit similar profiles in the cylindrical radial direction, while the vertical distribution of metals is highly dependent on the simulation setups. See Section 3.3 for more information on this figure. difference in how each simulation transports metals in and out of the galactic disk (e.g., different feedback strength, inherent differences in hydrodynamics), not because any one simulation harbors a larger/smaller amount of metals. Figure 3 displays the face-on and edge-on projections of metal density and the vertical velocity of gas outflows from the disk plane at 500 Myr after the simulation starts. In the face-on view, all simulations show a similar distribution of metals along the spiral arms. The edge-on view of metal distribution, however, varies substantially depending on the simulation setup. The first stark contrast is between the Enzo run and the Gad2 run (1st and 2nd column in Figure 3, respectively). As previously reported by the AGORA Collaboration (Kim et al. 2016), in a particle-based simulation with neither the gas halo nor an explicit diffusion scheme, metals are inevitably scarce in the halo away from the disk (i.e., 2nd and 3rd rows of the Gad2 run). This is because the halo is only populated with very few (metal-enriched) gas particles ejected from the disk -as discussed in Section 2.2.2 and will become more obvious in later sections. The enclosed metal mass in the halo region (vertical distance from the disk z > 5 kpc) is about 30 times lower in the Gad2 run than in the Enzo run (see also the left panel of Figure 5; to be discussed in details in Section 3.3). Despite having a very metal-poor halo, the Gad2 run shows the fastest gas outflows from the disk among all the runs, in the bottom row of Figure 3 that displays the density-weighted projection of vertical velocity. This counter-intuitive result is due to the unphysical nature of the Gad2 run's IC, and to the fact that the halo region contains only a few gas particles with an extremely high velocity (see also the right panel of Figure 4; to be further discussed in details in Section 3.3). In the absence of gas in the halo region, the high-velocity SN ejecta travels into the halo with-out experiencing any pressure that impedes its motion. Not suffering any deceleration, these high-velocity gas particles may reach hundreds of kiloparsecs away from the disk. Therefore, to rectify the unphysical results of the Gad2 run, another IC has been tested, which now includes the gas halo around the disk -i.e., the Gad2-GasHalo run (3rd column in Figure 3; see also Section 2.2.2). As we compare the Gad2 and the Gad2-GasHalo run, we first find that in terms of metals in the halo, the Gad2-GasHalo run is hardly different from the Gad2 run (3rd row). This is because the additional halo gas particles still cannot receive metals unless there is an explicit way for the metals to diffuse into the halo. It is also because the SN ejecta cannot easily penetrate the gas halo as it did in the Gad2 run. The halo gas applies ram pressure on the gas outflows at the disk-halo boundary and restricts the reach of the metal-enriched ejecta. 4 As a result, the galactic outflow becomes very weak in the Gad2-GasHalo run (bottom row of Figure 3), and the metal-enriched ejecta remains near the galactic disk (3rd row).

Overview: Metal Distribution and Outflow Velocity
In other words, our experiment suggests that particlebased codes may require more stellar feedback energy than mesh-based codes to launch galactic outflows into the gas halo. Indeed, the GADGET-2 simulation with twice the thermal feedback energy -i.e., the Gad2-GasHalo+TFB2 run (4th column in Figure 3) -shows a similar metal distribution and outflow velocity to the Enzo run. 5 Comparing the runs with varying thermal feedback energies -Gad2-GasHalo/+TFB2/+TFB3 runs -we find that the 4 As noted in Section 2.1.2, the original GADGET-2 tends to suppress fluid instabilities due to the shielding effect between the two media with extreme contrast in density (Agertz et al. 2007). Later variants of GADGET-2 have improved to capture such instabilities, which can affect the penetration of high-velocity winds into the halo (Hopkins 2013). 5 Although the outflow is with a larger opening angle than in the Enzo run.  Here, the cylindrical profile is made for the disk gas only, defined as the region where the vertical distance from the disk plane z is less than 5 kpc. The shapes of the metallicity profiles are affected by the existence of the halo gas and by the amount of stellar feedback energy. See Section 3.3 for more information on this figure. metal enrichment in the halo is highly sensitive to the feedback strength. The amount of thermal energy injected directly determines the momentum of SN ejecta, and consequently, the mass of metal-enriched gas in more turbulent bubbles. The increased turbulence enables more metalenriched gas to be coupled with large momentum, allowing the SN ejecta to escape from the disk easily. In particle-based simulations, this process may require more energy than in mesh-based ones, due to the inherent inter-code discrepancies in how the thermal feedback energy is distributed in the neighborhood of newly-born stars, and how the Riemann problem is solved at the disk-halo boundary. For fluids in vastly different phase -e.g., SN hot bubbles in the colddense gas clouds -in particle-based simulations, the density of the dilute fluid can be overestimated by the SPH kernels in insufficient resolution, which gives rise to overcooling and inhibiting the development of hot gas (Marri & White 2003;Creasey et al. 2011). 6 Lastly, we include the explicit metal diffusion scheme (see Section 2.3.3) to allow metals of highly-enriched gas particles to slowly diffuse into the less-enriched neighbors. Comparing the Gad2-GasHalo+TFB2 and Gad2-GasHalo+TFB2+diff1 run in Figure 3, we discover 6 Note that we have only tested the thermal feedback prescription based on Kim et al. (2016) Figure 7. Two-dimensional probability distribution functions (PDFs) of metal mass on the density-temperature plane for different simulation setups at 500 Myr. The color represents the mass of the metals in each two-dimensional bin. Only with the gas halo sufficiently resolved in the IC (i.e., Gad2-GasHalo/+TFB2/+TFB3/+TFB2+diff1 runs) can the particle-based simulations identify the hot diffuse gas surrounding the disk at [∼ 10 −29 g cm −3 , ∼ 10 6 K]. The metal mass distribution in the PDF, especially in the halo, heavily depends on the thermal stellar feedback energy as well. See Table 2 for more information on our list of simulations, and Section 3.4 for more information on this figure. that the metal distribution in space does not highly depend on the diffusion scheme. However, in Section 4.1 we will demonstrate why the diffusion scheme must be included.

Spatial Profiles of Metals and Its Evolution In Time
Thus far, using Figure 3 we have shown that the presence of the gas particles in the halo region, despite its negligible density, affects the metal distribution therein. We have also demonstrated that the amount of metals expelled from the disk depends on the stellar feedback energy. In this subsection, we further investigate these points quantitatively. Figure 4 illustrates the density-weighted metal density profiles and the velocity of gas outflows perpendicular to the disk plane 500 Myrs after the simulation starts. As observed in Figure 3, the metal density profiles in the disk's radial direction (left panel of Figure 4) are similar across different simulation setups. In contrast, the metal density distribution in the halo and the amount of metal transported to the halo are notably different between the runs. The middle and right panels of Figure 4 illustrate this difference. Here, the height range is chosen to be between 5 and 200 kpc from the galactic disk in order to avoid including the disk gas in our analysis. The metal density in the Enzo run decreases smoothly out to z ∼ 150 kpc -the edge of the metal-enriched halo gas -at which point the density drops sharply. The Gad2 run without a gas halo in the IC shows only a few discrete points, indicating that only a small number of gas particles have been ejected into and remained in the halo. These discrete points have high outflow velocities, nearly 500 − 1000 km s −1 , as they do not have to move through any medium that decelerates the outflow. In contrast, once a gas halo is included in the IC (Gad2-GasHalo run), the outflow velocity can reach only up to a few of ∼ 50 km s −1 . Finally, comparing the runs with a gas halo but with different thermal feedback energies -Gad2-GasHalo/+TFB2/+TFB3 runs -we find that the extent of metal-enriched gas is dictated by the feedback strength. The higher the feedback energy is, the faster the gas outflow becomes, enriching a larger volume of the halo. The inclusion of the diffusion scheme does not affect the spatially averaged distribution of metals. Figure 5 displays the enclosed metal masses as functions of the vertical height (left panel) and their time evolution (right panel). We can again observe that the simulation with more feedback energy transports more metals from the disk to the halo. In terms of the total metal mass in the halo, the mesh-based Enzo run is most compatible with the Gad2-GasHalo+TFB2(+diff1) run in both panels (as men- Figure 8. The probability distribution functions (PDFs) of metal mass in gas density (left panel) and temperature (right panel), for different simulation setups at 500 Myr. Only with the gas halo sufficiently resolved in the IC (i.e., Gad2-GasHalo/+TFB2/+TFB3/+TFB2+diff1 runs) can the particle-based simulations identify the hot diffuse gas in the halo at 10 −30 − 10 −27 g cm −3 and 10 4 − 10 7 K. The amount of metal mass in this hot diffuse medium is dictated by the strength of stellar feedback. See Section 3.4 for more information on this figure. tioned in Section 3.2 for Figure 3). Without a gas halo (Gad2 run), the few metal-enriched gas particles rarely stay in the halo due to their high velocity, yielding an unrealistically metal-poor halo throughout the simulation.
In the right panel of Figure 5, the role of a gas halo in containing the SN ejecta is again illustrated. The metal mass in the region of R 200 < r < 2R 200 in the Gad2 run (where R 200 = 205.5 kpc) is shown with a thin dotted line, and this indicates that a few high-velocity SN ejecta particles have escaped the virial radius. They occasionally -and unphysically -reach thousands of kiloparsecs away from the galactic center. On the contrary, due to the presence of the halo gas, no ejected particle escapes the virial radius in all other runs. The gas halo, even when its density is negligible, imposes pressure on the gas outflows and decelerates them. The confinement of metal-enriched outflows has been proposed by Ferrara et al. (2005), who suggested that the gas surrounding the galactic disk exerts ram pressure on to the outflows so that the ejected metals are in a hot-diffuse phase.
Finally, in Figure 6, we present the mass-weighted gas metallicity profiles in both the disk's radial direction and the vertical direction from the disk plane. In the cylindrical radial direction (left panel), the metallicity in most runs is near the initial disk metallicity Z disk = 0.02041 (Section 2.2.1), and drops sharply at r ∼ 25 kpc -the edge of the galactic disk. However, the (mass-weighted) metallicity is higher in all GADGET-2 runs in the galactic core than in the Enzo run. It is because, in the particle-based simulations, metals tend to be locked in the dense region before they slowly disperse or diffuse into less dense regions. Meanwhile, the vertical metallicity profiles (right panel of Figure 6) show a similar trend to the metal density profiles (the middle of Figure 4) with one exception, the Gad2 run. In the Gad2 run, the halo is insufficiently resolved with only a few high-velocity gas particles ejected from the metal-enriched star-forming regions; thus, the metal fraction in this region is not reliable.

Metal Distribution in The Density-Temperature Plane
We now investigate the metal distribution in the densitytemperature phase space. In Figure 7, we draw the twodimensional probability distribution functions (PDFs) of metal mass for various simulation setups with ENZO and GADGET-2. The one-dimensional projections along one of the axes -i.e., density PDF and temperature PDF -are shown in Figure 8 for more quantitative comparison.
In Figure 7, for all the runs considered, the majority of metal masses are on the thermal equilibrium curve stretching from ∼ 10 4 K to ∼ 10 2 K where cooling and heating rates are equal. Three distinct phases of gas -cold (< 10 3 K), warm (10 3−5 K), and hot phase (> 10 5 K) -are visible in all panels except Gad2 (see also the right panel in Figure 8). The gas surrounding the disk (missing in the Gad2 run) is fed with hot gas particles expelled by SNe, and in turn, exerts pressure on the galactic outflow. The dilute gas with varying entropy and pressure present in the Gad2 run, is now collapsed and confined to a constant pressure line between 10 −30 and 10 −27 g cm −3 in the Enzo and Gad2-GasHalo runs. This constant pressure phase develops via the pressure balance between the disk and the halo gas, and subsequently, a hot diluted halo is built. The absence of a radiation channel via line emission at ∼ 10 6 K thus creates a hot galactic halo in a thermodynamic equilibrium (Ferrara et al. 2005). The metals in this hot phase gas can be hard to detect, potentially presenting a solution for the missing metal problem.
As we have discussed previously, the metal distribution in the galactic halo is greatly affected by the strength of thermal Figure 9. The metallicity distribution functions (MDFs) for the disk gas (left panel) and the halo gas (right panel), in different simulation setups at 500 Myr. Here, the disk is defined as the region where the vertical distance from the disk plane z is less than 5 kpc. The halo is defined as the region where z > 5 kpc or the radial distance from the disk center r > 60 kpc. All GADGET-2 runs are with the gas halo and the thermal feedback energy of 2 × 10 51 ergs per SN (i.e., GasHalo+TFB2). Only with a sufficiently high metal diffusion coefficient C d (≥ 0.02, i.e., Gad2-(GasHalo+TFB2)+diff1/+diff3 runs) can the particle-based simulations fill the domain of low-metallicity gas between Z = 0 and Z disk = 0.02041. See Table 2 for the list of our simulations, and Section 4.1 for more information on this figure. stellar feedback. Comparing the runs with different thermal feedback energies -Gad2-GasHalo/+TFB2/+TFB3 runsin Figure 8, one can observe that the run with higher energy transports more metals to the hot-diffuse region. In terms of the density and temperature PDFs in Figure 8, the Gad2-GasHalo+TFB2(+diff1) run is the most compatible with the mesh-based Enzo run (as discussed in Section 3.2 for Figure 3, and in Section 3.3 for Figure 5). 7 The inclusion of the diffusion scheme does not significantly change the metal distribution in the density or temperature phase space.

METAL DISTRIBUTION IN HALO: DEPENDENCE ON EXPLICIT METAL DIFFUSION SCHEMES
In this section, we investigate how the explicit turbulent metal diffusion scheme changes the metal content in a galaxy simulated with particle-based codes. We test different values for the metal diffusion coefficient (C d in Section 2.3.3; see also Table 2) with a fixed stellar feedback model identical to Gad2-GasHalo+TFB2 that is shown to exhibit similar halo properties to the Enzo run in Section 3. Figure 9 presents the metallicity distribution functions (MDFs; metallicity PDFs) for the disk (left panel) and the halo gas (right panel). We compare simulations with four diffusion coefficients, C d = 0, 0.006, 0.02 and 0.06, labelled as Gad2-GasHalo+TFB2, Gad2-(GasHalo+TFB2)+diff0.3, Gad2-diff1, and Gad2-diff3 run, respectively (see Section 2.3.3 and Table 2 for more information). Since we have calibrated these runs so that the SFRs are similar (Figure 2), and their averaged metal profiles and metal masses are comparable (Figures 4, 5 and 8), we can conjecture that any discrepancy we see in the MDF is caused by varying the diffusion coefficient.

Metallicity Distribution Function
In the Gad2-GasHalo+TFB2 run that does not include a diffusion scheme, the MDF (for both the disk and the halo) shows one sharp peak at Z 0 and another broader peak starting at Z 0.02. The two values correspond to the initial metallicity values of the halo (Z halo = 10 −6 Z disk ) and the disk (Z disk = 0.02041), respectively. Unless an explicit diffusion scheme is used, gas particles residing only within the SN bubbles can acquire metals in a particle-based simulation. The gas in the outer region, away from the star-forming regions, never receives or loses metals, making the MDF overly inhomogeneous. Then, by increasing the diffusion strength, we can see that the gas metallicity is more evenly distributed. For example, in the Gad2-diff1/-diff3 runs, the gap between Z = 0 and Z disk is now filled, to a level of the Enzo run. The Gad2-diff1 run with C d = 0.02 shows the most similar MDF to the Enzo run, 8 although none of the GADGET-2 runs produces the high peak at Z = Z disk in Enzo's MDF for Figure 10. Projected volume of metal-poor gas (Z < 10 −3 ) from the disk's edge-on angle at 100, 300, and 500 Myr. We perform the projection in the cylindrical volume of 30 kpc in radius and 200 kpc in height, centered on the galactic disk in the z = 0 plane. The metal-enriched gas particles ejected by the stellar feedback gradually reduces the volume of metal-poor gas in time. By comparing the two particle-based simulations with and without the metal diffusion scheme, Gad2-GasHalo+TFB2 (left panel) and Gad2-GasHalo+TFB2+diff1 (middle panel), one can see that the scheme helps to enrich the halo with metals more effectively. See Section 4.2 for more information on this figure.
the halo. In the meantime, the run with a weaker diffusion strength, the Gad2-diff0.3 run, fails to fully populate the domain between Z = 0 and Z disk .
Therefore, one can conclude that an explicit metal diffusion scheme is essential in making the realistic lowmetallicity gas in particle-based simulations. Our results suggest that the shape of an MDF is highly sensitive to the diffusion strength, both in the disk and in the halo. In particular, the diffusion scheme can help efficiently enrich the pristine gas at rest in the halo when it is penetrated by high-velocity metal-enriched outflows. It is because, in the Smagorinsky-Lilly model, the flux of the metal field is proportional to its gradient and the velocity shear between the pockets of gas (see Section 2.3.3). We however note that even the diffusion scheme has hard time mitigating the discrepancy at the high metallicity end (Z > 0.02 in Figure 9) between the GADGET-2 and ENZO runs. This is related to another discrepancy seen in the left panel of Figure 6 at small cylindrical radii. The gas in the star-forming region remains metal-rich, and the metals therein are not readily dispersed even with the help of the explicit diffusion scheme (even with high C d ).

Pockets of Metal-poor Gas
We now look into the spatial distribution of pristine, metalpoor gas. Finding pockets of pristine gas -if any -that survived the metal contamination by its host galaxy has important implications in many studies in astrophysics: observations of the metal absorption lines in the CGM (e.g., Roca-Fàbrega et al. 2019;Strawn et al. 2020), the search for possible birthplaces for massive stars (e.g., Turk et al. 2009;Sarmento et al. 2017), and the search for massive black holes stemming (arguably) from the merging of massive stars (e.g., Belczynski et al. 2010), among others. Therefore, here we particularly focus on the metal-poor gas in the halo and investigate how its volume changes due to the inclusion of the explicit metal diffusion scheme. Figure 10 displays the projected volume of metal-poor gas (Z < 10 −3 ) from the disk's edge-on angle for different simulation setups at 100, 300, and 500 Myr. At 0 Myr, all metalpoor gas is by design only in the halo. As the galaxy evolves in time, the SN-driven winds make the metal-poor gas gradually disappear, starting from the regions closer to the disk. 9 Comparing the two GADGET-2 simulations with and without the metal diffusion scheme, Gad2-GasHalo+TFB2 and Gad2-GasHalo+TFB2+diff1 (left and middle panel), we find that the run with the scheme reduces the metal-poor gas more efficiently. As discussed in Section 4.1, the metal diffusion scheme in particle-based codes helps to redistribute the metals homogeneously in the halo and enriches a larger volume with metals.
Metal diffusion is a vital component in the process of galaxy formation. In Section 4, we have demonstrated that it has to be explicitly included in particle-based simulations to produce a realistic ISM and CGM in and around a simulated galaxy. Without considering the transport of metals via diffusion between gas particles, the MDF may become unreasonable (Figure 9), and pockets of unrealistically metal-poor gas may survive in the halo (Figure 10). Figure 11. 500 Myr simulation snapshots similar to Figure 3, but this time for the runs using the GIZMO-PSPH code and our proposed prescription that makes a particle-based simulation compatible with the mesh-based Enzo run -stellar feedback with boosted thermal energy and the metal diffusion coefficient of C d = 0.02. The GIZMO-PSPH simulation without the halo gas -the PSPH run (see Table 2); the same runtime condition as the particle-code runs in Kim et al. (2016) -substantially differs from the Enzo run in the metal distribution inside the halo. In contrast, another GIZMO-PSPH simulation, but this time with a gas halo, more feedback energy, and the explicit metal diffusion scheme -the PSPH-GasHalo+TFB1.8+diff1 run -is compatible with the Enzo run and the Gad2-GasHalo+TFB2+diff1 run in Figure 3. See Table 2 for the list of our simulations, and Section 5 for more information on this figure.

GENERALIZING OUR FINDINGS IN ANOTHER PARTICLE-BASED CODE
In Sections 3 and 4, using various metrics such as PDFs in density, temperature and metallicity, we have found that the Gad2-GasHalo+TFB2+diff1 run is the most compatible with the Enzo run in reproducing its metal properties. The Gad2-GasHalo+TFB2+diff1 run features a sufficientlyresolved gas halo in the IC, stellar feedback with boosted thermal energy (twice the value used in the Enzo run), and the metal diffusion with coefficient C d = 0.02. Before concluding our paper, in this section, we briefly test if our prescription for the GADGET-2 code is also applicable in the GIZMO-PSPH code (see Section 2.1.2), and if our findings in GADGET-2 can be generalized in other particle-based simulations. Readers should note that the authors never mean to imply that the Enzo run is a gold standard that all other simulations should match. We adopt the Enzo run only as a refer-ence while trying to find a setup that makes the mesh-based and particle-based codes behave in a similar fashion. Figure 11 is similar to Figure 3, but now the particle-based simulations are performed on the GIZMO-PSPH code. The PSPH run in Figure 11 (2nd column; see also Table 2) behaves very similarly as the Gad2 run in Figure 3, harboring an extremely metal-poor halo with only a few metalenriched particles of very high outflow velocity. Meanwhile, the PSPH-GasHalo+TFB1.8+diff1 run that utilizes our proposed prescription (3rd column; but with 10% less energy than Gad2-GasHalo+TFB2+diff1) shows a similar metal distribution and outflow velocity as the mesh-based Enzo run. The two-dimensional phase plot in Figure 12 verifies the same trend. The PSPH run in Figure 12 (2nd panel) is similar to the Gad2 run in Figure 7, lacking the hot diffuse medium around the disk, unlike Enzo. But with our proposed setup, the PSPH-GasHalo+TFB1.8+diff1 run (3rd panel) now features the hot metal-enriched medium just like the Figure 12. Two-dimensional probability distribution functions (PDFs) of metal mass on the density-temperature plane at 500 Myr. The figure is similar to Figure 7, but this time for simulations using GIZMO-PSPH and our proposed prescription that makes a particle-based simulation compatible with the mesh-based Enzo run. Only with a gas halo sufficiently resolved in the IC, more feedback energy, and the explicit metal diffusion (i.e., PSPH-GasHalo+TFB1.8+diff1) can the particle-based simulations identify the hot diffuse gas around the disk at [∼ 10 −29 g cm −3 , ∼ 10 6 K]. See Section 5 for more information on this figure.
Gad2-GasHalo+TFB2+diff1 run in Figure 7. As a note, we have found that in the GIZMO-PSPH simulation, slightly less (10%) thermal feedback energy is required than in GADGET-2, to best match the mesh-based Enzo run. This small difference could be attributed to an inherent inter-code discrepancy between the pressure-energy formulation of SPH in GIZMO-PSPH and the density-entropy formulation in GADGET-2.
Based on these experiments, we argue that our proposed setup help to alleviate the discrepancy between mesh-based and particle-based codes previously reported in, e.g., Kim et al. (2016). Because our prescription is straightforward and relies only on the fundamental properties of SPH (see Section 2.1.2 for more discussion), rather than a novel feature in any one code, we expect it to be widely applicable in many SPH codes. One may also argue that our criteria can be used to check if any particle-based simulation is robust and reproducible -especially by a mesh-based code. For example, in a cosmological zoom-in simulation using a particle-based code, one may check if a galactic halo is resolved with a sufficient number of gas particles before analyzing its metal content or performing simulated metal line observations. 6. DISCUSSION AND CONCLUSION Acquiring a realistic metal distribution in numericallyformed galaxies is vitally important, yet it is highly sensitive to the hydrodynamic schemes used and the diffusion model employed. Indeed, the AGORA code comparison project has previously reported a nontrivial discrepancy in the metal distribution of an idealized galaxy simulation between meshbased and particle-based codes (Kim et al. 2016). Following up on their observations, in this paper, we have investigated what causes the discrepancy and how it could be alleviated by changing the setup of a particle-based simulation. First, we have tested a modified IC for particle-based codes (Section 2.2.2) that contains a large number of gas particles in the galactic halo to match the initial gas distribution of a meshbased simulation. Then, we have examined the metal distributions in a suite of GADGET-2 simulations with different stellar feedback strengths and compare them with that of a ENZO simulation (see Section 3). We have also discussed the effect of an explicit metal diffusion scheme (Section 2.3.3) described in Hopkins et al. (2018) and Escala et al. (2018), and tested various coefficient values (Section 4).
We propose that, to alleviate the discrepancy in metal distributions between mesh-based and particle-based codes, the following three factors should be considered in a particlebased simulation: (1) Sufficiently-resolved gas halo: Our study finds that a gas halo with density n H = 10 −6 cm −3 can provide enough pressure to contain the galactic outflows within the virial radius. A sufficient number of gas particles is needed in the halo to describe a well-resolved medium into which the energy and metals of the SN-driven outflows could be transferred. Consequently, the existence of the gas haloor the lack thereof -heavily affects the metal distribution in it.
(2) Stellar feedback: Stellar feedback is the main source of energy that maintains the hot-diffuse medium around the galactic disk. We find that the amount of metal-enriched gas and the metallicity profiles in the halo are dictated by the strength of thermal stellar feedback. Particle-based codes require approximately twice the thermal feedback energy as the mesh-based ENZO code to produce compatible metal distributions in the halo. (3) Turbulent metal diffusion: We find that the explicit metal diffusion scheme based on turbulent mixing is essential to render a realistic low-metallicity gas in the galactic ISM and CGM. The shape of a metallicity PDF (or MDF) is highly sensitive to the strength of diffusion, both in the disk and in the halo. The diffusion coefficient C d = 0.02 in a particle-based simulation provides the best match to a mesh-based ENZO simulation. Our proposed prescription combining the three factors above has been tested with two particle-based codes, GADGET-2 (Sections 3 to 4) and GIZMO-PSPH (Section 5), and is generally applicable in many SPH codes.
Even though the experiments reported in this paper have been performed with an idealized, isolated galaxy, our study offers a useful reference point for cosmological (zoom-in) simulations as well. For example, one may check if a galactic halo is sufficiently resolved in a particle-based simulation to make sure that any metal-related properties in the halo are reproducible by a mesh-based code (Section 5). In the forthcoming paper, we will investigate the metal distribution inside the CGM in a full cosmological simulation with meshbased and particle-based codes. We aim to examine how the predicted metal lines in the CGM and the pockets of pristine gas change as we adopt different hydrodynamic schemes (e.g., AMR vs. SPH vs. SPH+diffusion scheme). In addition, we will study the possibility of producing an extended metalenriched CGM via a galaxy merger, inspired by the recent ob-servations of widely extended or confined [CII] lines in highz galaxies (e.g., Fujimoto et al. 2019;Ginolfi et al. 2020).