TY - JOUR
AU - Moa,, Belaid
AB - ABSTRACT Modelling the turbulent diffusion of thermal energy, momentum, and metals is required in all galaxy evolution simulations due to the ubiquity of turbulence in galactic environments. The most commonly employed diffusion model, the Smagorinsky model, is known to be overdiffusive due to its strong dependence on the fluid velocity shear. We present a method for dynamically calculating a more accurate, locally appropriate, turbulent diffusivity: the dynamic localized Smagorinsky model. We investigate a set of standard astrophysically relevant hydrodynamical tests, and demonstrate that the dynamic model curbs overdiffusion in non-turbulent shear flows and improves the density contrast in our driven turbulence experiments. In galactic discs, we find that the dynamic model maintains the stability of the disc by preventing excessive angular momentum transport, and increases the metal-mixing time-scale in the interstellar medium. In both our isolated Milky Way-like galaxies and cosmological simulations, we find that the interstellar and circumgalactic media are particularly sensitive to the treatment of turbulent diffusion. We also examined the global gas enrichment fractions in our cosmological simulations, to gauge the potential effect on the formation sites and population statistics of Population III stars and supermassive black holes, since they are theorized to be sensitive to the metallicity of the gas out of which they form. The dynamic model is, however, not for galaxy evolution studies only. It can be applied to all astrophysical hydrodynamics simulations, including those modelling stellar interiors, planetary formation, and star formation. diffusion, hydrodynamics, turbulence, methods: numerical, galaxies: intergalactic medium, galaxies: ISM 1 INTRODUCTION Galaxies form at the confluence of gas streams and cooling flows at the centres of virialized haloes and evolve via a constant exchange of baryons with their environments. Despite significant recent progress on understanding the details of this picture, developing predictive models for the evolution and observed properties of galaxies has proven to be an immense challenge (see, for example, Guedes et al. 2011; Hopkins et al. 2014; Vogelsberger et al. 2014; Schaye et al. 2015; Genel 2016; Davé et al. 2017; – see also Somerville & Davé 2015 and Naab & Ostriker 2017 for a recent review and additional references). The problem lies in the large number of complex interconnected processes involved, and the huge dynamic range in spatio-temporal scales over which they operate. One important interplay involves the interstellar medium (ISM), the gas that permeates a galaxy and provides fuel for star formation, and circumgalactic medium (CGM), the gas that cocoons the galaxy. The amount of gas in the CGM and the efficiency with which it can cool, fall into the galaxy and replenish the ISM, are important variables in setting the duration and the rate of star formation (Somerville & Davé 2015). Stellar winds and supernova explosions (SNe) – processes directly related to star formation – provide competition for gas cooling (Springel & Hernquist 2003b; Oppenheimer & Davé 2008). These deposit energy and momentum into the ISM, engendering outflows of gas. If these outflows are sufficiently powerful, they not only heat the CGM, but also expel most of the CGM from the galaxy’s halo. This ongoing competition between gas into and out of galaxies provides a basic framework for understanding a number of observed properties (Shen et al. 2012; Crain et al. 2013; Christensen et al. 2016; Oppenheimer et al. 2016; Sokołowska et al. 2016). The cooling efficiency and ionization state of the gas in the CGM depend sensitively not only on the spatial injection and redistribution of thermal energy and momentum (Suresh et al. 2017), but also on the metals that are transported from the galaxy via the galactic outflows (Davé, Finlator & Oppenheimer 2006; Oppenheimer & Davé 2006; Finlator & Davé 2008; Hani et al. 2018). Metals, although negligible in terms of mass fraction, play an out-sized role in galaxy evolution because they can dramatically alter the CGM’s cooling profile (van de Voort et al. 2012; Oppenheimer & Schaye 2013; Sokołowska et al. 2017) and, hence, the delicate balance between gas in- and outflow. In addition to cooling, the redistribution of metals can impact other processes, such as the sites and formation history of the putative Population III (Pop III) stars and supermassive black holes (SMBHs). Pop III stars are associated with star formation involving near-pristine gas with an upper limit on metallicities somewhere in the range [Z] ∼ −6 to [Z] ∼ −3 (Sarmento, Scannapieco & Pan 2016). As for SMBHs, a number of authors (e.g. Volonteri 2010 and for use in recent simulations Tremmel et al. 2017) postulate that they form via direct collapse of gas clouds with metallicity [Z] ≲ −4.0, with the idea that very low metallicity would prevent the gas from cooling rapidly and fragmenting into Pop III stars during collapse. It is therefore crucial to identify which physical processes redistribute thermal energy, momentum, and metals in galactic environments spatially, and to include these accurately in the numerical galaxy evolution experiments. One critical, often overlooked, redistribution mechanism is gas turbulence. Turbulence occurs when inertial forces dominate viscous forces in a gaseous environment, and kinetic energy injected on large scales cannot immediately dissipate as heat. This leads to the formation of a kinetic energy cascade, as coherent turbulent eddies on large scales spawn eddies on successively smaller scales, until the energy thermalizes. Galactic environments, for example, are expected to be highly turbulent (Evoli & Ferrara 2011; Iapichino, Viel & Borgani 2013). In the case of the CGM, this is strongly suggested by the kinematic complexity revealed by absorption and emission line measurements (Tumlinson, Peeples & Werk 2017) and it has long been recognized that the cold ISM is also highly turbulent. The susceptibility of a medium to turbulence is quantified by its dimensionless Reynolds number,1Re. The Reynolds number for the cold ISM has been estimated to be as high as Re ∼ 107 (Elmegreen & Scalo 2004), whereas the onset of turbulence usually occurs at Re ∼ 103. The Reynolds number is also a measure of the separation of scales in the energy cascade and, in the case of incompressible turbulence, L/η ∼ Re3/4, where L is the kinetic energy injection scale, and η is the dissipation scale.2 Therefore the degrees of freedom in a 3D simulation scales as Re9/4, and simulating a Re ∼ 107 flow would require 1015 fluid elements! Contemporary cosmological simulations have reached ∼1012 fluid elements and a dynamic range typically of the order ∼106 (Somerville & Davé 2015), with the smallest scale being the resolution limit h. Therefore all cosmological simulations that involve turbulence – independent of hydrodynamical method – have a natural cut-off scale h, in the range η ≪ h ≪ L, where discretization truncates the turbulent cascade. Physically, small-scale turbulent fluctuations, by promoting mixing, provide a transport mechanism for the fluid properties such as momentum, thermal energy, and metals. In numerical simulations, this implies that turbulence on scales smaller than h can potentially impact the resolved properties of the flow and consequently must be accounted for (Germano et al. 1991). The crux of the issue is that, as stated above, the kinetic energy flux down the turbulent cascade is truncated at h in the simulations whereas, in reality, the kinetic energy cascade should continue to smaller scales until it is dissipated. Numerical simulations break the physical coupling between the scales that are resolved (>h) and unresolved (<h), and, therefore, they require models of (i) the kinetic energy flux from the resolved to the unresolved scales, (ii) the effect of unresolved eddies on the resolved scales, and (iii) the transport properties of kinetic energy on unresolved scales. Typically, this is done using the turbulent eddy-diffusion models that treat sub-grid turbulent eddy motion as a diffusive process. Several implementations of this approach have been proposed (cf. Schmidt & Federrath 2011; Di Mascio et al. 2017; Wadsley, Keller & Quinn 2017; – see Sagaut 2006 and Garnier, Adams & Sagaut 2009 for extensive lists); none the less, many cosmological hydrodynamical simulation studies assume that numerical diffusion adequately accounts for sub-grid turbulent transport (Schmidt 2015). Numerical diffusion, however, can lead to diffusive behaviour that poorly represents turbulent flow statistics (Sagaut 2006). Extensive effort has been devoted towards developing models for treating turbulent diffusion in Eulerian cosmological simulations, which employ grids to discretize the fluid. Here we refer the reader to Scannapieco & Brüggen (2008), Pan, Scannapieco & Scalo (2013), Federrath (2013), Schmidt et al. (2014), Schmidt (2015), Semenov, Kravtsov & Gnedin (2016), and Sarmento et al. (2016). In this paper we focus on the Lagrangian hydrodynamical approach. Within the Lagrangian framework, the fluid equations of motion are approximated by tracking individual fluid elements as they move with the flow. Commonly used Lagrangian methods in computational cosmology include smoothed particle hydrodynamics (SPH) (Gingold & Monaghan 1977; Lucy 1977; Hernquist & Katz 1989), moving-meshes (MM) (Springel 2010), and higher order mesh-free methods (MF) (Lanson & Vila 2008a,b; Gaburov & Nitadori 2011; Hopkins 2015) – see Springel (2010) and Hopkins (2015) for extensive discussion. SPH has no inherent diffusion and therefore, by construction, explicitly requires additional transport terms, one benefit of which is full control over the strength of mixing. In MM and MF methods, numerical diffusion arises as a by-product of the numerical scheme used to solve the Riemann problem between fluid elements. However, as noted above, inherent numerical diffusion does not, in general, reproduce the correct turbulent properties of a simulated flow. The key problem here is: how can we address this? Or, more precisely, how do we model the interaction between the unresolved and the resolved scales? A common approach is to assume that the interaction between the resolved and the unresolved scales reduces to a local transfer of kinetic energy from the large to the small scales. If we additionally assume that the kinetic energy transfer mechanism is analogous to a diffusive process, where kinetic energy is redistributed on progressively smaller scales via momentum diffusion, then we can treat the action of unresolved eddies in a similar fashion to molecular viscosity. Under this equivalency, on scales close to the resolution scale h, the unresolved eddies extract kinetic energy from the resolved flow via momentum diffusion, and simultaneously allow for the dissipation of kinetic energy. This is the eddy-viscosity hypothesis, and is fundamentally a model assumption that describes the simplest view of the interactions between the scales. The action of the unresolved eddies can then be modelled by a viscous term in the fluid equations of motion – diffusing momentum and dissipating kinetic energy (Pope 2000). The above assumptions lead to a simple model but it does not provide guidance about the appropriate choice for the effective diffusivity/viscosity. In a turbulent cascade, the diffusive action of the eddies depends on their velocity and length-scale. On the sub-grid level, this involves estimating the velocity scale that transports fluid properties over the resolution scale h (Wadsley, Veeravalli & Couchman 2008; Greif et al. 2009). In the simplest model, the Smagorinsky model (Smagorinsky 1963), this velocity is assumed to be proportional to the gradients of the velocity field. In this model, the sub-grid eddies both diffuse (characterized by diffusivity D) and dissipate via eddy viscosity, νsgs, given by νsgs = D = (Csh)2|S*|, where |S*| is the norm of the trace-free shear tensor, and Cs is the Smagorinsky model constant. Contemporary fluid mechanics literature notes that in general Cs needs to be tuned to a value between 0.1 and 0.2 for optimal results under different flow conditions (Garnier et al. 2009). The main advantage of the Smagorinsky model is its simplicity and as a result, many researchers have started to incorporate this model into their cosmological simulation codes – specifically as a model for the diffusivity when treating thermal energy and metal mixing (Shen, Wadsley & Stinson 2010; Brook et al. 2014; Williamson, Martel & Kawata 2016; Sokołowska et al. 2017; Tremmel et al. 2017; Escala et al. 2018). The model, however, has some drawbacks: (i) A single valued Cs is incapable of correctly describing different types of turbulent flows. Studies show that the Smagorinsky model introduces too much diffusion into the flow in almost all cases except for homogeneous, isotropic turbulence (Garnier et al. 2009). And, (ii) the sub-grid eddy viscosity does not vanish for laminar shear flows where there ought to be no diffusion due to turbulence. Overdiffusion is especially worrisome given the push to resolve the multiphase structure in the ISM and CGM at greater levels, and the recent results that differential, localized metal mixing could change our understanding of galactic chemical evolution (Emerick et al. 2018). As noted above, most implementations of the Smagorinsky model only consider thermal energy and metal diffusion, but not momentum diffusion because the latter drawback above is a concern for differentially rotating structures, such as galactic discs. Specifically, it results in undesired angular momentum transport and, consequently, unphysical flows in the discs. Most galaxy evolution studies ignore momentum diffusion in order to avoid this viscous instability. The crux of the problem is that a constant Cs cannot automatically adjust to either the local conditions or the changing character of the flow with time. This led Germano et al. (1991) to propose the dynamic Smagorinsky model where Cs is a function of space and time, i.e. Cs = Cs|$(\boldsymbol{x}, t)$|⁠. The dynamic Smagorinsky model has, to our knowledge, not yet been implemented and investigated in cosmological simulations. The model has, however, been validated extensively in the fluid simulation community by comparing to the results of standard numerical tests and experiment data, and has been shown to improve upon the constant-coefficient model (Kirkpatrick et al. 2006; Kleissl et al. 2006; Benhamadouche, Arenas & Malouf 2017; Lee & Cant 2017; Kara & Çağlar 2018; Taghinia, Rahman & Lu 2018). In the dynamic model, the sub-grid properties of a turbulent fluid are computed under two assumptions3: (i) the behaviour of the largest unresolved eddies is entirely determined by their interactions with the eddies on the smallest resolved scales, and (ii) these interactions are analogous to those between the fluid motions on the smallest resolved scale and the motions on larger scales. In practice, determining the characteristics of these interactions involves filtering (or smoothing) the resolved velocity field on two different scales. When the sub-grid turbulent properties are calculated based on the local fluid properties, Cs|$(\boldsymbol{x}, t)$| reduces to zero (i.e. the eddy viscosity/diffusivity vanishes) in non-turbulent (or laminar) shear flows (Piomelli & Liu 1995). Consequently, this allows for the self-consistent treatment of momentum diffusion, along with thermal energy and metal diffusion, in numerical studies of cosmic baryons and galaxy evolution. In this study, we introduce an implementation of the dynamic Smagorinsky model for the first time in Lagrangian astrophysical simulations, focusing on the higher order MF approach. We discuss the implementation details in Section 2. In order to test the impact of dynamic localized turbulent mixing on the galactic ecosystem we run a series of hydrodynamical and physical experiments relevant to galaxy evolution. In Section 3 we show the results for a set of standard hydrodynamic tests and explicitly check on the extent of sub-grid diffusion in laminar shear flows. We then go on to examine the effects of dynamic mixing in an isolated disc galaxy in Section 4, followed by a set of cosmological simulations in Section 5. We investigated SPH and performed all of the experiments presented below, but we do not include them in this paper because the results are qualitatively similar to the MF results, as we mention in Section 6. 2 METHODS 2.1 Hydrodynamics In order to test the impact and robustness of the localized turbulent diffusion model, we employ a modified version of the gizmo gravity plus hydrodynamics solver code (Hopkins 2015). gizmo builds on the gadget-3 code base (Springel 2005), with improvements in numerical accuracy and includes an implementation of the novel mesh-free finite mass (MFM) method, in addition to various implementations of SPH methods. The MFM method evolves the fluid equations of motion in a Lagrangian manner similar to SPH. However, while the fluid mass elements in SPH are discretized into particles and their motions are determined by fluid properties smoothed over neighbouring particles, the conservation laws in MFM are evolved by calculating the flux of basic variables between neighbouring particles.4 These fluxes depend on the effective face area between the two particles and are determined by solving the Riemann problem along the line connecting them. This removes the need for additional terms, such as artificial viscosity and conductivity as is necessary in SPH to ensure proper treatment of shocks. The numerical Riemann solvers have inherent numerical dissipation,5 hence improved shock capturing capabilities. We, therefore, focus on the use of the MFM method in our investigation of the dynamic diffusion model. For a thorough exposition of gizmo and an extensive comparison between MFM and SPH, see Hopkins (2015). 2.2 Sub-grid turbulent diffusion terms gizmo solves the conservation equations for momentum, energy, and mass using the MFM method, and like other hydrodynamic methods, it is limited in resolution down to a scale h. The minimum resolution limits the ability to resolve high Reynolds number flows down to the viscous dissipation scale, impacting its ability to resolve the turbulent cascade. As mentioned in Section 1, the interaction of the resolved and unresolved scales must be modelled. These models apply to the additional sub-grid scale terms that appear in the equations of motion when treating discretization as a filtering process (Sagaut 2006; Garnier et al. 2009; Schmidt 2015). In this section, we detail the origin of the sub-grid scale terms and discuss which terms we include in the gizmo code. Discretization of the conservation equations is equivalent to applying a low-pass filter, damping out high-frequency turbulent fluctuations. When we discuss filtering, we refer to the definition of a general filtered scalar field |$f(\boldsymbol{x})$|⁠, \begin{eqnarray*} \overline{f}(\boldsymbol{x}) \equiv \int _D f(\boldsymbol{x}^{\prime })G(|\boldsymbol{x}^{\prime } - \boldsymbol{x}|, \overline{h})\, \mathrm{ d}\boldsymbol{x}^{\prime }, \end{eqnarray*} (1) where |$G(|\boldsymbol{x}^{\prime }-\boldsymbol{x}|, \overline{h})$| is a filter function. |$\overline{h}$| is the characteristic size of the filtering operation below which fluctuations are damped (essentially the resolution scale, in the present context), |$\mathrm{ d}\boldsymbol{x}^{\prime }$| is a volume element, and the integral is evaluated over the entire domain. We discuss our filtering implementation in more detail, in Section 2.4. We now apply this equation to the conservation equations in order to see that additional terms in the hydrodynamical equations are required. The momentum conservation equation for a compressible fluid follows6 (Landau & Lifshitz 1987), \begin{eqnarray*} \frac{\partial }{\partial t}(\rho u_i) + \frac{\partial }{\partial x_j}[\rho u_i u_j + p\delta _{ij}] = 0, \end{eqnarray*} (2) where |$u_i = u_i(\boldsymbol{x},t)$| is the fluid velocity vector in the i = {x, y, z} direction, |$\rho = \rho (\boldsymbol{x}, t)$| is the fluid density, and |$p = p(\boldsymbol{x}, t)$| is the pressure. When we filter the momentum equation, assuming the filtering operation and derivatives commute, we end up with an extra term τij, \begin{eqnarray*} \frac{\partial }{\partial t}(\overline{\rho } \tilde{u_i}) + \frac{\partial }{\partial x_j}[\tau _{ij} + \overline{p}\delta _{ij} + \overline{\rho }\tilde{u}_i\tilde{u}_j] = 0, \end{eqnarray*} (3) where we have also switched to density weighted variables such that |$\tilde{u}_i = \overline{\rho u_i} / \overline{\rho }$|⁠. τij is the subgrid-scale turbulent stress tensor or residual stress tensor and is defined as, \begin{eqnarray*} \tau _{ij} \equiv \overline{\rho } (\widetilde{u_i u_j} - \tilde{u}_i\tilde{u}_j). \end{eqnarray*} (4) This term must be modelled because |$\widetilde{u_i u_j}$| is unknown at the time of simulation, i.e. the system of equations is not closed. A common model, or closure, involves the eddy-viscosity assumption where the sub-grid scales impart a momentum flux on the resolved scales that is linearly dependent on the rate of strain of the resolved scale, \begin{eqnarray*} \tau _{ij} = -2 \overline{\rho } \nu _\mathrm{sgs} \widetilde{S^*}_{ij}, \end{eqnarray*} (5) where νsgs is the sub-grid eddy viscosity, and |$\widetilde{S^*}_{ij}$| is the trace-free resolved rate of strain tensor. The same logic can be applied to any of the conservation equations and any filtered multiplicative terms require modelling. In the derivations below, we model (as in equation 5 above) the unknown terms under the assumption that they act as diffusive processes. There is a similar term when filtering the total energy equation (e is the specific total energy), \begin{eqnarray*} \frac{\partial }{\partial t}(\rho e) + \frac{\partial }{\partial x_j}[\rho u_j e + u_j p] = 0. \end{eqnarray*} (6) leads to additional terms, \begin{eqnarray*} \frac{\partial }{\partial t}(\overline{\rho } \tilde{e}) + \frac{\partial }{\partial x_j}[Q_j + P_j + \overline{\rho } \tilde{u}_j \tilde{e} + \overline{u}_j \overline{p}] = 0, \end{eqnarray*} (7) where |$e = \theta + \frac{1}{2} |\boldsymbol{u}|^2$| (θ is the specific internal energy). We also have defined: \begin{eqnarray*} Q_j \equiv \overline{\rho } (\widetilde{u_j e} - \tilde{u}_j \tilde{e}),\quad \quad P_j \equiv \overline{u_j p} - \overline{u}_j \overline{p}. \end{eqnarray*} (8) In this study, we ignore the term associated with pressure, Pj, and the term in Qj that arises from |$\frac{1}{2}|\boldsymbol{u}|^2$| in e, and focus on the application of the dynamic model to the sub-grid momentum term working in concert with the widely employed thermal energy term (Shen et al. 2013; Brook et al. 2014; Tremmel et al. 2017; Wadsley et al. 2017), \begin{eqnarray*} q_j \equiv \overline{\rho } (\widetilde{u_j \theta } - \tilde{u}_j \tilde{\theta }) = -\overline{\rho } \nu _\mathrm{sgs} \frac{\partial \tilde{\theta }}{\partial x_j}. \end{eqnarray*} (9) Not only are there additional terms for momentum and thermal energy, but any scalar quantities, such as the concentration of different metal species, in the gas should be transported in a turbulent flow. In order to model this, we treat metal concentrations, ϕz, as passive scalars that obey a diffusion equation (Pope 2000; Shen et al. 2010), \begin{eqnarray*} \frac{\partial \phi _z}{\partial t} = \frac{\partial }{\partial x_j} \bigg (\overline{\rho } \nu _\mathrm{sgs} \frac{\partial \phi _z}{\partial x_j} \bigg). \end{eqnarray*} (10) For a detailed description of the incorporation of these fluxes into gizmo, see Hopkins (2017). Throughout this paper, we refer to the action of the terms in equations (5), (9), and (10) as turbulent diffusion because they contribute to the conservation equations as ∇2f, where f is the flux quantity. Additionally, when we mention energy diffusion, we are referring to the term |$\boldsymbol{\nabla }\cdot \boldsymbol{q}$| (thermal energy diffusion), and, similarly, when we mention momentum diffusion we are referring to the action of the stress tensor through |$\boldsymbol{\nabla }\cdot \boldsymbol{\tau }$|⁠, along with the corresponding kinetic energy dissipation. 2.3 Diffusivity Physically, turbulent mixing can be modelled as a diffusive process with diffusivity D and, in the simplest model, the fluid properties are assumed to mix over the resolution scale h with a velocity h|S*|, where |S*| is the norm of the trace-free shear tensor. This is, as we mentioned in Section 1 the Smagorinsky model and the corresponding diffusivity is parametrized as D = νsgs = (Csh)2|S*|. The Smagorinsky model inherently assumes that the kinetic energy transfer rate down the turbulent cascade is equal on all scales, and is equal to the physical dissipation rate (i.e. the flow is in local equilibrium). In simulations, the resolution scale h inhibits kinetic energy from moving to progressively smaller scales, and results in a build-up of kinetic energy at the resolution scale – so long as numerical dissipation cannot extract kinetic energy sufficiently rapidly. The Smagorinsky model combined with the local equilibrium assumption only provides a model for the turbulent stress, τij, and dissipation in the flow, |$\Sigma = \tau _{ij} \tilde{S}^*_{ij}$|⁠, and ignores the additional terms we discuss in Section 2.2. In order to consistently model all of the energetic terms, one must relax the local equilibrium assumption and follow the sub-grid kinetic energy, K, directly. It is possible to derive a one-equation model for K that includes a third-order term for the transport of K on sub-grid scales (Schmidt 2015), and self-consistently follows all of the sub-grid quantities. Each sub-grid term can then be modelled using diffusive terms similar to the models in Section 2.2, in order to close the system of equations. If the Smagorinsky model only considers the turbulent stress, is it then valid to apply this model (as we have done following Shen et al. 2010) to the thermal energy and metal sub-grid terms? In order for this to be possible, the local equilibrium condition must be approximately true in the regime of interest. We are specifically interested in cosmological-scale gas, and Schmidt et al. (2016) show that the local equilibrium condition holds – on average – in a cosmological-scale volume. Introducing the dynamic model on top of these approximations further supports our model assumption, because the dynamic model inherently accounts for the deviations from local equilibrium. The Smagorinsky model diffusivity is parametrized in gizmo for a particle ‘a’ as, \begin{eqnarray*} D_{\mathrm{ a}} = \rho _{\mathrm{ a}} (C_\mathrm{s} h_{\mathrm{ a}})^2 |S^*|_{\mathrm{ a}}, \end{eqnarray*} (11) where Cs is the Smagorinsky constant, ha is the mean inter-particle spacing in the kernel, and |S*|a is the magnitude of the trace-free symmetric shear tensor. Note that we absorb the densities |$\overline{\rho }$| from Section 2.2 into νsgs via ρa. Da is used in the diffusion equations for thermal energy, momentum, and metal mass fractions as described in Hopkins (2017). There are a myriad of values quoted for Cs in the literature (see section 5 of Sagaut 2006 for an extensive list), but we choose the value calibrated for fully eveloped isotropic turbulence, Cs = 0.2 (Clark, Ferziger & Reynolds 1979), because the Smagorinsky model was developed for this specific regime. In order to compare our Cs with other values in the literature, it is important to consider the definitions of the quantities in equation (11). In many SPH studies, the length-scale in the diffusion coefficient definition is taken as the kernel-support radius, hSPH, which is the maximum extent from a particle that gives a non-zero weight. In contrast, as we mentioned above, we employ the mean inter-particle spacing within the kernel or ha ≈ 0.5hSPH. Additionally, some studies use |$C \equiv \sqrt{2}$|Cs2 whereas for this study |$\sqrt{2}$| is absorbed into our definition of the norm of the shear tensor (see equation 13 below), closely following the fluid simulation literature (Piomelli & Liu 1995). Using these definitions, our adopted value of Cs is lower than those quoted in the astrophysics literature. For example, the value C = 0.05 (see Shen et al. 2010, 2013; Brook et al. 2014), corresponds to Cs = 0.37 whereas C = 0.03 corresponds to Cs =0.29 (see Wadsley et al. 2017). However, our Cs is higher than the value recently calibrating from studying metal mixing in dwarf galaxies, where Escala et al. (2018) found Cs = 0.046 reproduced more realistic stellar metal distribution functions (MDFs) via supersonic mixing in the ISM. We compute the trace-free symmetric shear tensor via the high-order accurate gradient estimators in gizmo, \begin{eqnarray*} S_{ij}^* = \frac{1}{2}\bigg (\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\bigg) - \frac{1}{3}\delta _{ij}\frac{\partial u_k}{\partial x_k}, \end{eqnarray*} (12) where ui is the fluid velocity vector and xi the spatial coordinate, and i, j = {x, y, z}. The magnitude of equation (12) is implemented using the Frobenius norm (Piomelli & Liu 1995), \begin{eqnarray*} |S^*| = \sqrt{2 S_{ij}^* S_{ij}^*}. \end{eqnarray*} (13) We note that in gizmo, Escala et al. (2018) impose an ad hoc cap on the diffusivity based on the expected maximum mass flux between resolution elements. The cap does not, however, mitigate the fact that the constant Smagorinsky model induces diffusion whenever there is shear, regardless of whether the fluid is laminar or turbulent. It is only intended to prevent unphysical mass/energy transport that can potentially arise due to the excessive diffusivity of the Smagorinsky model. We adopt the same limiter in this study, but note that the diffusivity rarely reaches the maximum limit. 2.4 Dynamic model The Smagorinsky model provides an approximate model of sub-grid mixing for fully developed, homogeneous turbulence but it is far too simple for complex flows. In fact, in laminar shear flows, the constant-coefficient Smagorinsky model predicts a non-zero diffusivity through its dependence on the shear strength (equation 11). However, in this situation, the diffusivity should be zero since the fluid is not turbulent. In more complex flows, such as those in astrophysical contexts, the value of the constant ought to depend on the spatio-temporal coordinates |$C_\mathrm{ s} = C_\mathrm{ s}(\boldsymbol{x}, t)$| (Germano et al. 1991). Piomelli & Liu (1995) showed that by assuming scale-similarity (cf. Section 1), the local Smagorinsky constant in a neighbourhood can be calculated at each point (at a fixed simulation time-step) as follows: \begin{eqnarray*} C_{\mathrm{dyn}}(\boldsymbol{x}) = C_\mathrm{ s}^2 = -\dfrac{1}{2} \dfrac{(L_{ij} - 2 \widehat{C_{\mathrm{dyn}}^{p}\beta _{ij}})\alpha _{ij}}{\alpha _{mn}\alpha _{mn}}. \end{eqnarray*} (14) Here Lij is the Leonard tensor, \begin{eqnarray*} L_{ij} = \widehat{\overline{u}_i \overline{u}_j} - \widehat{\overline{u}}_i \widehat{\overline{u}}_j, \end{eqnarray*} (15)|$C_{\mathrm{dyn}}^{p}$| is the value of Cdyn at the previous time-step, and αij and βij are defined as, \begin{eqnarray*} \begin{split} \alpha _{ij} \equiv \widehat{h}^2 |\widehat{\overline{S^*}}| \widehat{\overline{S^*}}_{ij},\\ \beta _{ij} \equiv \overline{h}^2 |\overline{S^*}| \overline{S^*}_{ij}. \end{split} \end{eqnarray*} (16) Here |$\overline{f}$| represents a filtering (or smoothing) operation on f over a length-scale |$\overline{h}$|⁠, and |$\widehat{f}$| represents smoothing on a scale |$\widehat{h}$|⁠. |$\overline{h}$| usually is equated with the lowest resolvable scale, and extensive work has been done in the fluid simulation community to show that the optimal value for |$\widehat{h}$| is |$\widehat{h} = 2 \overline{h}$| (Germano et al. 1991; Piomelli & Liu 1995; Spyropoulos & Blaisdell 1996; Schmidt, Niemeyer & Hillebrandt 2006; Grete et al. 2017). For this work, we choose |$\overline{h}$| to be the compact support radius of the kernel, and |$\widehat{h} = 2 \overline{h}$|⁠. More precisely, we smooth a scalar field |$f(\boldsymbol{x})$| by convolving it with the filter function |$G(|\boldsymbol{x}^{\prime } - \boldsymbol{x}|, \overline{h})$| over the domain,7 \begin{eqnarray*} \overline{f}(\boldsymbol{x}) = \int _D f(\boldsymbol{x}^{\prime })G(|\boldsymbol{x}^{\prime } - \boldsymbol{x}|, \overline{h})\,\mathrm{ d} \boldsymbol{x}^{\prime }. \end{eqnarray*} (17) This is similar to the SPH method of interpolating a scalar function, with |$G(|\boldsymbol{x}^{\prime } - \boldsymbol{x}|, \overline{h})$| sharing the same properties as the smoothing kernel |$W(|\boldsymbol{x}^{\prime } - \boldsymbol{x}|, \overline{h})$|⁠, \begin{eqnarray*} \begin{split} \int G(|\boldsymbol{x}^{\prime } - \boldsymbol{x}|, \overline{h})\,\mathrm{ d}\boldsymbol{x}^{\prime } = 1,\\ \lim _{\overline{h}\rightarrow 0} G(|\boldsymbol{x}^{\prime } - \boldsymbol{x}|, \overline{h}) = \delta (|\boldsymbol{x}^{\prime } - \boldsymbol{x}|). \end{split} \end{eqnarray*} (18) The MFM method employs a similar technique for evaluating integrals and in order to be consistent, we choose G = W. The integral in equation (17) is expensive, but we simplify the computation using XSPH smoothing (Monaghan 1989, 2005, 2011), \begin{eqnarray*} \overline{f}(\boldsymbol{x}) = f(\boldsymbol{x}) + \epsilon \int _D (f(\boldsymbol{x}^{\prime }) - f(\boldsymbol{x}))W(|\boldsymbol{x}^{\prime } - \boldsymbol{x}|, \overline{h})\,\mathrm{ d}\boldsymbol{x}^{\prime }. \end{eqnarray*} (19) Following Monaghan (2011), the Fourier coefficients ak of the velocity satisfy |$\overline{a}_\mathrm{k} = a_\mathrm{k}[1 + \epsilon (\widetilde{G}(k)-1)]$| where |$\overline{a}_\mathrm{k}$| are the coefficients of the smoothed field, |$\widetilde{G}(k)$| is the Fourier transform of the filter function, and k is the spatial frequency k = 2π/x. In the limit k → ∞, the coefficients satisfy |$\overline{a}_\mathrm{k} \rightarrow (1 - \epsilon) a_\mathrm{k}$|⁠. The value of ε controls the magnitude of the smoothing on a scale of |$\le \overline{h}$|⁠, and is constrained to 0 < ε ≤ 1. We choose ε = 0.8 to be consistent with the tests in Monaghan (2011). We discretize equation (19) as, \begin{eqnarray*} \overline{f}_a = f_a + \epsilon \sum _b \frac{m_b}{\langle \rho _{ab}\rangle _{\overline{h}}} (f_b - f_a) W(|\boldsymbol{x_a} - \boldsymbol{x_b}|, \overline{h}_{ab}), \end{eqnarray*} (20) where fa represents the quantity at particle a, |$\overline{h}_{ab}$| is the arithmetic mean of |$\overline{h}_a$| and |$\overline{h}_b$|⁠,8|$\langle \rho _{ab}\rangle _{\overline{h}}$| is the harmonic mean of the densities |$\overline{\rho }_a$| and |$\overline{\rho }_b$|⁠,9 and the sum is taken over b nearest neighbours. We also require doubly filtered quantities, which involves another application of equation (20) to the singly filtered quantities, \begin{eqnarray*} \widehat{\overline{f}}_a = \overline{f}_a + \epsilon \sum _b \frac{m_b}{\langle \rho _{ab}\rangle _{\widehat{h}}} (\overline{f}_b - \overline{f}_a) W(|\boldsymbol{x_a} - \boldsymbol{x_b}|, \widehat{h}_{ab}). \end{eqnarray*} (21) In order to calculate the average densities, or weights, in equations (20) and (21), we require the density at each particle for a given scale, \begin{eqnarray*} \begin{split} \overline{\rho }_a = \sum _b m_b W(|\boldsymbol{x}_a - \boldsymbol{x}_b|, \overline{h}_{a}), \\ \widehat{\rho }_a = \sum _b m_b W(|\boldsymbol{x}_a - \boldsymbol{x}_b|, \widehat{h}_{a}). \end{split} \end{eqnarray*} (22) It is important to note that the values of |$\overline{S^*}_{ij}$| and |$\widehat{\overline{S^*}}_{ij}$| are built from the smoothed velocity field and are not smoothed versions of the trace-free symmetric shear tensor in equation (12) (Schmidt 2015). They have the following corresponding equations: \begin{eqnarray*} \begin{split} \overline{S^*}_{ij} = \frac{1}{2}\bigg (\frac{\partial \overline{u}_i}{\partial x_j} + \frac{\partial \overline{u}_j}{\partial x_i}\bigg) - \frac{1}{3}\delta _{ij}\frac{\partial \overline{u}_k}{\partial x_k}, \\ \widehat{\overline{S^*}}_{ij} = \frac{1}{2}\bigg (\frac{\partial \widehat{\overline{u}}_i}{\partial x_j} + \frac{\partial \widehat{\overline{u}}_j}{\partial x_i}\bigg) - \frac{1}{3}\delta _{ij}\frac{\partial \widehat{\overline{u}}_k}{\partial x_k}. \end{split} \end{eqnarray*} (23) The magnitudes |$|\overline{S^*}|$| and |$|\widehat{\overline{S^*}}|$| are given via equation (13). The dynamic method relies on local scale-similarity in the neighbourhood of a point |$\boldsymbol{x}$|⁠, which in turn implies that the Smagorinsky model is an accurate description of the flow, albeit with a variable constant. This assumption breaks down in highly complex flows and in some cases, the dynamic model predicts negative values (Piomelli & Liu 1995; Urzay et al. 2013). Negative values of Cdyn are usually explained as backscatter in the cascade (Piomelli et al. 1991; Piomelli & Liu 1995; Meneveau & Katz 2000; Vreman 2004; Urzay et al. 2013) where in some circumstances a fraction of the energy cascading to small scales can return to large scales as smaller eddies unite to form larger eddies. An alternate explanation for Cs <0 is that the Smagorinsky model fails, and a more appropriate model should be employed. For now, we adopt the latter view and follow the usual approach in restricting Cdyn to the range Cdyn ∈ [0, Cs] with values Cdyn < 0 set to Cdyn = 0 (Garnier et al. 2009; Schmidt 2015). The upper limit is imposed since large values are thought to be due to numerical instability. For the remainder of this paper, we identify and discuss Cs|$=\sqrt{C_\mathrm{dyn}}$|⁠. The distinction between the simple and dynamic model is distinguishable based on context. 3 HYDRODYNAMICAL TESTS As indicated previously, a number of studies have carried out robust validation of the dynamic Smagorinsky model against, for example, experimental results within the fluid mechanics community (Kleissl et al. 2006; Benhamadouche et al. 2017; Lee & Cant 2017; Kara & Çağlar 2018; Taghinia et al. 2018) and the model has also been adopted by other users, including researchers studying atmospheric phenomena (e.g. Kirkpatrick et al. 2006). To motivate its use in cosmological and astrophysical simulations, we start by discussing the model within the context of three hydrodynamical tests. First, we investigate the distribution of predicted Cs values in homogeneous driven turbulence, and the sensitivity of the distributions to variations in the smoothing parameter, ε. Next, we examine the evolution of a Keplerian disc where turbulence is not expected to develop a priori yet numerical instabilities and noise lead to disorder in the velocity fields, and subsequent artificially enhanced diffusivities via the trace-free shear strength, |S*|. Last, we consider the linear regime of the Kelvin–Helmholtz instability which suffers from similar challenges as the Keplerian disc. In both cases, numerical instability causes the constant-coefficient Smagorinsky model to fail, and we investigate whether the dynamic model can mitigate spurious sub-grid turbulent mixing. 3.1 Homogeneous turbulence We investigate homogeneous, isotropic, driven turbulence to determine to what degree the dynamic model predicts different diffusivities on a per particle basis. We initialize periodic boxes of side length L = 1 with 643, 1283, and 2563 equal mass particles of an ideal isothermal gas, initial density ρ = 1, and energy per unit mass u = 1000, placed on uniform Cartesian grids.10 Following the methods of Bauer & Springel (2012) and Hopkins (2015), as the system is driven, the thermal energy of the gas particles are reset to the initial value in order to simulate isothermal turbulence. We investigate the five combinations of mixing models described in Table 1, and examine subsonic (⁠|$\mathcal {M} \approx 0.3$|⁠) and supersonic (⁠|$\mathcal {M} \approx 8.4$|⁠) test cases. Table 1. We compare five mixing models involving combinations of the standard implementation, dynamic implementation, as well as the mixing of energy, velocity, and metals. We prefix models using the standard implementation by S-, and those involving the dynamic model by D-. Models which mix thermal energy, velocity (momentum), or metals have combinations of the suffixes u, v, or z, respectively. Name Dynamic Thermal energy Velocity Metals None N/A × × × S-uz × |$\checkmark$| × |$\checkmark$| D-uz |$\checkmark$| |$\checkmark$| × |$\checkmark$| S-uvz × |$\checkmark$| |$\checkmark$| |$\checkmark$| D-uvz |$\checkmark$| |$\checkmark$| |$\checkmark$| |$\checkmark$| Name Dynamic Thermal energy Velocity Metals None N/A × × × S-uz × |$\checkmark$| × |$\checkmark$| D-uz |$\checkmark$| |$\checkmark$| × |$\checkmark$| S-uvz × |$\checkmark$| |$\checkmark$| |$\checkmark$| D-uvz |$\checkmark$| |$\checkmark$| |$\checkmark$| |$\checkmark$| View Large Table 1. We compare five mixing models involving combinations of the standard implementation, dynamic implementation, as well as the mixing of energy, velocity, and metals. We prefix models using the standard implementation by S-, and those involving the dynamic model by D-. Models which mix thermal energy, velocity (momentum), or metals have combinations of the suffixes u, v, or z, respectively. Name Dynamic Thermal energy Velocity Metals None N/A × × × S-uz × |$\checkmark$| × |$\checkmark$| D-uz |$\checkmark$| |$\checkmark$| × |$\checkmark$| S-uvz × |$\checkmark$| |$\checkmark$| |$\checkmark$| D-uvz |$\checkmark$| |$\checkmark$| |$\checkmark$| |$\checkmark$| Name Dynamic Thermal energy Velocity Metals None N/A × × × S-uz × |$\checkmark$| × |$\checkmark$| D-uz |$\checkmark$| |$\checkmark$| × |$\checkmark$| S-uvz × |$\checkmark$| |$\checkmark$| |$\checkmark$| D-uvz |$\checkmark$| |$\checkmark$| |$\checkmark$| |$\checkmark$| View Large In order to mix the fluid over time, we use an identical forcing routine as in Bauer & Springel (2012). The accelerations are calculated in Fourier space and only contain power over a small range of modes corresponding to a spatial range ℓ ∈ [L/2, L] (i.e. the largest scales), and the Fourier mode phases are drawn from an Ornstein–Uhlenbeck process. In the subsonic case, the forcing is purely solenoidal (or incompressible) since the compressive part of the acceleration is removed via a Helmholtz decomposition in Fourier space. It is important to note that Grete, O’Shea & Beckwith (2018) showed that this is not completely correct, and that compressive modes still exist even with purely solenoidal forcing. However, we are comparing the effects of the dynamic model using the same forcing methodology across all of our test cases, and additionally we construct the shear and sub-grid scale stress tensor to be trace-free, removing any contributions from compression of the fluid. We use the exact parameters in table 1 of Bauer & Springel (2012), and point the interested reader to their section 2.2 for the precise details of the driving routine. The systems enter an approximate steady state after t ≳ 5. We measure the probability density functions (PDFs) of Cs in each test in order to determine its sensitivity to the smoothing parameter ε in equation (20). In addition, we measure the PDF of the density field in order to gauge the ability of each model to resolve different density ranges in the turbulent flow. The left-hand panel in Fig. 1 shows the distribution of Cs in the subsonic case as predicted with the dynamic model for three separate resolutions: 643, 1283, and 2563. The median value and the shape of the distributions do not change much with resolution, indicating excellent convergence. At 643 resolution we find a median value Cs = 0.1 and approximately 9.61 per cent of the particles have Cs = 0. In the fluid mechanics literature, as many as 50 per cent of the fluid elements have been reported to have Cs = 0 (Piomelli et al. 1991; Urzay et al. 2013). We also test for convergence in supersonic turbulence (see the middle plot in Fig. 1). Compared to the subsonic case, the dynamic model predicts more particles at Cs = 0, with a total fraction below 50 per cent. The median value is much lower than in the subsonic case, Cs = 0.066. The lower median agrees with calibration results from Colbrook et al. (2017) who found that Cs ≈ 0.05 reproduces the turbulent scaling relationships in supersonic turbulence. Figure 1. View largeDownload slide (Left) Probability density function of Cs, calculated with the dynamic model, in homogeneous subsonic turbulence at three resolutions – 643, 1283, and 2563 particles. The MFM method with a cubic spline kernel is employed, with Nngb = 32. Median values of the predicted Smagorinsky constant are Cs,64 = 0.1005, Cs,128 = 0.1017, Cs,256 = 0.1009 and are well within 2 per cent at maximum difference. (Middle) Supersonic, homogeneous turbulence at three resolutions, with Cs,64 = 0.0668, Cs,128 = 0.0632, Cs,256 = 0.0679. (Right) Probability density function of Cs as we vary the smoothing parameter, ε, in homogeneous, subsonic turbulence using the MFM method, with a quintic spline kernel and Nngb = 64. Figure 1. View largeDownload slide (Left) Probability density function of Cs, calculated with the dynamic model, in homogeneous subsonic turbulence at three resolutions – 643, 1283, and 2563 particles. The MFM method with a cubic spline kernel is employed, with Nngb = 32. Median values of the predicted Smagorinsky constant are Cs,64 = 0.1005, Cs,128 = 0.1017, Cs,256 = 0.1009 and are well within 2 per cent at maximum difference. (Middle) Supersonic, homogeneous turbulence at three resolutions, with Cs,64 = 0.0668, Cs,128 = 0.0632, Cs,256 = 0.0679. (Right) Probability density function of Cs as we vary the smoothing parameter, ε, in homogeneous, subsonic turbulence using the MFM method, with a quintic spline kernel and Nngb = 64. We also investigate the sensitivity of the Cs to the smoothing parameter, ε, from equation (20). The right plot in Fig. 1 shows the PDFs of Cs in homogeneous, subsonic turbulence using the MFM method. We vary ε between 0.2 and 1.0 since ε cannot be greater than 1.0, as it is derived from a series expansion, and should be ≥0 in order to have positive kinetic energies in the smoothed fields (Monaghan 2011). Values in the range 0.7 ≤ ε < 1.0 produce comparable distributions with medians Cs ≈ 0.1. Monaghan (2011) found, using a version of SPH with smoothed velocities, that ε = 0.8 reproduced turbulent flow trends in decaying wall-bounded turbulence. For this reason we employ ε = 0.8 in all of our tests. Turning to the gas properties, Fig. 2 shows histograms of the gas densities in the 643 homogeneous subsonic turbulence simulations, averaged at 150 equally spaced times from t = 10 to t = 25 (inclusive), normalized such that each maximum is Nmax = 1. After t = 10, each turbulent field across all models is in an approximate statistically stationary state with |$\overline{\rho } \approx 1$|⁠. Although the model labels include the -z flag, we do not include metal mixing in our driven turbulence simulations. We retain the suffix for cross-comparison across the different cases in paper, and point the reader to Section 3.3 for a discussion of turbulent metal mixing in an idealized experiment. We note that the None and S-uz cases coincide and the lines in the figure overlap, indicating that internal energy diffusion with a global Cs value has no effect on the density distribution. Figure 2. View largeDownload slide Histograms of the density field in homogeneous subsonic turbulence (643 case), normalized to the bin with the maximum particle count. Density measurements taken from 150 snapshots between t = 10 and t = 25. Here, the None and S-uz cases coincide. The particle density contrast in the S-uvz case is much tighter than the other cases with particles closer to the mean density (ρ = 1). The D-uvz case is able to represent the widest range of densities, including much higher density regions, with the same number of fluid elements. Figure 2. View largeDownload slide Histograms of the density field in homogeneous subsonic turbulence (643 case), normalized to the bin with the maximum particle count. Density measurements taken from 150 snapshots between t = 10 and t = 25. Here, the None and S-uz cases coincide. The particle density contrast in the S-uvz case is much tighter than the other cases with particles closer to the mean density (ρ = 1). The D-uvz case is able to represent the widest range of densities, including much higher density regions, with the same number of fluid elements. First, using the None case as a reference curve, the D-uz case shows a narrower distribution between 0.75 < ρ < 1.15 with a prominent extended tail towards higher densities. This indicates that the D-uz case can represent a wider range of densities in the turbulent cascade. Since the S-uz case follows the None case exactly, this suggests that the localization of the diffusivity D impacts the density resolution much more than the dependence on velocity shear. Introducing momentum diffusion alters the density distributions significantly compared to the None case. In the S-uvz case, the density distribution is tighter and exhibits no apparent wings, with the majority of densities falling in the range 0.75 < ρ < 1.2. Here the increased diffusivity destroys any small-scale structure by causing densities to remain closer to the mean. However, when we employ the D-uvz model, we find the opposite effect – localizing momentum diffusion leads to a wider range of densities in the turbulent gas compared to the all other cases. The effect is strongest at higher densities and therefore we can conclude that the dynamic model is able to resolve higher densities in a turbulent flow at the same mass resolution. 3.2 Keplerian disc In numerical studies of galaxy formation, inherent or artificial dissipation in the hydrodynamical method can cause gas to lose angular momentum and flow radially inwards (Hosono, Saitoh & Makino 2016), i.e. the viscous instability. Numerical simulations require sub-grid diffusive terms as they cannot resolve the viscous scale, however, additional momentum diffusion enhances the viscous instability. This is an important consideration for the constant-coefficient Smagorinsky model: in a simulation of a gaseous disc where the rotational velocity curve depends on the radius, including momentum diffusion will trigger the viscous instability even if there is no turbulence as the radial velocity gradient contributes to |S*|. A quick analytic calculation demonstrates this. Let us consider a 2D idealized rotating gaseous disc with constant surface density that follows a Keplerian velocity profile vϕ ∝ r−1/2. For this disc, |S*| ∝ r−3/2 and inserting this into equation (11) with h = const., we find D ∝ r−3/2. Generally, for any non-constant velocity profile vϕ = vϕ(r), D ∝ ∂rvϕ in the constant density case. In a Keplerian disc, particles near the inner radii of a rotating disc will diffuse the strongest in the standard Smagorinsky model as the difference in velocity between each concentric ring is much higher in this region; leading to the rapid break-up of the disc. One could mitigate the overdiffusion by using a smaller value of Cs in equation (11) but then the model would lose its advantages in turbulent flows. The dynamic model provides a solution to this problem. We use the 2D idealized Keplerian disc as a representative case of an astrophysical laminar shearing flow to illustrate the aforementioned overdiffusion and compare to the results of the dynamic model. We simulate a gas annulus of constant surface density using the MFM method, with particles initialized on circular orbits about the centre. We smooth the inner and outer edges of the annulus in order to suppress numerical instabilities that occur at sharp boundaries. The particles are subject to an external softened gravitational acceleration [|$\boldsymbol{a} = -\boldsymbol{r}(r^2 + \epsilon ^2)^{-3/2}$|] directed towards the centre of the annulus, and follow a corresponding Keplerian velocity profile. This initial condition (IC) is identical to that in section 4.2.4 of Hopkins (2015), with surface density as a function of radial distance, r, \begin{eqnarray*} \Sigma (r) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}(2r)^3 & r \lt 0.5, \\ 1 & 0.5 \leqslant r \leqslant 2, \\ (1 + 10(r - 2))^{-3} & r \gt 2. \end{array}\right. \end{eqnarray*} (24) In the ideal case, the disc should remain intact at any time t > 0. We study the surface densities of the test cases involving thermal energy and momentum diffusion, described in Table 1. As with the driven turbulence experiments, we do not include metals despite the model suffix -z. We include the suffix to allow for cross-comparison across the various physical tests. The leftmost plot in Fig. 3 shows the surface density of the disc in the None model at t ≈ 2torb,r=1. We focus on relatively early times to decouple the effects of inherent numerical diffusion in the MFM method with those of the turbulent mixing models. In the None case, we see that the inner half of the disc is noisy and in the outer region, there are density waves propagating outward, similar to the results in Hopkins (2015). The noise in the inner region, where the orbital time is short, is due to numerical diffusion randomizing the particle motions; short of altering the hydrodynamical solver, this effect is unavoidable and is present in all of the tests we investigate here, with or without mixing. We therefore use the None model as a baseline experiment to compare the four mixing models. Figure 3. View largeDownload slide Normalized surface density differences between each mixing model. The leftmost plot shows the surface density profile of the None case for comparison. The S- models lead to a more rapid break-up of the disc, especially in the case including momentum diffusion (S-uvz). The D-uz model minimizes the difference to the None case. The dynamic diffusion of thermal energy and momentum (D-uvz) leads to an equivalent amount of differences to the S-uz case but the break-up of the disc, due to the overdiffusive S-uvz case, is avoided. Figure 3. View largeDownload slide Normalized surface density differences between each mixing model. The leftmost plot shows the surface density profile of the None case for comparison. The S- models lead to a more rapid break-up of the disc, especially in the case including momentum diffusion (S-uvz). The D-uz model minimizes the difference to the None case. The dynamic diffusion of thermal energy and momentum (D-uvz) leads to an equivalent amount of differences to the S-uz case but the break-up of the disc, due to the overdiffusive S-uvz case, is avoided. In the rest of the four panels in Fig. 3, we show the point-wise difference in surface density between the model in question, Σi(r), and the None case, ΣNone(r), normalized to the mean surface density in the None case; i.e. ΔΣ(r)/〈ΣNone〉 = (Σi(r) − ΣNone(r))/〈ΣNone〉. This allows us to compare the diffusion of energy and momentum spatially, by observing the differences directly on the surface of the disc. In the mixing tests without momentum diffusion, S-uz (second panel in Fig. 3) and D-uz (third panel, Fig. 3), the inner region (0.5 < r < 1.0) of each annulus shows differences compared to the None case. These are due to the false identification of turbulence caused by two effects: (i) both models identify random particle motions, like those in the central region, with turbulence and (ii) the diffusivity scales as D ∝ r−3/2 in the S-uz case. The advantage of the dynamic model is that the radial extent of the affected region is significantly smaller compared to the S-uz case. D-uz predicts much smaller values of Cs (median of Cs ≈ 0.026 in the D-uz case compared to Cs = 0.2 in S-uz) and, therefore, the differences with the None case are mostly limited to the relatively noisy central annulus. From Fig. 3 we see that the addition of momentum diffusion (the S-uvz and D-uvz cases, two rightmost panels) causes increased noise throughout the disc. Specifically for the S-uvz case, the noise in the disc extends radially to the outer boundary, and the material spreads into the central region, compared to the None case (notice the faint pink annulus). There is also a corresponding deficit of gas at R ≈ 0.5 (blue ring) showing that gas has collapsed due to the viscous instability. Comparing D-uvz to the None case, we see that the differences do not extend to the boundary, but rather follow the density waves in the disc caused by natural dissipation in MFM. The central noisy region has the same extent as in the S-uz case but with no apparent in-fall of material into the central region. These results indicate that the dynamic model allows momentum diffusion in laminar shear flows without instigating the viscous instability. In Fig. 4, we show the azimuthally averaged difference between the surface density in each mixing model, Σi, and the None case, ΣNone, in order to get a more quantitative estimate of the extent of overdiffusion. Both S-uz and D-uz are nearly identical to the None case, except in the innermost region (R ≲ 1 for S-uz and R ≲ 0.7 for D-uz) where there are slight fluctuations about the None value. The differences compared to the None case are reduced because of the azimuthal averaging. Figure 4. View largeDownload slide The azimuthally averaged differences in surface density for the Keplerian disc experiments, between each mixing model i and the None case at t ≈ 2torb. The D-uz, D-uvz, and S-uz cases show small fluctuations around the None case whereas the S-uvz model causes large differences due to overdiffusing momentum. Figure 4. View largeDownload slide The azimuthally averaged differences in surface density for the Keplerian disc experiments, between each mixing model i and the None case at t ≈ 2torb. The D-uz, D-uvz, and S-uz cases show small fluctuations around the None case whereas the S-uvz model causes large differences due to overdiffusing momentum. Introducing momentum diffusion results in greater quantitative differences in the disc. The S-uvz model results in significant flows (inwards and outwards) in the annulus relative to the None case. For example, angular momentum transport causes some of the material at R ≈ 0.5 − 0.6 to flow inwards due to loss of angular momentum, and some to flow outward. This occurs throughout the disc, resulting in fluctuations extending to the outer edge of the disc. These flows lead to the S-uvz case showing regions of higher densities between 0.6 ≲ R ≲ 1.8. Comparatively, the small fluctuations in the D-uvz case are similar to the models without momentum diffusion. The dynamic model clearly reduces the impact of the viscous instability, but a question then arises: why does the dynamic model not eliminate the viscous instability completely? The dynamic model formally predicts Cs = 0 for all of the gas particles. Prior to t ≈ 0.4torb, virtually all of the particles have a near zero value of Cs. However, by t ≈ 2torb, turbulent momentum diffusion and the inherent noise in the inner regions lead to a small non-zero distribution centred at approximately Cs ≈ 0.04 (ignoring particles with Cs = 0). Overall, the particles have a median value of Cs ≈ 0.026 when including the 21.4 per cent of particles at Cs = 0. Although we cannot avoid the inherent numerical noise, due to MFM’s Riemann solver, the dynamic method does minimize the damage: in standard implementations, Cs is in the range ≈0.1−0.2 (Garnier et al. 2009) while in the dynamic model, only a small fraction (≈7 per cent) of the particles attain such values. These tests show that in the case of a rotating, laminar shear flow, the dynamic model (D-uvz) can indeed minimize turbulent mixing of thermal energy and momentum – preventing unphysical viscous flows within the disc. Therefore, we recommend incorporating the dynamic Smagorinsky model in all numerical simulations involving rotating galactic discs which simultaneously include turbulent mixing models. 3.3 The Kelvin–Helmholtz instability In a fluid with high-velocity shear, or at the shearing interface between two fluids, rapidly growing perturbations cause mixing within the fluid. In the case of a shear interface, the perturbations cause the two fluids to encroach the boundary, transporting and mixing fluid properties such as thermal energy, momentum, and metals. This is the Kelvin–Helmholtz (KH) instability, and the time-scale characterizing the growth of perturbations is given by, \begin{eqnarray*} \tau _\mathrm{KH} = \frac{(\rho _1 + \rho _2)}{\sqrt{\rho _1\rho _2}}\frac{\lambda }{\Delta v}, \end{eqnarray*} (25) where ρ1, ρ2 are the densities of the two fluids, Δv is the velocity difference, and λ is the wavelength of the perturbation (Chandrasekhar 1961). In the t < τKH regime (the linear regime), the flow has not completely transitioned to turbulence and sub-grid turbulent mixing does not dominate the resolved mixing. Shear flows, and hence KH instabilities, are ubiquitous in galactic environments: ram pressure stripping of galaxies falling into groups and clusters, galaxy mergers, galactic winds streaming into the CGM, and a myriad of other processes involve the KH instability. The constant-coefficient Smagorinsky model overdiffuses in such situations because it identifies shearing motion with turbulence, as discussed in Section 3.2, and the presence of high shear increases the diffusivity to the maximum at the interface. Ideally, sub-grid turbulent mixing models that better reflect physical reality are preferable. More precisely, models that capture unresolved mixing in turbulent situations and avoid diffusion in laminar shear flows. We investigate a simple 2D KH test in order to demonstrate the overdiffusive nature of the standard Smagorinsky model and determine if the dynamic model mitigates the problem. We set up a 2D configuration of 2562 ideal gas particles in a square (L = 1) domain, with constant pressure, and with an initial density profile, \begin{eqnarray*} \rho (y) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\rho _2 - (\Delta \rho / 2)\exp {[(y - 1/4)/\Delta y]}, & 0\leqslant y \lt 1/4\nonumber \\ \rho _1 + (\Delta \rho / 2)\exp {[(1/4 - y)/\Delta y]}, & 1/4 \leqslant y \lt 1/2\nonumber \\ \rho _1 + (\Delta \rho / 2) \exp {[(y - 3/4)/\Delta y]}, & 1/2 \leqslant y \lt 3/4\nonumber \\ \rho _2 - (\Delta \rho / 2) \exp {[(3/4 - y)/\Delta y]}, & 3/4 \leqslant y \lt 1\nonumber \\ \end{array}\right. \nonumber \\ \end{eqnarray*} (26) and initial velocity profile, \begin{eqnarray*} v_x(y) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}-1/2 + (1/2) \exp {[(y - 1/4)/\Delta y]}, & 0\leqslant y \lt 1/4\nonumber \\ 1/2 - (1/2) \exp {[(1/4 - y)/\Delta y]}, & 1/4 \leqslant y \lt 1/2\nonumber \\ 1/2 - (1/2) \exp {[(y - 3/4)/\Delta y]}, & 1/2 \leqslant y \lt 3/4\nonumber \\ -1/2 + (1/2) \exp {[(3/4 - y)/\Delta y]}, & 3/4 \leqslant y \lt 1.\nonumber \\ \end{array}\right. \nonumber \\ \end{eqnarray*} (27) We choose ρ2 = 2, ρ1 = 1, and Δy = 0.025, and introduce a sine wave velocity perturbation, with period T = 2 and amplitude A = 0.01, at t = 0. This gives a perturbation wavelength λ = 1/2, and therefore τKH ≈ 0.71. In addition, we add a uniform passive scalar tracer of concentration Q = 1 to all gas particles in the range 0 ≤ y < 1/4 and 3/4 ≤ y < 1, which is the higher density gas. We focus our analysis on the evolution of the tracer concentration Q. Fig. 5 shows the tracer concentration at t = 0.28τKH (i.e. in the linear regime). The five panels show the five mixing models described in Table 1. We first consider the None case. The tracer concentration Q follows the high-density regions of the experiment and we see individual particles advecting across the shear interface. Although there appears to be less tracer on this interface (orange-red line), it is impossible for particles to exchange tracer in the None case. The interpolation method we employ causes this effect as it smooths the particle properties over the resolution scale h, leading to a minuscule amount of artificial mixing. Figure 5. View largeDownload slide A comparison of tracer concentrations Q, in the 2D Kelvin–Helmholtz instability test at t = 0.28τKH, simulated with the MFM method. The columns represent the five turbulent mixing models. In the None case, gas particles cannot exchange the tracer concentrations and the interface of the instability remain unsmoothed. In this case, inter-fluid mixing only occurs when particles move across the boundary. The S- cases diffuse rapidly due to the presence of strong shear at the boundary, and the tracer engulfs the whorls during the early evolution of the instability. The D- cases provide a compromise between the two situations – they limit diffusion strictly to the interface between the two fluids, and the internal structure of the whorls are distinguishable. Figure 5. View largeDownload slide A comparison of tracer concentrations Q, in the 2D Kelvin–Helmholtz instability test at t = 0.28τKH, simulated with the MFM method. The columns represent the five turbulent mixing models. In the None case, gas particles cannot exchange the tracer concentrations and the interface of the instability remain unsmoothed. In this case, inter-fluid mixing only occurs when particles move across the boundary. The S- cases diffuse rapidly due to the presence of strong shear at the boundary, and the tracer engulfs the whorls during the early evolution of the instability. The D- cases provide a compromise between the two situations – they limit diffusion strictly to the interface between the two fluids, and the internal structure of the whorls are distinguishable. When we allow for the turbulent mixing of the tracer in the S- and D- cases, we see diffusion occurring along the shear interface. In both S- models, Q ∼ 10−2 at the interface and the tracer engulfs the initially pristine, lower density gas in the whorls. This is in contrast to the two D- cases where the majority of diffusion occurs in the whorls themselves with comparatively little along the rest of the interface. The constant-coefficient Smagorinsky model diffuses the most because of the false identification of strong turbulence through the norm of the shear tensor |S*|. The velocity profile in equation (27) shows that although the flow is laminar, there is a gradient, |∂vx/∂y| > 0, across the entire domain. Based on our arguments for the Keplerian disc in Section 3.2, the fact that D ∝ |∂vx/∂y| directly leads to the overdiffusion of the tracer. The diffusion coefficient D also depends on |S*| in the D- models but in the dynamic model most of the values of Cs are near zero. Only in the whorls do we find values of Cs as large as Cs = 0.2 but these are limited to this region – where the transition to turbulence begins. Overdiffusion in non-turbulent shear flows can have important consequences for the gas in numerical studies of galaxy formation. If the KH time-scale is longer than the gas cooling time, then the gas is susceptible to overcooling in the whorls due to metals transferred to the region. This does not accurately capture what physically occurs at the sub-grid level; in the linear regime, the fluids do not mix completely. The dynamic model solves this issue by minimizing mixing along the interfaces while allowing it to proceed in the whorls where, in principle, the KH instability continues down to unresolved scales. 4 ISOLATED DISC GALAXY Here we investigate an isolated galaxy in order to test the effects of localized diffusion in a more realistic, physical environment. 4.1 Initial conditions We follow the method outlined in Springel, Di Matteo & Hernquist (2005) to set up the ICs, using the galstep package11 (Ruggiero & Lima Neto 2017) in combination with dice12 (Perret et al. 2014). The galaxy is a Milky Way-like system (Sokołowska et al. 2016) consisting of a dark matter halo of mass Mh = 1012 M⊙, a gaseous halo of mass Mgh = 3 × 1010 M⊙, a stellar bulge of mass Mb = 1010 M⊙, and gas and stellar discs of masses Mg = 1010 M⊙ and Ms = 5 × 1010 M⊙, respectively. The dark matter and bulge components follow a Hernquist density profile with scale factors a = 47 kpc and a = 1.5 kpc, respectively. The stellar and gaseous discs follow an exponential density profile with a radial scale Rd = 3.5 kpc, and the scale heights for these components are z0 = 0.7 kpc and z0 = 0.0175 kpc, respectively. We initialize the gaseous halo metallicity at Zgh = 10−3 Z⊙, and the gaseous disc metallicity at Zg = Z⊙/3, with Z⊙ = 0.02 (Anders & Grevesse 1989). Our fiducial run is carried out at a gas mass resolution of Mg,res = 5 × 103 M⊙, along with the softening values specified in Table 2. We evolve the disc for 2 Gyr in an isolated (non-cosmological) setting. Table 2. Parameters for the isolated galaxy. Component Particle mass (M⊙) Min. softening (pc) Npart Gas 5.0 × 103 1.4 8 × 106 Halo 5.0 × 105 12.0 2 × 106 Disc 5.0 × 105 3.2 1 × 105 Bulge 2.5 × 105 1.4 4 × 104 Component Particle mass (M⊙) Min. softening (pc) Npart Gas 5.0 × 103 1.4 8 × 106 Halo 5.0 × 105 12.0 2 × 106 Disc 5.0 × 105 3.2 1 × 105 Bulge 2.5 × 105 1.4 4 × 104 View Large Table 2. Parameters for the isolated galaxy. Component Particle mass (M⊙) Min. softening (pc) Npart Gas 5.0 × 103 1.4 8 × 106 Halo 5.0 × 105 12.0 2 × 106 Disc 5.0 × 105 3.2 1 × 105 Bulge 2.5 × 105 1.4 4 × 104 Component Particle mass (M⊙) Min. softening (pc) Npart Gas 5.0 × 103 1.4 8 × 106 Halo 5.0 × 105 12.0 2 × 106 Disc 5.0 × 105 3.2 1 × 105 Bulge 2.5 × 105 1.4 4 × 104 View Large 4.2 Galactic physics 4.2.1 Star formation and cooling We employ the sub-grid multiphase ISM model of Springel & Hernquist (2003a), which places the ISM gas (nH > n*,crit, where n*,crit is the star formation density threshold) on an effective equation of state (EoS). This model has been used extensively in numerical galaxy formation studies (for recent examples, see Genel et al. 2014; Schaye et al. 2015; Grand et al. 2017), and provides well-converged results (Springel et al. 2005). For primordial and metal-line cooling, we use the grackle-2.1 cooling library (Smith et al. 2017) in combination with the UV background from Faucher-Giguère et al. (2009). In this sub-grid model, stars form stochastically, on a star formation time-scale tsfr, from gas that reaches a density n*,crit. Here we take n*,crit = 0.1 cm−3 and tsfr = 2.1 Gyr, which give a good fit to the Kennicutt law (Kennicutt 1998; Springel & Hernquist 2003a). We differ from the original model in the choice of the initial mass function (IMF). We assume the Chabrier IMF (Chabrier 2003) instead of the Salpeter IMF (Salpeter 1955). 4.2.2 Feedback Due to lack of resolution, it is necessary to include a sub-grid model for feedback from stars, including the effects of supernovae, stellar radiation, and stellar winds. For massive stars, we adopt the scheme used in the MUFASA simulations (Davé, Thompson & Hopkins 2016; Davé et al. 2017) and refer the reader to these references for a detailed description. In brief, stellar feedback is expected to drive galactic outflows and, in the MUFASA approach, stellar feedback directly launches a kinetic wind via a two-parameter model that characterizes the net effect of stellar feedback into a mass loading factor η, and the wind speed vw. These parameters are calibrated to the FIRE wind scalings (Muratov et al. 2015), where η scales with the stellar mass of the host galaxy and vw scales with the galaxy’s circular velocity. We fix η and vw to the values for our isolated system based on equations 6 and 7 in Davé et al. (2016). At launch, the outflow hydro-dynamically decouples, and only recouples if the wind speed drops to 50 per cent of the local sound speed, the density of the ambient medium is 1 per cent of the ISM density, or the outflow has travelled for 2 per cent of the Hubble time at launch (Davé et al. 2016). For a more detailed description of the decoupled outflow model, see Springel et al. (2005), Oppenheimer & Davé (2008), Liang et al. (2016), and Davé et al. (2016). The contribution from supernovae type Ia (SNIa) are modelled following Scannapieco & Bildsten (2005) as a prompt and delayed component, where the prompt component occurs simultaneously with SNII, and the delayed component begins 0.7 Gyr after the star formation time. The prompt component is assumed to release 1051 erg of thermal energy to the star-forming gas, whereas the delayed component is added in a kernel-weighted manner to the nearest 16 gas particles. 4.2.3 Chemical enrichment The chemical enrichment of gas is paramount to the study of turbulent mixing as the metallicity follows the mixing of energy and momentum, and, therefore, provides a tracer for the diffusion equation. We track 11 chemical elements in our isolated and cosmological simulations: H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe. These elements are produced from three sources in the simulations: SNIa, SNII, and the winds from AGB stars. For SNIa, the prompt component returns mass to the ISM and enriches the star-forming gas instantaneously. Each SNIa is assumed to release 1.4 M⊙ of metals, with yields from Iwamoto et al. (1999). For the delayed component, stars deposit metals over their nearest 16 neighbouring gas particles in a kernel-weighted fashion. SNII return mass and enrich the gas via the instantaneous recycling approximation (Springel & Hernquist 2003a; Oppenheimer & Davé 2008; Davé et al. 2016) following, \begin{eqnarray*} \Delta Z_\mathrm{\mathit{ i}} = (1 - f_\mathrm{SN})\cdot y_\mathrm{\mathit{ i}}(Z)\cdot \frac{\Delta t}{t_\mathrm{sfr}} , \end{eqnarray*} (28) where fSN is the fraction of stars in the Chabrier IMF expected to go supernova, yi(Z) is the metallicity dependent yield of species i, Δt is the time-step, and tsfr is the aforementioned star formation time-scale. The SNII yields follow Nomoto et al. (2006) and are a function of the metallicity of the gas receiving the metals. Following Davé et al. (2016), the SNII yields are reduced by a factor of 0.5 in order to match the mass–metallicity relationship. SNII also return mass into the gas via the instantaneous recycling approximation. For AGB stars, chemical enrichment is done in a kernel-weighted fashion over the nearest 16 neighbours. AGB yields are obtained from a lookup table as a function of age and metallicity based on the study in Oppenheimer & Davé (2008). The mass-loss rates of the AGB stars are calculated from a lookup table based on Bruzual & Charlot (2003) stellar models. 4.3 Results: disc stability Fig. 6 shows radial disc profiles of the surface density (left-hand panel) and entropy13 (right-hand panel) in our isolated galaxies after 1.27 Gyr of evolution. The dotted line represents the IC for all models. The differences between models only appear after 1 Gyr (≈4 rotations in the mid-disc) and continue until star formation consumes the bulk of the gas after 2 Gyr. Figure 6. View largeDownload slide Isolated galaxy radial disc profiles at t = 1.27 Gyr, averaged azimuthally and vertically between ±0.5 kpc from the plane of the disc. (Left) The gas surface density. (Right) A monotonic measure of the gas entropy. The dotted line in both plots represents the corresponding profile in the IC, which is identical for all models. Over time, the normalization of the surface density decreases due to gas consumption and stellar feedback. The dynamic model with thermal energy and momentum diffusion (D-uvz) produces a more stable disc whereas the constant-coefficient model (S-uvz) shows a more concentrated central region. Figure 6. View largeDownload slide Isolated galaxy radial disc profiles at t = 1.27 Gyr, averaged azimuthally and vertically between ±0.5 kpc from the plane of the disc. (Left) The gas surface density. (Right) A monotonic measure of the gas entropy. The dotted line in both plots represents the corresponding profile in the IC, which is identical for all models. Over time, the normalization of the surface density decreases due to gas consumption and stellar feedback. The dynamic model with thermal energy and momentum diffusion (D-uvz) produces a more stable disc whereas the constant-coefficient model (S-uvz) shows a more concentrated central region. First, we consider the None case. The radial surface density (left-hand panel, Fig. 6) gives a measure of the stability of the disc. We see that by t = 1.27 Gyr, compared to the shape of the IC density profile, the gas has moved inwards towards the centre, especially within R < 3 kpc. A combination of the bar instability and inherent numerical dissipation causes the gas to concentrate inside R ≈ 2 kpc, whereas gas consumption and galactic winds due to supernova feedback cause the difference in normalization compared to the IC. Correspondingly in the right-hand panel, there is an order-of-magnitude increase in entropy from R = 1 kpc to R = 3 kpc. In the S-uz and D-uz models, we see only minor difference in the surface density and entropy compared to the None case. Evidently, thermal energy diffusion, combined with metal diffusion, has negligible impact on the structure of the disc. When we introduce momentum diffusion in the S-uvz and D-uvz cases, the differences compared to the None case at t = 1.27 Gyr are more significant. Fig. 6 shows that in the S-uvz case there is an order-of-magnitude deficit of gas surface density (left plot) between R = 1 kpc and R = 3 kpc with a corresponding jump in entropy (right plot). Overdiffusion of momentum due to the diffusivity scaling strongly with the shear causes the inward flowing gas to be more centrally concentrated compared to the None case, and also engenders an outward flow leading to slightly higher density (again, compared to None) at R ≈ 4–7 kpc. Once the instability occurs, the effect accelerates and the trend remains throughout the evolution of the disc. However, with the dynamic model (D-uvz), less gas gets redistributed. From the surface density, we see that the disc stabilizes when momentum is diffused locally, based on the turbulent character of the flow. 4.4 Results: metal distribution functions Supernovae and stellar winds inject energy and metals into the ISM, engendering turbulent motion that mixes and spreads thermal energy, momentum, and metals throughout the medium. These processes also drive a galactic-scale wind that deposits metals and drives turbulence in CGM (Evoli & Ferrara 2011). In order to obtain a measure of energy and metal mixing in the five models under consideration (see Table 1), we examine MDFs in three phases of gas: (1) gas with density above the star formation critical density (nH > 0.1 cm−3) – i.e. the ISM gas, (2) warm CGM gas in the halo in the range 105 K <T < 106 K, and (3) hot CGM gas with T > 106 K. The gas density in the latter two phases is nH ≤ 0.1 cm−3. In our isolated system, we find that the cool, non-star-forming gas (T < 104 K) is mostly at the outskirts of the halo and that there are no significant variation in its MDFs between models. Therefore, we do not discuss this phase further in the isolated case. Fig. 7 shows the MDFs for the ISM, the warm CGM, and the hot CGM, at three separate evolutionary times: t = 0.25 Gyr (top row), 0.50 Gyr (middle row), and 0.75 Gyr (bottom row) after the ICs for the MFM method. Henceforth, we use [Z] ≡ log (Z/Z⊙) as a proxy for metallicity. Figure 7. View largeDownload slide The MDFs for the isolated disc galaxy in solar units (Z⊙ = 0.02). Each column represents a specific region of T−ρ phase space. The leftmost column represents the ISM, whereas the middle and rightmost columns show the MDFs for the warm (105 K <T < 106 K) and hot (T > 106 K) CGM, respectively. Each row represents a different time in Gyr: t = 0.25 Gyr, 0.50 Gyr, and 0.75 Gyr from top to bottom, respectively. Dynamic diffusion slows metal mixing in the ISM and causes increased mixing in the CGM, at approximately the same levels as the constant-coefficient Smagorinsky model. Figure 7. View largeDownload slide The MDFs for the isolated disc galaxy in solar units (Z⊙ = 0.02). Each column represents a specific region of T−ρ phase space. The leftmost column represents the ISM, whereas the middle and rightmost columns show the MDFs for the warm (105 K <T < 106 K) and hot (T > 106 K) CGM, respectively. Each row represents a different time in Gyr: t = 0.25 Gyr, 0.50 Gyr, and 0.75 Gyr from top to bottom, respectively. Dynamic diffusion slows metal mixing in the ISM and causes increased mixing in the CGM, at approximately the same levels as the constant-coefficient Smagorinsky model. First, we focus on the ISM shown in the left column of Fig. 7. At t = 0.25 Gyr, the MDF of the None case is bimodal with a narrow peak at [Z] ≈ −0.3, a broad distribution at [Z] < −0.5 with a peak at [Z] = −1.2, and a dearth at [Z] = −0.5. Gas in the ISM with metallicity lower than the initial [Z] ≈ −0.5 comes from the lower metallicity (initial [Z] = −3.0) halo gas that cools on to the disc and is steadily enriched. Recall that in the None case, gas cannot exchange metals between particles and therefore it gives an upper bound on the mixing time-scale for a given metal injection rate from stars. The D- cases follow the None case with minimal differences, therefore we conclude that the D- cases provide minimal turbulent mixing in the ISM. The S- cases share the basic shape as the other cases but the lower metallicity component is narrower (there is very little gas with [Z] < −2.0), and is shifted to the right, with peak [Z] ≈ −0.9. The overdiffusivity of the S- cases is apparent and is caused by the high levels of shearing motion in the supersonic ISM, leading to high values of the diffusivity D, similar to the Keplerian disc in Section 3.2. At t = 0.5 Gyr in the ISM (left column, middle row), the low-metallicity component of the MDF for all cases has narrowed and shifted towards higher metallicity. In the None case, the tightening is due to stars depositing metals into the medium and driving all of the gas towards highly metallicities. Here, the D- cases follow the None case closely – as at t = 0.25 Gyr – while the S- cases have a smaller spread. We attribute the differences between the D- and S- cases to the dynamic model predicting median diffusivities orders-of-magnitude lower than the constant-coefficient Smagorinsky model in the ISM at t > 0.25 Gyr. By t = 0.75 Gyr, all cases have tight distributions with very few particles having [Z] < −1.5. None and D- cases continue to exhibit very similar distributions with a larger spread compared to the S- cases. Now we consider the MDFs of the warm CGM, in the middle column of Fig. 7. In the None case at t = 0.25 Gyr, the distribution is bimodal with peaks at [Z] ≈ −0.5 and [Z] = −3.0. The peak at [Z] = −3.0 is the initial halo metallicity whereas the peak at [Z] ≈ −0.5 is due to enriched gas either pushed out of the ISM by stellar winds or at the interface between the CGM and ISM. The D- and S- cases closely follow the None case, but with a slight spread towards [Z] > −3.0 for the low-metallicity gas. Next we consider the t = 0.5 Gyr and t = 0.75 Gyr cases. In all of the experiments, the flow of the enriched ISM gas into the CGM results in an increase in the fraction of particles near [Z] ≈ −0.5. The differences near [Z] = −3.0 between the None and S-/D- cases are due to the inability of the gas in former case to mix metals. For the None case, a spread only occurs if the gas enters the CGM at [Z] > −3.0 from the ISM, or the CGM gas at the ISM–CGM boundary becomes enriched via the kernel-weighting approach associated with delayed SNIa and AGB wind feedback. Recall that in this implementation, metals are deposited in a kernel-weighted fashion over the nearest 16 neighbouring particles. Therefore, gas classified as CGM, that is spatially adjacent to the star-forming gas in the ISM, can be enriched at low levels. The gas in the turbulent mixing cases consistently exchange metals if the value of Cs is non-zero, and the gas-phase metals in the halo further mix – leading to the greater spread in the lower peak at [Z] = −3.0, corresponding to the bulk of the initial gas. The D- cases show the aforementioned spread in the metallicity in the range −3.0 < [Z] < −1.75, indicating increased mixing due to sub-grid turbulence, but then coincides with the None case for metallicities [Z] > −2.0. The S- and D- cases share a similar distribution below [Z] < −2.25 but the S- cases show more particles with metallicities in the range −2.25 < [Z] < −0.75, indicating even greater mixing than the D- cases. The MDFs in the hot CGM, the rightmost column of Fig. 7, follow similar trends to the those in the warm CGM. At t = 0.25 Gyr for the None case, the peak due to enriched ([Z] ≈ −0.3) ISM gas entering the halo is prominent. There is also a small spread at [Z] = −3.0 for reasons already noted. At t ≥ 0.5 Gyr, the None case undergoes slight evolution and by t = 0.75 Gyr, a mild positive slope develops in the range −2.5 ≤ [Z] ≤ −0.75. The D- and S- cases follow near identical evolution at these later stages of the simulation, each with a spread in the distribution at [Z] = −3.0. Here the differences between the S- and D- cases are minimal since the hot gaseous halo is turbulent and dominated by random motions, driven by the galactic outflows. Our isolated experiment demonstrates that by endowing particles with diffusivities based on the local fluid properties, we obtain significant differences in the ISM. The constant-coefficient Smagorinsky model causes the MDFs in the ISM to rapidly tighten towards the mean value, whereas the dynamic model predicts diffusivities orders-of-magnitude smaller and the corresponding MDFs closely follow the None case. Simultaneously, in the hot turbulent halo, the constant-coefficient and dynamic Smagorinsky models produce similar distributions by t = 0.5 Gyr due to the latter having higher values of Cs. Overall, for the constant-coefficient model, non-negligible shear in all gas phases causes the rapid diffusion of fluid properties, whereas the dynamic model allows different regions of phase space to undergo unique evolution in terms of the MDFs. We stress the importance of the unique evolution of both phases of gas: with a constant Cs, it is impossible to capture the decreased mixing in the ISM while simultaneously capturing the high level of mixing in the hot turbulent halo. The dynamic model provides an interesting avenue for follow-up study with zoom-in simulations, in order to gauge the impact on the ISM and CGM in a cosmological context. 5 COSMOLOGICAL VOLUMES As we demonstrate in the previous section, turbulent mixing can alter the distribution of metals in various gas phases of an isolated disc experiment depending on the localization of Cs. In a realistic cosmological environment, the evolution of galaxies is much more complex due to interactions between the galaxies and their environments. These interactions include galaxy mergers, gas inflows and outflows, tidal interactions, ram pressure stripping, etc. – all of which contribute to the production of turbulence in the galactic environments (Iapichino et al. 2013; Schmidt et al. 2016). The resulting turbulence redistributes thermal energy, momentum, and metals, and must be included in numerical studies of galaxy evolution in order to have a self-consistent treatment of the physical models. In contemporary Lagrangian-based numerical cosmological experiments, SPH has been employed in simulation programs such as EAGLE (Schaye et al. 2015), OWLS (Schaye et al. 2010), and Romulus (Tremmel et al. 2017), whereas recently MFM method has been employed in the MUFASA (Davé et al. 2016, 2017) simulations. These experiments have produced a wealth of results for understanding galaxy evolution and gas properties (see Somerville & Davé 2015; Naab & Ostriker 2017 for a summary), despite the fact that it is not possible for contemporary models to include all of the relevant physics.14 Including sub-grid turbulent mixing could alter the results of such large-scale simulations. Indeed, Tremmel et al. (2018) argue that turbulent mixing is critical for efficient redistribution of thermal energy released during active galactic nuclei episodes, and in previous sections, we showed that metal redistribution is also affected by turbulent mixing. In models where the overdiffusive Smagorinsky model is used, the dynamic model can lead to differences in, for example, gas-phase metal abundances and hence, stellar abundances. In principle, the differences in the manner and the rate at which metals are distributed could also affect the formation sites and population statistics of Population III stars and direct collapse seed SMBHs (see Section 1). In this section, we examine a set of cosmological simulations in order to test the effects of the dynamic model on the global gas enrichment levels and distributions. In what follows, we adopt the MUFASA model (Davé et al. 2016) in combination with the diffusion models we describe in Table 1. We choose to use the MUFASA model partly because it is the only cosmological model that has been implemented using the MFM method at the present. The ICs were created using a modified version of grafic-215 (Bertschinger 2011) and the parameters describing our simulations are listed in Table 3. In the following subsection, we briefly describe the MUFASA models, and point interested readers to Davé et al. (2016) and Davé et al. (2017) for a more detailed explanation of the physical models. In Sections 5.3 and 5.4 we examine the global gas-phase metallicity fractions and the MDFs, respectively, from the simulation suite. Table 3. Parameters for the cosmological simulations. We use the Planck Collaboration XIII (2015) cosmological model. Simulation parameters L 25 Mpc h−1 N 2 × 2563 mg 1.26 × 107 M⊙ h−1 mdm 6.88 × 107 M⊙ h−1 εsoft,min 0.5 kpc h−1 zinit 70 Tinit 59 K Cosmological model Ωm,0 0.308 ΩΛ,0 0.692 Ωb,0 0.048 h 0.678 σ8 0.815 ns 0.968 Simulation parameters L 25 Mpc h−1 N 2 × 2563 mg 1.26 × 107 M⊙ h−1 mdm 6.88 × 107 M⊙ h−1 εsoft,min 0.5 kpc h−1 zinit 70 Tinit 59 K Cosmological model Ωm,0 0.308 ΩΛ,0 0.692 Ωb,0 0.048 h 0.678 σ8 0.815 ns 0.968 View Large Table 3. Parameters for the cosmological simulations. We use the Planck Collaboration XIII (2015) cosmological model. Simulation parameters L 25 Mpc h−1 N 2 × 2563 mg 1.26 × 107 M⊙ h−1 mdm 6.88 × 107 M⊙ h−1 εsoft,min 0.5 kpc h−1 zinit 70 Tinit 59 K Cosmological model Ωm,0 0.308 ΩΛ,0 0.692 Ωb,0 0.048 h 0.678 σ8 0.815 ns 0.968 Simulation parameters L 25 Mpc h−1 N 2 × 2563 mg 1.26 × 107 M⊙ h−1 mdm 6.88 × 107 M⊙ h−1 εsoft,min 0.5 kpc h−1 zinit 70 Tinit 59 K Cosmological model Ωm,0 0.308 ΩΛ,0 0.692 Ωb,0 0.048 h 0.678 σ8 0.815 ns 0.968 View Large 5.1 MUFASA The MUFASA simulations include the sub-grid models we described in Section 4.2 and sub-sections therein, with some modifications to the star formation recipe and feedback described in the following sub-sections. 5.1.1 Star formation Star formation is based on the molecular gas model of Krumholz, McKee & Tumlinson (2009), and the implementation dynamically calculates the fraction of molecular hydrogen in gas particles, |$f_\mathrm{H_2}$|⁠, based on the gas surface density and the metallicity. For more precise details, see Davé et al. (2016) and Krumholz et al. (2009). The star formation rate follows, \begin{eqnarray*} \frac{\mathrm{d}\rho _*}{\mathrm{d}t} = \epsilon _* \frac{\rho f_\mathrm{H_2}}{t_\mathrm{dyn}}, \end{eqnarray*} (29) where tdyn = (Gρ)−1/2 is the local dynamical time, ρ is the density of the gas, and ε* = 0.02 is the efficiency of star formation (Kennicutt 1998). The critical density of star formation is taken at n*,crit = 0.2 cm−3 following Davé et al. (2016). For our cosmological simulations, we do not employ the sub-grid ISM model of Springel & Hernquist (2003a), however we still require that the Jeans mass is resolved. Therefore, in order to prevent numerical fragmentation at high densities, an artificial pressure is applied above a density \begin{eqnarray*} n_\mathrm{th} = \frac{3}{4\pi \mu m_\mathrm{p}} \bigg (\frac{5k_\mathrm{b} T_\mathrm{0}}{G\mu m_\mathrm{p}}\bigg)^3 \bigg (\frac{1}{N_\mathrm{ngb} m_\mathrm{g}}\bigg)^2, \end{eqnarray*} (30) where mg is the gas particle mass, μ = 1.22, T0 = 104 K, and Nngb = 64 is the number of neighbours. The pressure is applied in the form of a minimum temperature (Teyssier et al. 2011; Davé et al. 2016), \begin{eqnarray*} T_\mathrm{JMT} = T_\mathrm{0} \bigg (\frac{n}{n_\mathrm{th}}\bigg)^{1/3}. \end{eqnarray*} (31) 5.1.2 Feedback In Section 4.2.2, we described the decoupled-wind model for massive star feedback in addition to feedback from supernova type-Ia (SNIa) and asymptotic giant branch (AGB) stars. We use the same in the following experiments. The simulations we present here do not include explicit active galactic nucleus (AGN) feedback. AGN feedback is thought to be necessary to prevent excessive cooling and quench star formation in massive systems (see King & Pounds 2015 for a recent review), but our smaller simulation volumes do not include many such systems. We do, however, include an effective AGN feedback model from the MUFASA simulations that mimics the action of AGN feedback, and suppresses cooling of the diffuse halo gas. Specifically, gas that is not self-shielded in haloes with Mhalo > Mq [where Mq = (0.96 + 0.48z) × 1012 M⊙] is heated to 20 per cent above the virial temperature of the halo (Mitra, Davé & Finlator 2015). The virial temperature follows (Balogh, Babul & Patton 1999; Voit 2005), \begin{eqnarray*} T_\mathrm{vir} = 9.52\times 10^7 \bigg (\frac{M_\mathrm{halo}}{10^{15} \mathrm{M}_\odot \mathrm{h}^{-1}}\bigg)^{2/3}\,\, \mathrm{K}. \end{eqnarray*} (32) 5.2 Gas phases We define five separate gas phases for the following subsections, and examine their properties in a global sense over the entire simulation volume. We give a summary in Table 4. The definitions are largely from Davé et al. (2010), except for the definition of the gas associated with the CGM of galaxies. For density, we cut the gas phase space using two thresholds ρbound and ρ*, with ρbound = ρbound(z) following, \begin{eqnarray*} \frac{\rho _\mathrm{bound}(z)}{\Omega _\mathrm{b}(z)\rho _\mathrm{c}(z)} = 6\pi ^2 \bigg (1 + 0.4093 \bigg (\frac{1}{\Omega _\mathrm{m}(z)} - 1\bigg)^{0.9052}\bigg) - 1, \end{eqnarray*} (33) where \begin{eqnarray*} \Omega _\mathrm{m}(z) = \frac{\Omega _\mathrm{m,0} (1 + z)^3}{\Omega _\mathrm{m,0}(1 + z)^3 + \Omega _\mathrm{\Lambda ,0}}\,\,, \end{eqnarray*} (34) \begin{eqnarray*} \Omega _\mathrm{b}(z) = \frac{\Omega _\mathrm{b,0} (1 + z)^3}{\Omega _\mathrm{m,0}(1 + z)^3 + \Omega _\mathrm{\Lambda ,0}}\,\,, \end{eqnarray*} (35) ρc(z) = 3(H(z))2/(8πG), and |$H(z) = H_\mathrm{0} \sqrt{\Omega _\mathrm{m,0}(1 + z)^3 + \Omega _\mathrm{\Lambda , 0}}$|⁠. For the second density cut, we adopt ρ* = 4.4 × 10−25 g cm−3, the star formation density threshold. We have also applied a single temperature cut at T5 = 105 K to separate the warm-hot intergalactic medium (WHIM) from the cold diffuse gas (DIFF). Following Torrey et al. (2017), we apply a cut to distinguish the hot halo gas (HHG) from the cool CGM, \begin{eqnarray*} \log {\bigg (\frac{T}{10^6\,\, \mathrm{K}}\bigg)} = 0.25 \log {\bigg (\frac{n}{405\,\, \mathrm{cm}^{-3}}\bigg)}, \end{eqnarray*} (36) where T is the gas temperature, and n is the gas density. We define the HHG to be above the temperature threshold in equation (36), and the CGM to be below. Table 4. We separate gas in our cosmological simulations into five phases: the ISM, cool CGM, HHG, WHIM, and cool diffuse gas (DIFF). ρ* is the star formation threshold, n*,crit = 0.2 cm−3 and we give ρbound in equation (33). Name Density range Temperature range ISM ρ > ρ* Any CGM ρ* > ρ > ρbound Below equation (36) HHG ρ > ρbound Above equation (36) WHIM ρ < ρbound T > 105 K DIFF ρ < ρbound T < 105 K Name Density range Temperature range ISM ρ > ρ* Any CGM ρ* > ρ > ρbound Below equation (36) HHG ρ > ρbound Above equation (36) WHIM ρ < ρbound T > 105 K DIFF ρ < ρbound T < 105 K View Large Table 4. We separate gas in our cosmological simulations into five phases: the ISM, cool CGM, HHG, WHIM, and cool diffuse gas (DIFF). ρ* is the star formation threshold, n*,crit = 0.2 cm−3 and we give ρbound in equation (33). Name Density range Temperature range ISM ρ > ρ* Any CGM ρ* > ρ > ρbound Below equation (36) HHG ρ > ρbound Above equation (36) WHIM ρ < ρbound T > 105 K DIFF ρ < ρbound T < 105 K Name Density range Temperature range ISM ρ > ρ* Any CGM ρ* > ρ > ρbound Below equation (36) HHG ρ > ρbound Above equation (36) WHIM ρ < ρbound T > 105 K DIFF ρ < ρbound T < 105 K View Large It is important to note that these density and temperature cuts do not distinguish gas domains that precisely correspond to their associated galactic or inter-galactic regions. At high redshift (z ≳ 5), most of the gas classified as the WHIM phase is spatially located in the region that mostly corresponds to the CGM/HHG, and corresponds to outflowing galactic winds from early star formation, which gives rise to low-density, high-temperature gas. Similarly, a fraction of the gas classified as the CGM at z ≳ 1 is in the cores of the cosmic filaments. However, by z = 0, the majority of what we consider the CGM is indeed inside the haloes, and the cores of the filaments eventually end up in the WHIM phase. This introduces a transition from CGM to WHIM gas that might not be obvious at first glance. We have not tried to address these trends or optimize the phase cuts because the present study does not pertain to the evolution of the phases per se, but rather the effect of differences in mixing strength on approximately physical phase-space cuts. 5.3 Global gas evolution Fig. 8 shows the enriched fraction as a function of redshift, i.e. the ratio of the enriched gas mass to the total gas mass in each phase for two metallicity cuts: [Z] > −5.0 (left) and [Z] > −3.0 (right). We use [Z] ≡ log Z/Z⊙ as a proxy for metallicity, where Z is the mass fraction of metals in a gas particle, and Z⊙ = 0.02 is the solar metal mass fraction (Anders & Grevesse 1989). The rows represent the five phases defined in Table 4. Figure 8. View largeDownload slide The enriched fraction as a function of redshift, defined as the mass fraction of gas with Z > 10−5 Z⊙ (left) and Z > 10−3 Z⊙ (right) in each phase. The phases are, from top to bottom: ISM, cool CGM, HHG, DIFF, and the WHIM. Dynamic diffusion results in higher enriched fraction at the [Z] > −5 level, while maintaining a similar fraction to the no-mixing case above [Z] > −3. The constant-coefficient Smagorinsky model increases the enriched fractions at all metallicity cuts. Figure 8. View largeDownload slide The enriched fraction as a function of redshift, defined as the mass fraction of gas with Z > 10−5 Z⊙ (left) and Z > 10−3 Z⊙ (right) in each phase. The phases are, from top to bottom: ISM, cool CGM, HHG, DIFF, and the WHIM. Dynamic diffusion results in higher enriched fraction at the [Z] > −5 level, while maintaining a similar fraction to the no-mixing case above [Z] > −3. The constant-coefficient Smagorinsky model increases the enriched fractions at all metallicity cuts. We first focus on the [Z] > −5.0 cut and start by examining the ISM results in the top row of Fig. 8. In the None case, the fraction of gas enriched to [Z] > −5.0 exceeds 90 per cent at z = 4, reaches a peak at z = 0.5, and then very slightly downturns by z = 0 due to accretion of low-metallicity gas. Note that, in the following discussion, the differences between models are more important than the absolute values. In the S- models, gas is enriched much earlier and we see 95 per cent of gas above [Z] > −5.0 by z = 9. We do not observe the same slight downturn as in the None case. This is not surprising. The overdiffusive nature of the S- model leads to a reduced fraction of low-metallicity gas. By z = 0, 100 per cent of the gas is above [Z] > −5.0 in the S- models. Until z ∼ 1, the D- models show enrichment levels intermediate between the None and S- cases. At z < 1, the D-uz case continues this trend while the D-uvz model exhibits the biggest downturn. In the CGM, the second row in Fig. 8, we notice similar trends to the ISM, at enrichment levels [Z] > −5.0. For the None case, the gas is 30 per cent enriched at z = 385 per cent (the peak) at z = 0.5, followed by a slight downturn. The S- models follow the same qualitative trend but are at a constant 10 per cent above the None case, while the D- models remain in between the S- and None curves at all times. We do not discuss the details of the trends in the following three rows: the HHG, DIFF, and WHIM phases, respectively, but include them for completeness. In these gas phases, the differences between mixing models are qualitatively the same as the CGM, whereas the S- and D- cases show significant increased enrichment above [Z] > −5.0. Now we focus on the higher metallicity cut, [Z] > −3.0, in the ISM (top, right panel in Fig. 8). In the None case, gas is enriched over 90 per cent starting at redshift z = 2, reaches a peak at z ≈ 0.5, and turns down by z = 0. The trend for the [Z] > −3.0 gas is similar to that for the [Z] > −5.0 gas. The normalization of the curve is lower than at the [Z] > −5.0 threshold, as expected, since gas is enriched at higher metallicities later in cosmic evolution. In the S- cases, gas is enriched above 90 per cent (for the [Z] > −3.0 cut) earlier compared to the None case, starting at z = 3, and remains above the None case at all times. The S- cases show a similar downturn to the None case near z ≈ 0.5. Enrichment levels in the D-uz case remain slightly above the None case, but in D-uvz, there is a sharp downturn at z ≈ 0.75 and the final enrichment level is below the None case. The enrichment downturn in all of the mixing models is due a fresh supply of lower metallicity gas entering the medium, and the D-uvz model appears to amplify this effect. We now discuss the enrichment levels above [Z] > −3.0 for the CGM in the right column of Fig. 8, in the second row from the top. In the None case, the gas is enriched above 30 per cent by z ≈ 3, rises to a maximum of 80 per cent at z ≈ 0.4, and slightly decreases to 75 per cent at z = 0. Compared to the None case at the [Z] > −5.0 metallicity cut, the qualitative trend remains the same while the normalization has decreased, for the same reason we describe above. In the S- cases, the gas reaches 30 per cent enrichment levels somewhat earlier, by z ≈ 3.5, and reaches a maximum enrichment level of 90 per cent at z ≈ 0.4. The maximum enrichment is 10 per cent higher in the S- cases compared to the None case showing that the constant-coefficient Smagorinsky model affects higher and lower metallicities equally. The enrichment levels of the gas in the D- cases closely follow the None case for metallicities [Z] > −3.0, with only a slight divergence at z = 0. The HHG, DIFF, and WHIM phases follow the same qualitative trends as the CGM and, as noted above, we do not examine them in detail here. However, these phases show that the dynamic model increases the gas enrichment levels above [Z] > −3.0 more so than in the CGM, but only slightly. As in the isolated galaxy experiment (see Section 4), the diffusivity plays a significant role in the enrichment levels of cosmic gas on a global scale. The overdiffusive character of the constant-coefficient Smagorinsky model consistently results in the highest gas enrichment levels throughout the evolution of the simulations. Interestingly, these same results also indicate that the dynamic model has the most effect on the metallicities in the range −5.0 < [Z] < −3.0 while leaving those above [Z] > −3.0 near the no-mixing level. Early star formation and SMBH formation are very sensitive to metallicities in this range (see Section 1), and we have demonstrated that the dynamic model maximally affects those metallicities. Although we do not resolve early star formation or include SMBHs in our simulations, our results show that the dynamic model ought to be investigated further in cosmological simulations. 5.4 Global gas-phase metallicity MDFs provide additional insight into the spatial redistribution of metals when compared with global enriched fractions, and we therefore investigate the MDFs of galactic gas at three separate evolutionary times: z = 2, 1, 0. Our MDFs are probability density functions (dn/d[Z], where n is the number of particles) and are normalized such that the area under each curve is unity. To facilitate plotting, we set a metallicity minimum of [Z] = −10 for all gas particles. Fig. 9 shows the MDFs in the ISM, CGM, and HHG in the columns, left to right, and at z = 2, 1, and 0, from top to bottom. The narrow, leftmost spike in the panels corresponds to the metallicity minimum. Figure 9. View largeDownload slide The MDFs across three global gas phases (columns) in the simulation volume, at three separate evolutionary times (rows from top to bottom): z = 2, 1, and 0. We set a minimum metallicity of [Z] = −10 in our simulations, in order to show the abundance of unenriched gas. The overall trends in each phase do not change significantly over time, yet comparing between mixing models reveal slight differences in the distributions of gas-phase metals. Without metal mixing (None), galactic winds build MDFs in the non-star-forming gas that are similar to the distribution of the ISM. Turbulent mixing allows for further redistribution of metals once they arrive in the exterior media, altering the original distribution. The diffusivity has a larger role in this method rather than the quantities diffused. Figure 9. View largeDownload slide The MDFs across three global gas phases (columns) in the simulation volume, at three separate evolutionary times (rows from top to bottom): z = 2, 1, and 0. We set a minimum metallicity of [Z] = −10 in our simulations, in order to show the abundance of unenriched gas. The overall trends in each phase do not change significantly over time, yet comparing between mixing models reveal slight differences in the distributions of gas-phase metals. Without metal mixing (None), galactic winds build MDFs in the non-star-forming gas that are similar to the distribution of the ISM. Turbulent mixing allows for further redistribution of metals once they arrive in the exterior media, altering the original distribution. The diffusivity has a larger role in this method rather than the quantities diffused. First, as we mention at the end of Section 5.2, it is important to note that the differences in the MDFs are not solely due to turbulent mixing. There are transitions that occur across the sharp boundaries of the phase cuts we use. In all of our mixing models, the cuts mostly affect the gas that falls within the CGM region of phase space and, in examining its spatial distribution at z ≳ 1, we find that the majority of the low-metallicity ([Z] < −3) gas in the CGM is in the cosmic filaments. This gas is not spatially associated with dark matter haloes, but has the correct density and temperature to belong to the CGM phase. Conversely, the gas in the ISM and HHG phase-space regions do correspond to what we consider their spatial counterparts, and the differences in their MDFs across the evolution of the simulation are dominated by turbulent mixing. Turning to the distributions in Fig. 9, we first investigate the ISM. The distributions across all models appear qualitatively similar, with slight differences at lower metallicities. At z = 2, there is a peak at [Z] ≈ −1 in all models, and in the None case, there is also near-pristine gas in the ISM. The D- cases show a slightly extended tail covering the range −6 ≲ [Z] ≲ −4, compared to the S- and None cases, corresponding to the slight enrichment of the low-metallicity gas in the None case. In the S- cases, the enrichment process is more efficient. By z = 1, the tail tightens and all the distributions have negligible differences. At z = 0, however, the None and D- cases share the same distribution whereas the S- models show a tighter distribution. Additionally, there is much more gas at the minimum metallicity in the None case compared to both the D- and S- cases because metal enrichment of pristine particles in the None case only occurs when they are spatially adjacent to the stellar feedback sources. In the D- and S- cases, any particle that is enriched acts as a local source of metals for neighbouring pristine particles, driving down the number of particles at the metallicity minimum. The trends here are similar to, and caused by, the same effect we see in the ISM of the isolated galaxy in Section 4, where the spread in the ISM MDF strongly depends on the diffusivity. Next, we examine the MDFs of the CGM. In all models at z = 2, the gas in the CGM is a combination of the dense cores of the cosmic filaments and the cool, dense gas within dark matter haloes. In the None case, the MDF resembles the ISM distribution albeit with a more extended tail towards lower metallicities and more gas at the minimum metallicity. None the less, most of the enriched gas is above [Z] ≥ −6.0. The CGM distribution in the None case builds from the galactic winds transporting metals into the medium, sampling the underlying ISM distribution. Additionally, SNIa and AGB stars contribute to varying metal distributions via the kernel-weighting procedure as gas spatially adjacent to the ISM is enriched at low levels. In the None case, the gas that is spatially in the cosmic filaments is at the metallicity minimum. The z = 2 D- models show a bi-modality with peaks at [Z] ≈ −1.0 and [Z] ≈ −5.5, and in the S- models we also see a bi-modality but with peaks at [Z] ≈ −1.0 and [Z] ≈ −3.0. The higher metallicity peaks in the D- and S- cases correspond to the peak metallicity in the ISM whereas the lower peak in each case is due to the gas that is classified as the CGM, yet is spatially in the cosmic filaments. The spatial location of the gas does not change the effect of varying turbulent mixing strength; at redshift z = 2, the increased diffusivity in the S- models leads to a 2 - 3 order of magnitude shift compared to the D- models in the secondary low-metallicity peak, indicating that the metals are much more dispersed in the S- cases. Comparing to the None case, the question arises as to why the D- and S- cases have a broad lower metallicity peak. At high redshift (z ≳ 2), in the D- and S- cases, the metal-enriched galactic winds escape the galactic haloes and under the action of turbulent mixing, contaminate the gas in the filaments. This is evidenced by the lower amplitude spike at [Z] = −10 in these models. This does not occur in the None case, because the particles are unable to exchange metals directly. By z = 1 the fraction of CGM gas in the filaments has dropped as the dense filament cores are heated and enter the WHIM. However, there is still a small fraction of gas associated with the filaments. The None case shows a qualitatively similar distribution compared to z = 2 because the gas leaving the cosmic filaments is at the metallicity minimum. In the D- cases, the previous peak at [Z] ≈ −5.5 becomes a broad shelf between −5.0 < [Z] < −2.0. In the S- cases there is also a shelf of gas at [Z] ≈ −2.0 but the main peak dominates. The strong bi-modality from z = 2 has disappeared by z = 1 in the D- and S- models partly because the filamentary structure is increasingly classified as the WHIM and partly because the metallicity of the gas particles continues increasing due to mixing. We turn now to the MDF of the CGM at z = 0. At this redshift, most of the gas in the CGM region of phase space is associated with dark matter haloes. The None case has a similar distribution to z = 1 and z = 2. The D- MDFs have narrowed further, with the low extended metallicity shelf at z = 1 transforming into a tail that extends from [Z] ≈ −7.5 to [Z] ≈ −3.0, and a small shelf at [Z] ≈ −3.5. Above [Z] > −2.0, the D- MDFs coincide with the None case. In the S- cases, the distribution also tightens and the tail is approximately a power law from [Z] ≈ −6.0 to [Z] ≈ −1.0, where the latter value is the peak of the distribution. The variations in the MDFs between the None case and both the D- and S- cases represent enriched fractions that have been altered via turbulent mixing rather than the aforementioned transitions between phases. Shifting focus to the HHG at z = 2, we find that the gas in this phase is spatially associated with the dark matter haloes. We emphasize again that in the None case it is not possible for gas particles to exchange metals directly and consequently, the metal distributions only change via direct enrichment, or if the gas from the ISM reaches the phase under consideration via winds or via gas cooling from the intergalactic medium. We see that the None case shows a similar distribution to both the ISM and the CGM, albeit with a slightly broader low-metallicity tail. Like the z = 2 CGM distributions, the D- MDFs for the HHG show a bi-modality with peaks at [Z] ≈ −5.0 and [Z] ≈ −1.0 and a valley between −4.0 < [Z] < −2.0. Unlike the CGM MDF, the HHG low-metallicity peak is broader and it is not due to gas transitioning from the HHG phase, but rather it is due to inflowing, low-metallicity gas mixing with the already enriched gas in the haloes. The S- cases share the peak at [Z] ≈ −1.0, but the distribution is flat between −6.0 < [Z] < −2.0, before dropping-off towards low metallicities, in lockstep with the D- results. The gap that is apparent in the D- cases has disappeared, indicating that in the S- cases there is much more gas with metallicities in the range −4.0 < [Z] < −2.0 than in the D- cases, which is not surprising given the overdiffusive nature of the constant-coefficient Smagorinsky model. At z = 1 in the HHG, the None case remains unchanged except for a slight increase in the amount of gas near [Z] ≈ −3.0, and the tail has extended slightly towards lower metallicity. The latter is due to less enriched gas from the WHIM and DIFF phases accreting on to haloes, diluting the metal distribution, and this dilution continues through to z = 0. In the D- and S- cases, the low-metallicity feature (peak/shelf and extended tail) has shifted slightly towards higher metallicity as turbulent mixing redistributes the metals from the highly enriched particles. Next, we examine the HHG at z = 0. The enrichment level, in the None case, is decreasing as evidenced by the tail of the MDF due to less enriched gas from the WHIM and DIFF phases accreting on to haloes, diluting the metal distribution. In the D- and S- cases, there are coincident peaks in the MDFs at [Z] ≈ −6.0, although the D- cases show more gas at lower metallicities. Comparing across phases, the MDFs of the CGM gas are more sensitive to non-zero turbulent mixing strength than the ISM or the HHG. For instance, even though the bi-modalities at z = 2 in the CGM MDFs, in the S- and D- cases, disappear by z = 0 due to the CGM to WHIM transition, the dynamic model results in a residual extended tail in the CGM MDF as metals mix throughout the medium, whereas the constant-coefficient Smagorinsky model tightens the distribution in this phase and the gas metallicities rapidly approach the mean. The peak in the CGM moves towards higher metallicities as gas flows between phases (CGM to WHIM). The consequences of not including sub-grid turbulent mixing, regardless of the diffusivity, are clear – complex structure in the MDFs is highly dependent on the mixing strength. This is in contrast to Su et al. (2017) and Escala et al. (2018), who found that turbulent metal mixing strength had low-level effects in their simulated ISM. We posit that the low-level effects were due the authours** use of a constant value of Cs, rather than localizing the mixing strength to the appropriate regions. We do not see the aforementioned trends in the None case because it is not possible for gas particles to exchange metals in our simulations. Our stellar feedback model drives a decoupled wind from the ISM and that wind samples the MDF in the ISM, building up a similar metal distribution in the other phases that cannot change over time, unless material flows between the phases or due to delayed SNIa and AGB stars, via the kernel-weighted enrichment procedure. One important difference between the model we present and realistic environments is that our winds do not mix as they free-stream out of the galaxies. While the coupling and mixing strength of galactic winds is uncertain, the wind could redistribute thermal energy, momentum, and metals internally as a cohesive unit, even if they do not couple strongly to the surrounding medium (Huang & Katz, private communication). In our model, once the winds reach the criteria for recoupling in our simulations, they are free to mix their fluid properties and this, subsequently, allows for diverse MDFs in the gas phases exterior to the ISM. We do not investigate the details of individual galaxies here but note that the differences outlined above will impact the gas-phase and stellar metallicities of those systems. Therefore, including the dynamic model with a more accurate estimation of Cs is necessary, moving forward, in order to capture the physical redistribution of metals in galaxy evolution. 6 CONCLUSIONS All hydrodynamical methods that are used to investigate galaxy evolution, whether Lagrangian or Eulerian, require additional sub-grid thermal energy and momentum diffusion terms in order to account for sub-grid turbulence. In Lagrangian methods, such as MFM and SPH, metal diffusion is also required due to the inability of fluid elements to exchange metals by construction. Most implementations use the constant-coefficient Smagorinsky model – one that has been shown to be overdiffusive in almost all cases, especially laminar shear flows. We implemented and investigated the impact of the localized dynamic Smagorinsky model on global gas-phase properties in a series of numerical experiments using the gizmo code. In the dynamic case, the model coefficient depends on the local turbulent flow conditions, hence on the spatio-temporal coordinates. This is in contrast to the constant-coefficient Smagorinsky model where diffusivities depend directly on the magnitude of the velocity shear in the fluid. Compared to the constant-coefficient Smagorinsky model, the dynamic model has been shown to produce more accurate representations of fluid mechanical experiments (Kirkpatrick et al. 2006; Kleissl et al. 2006; Khani & Waite 2015; Benhamadouche et al. 2017; Lee & Cant 2017; Kara & Çağlar 2018; Taghinia et al. 2018). While we focused on cosmological experiments, the dynamic model has applications to any numerical experiment involving turbulent astrophysical flows, including stellar interiors, planetary formation, and star formation. Moreover, the method we describe in this paper, following Germano et al. (1991) and Piomelli & Liu (1995), is general and not only limited to Lagrangian hydrodynamics but also applicable to the Eulerian cases (see Schmidt 2015). For the MFM method, we showed that the dynamic model improves the density contrast in subsonic turbulence, allowing higher and lower density regions at fixed mass resolution. In an idealized Keplerian disc, an example of a laminar shear flow where the Smagorinsky model is known to be overdiffusive from basic analytic arguments, the dynamic model produced near-zero values of turbulent diffusivity. When we included thermal energy and momentum diffusion, the lower diffusivities prevented the rapid break-up of the disc due to excessive angular momentum transport. We observed similar minimized diffusivities in a Kelvin–Helmholtz instability experiment, where the constant-coefficient Smagorinsky model smoothed, and rapidly diffused, our metal tracer, whereas the dynamic model captured the fine level of mixing at the interface of the two fluids. We also investigated an isolated, Milky Way-like galaxy in order to test the dynamic model in a more complex, but still controlled, environment. The dynamic model in combination with momentum diffusion improved the stability of the gaseous disc compared to the constant-coefficient Smagorinsky model, and affected the spatial metal distributions, as indicated by the MDFs, as shown in Fig. 7. Rapid star formation early in the evolution of the disc leads to a higher diffusivity of thermal energy, momentum, and metals, and the subsequent exponential decay of star formation lowers the diffusivity in the dynamic model significantly. This results in a broader MDF in the ISM in the dynamic case, pointing towards less mixing in the ISM. When using the constant-coefficient Smagorinsky model, the diffusivity remained high throughout the evolution of the disc because of its strong dependence on the fluid velocity shear. We found similar variations in the CGM for both the dynamic and constant-coefficient Smagorinsky models that we attribute to the turbulence generated from stellar feedback increasing the diffusivity in both cases. We also examined the global gas enrichment fractions in a set of cosmological simulations. Global gas enrichment fractions are important for the formation of Population III stars and SMBHs because they are theorized to be sensitive to the metal content in the gas out of which they form (see Section 1 for more details). We found that the dynamic model lowers overall enrichment compared to the standard Smagorinsky model, and that it maximally impacts metallicities in the range −5.0 < [Z] < −3.0. This is precisely the metallicity regime that constrains the formation sites of SMBHs and Population III stars (Volonteri 2010; Sarmento et al. 2016). Specifically, the dynamic model increases the amount of gas above [Z] > −5.0 while maintaining the same enriched fraction of gas above [Z] > −3.0, compared to the no-mixing case. The standard Smagorinsky model increased the enriched fraction at all metal thresholds and in all gas phases. In our cosmological simulations without turbulent mixing, we found that each gas phase external to the ISM has a qualitatively similar MDF to the ISM itself. Turbulent mixing allows for regions to mix their metals, and introduces additional structure in MDFs of each phase. We found that the diffusivity had a significant impact on the MDFs of the ISM and CGM – the dynamic model shows broader MDFs in both phases at z = 0. In these regions, we found a bi-modality in the CGM at z = 2 which disappeared by z = 0 in both cases, yet more lower metallicity gas remained in the dynamic case. Our broad density and temperature phase-space criteria led to the bi-modality, as we found low-temperature and dense gas in the cores of the cosmic filaments at z ∼ 2. These spatial regions were of lower metallicity, and eventually return to the WHIM phase by z = 0. The peaks of the bi-modality, however, depend on the diffusivity: the dynamic model produced more metal poor (by several orders of magnitude) gas than the constant-coefficient Smagorinsky model. Finally, we briefly touch on our conclusions for SPH. Most authors apply the constant-coefficient Smagorinsky model to SPH (Wadsley et al. 2008; Shen et al. 2010, 2013; Williamson et al. 2016; Tremmel et al. 2017; Wadsley et al. 2017) and only include thermal energy and metal mixing. In reality, there are additional turbulent transport terms that are unique to SPH (see Di Mascio et al. 2017 for an introduction and derivation) that must be included.16 Introducing momentum diffusion (e.g. via turbulent mixing as in D-uvz and S-uvz) in SPH is problematic because of unknown interactions with artificial viscosity. Additionally, by construction, the smoothing kernel in SPH acts to produce coherent flows rather than fine structure observed in mesh-free or grid methods. When we introduced momentum diffusion into SPH, the results from all of the experiments in this study were amplified when compared to the MFM method, but the qualitative trends remained. Specifically, in our cosmological experiments, we found that momentum diffusion with the constant-coefficient model causes a delay in the formation of the ISM by ∼1 Gyr, compared to the None case. The dynamic model reduces the delay, but not to the no-mixing case. We attribute this to momentum diffusion and dissipation causing the Jeans mass to increase resulting in the damping of mass fluctuations. While we note that the dynamic turbulent mixing model introduced here is a step forward in understanding the redistribution of fluid properties in Lagrangian codes, there are caveats that must be explored. The dynamic model predicts the correct behaviour in supersonic flows, but we did not include compressive mixing terms into our equations of motion. It may be that compressive sub-grid mixing models further improve supersonic turbulence in the MFM method. Also, we justified using the Smagorinsky model by assuming that the local equilibrium condition holds, where the kinetic energy transfer rate down the turbulent cascade is equal on all scales. While the assumption is approximately true on average in the regimes we investigated (Schmidt et al. 2016), a fully consistent turbulence model involves tracking the sub-grid kinetic energy via an additional transport equation that includes all of the necessary, higher order, sub-grid scale terms (Schmidt 2015). However, the dynamic model mitigates the issue by inherently calculating the deviations from local equilibrium. Furthermore, the approximations for filtering the fluid fields require care and attention. While Monaghan’s filtering approximation (equation 19) holds on the singly filtered quantities for variations on scales larger than the resolution scale h, doubly filtered quantities may be over- or undersmoothed. A more robust, efficient, filtering procedure will need to be derived specifically for Lagrangian methods in the highly compressible case, for filtering scales larger than h. In summary, the dynamic Smagorinsky model localizes the strength of turbulent mixing to only turbulent regions of the flow. This provides a turbulent mixing model that does not rely on pre-calibrated parameters – therefore simultaneously allowing near-zero diffusion in laminar shear flows and the expected diffusion in turbulent flows. The physical experiments to which we subjected the model show that the dynamic model significantly alters the MDFs of the ISM and CGM in a global sense. In future work we will examine the extent of small-scale differences associated with dynamic diffusion, and its impact on galaxy properties. ACKNOWLEDGEMENTS This research was enabled in part by support provided by WestGrid and Compute/Calcul Canada. DR and AB acknowledge support from NSERC (Canada) through the Discovery Grant program. DR thanks the organizers of the Computing the Universe: At the Intersection of Computer Science and Cosmology conference in Oaxaca, Mexico for an invited talk, and also James Wadsley and Andrey Kravtsov for their recommendations at the conference that led to the further refinement of this research project. DR also thanks Valentin Perret for the DICE code, and Fabrice Durier, Ondrea Clarkson, Austin Davis, and Maan Hani for many useful discussions during the course of this research. Support for PFH was provided by an Alfred P. Sloan Research Fellowship, NSF Collaborative Research Grant #1715847 and CAREER grant #1455342, and NASA grants NNX15AT06G, JPL 1589742, 17-ATP17-0214. We would also especially like to thank our referee, Wolfram Schmidt, for his contributions in improving the final version of this study. Footnotes 1 Defined as the ratio of inertial forces to viscous (dissipative) forces. 2 In compressible turbulence, the scaling is much steeper (see Kritsuk et al. 2007; Federrath 2013). 3 The two assumptions are often combined together and referred to as the scale-similarity hypothesis. 4 We refer to any fluid element as a particle, for simplicity. Fluid elements in the MFM method are not particles in an SPH sense, and are defined by the effective geometrical faces moving along lines connecting cells enclosing a finite mass. 5 Riemann solver dissipation arises from the high, even-order, truncated terms in the Taylor series expansion of the basic variables in the conservation equations. 6 We follow this Einstein notation throughout this paper. 7 Our filtering implementation naturally density-weights quantities because we follow the hydrodynamical weighting scheme. 8 Taken as |$\overline{h}_{ab} = \frac{1}{2}(\overline{h}_a + \overline{h}_b)$| in order to equally weight each smoothing scale. 9 The harmonic mean is of the form |$\langle \rho _{ab}\rangle _{\overline{h}} = 2\overline{\rho }_a\overline{\rho }_b / (\overline{\rho }_a + \overline{\rho }_b)$| and weights towards the lowest value. This allows high-density particles to have a fair contribution to the differences within the kernel. 10 Units are arbitrary code units. 11 https://github.com/ruggiero/galstep 12 https://bitbucket.org/vperret/dice 13 The entropy scales as S ∝ ln(P/ργ), and therefore our measure is off by a multiplicative constant. 14 For several examples see Naab & Ostriker (2017). 15 http://web.mit.edu/edbert/ 16 It is possible to apply the dynamic model to the transport terms in Di Mascio et al. (2017), further improving upon their work. REFERENCES Anders E. , Grevesse N. , 1989 , Geochim. Cosmochim. Acta , 53 , 197 https://doi.org/10.1016/0016-7037(89)90286-X Crossref Search ADS Balogh M. L. , Babul A. , Patton D. R. , 1999 , MNRAS , 307 , 463 https://doi.org/10.1046/j.1365-8711.1999.02608.x Crossref Search ADS Bauer A. , Springel V. , 2012 , MNRAS , 423 , 2558 https://doi.org/10.1111/j.1365-2966.2012.21058.x Crossref Search ADS Benhamadouche S. , Arenas M. , Malouf W. J. , 2017 , Nucl. Eng. Des. , 312 , 128 https://doi.org/10.1016/j.nucengdes.2016.09.010 Crossref Search ADS Bertschinger E. , 2011 , Astrophysics Source Code Library . ascl:1106.008 Brook C. B. , Stinson G. , Gibson B. K. , Shen S. , Macciò A. V. , Obreja A. , Wadsley J. , Quinn T. , 2014 , MNRAS , 443 , 3809 https://doi.org/10.1093/mnras/stu1406 Crossref Search ADS Bruzual G. , Charlot S. , 2003 , MNRAS , 344 , 1000 https://doi.org/10.1046/j.1365-8711.2003.06897.x Crossref Search ADS Chabrier G. , 2003 , PASP , 115 , 763 https://doi.org/10.1086/376392 Crossref Search ADS Chandrasekhar S. , 1961 , Hydrodynamic and Hydromagnetic Stability . Oxford Univ. Press , London and New York Christensen C. R. , Davé R. , Governato F. , Pontzen A. , Brooks A. , Munshi F. , Quinn T. , Wadsley J. , 2016 , ApJ , 824 , 57 https://doi.org/10.3847/0004-637X/824/1/57 Crossref Search ADS Clark R. A. , Ferziger J. H. , Reynolds W. C. , 1979 , J. Fluid Mech. , 91 , 1 https://doi.org/10.1017/S002211207900001X Crossref Search ADS Colbrook M. J. , Ma X. , Hopkins P. F. , Squire J. , 2017 , MNRAS , 467 , 2421 https://doi.org/10.1093/mnras/stx261 Crossref Search ADS Crain R. A. , McCarthy I. G. , Schaye J. , Theuns T. , Frenk C. S. , 2013 , MNRAS , 432 , 3005 https://doi.org/10.1093/mnras/stt649 Crossref Search ADS Davé R. , Finlator K. , Oppenheimer B. D. , 2006 , MNRAS , 370 , 273 https://doi.org/10.1111/j.1365-2966.2006.10464.x Crossref Search ADS Davé R. , Oppenheimer B. D. , Katz N. , Kollmeier J. A. , Weinberg D. H. , 2010 , MNRAS , 408 , 2051 https://doi.org/10.1111/j.1365-2966.2010.17279.x Crossref Search ADS Davé R. , Thompson R. , Hopkins P. F. , 2016 , MNRAS , 462 , 3265 https://doi.org/10.1093/mnras/stw1862 Crossref Search ADS Davé R. , Rafieferantsoa M. H. , Thompson R. J. , Hopkins P. F. , 2017 , MNRAS , 467 , 115 https://doi.org/10.1093/mnras/stx108 Di Mascio A. , Antuono M. , Colagrossi A. , Marrone S. , 2017 , Phys. Fluids , 29 , 035102 https://doi.org/10.1063/1.4978274 Crossref Search ADS Elmegreen B. G. , Scalo J. , 2004 , ARA&A , 42 , 211 https://doi.org/10.1146/annurev.astro.41.011802.094859 Crossref Search ADS Emerick A. , Bryan G. L. , Mac Low M.-M. , Côté B. , Johnston K. V. , O’Shea B. W. , 2018 , ApJ , 869 , 94 Crossref Search ADS Escala I. et al. , 2018 , MNRAS , 474 , 2194 https://doi.org/10.1093/mnras/stx2858 Crossref Search ADS Evoli C. , Ferrara A. , 2011 , MNRAS , 413 , 2721 https://doi.org/10.1111/j.1365-2966.2011.18343.x Crossref Search ADS Faucher-Giguère C.-A. , Lidz A. , Zaldarriaga M. , Hernquist L. , 2009 , ApJ , 703 , 1416 https://doi.org/10.1088/0004-637X/703/2/1416 Crossref Search ADS Federrath C. , 2013 , MNRAS , 436 , 1245 https://doi.org/10.1093/mnras/stt1644 Crossref Search ADS Finlator K. , Davé R. , 2008 , MNRAS , 385 , 2181 https://doi.org/10.1111/j.1365-2966.2008.12991.x Crossref Search ADS Gaburov E. , Nitadori K. , 2011 , MNRAS , 414 , 129 https://doi.org/10.1111/j.1365-2966.2011.18313.x Crossref Search ADS Garnier E. , Adams N. , Sagaut P. , 2009 , Large Eddy Simulation for Compressible Flows . Scientific Computation, Springer Netherlands , Dordrecht https://doi.org/10.1007/978-90-481-2819-8 Genel S. , 2016 , ApJ , 822 , 107 https://doi.org/10.3847/0004-637X/822/2/107 Crossref Search ADS Genel S. et al. , 2014 , MNRAS , 445 , 175 https://doi.org/10.1093/mnras/stu1654 Crossref Search ADS Germano M. , Piomelli U. , Moin P. , Cabot W. H. , 1991 , Phys. Fluids A , 3 , 1760 https://doi.org/10.1063/1.857955 Crossref Search ADS Gingold R. A. , Monaghan J. J. , 1977 , MNRAS , 181 , 375 https://doi.org/10.1093/mnras/181.3.375 Crossref Search ADS Grand R. J. J. et al. , 2017 , MNRAS , 207 , 179 https://doi.org/10.1093/mnras/stx071 Greif T. H. , Glover S. C. O. , Bromm V. , Klessen R. S. , 2009 , MNRAS , 392 , 1381 https://doi.org/10.1111/j.1365-2966.2008.14169.x Crossref Search ADS Grete P. , Vlaykov D. G. , Schmidt W. , Schleicher D. R. , 2017 , Phys. Rev. E , 95 , 033206 https://doi.org/10.1103/PhysRevE.95.033206 Crossref Search ADS PubMed Grete P. , O’Shea B. W. , Beckwith K. , 2018 , ApJ , 858 , L19 https://doi.org/10.3847/2041-8213/aac0f5 Crossref Search ADS Guedes J. , Callegari S. , Madau P. , Mayer L. , 2011 , ApJ , 742 , 76 https://doi.org/10.1088/0004-637X/742/2/76 Crossref Search ADS Hani M. H. , Sparre M. , Ellison S. L. , Torrey P. , Vogelsberger M. , 2018 , MNRAS , 475 , 1160 https://doi.org/10.1093/mnras/stx3252 Crossref Search ADS Hernquist L. , Katz N. , 1989 , ApJS , 70 , 419 https://doi.org/10.1086/191344 Crossref Search ADS Hopkins P. F. , 2015 , MNRAS , 450 , 53 https://doi.org/10.1093/mnras/stv195 Crossref Search ADS Hopkins P. F. , 2017 , MNRAS , 466 , 3387 https://doi.org/10.1093/mnras/stw3306 Crossref Search ADS Hopkins P. F. , Kereš D. , Oñorbe J. , Faucher-Giguère C.-A. , Quataert E. , Murray N. , Bullock J. S. , 2014 , MNRAS , 445 , 581 https://doi.org/10.1093/mnras/stu1738 Crossref Search ADS Hosono N. , Saitoh T. R. , Makino J. , 2016 , ApJS , 224 , 32 https://doi.org/10.3847/0067-0049/224/2/32 Crossref Search ADS Iapichino L. , Viel M. , Borgani S. , 2013 , MNRAS , 432 , 2529 https://doi.org/10.1093/mnras/stt611 Crossref Search ADS Iwamoto K. , Brachwitz F. , Nomoto K. , Kishimoto N. , Umeda H. , Hix W. R. , Thielemann F.-K. , 1999 , ApJS , 125 , 439 https://doi.org/10.1086/313278 Crossref Search ADS Kara R. , Çağlar M. , 2018 , Appl. Math. Comput. , 322 , 89 https://doi.org/10.1016/j.amc.2017.11.033 Kennicutt R. C. , 1998 , ApJ , 498 , 541 https://doi.org/10.1086/305588 Crossref Search ADS Khani S. , Waite M. L. , 2015 , J. Fluid Mech. , 773 , 327 https://doi.org/10.1017/jfm.2015.249 Crossref Search ADS King A. , Pounds K. , 2015 , ARA&A , 53 , 115 https://doi.org/10.1146/annurev-astro-082214-122316 Crossref Search ADS Kirkpatrick M. P. , Ackerman A. S. , Stevens D. E. , Mansour N. N. , 2006 , J. Atmos. Sci. , 63 , 526 https://doi.org/10.1175/JAS3651.1 Crossref Search ADS Kleissl J. , Kumar V. , Meneveau C. , Parlange M. B. , 2006 , Water Resources Res. , 42 , 1 https://doi.org/10.1029/2005WR004685 Crossref Search ADS Kritsuk A. G. , Norman M. L. , Padoan P. , Wagner R. , 2007 , ApJ , 665 , 416 https://doi.org/10.1086/519443 Crossref Search ADS Krumholz M. R. , McKee C. F. , Tumlinson J. , 2009 , ApJ , 699 , 850 https://doi.org/10.1088/0004-637X/699/1/850 Crossref Search ADS Landau L. , Lifshitz E. , 1987 , Course of Theoretical Physics. Vol. 6: Fluid Mechanics . Pergamon Press , Oxford Lanson N. , Vila J.-P. , 2008a , SIAM J. Numer. Anal. , 46 , 1912 https://doi.org/10.1137/S0036142903427718 Crossref Search ADS Lanson N. , Vila J.-P. , 2008b , SIAM J. Numer. Anal. , 46 , 1935 https://doi.org/10.1137/S003614290444739X Crossref Search ADS Lee C. Y. , Cant S. , 2017 , Flow Turbul. Combust. , 98 , 155 https://doi.org/10.1007/s10494-016-9751-4 Crossref Search ADS PubMed Liang L. , Durier F. , Babul A. , Davé R. , Oppenheimer B. D. , Katz N. , Fardal M. , Quinn T. , 2016 , MNRAS , 456 , 4266 https://doi.org/10.1093/mnras/stv2840 Crossref Search ADS Lucy L. B. , 1977 , AJ , 82 , 1013 https://doi.org/10.1086/112164 Crossref Search ADS Meneveau C. , Katz J. , 2000 , Annu. Rev. Fluid Mech. , 32 , 1 https://doi.org/10.1146/annurev.fluid.32.1.1 Crossref Search ADS Mitra S. , Davé R. , Finlator K. , 2015 , MNRAS , 452 , 1184 https://doi.org/10.1093/mnras/stv1387 Crossref Search ADS Monaghan J. J. , 1989 , J. Computat. Phys. , 82 , 1 https://doi.org/10.1016/0021-9991(89)90032-6 Crossref Search ADS Monaghan J. J. , 2005 , Rep. Prog. Phys. , 68 , 1703 https://doi.org/10.1088/0034-4885/68/8/R01 Crossref Search ADS Monaghan J. J. , 2011 , Eur. J. Mech. B , 30 , 360 https://doi.org/10.1016/j.euromechflu.2011.04.002 Crossref Search ADS Muratov A. L. , Kereš D. , Faucher-Giguère C.-A. , Hopkins P. F. , Quataert E. , Murray N. , 2015 , MNRAS , 454 , 2691 https://doi.org/10.1093/mnras/stv2126 Crossref Search ADS Naab T. , Ostriker J. P. , 2017 , ARA&A , 55 , 59 https://doi.org/10.1146/annurev-astro-081913-040019 Crossref Search ADS Nomoto K. , Tominaga N. , Umeda H. , Kobayashi C. , Maeda K. , 2006 , Nucl. Phys. A , 777 , 424 https://doi.org/10.1016/j.nuclphysa.2006.05.008 Crossref Search ADS Oppenheimer B. D. , Davé R. , 2006 , MNRAS , 373 , 1265 https://doi.org/10.1111/j.1365-2966.2006.10989.x Crossref Search ADS Oppenheimer B. D. , Davé R. , 2008 , MNRAS , 387 , 577 https://doi.org/10.1111/j.1365-2966.2008.13280.x Crossref Search ADS Oppenheimer B. D. , Schaye J. , 2013 , MNRAS , 434 , 1043 https://doi.org/10.1093/mnras/stt1043 Crossref Search ADS Oppenheimer B. D. et al. , 2016 , MNRAS , 460 , 2157 https://doi.org/10.1093/mnras/stw1066 Crossref Search ADS Pan L. , Scannapieco E. , Scalo J. , 2013 , ApJ , 775 , 111 https://doi.org/10.1088/0004-637X/775/2/111 Crossref Search ADS Perret V. , Renaud F. , Epinat B. , Amram P. , Bournaud F. , Contini T. , Teyssier R. , Lambert J.-C. , 2014 , A&A , 562 , A1 https://doi.org/10.1051/0004-6361/201322395 Crossref Search ADS Piomelli U. , Liu J. , 1995 , Phys. Fluids , 7 , 839 https://doi.org/10.1063/1.868607 Crossref Search ADS Piomelli U. , Cabot W. H. , Moin P. , Lee S. , 1991 , Phys. Fluids A , 3 , 1766 https://doi.org/10.1063/1.857956 Crossref Search ADS Planck Collaboration XIII , 2015 , A&A , 594 , A13 https://doi.org/10.1051/0004-6361/201525830 Pope S. , 2000 , Turbulent Flows. Cambridge Univ. Press , Cambridge , UK Ruggiero R. , Lima Neto G. B. , 2017 , MNRAS , 468 , 4107 https://doi.org/10.1093/mnras/stx744 Crossref Search ADS Sagaut P. , 2006 , Large Eddy Simulation for Incompressible Flows . Scientific Computation, Springer-Verlag , Berlin/Heidelberg Salpeter E. E. , 1955 , ApJ , 121 , 161 https://doi.org/10.1086/145971 Crossref Search ADS Sarmento R. , Scannapieco E. , Pan L. , 2016 , ApJ , 834 , 23 https://doi.org/10.3847/1538-4357/834/1/23 Crossref Search ADS Scannapieco E. , Bildsten L. , 2005 , ApJ , 629 , L85 https://doi.org/10.1086/452632 Crossref Search ADS Scannapieco E. , Brüggen M. , 2008 , ApJ , 686 , 927 https://doi.org/10.1086/591228 Crossref Search ADS Schaye J. et al. , 2010 , MNRAS , 402 , 1536 https://doi.org/10.1111/j.1365-2966.2009.16029.x Crossref Search ADS Schaye J. et al. , 2015 , MNRAS , 446 , 521 https://doi.org/10.1093/mnras/stu2058 Crossref Search ADS Schmidt W. , 2015 , Liv. Rev. Comput. Astrophys. , 1 , 64 https://doi.org/10.1007/lrca-2015-2 Schmidt W. , Federrath C. , 2011 , A&A , 528 , A106 https://doi.org/10.1051/0004-6361/201015630 Crossref Search ADS Schmidt W. , Niemeyer J. C. , Hillebrandt W. , 2006 , A&A , 450 , 265 https://doi.org/10.1051/0004-6361:20053617 Crossref Search ADS Schmidt W. et al. , 2014 , MNRAS , 440 , 3051 https://doi.org/10.1093/mnras/stu501 Crossref Search ADS Schmidt W. , Engels J. F. , Niemeyer J. C. , Almgren A. S. , 2016 , MNRAS , 459 , 701 https://doi.org/10.1093/mnras/stw632 Crossref Search ADS Semenov V. A. , Kravtsov A. V. , Gnedin N. Y. , 2016 , ApJ , 826 , 200 https://doi.org/10.3847/0004-637X/826/2/200 Crossref Search ADS Shen S. , Wadsley J. , Stinson G. , 2010 , MNRAS , 407 , 1581 https://doi.org/10.1111/j.1365-2966.2010.17047.x Crossref Search ADS Shen S. , Madau P. , Aguirre A. , Guedes J. , Mayer L. , Wadsley J. , 2012 , ApJ , 760 , 50 https://doi.org/10.1088/0004-637X/760/1/50 Crossref Search ADS Shen S. , Madau P. , Guedes J. , Mayer L. , Prochaska J. X. , Wadsley J. , 2013 , ApJ , 765 , 89 https://doi.org/10.1088/0004-637X/765/2/89 Crossref Search ADS Smagorinsky J. , 1963 , Mon. Weather Rev. , 91 , 99 https://doi.org/10.1175/1520-0493(1963)091 < 0099:GCEWTP>2.3.CO;2 Crossref Search ADS Smith B. D. et al. , 2017 , MNRAS , 466 , 2217 https://doi.org/10.1093/mnras/stw3291 Crossref Search ADS Sokołowska A. , Mayer L. , Babul A. , Madau P. , Shen S. , 2016 , ApJ , 819 , 21 https://doi.org/10.3847/0004-637X/819/1/21 Crossref Search ADS Sokołowska A. , Babul A. , Mayer L. , Shen S. , Madau P. , 2018 , ApJ , 867 , 73 Crossref Search ADS Somerville R. S. , Davé R. , 2015 , ARA&A , 53 , 51 https://doi.org/10.1146/annurev-astro-082812-140951 Crossref Search ADS Springel V. , 2005 , MNRAS , 364 , 1105 https://doi.org/10.1111/j.1365-2966.2005.09655.x Crossref Search ADS Springel V. , 2010 , MNRAS , 401 , 791 https://doi.org/10.1111/j.1365-2966.2009.15715.x Crossref Search ADS Springel V. , Hernquist L. , 2003a , MNRAS , 339 , 289 https://doi.org/10.1046/j.1365-8711.2003.06206.x Crossref Search ADS Springel V. , Hernquist L. , 2003b , MNRAS , 339 , 312 https://doi.org/10.1046/j.1365-8711.2003.06207.x Crossref Search ADS Springel V. , Di Matteo T. , Hernquist L. , 2005 , MNRAS , 361 , 776 https://doi.org/10.1111/j.1365-2966.2005.09238.x Crossref Search ADS Spyropoulos E. T. , Blaisdell G. A. , 1996 , AIAA J. , 34 , 990 https://doi.org/10.2514/3.13178 Crossref Search ADS Su K.-Y. , Hopkins P. F. , Hayward C. C. , Faucher-Giguère C.-A. , Kereš D. , Ma X. , Robles V. H. , 2017 , MNRAS , 471 , 144 https://doi.org/10.1093/mnras/stx1463 Crossref Search ADS Suresh J. , Rubin K. H. , Kannan R. , Werk J. K. , Hernquist L. , Vogelsberger M. , 2017 , MNRAS , 465 , 2966 https://doi.org/10.1093/mnras/stw2499 Crossref Search ADS Taghinia J. H. , Rahman M. M. , Lu X. , 2018 , Energy Build. , 170 , 47 https://doi.org/10.1016/j.enbuild.2018.03.075 Crossref Search ADS Teyssier R. , Moore B. , Martizzi D. , Dubois Y. , Mayer L. , 2011 , MNRAS , 414 , 195 https://doi.org/10.1111/j.1365-2966.2011.18399.x Crossref Search ADS Torrey P. et al. , 2017 , preprint (arXiv:1711.05261) Tremmel M. et al. , 2018 , preprint (arXiv:1806.01282) Tremmel M. , Karcher M. , Governato F. , Volonteri M. , Quinn T. R. , Pontzen A. , Anderson L. , Bellovary J. , 2017 , MNRAS , 470 , 1121 https://doi.org/10.1093/mnras/stx1160 Crossref Search ADS Tumlinson J. , Peeples M. S. , Werk J. K. , 2017 , ARA&A , 55 , 389 https://doi.org/10.1146/annurev-astro-091916-055240 Crossref Search ADS Urzay J. , O’Brien J. , Ihme M. , Moin P. , Saghafian a. , 2013 , Center for Turbulence Research , Annual Research Briefs , Stanford , p. 123 van de Voort F. , Schaye J. , Altay G. , Theuns T. , 2012 , MNRAS , 421 , 2809 https://doi.org/10.1111/j.1365-2966.2012.20487.x Crossref Search ADS Vogelsberger M. et al. , 2014 , MNRAS , 444 , 1518 https://doi.org/10.1093/mnras/stu1536 Crossref Search ADS Voit G. M. , 2005 , Adv. Space Res. , 36 , 701 https://doi.org/10.1016/j.asr.2005.02.042 Crossref Search ADS Volonteri M. , 2010 , A&AR , 18 , 279 https://doi.org/10.1007/s00159-010-0029-x Crossref Search ADS Vreman A. W. , 2004 , Phys. Fluids , 16 , 3670 https://doi.org/10.1063/1.1785131 Crossref Search ADS Wadsley J. W. , Veeravalli G. , Couchman H. M. P. , 2008 , MNRAS , 387 , 427 https://doi.org/10.1111/j.1365-2966.2008.13260.x Crossref Search ADS Wadsley J. W. , Keller B. W. , Quinn T. R. , 2017 , MNRAS , 471 , 2357 https://doi.org/10.1093/mnras/stx1643 Crossref Search ADS Williamson D. , Martel H. , Kawata D. , 2016 , ApJ , 822 , 91 https://doi.org/10.3847/0004-637X/822/2/91 Crossref Search ADS © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Dynamic localized turbulent diffusion and its impact on the galactic ecosystem
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/sty3376
DA - 2019-03-01
UR - https://www.deepdyve.com/lp/oxford-university-press/dynamic-localized-turbulent-diffusion-and-its-impact-on-the-galactic-THSsrK0ZU2
SP - 3810
VL - 483
IS - 3
DP - DeepDyve
ER -