# Mixing matters

Douglas Rennehan1
1Department of Physics & Astronomy, University of Victoria, BC V8P 5C2, Canada
E-mail<EMAIL_ADDRESS>(DR)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

All hydrodynamical simulations of turbulent astrophysical phenomena require
sub-grid scale models to properly treat energy dissipation and metal mixing.
We present the first implementation and application of an anisotropic eddy
viscosity and metal mixing model in Lagrangian astrophysical simulations,
including a dynamic procedure for the model parameter. We compare these two
models directly to the common Smagorinsky and dynamic variant. Using the mesh-
free finite mass method as an example, we show that the anisotropic model is
best able to reproduce the proper Kolmogorov inertial range scaling in
homogeneous, isotropic turbulence. Additionally, we provide a method to
calibrate the metal mixing rate that ensures numerical convergence. In our
first application to cosmological simulations, we find that all models
strongly impact the early evolution of galaxies leading to differences in
enrichment and thermodynamic histories. The anisotropic model has the
strongest impact, with little difference between the dynamic and the constant-
coefficient variant. We also find that the metal distribution functions in the
circumgalactic gas are significantly tighter at all redshifts, with the
anisotropic model providing the tightest distributions. This is contrary to a
recent study that found metal mixing to be relatively unimportant on
cosmological scales. In all of our experiments the constant-coefficient
Smagorinsky and anisotropic models rivaled their dynamic counterparts,
suggesting that the computationally inexpensive constant-coefficient models
are viable alternatives in cosmological contexts.

###### keywords:

methods: numerical – turbulence – hydrodynamics – galaxies: evolution –
galaxies: general

††pubyear: 2021††pagerange: Mixing matters–B

##  1 Introduction

Galaxies form and evolve in tempestuous gaseous environments where
hydrodynamics, radiative cooling, and gravity synergize to produce rich
emergent phenomena on a myriad of spatial scales. The immense dynamic range of
scales involved and their interconnectedness prove to be limiting factors in
advancing our understanding of the complete picture of galaxy evolution (see
Naab & Ostriker 2017 for an excellent review). At the forefront of the issue
is hydrodynamical turbulence as it is a multi-scale, non-linear phenomenon
that occurs in almost all galactic environments — directly impacting our
theoretical understanding of galaxy evolution.

While the importance of turbulence in the interstellar medium of galaxies has
been long recognized (see Elmegreen & Scalo 2004 and Scalo & Elmegreen 2004
for reviews), only recently has the role of turbulence in halo gas come under
careful consideration. Indeed, both the circumgalactic medium (CGM) of $L^{*}$
galaxies and the intracluster medium (ICM) of groups and clusters of galaxies
show signs of turbulence playing an important role in their evolution (see,
for example, Prasad et al. 2018 and Wang et al. 2020). Observationally, there
is evidence of complex kinematic structure in the CGM of $L^{*}$-galaxies that
emerged through the revolutionary Cosmic Origins Spectrograph halo survey on
the Hubble space telescope (COS-halos; Tumlinson et al. 2013, 2017). For
instance, Werk et al. (2016) found that turbulent velocities of $50$-$75$
$\mathrm{km}\,\mathrm{s}^{-1}$ explains the broadening of absorption lines in
the CGM that is not otherwise explainable that has subsequently been confirmed
using numerical studies (Buie et al., 2020). Moving up in mass scale, the ICM
also shows evidence of turbulence through indirect observational methods such
as X-ray surface brightness and Sunyaev-Zeldovich fluctuations (Zhuravleva et
al., 2014; Pinto et al., 2015; Zhuravleva et al., 2015; Khatri & Gaspari,
2016; Zhuravleva et al., 2018).

There are two main drivers of turbulence in galactic environments: (a) global
outflows that emerge from star formation processes and supermassive black
holes (SMBHs) within galaxies (Prasad et al., 2015, 2018; Karen Yang &
Reynolds, 2016; Bourne & Sijacki, 2017; Fielding et al., 2017; Fielding et
al., 2018; Sokołowska et al., 2018; Li et al., 2020) and (b) shearing motion
driven by gas in-fall during structure formation (Dekel et al., 2009; Vazza et
al., 2010, 2012; Vazza et al., 2017; Wittor et al., 2017; Bennett & Sijacki,
2020), mergers (ZuHone et al., 2013), and ram-pressure stripping of galaxies
moving through the ICM (Ruggiero & Lima Neto, 2017; Simons et al., 2020). In
both the CGM and ICM, turbulence could provide additional pressure support
(Poole et al., 2006; Vazza et al., 2018; Lochhaas et al., 2020) that prevents
the gas from collapsing and rapidly converting into stars as well as a
physical mechanism to transport energy and metals throughout gas directly –
impacting the cooling profile, star formation cycle, and metal distribution
functions (Shen et al., 2010; Shen et al., 2012, 2013; Brook et al., 2014;
Sokołowska et al., 2018; Escala et al., 2018; Tremmel et al., 2019; Rennehan
et al., 2019; Hafen et al., 2019, 2020). Therefore, understanding the nature
of turbulence is imperative to understand the complete picture of galaxy
evolution.

In a general sense, turbulence results from a large-scale injection of kinetic
energy into a gas, followed by the eventual transformation of that kinetic
energy into thermal energy through a small-scale viscous dissipation process.
Physically, the injected kinetic energy on large scales breaks the gas into
eddies of smaller and smaller sizes until the energy thermalises at the
dissipation scale. This transfer of kinetic energy from the large scale to the
small scale is the turbulent cascade and the ratio between the energy
injection scale $L$ and the viscous dissipation scale, $\eta$, determines the
dynamic range of scales within the cascade. The dissipation scale $\eta$ is
proportional to viscosity of the gas; a more viscous fluid has a larger $\eta$
and a smaller range over which the turbulent cascade may exist (Landau &
Lifshitz, 1987). It is in the cascade where the redistribution of energy and
metals occurs as the eddies transport the fluid properties as mixing
approaches a homogeneous state.

Understanding of galaxy evolution requires treating turbulence simultaneously
with the complexity of the astrophysics — a fact that demands the use of
numerical simulations. Contemporary computational power limits both the
spatial resolution, $h$, of cosmological simulations and the volume of the
Universe, $V$, that we may simulate. To obtain a statistically sound sample of
galaxies requires a representative (large) volume, while capturing the
relevant astrophysics, such as star formation, stellar and supermassive black
hole feedback, requires a small spatial resolution. Therefore, the dynamic
range of a simulation, $L/h$ (where $L\sim V^{1/3}$), is the limiting factor
for cosmological simulations. As a result, those important processes related
to star formation and supermassive black holes usually occur on scales much
smaller than $h$ and are, therefore, unresolved in most simulations. That fact
necessitates sub-grid astrophysical models in cosmological simulations that
reduce the relevant small-scale physics into a coarse-grained representation,
allowing the physics below a scale $h$ to interact with those on scales
$\ell$, where $L>\ell>h$ (for a more in-depth overview, see Somerville & Davé
2015). While there are many successful simulations that use a variety of sub-
grid assumptions to broadly reproduce galaxy populations (e.g. Guedes et al.
2011; Hopkins et al. 2014; Vogelsberger et al. 2014; Schaye et al. 2015; Genel
et al. 2014; Davé et al. 2017; Davé et al. 2019; Tremmel et al. 2019;
Pillepich et al. 2018), one aspect that is often overlooked is the treatment
of sub-grid turbulence.

Therein lies a problem: in hydrodynamical simulations, the physical
dissipation scale $\eta$ is almost always much smaller than the resolution
scale, $\eta\ll h$ (Pope, 2000). For that reason, the kinetic energy that is
flowing in the turbulent cascade reaches some scale $H\gtrsim h$ where it may
no longer progress. If the numerical viscosity of the hydrodynamical method is
not sufficiently strong to thermalise the kinetic energy, there will be a
build-up of kinetic energy at that scale $H\gtrsim h\gg\eta$ (Garnier et al.,
2009). The kinetic energy build-up is a completely unphysical representation
of turbulence and not only impacts the energetics, but the large scale flow
properties such as the redistribution of metals.

Although not explicitly stated, many cosmological hydrodynamical simulation
studies implicitly assume that numerical dissipation is sufficient to mimic
sub-grid turbulence. However, numerical dissipation is not be sufficient to
reduce the kinetic energy build-up and may not represent turbulent flow
statistics in all cases (Sagaut, 2006). This has led to the development of
sub-grid turbulence models whose purpose is to eliminate the build-up of
kinetic energy at the resolution scale by providing a new viscous dissipation
mechanism (or transport mechanism, in the case of metals).

In simulations that use Eulerian (i.e. grid-based) hydrodynamics, there has
been extensive effort to developing models for sub-grid turbulence (see Sagaut
2006; Garnier et al. 2009 and Schmidt 2015 for extensive reference lists). For
astrophysically-relevant work, we refer the reader to Scannapieco & Brüggen
(2008), Pan et al. (2013), Federrath (2013), Schmidt et al. (2014), and
Semenov et al. (2016) (for a review, see Schmidt 2015). In this paper, we
focus on the Lagrangian hydrodynamical methods.

In Lagrangian hydrodynamics, the fluid equations are approximated via fluid
elements that move with the flow. There are three main approaches that are
most common in cosmological simulations: smoothed particle hydrodynamics (SPH)
(Gingold & Monaghan, 1977; Lucy, 1977; Hernquist & Katz, 1989; Hopkins, 2013),
the moving-mesh method (MM) (Springel, 2010), and mesh-free methods (MF)
(Lanson & Vila, 2008a, b; Gaburov & Nitadori, 2011; Hopkins, 2015). Each of
these methods track fluid elements using different discretisation techniques
that lead to different levels of numerical dissipation. In SPH, there is no
inherent numerical dissipation and it, counter-intuitively, produces a deficit
of kinetic energy near the resolution scale rather than a build-up (Bauer &
Springel, 2012; Price, 2012). The MM and MF methods use Riemann solvers to
approximate the fluid equations of motion between neighbouring fluid elements.
Riemann solvers are generally diffusive due to their approximate nature, and
the build-up of kinetic energy is present in both methods on scales up to
$\sim 10$ times the resolution scale (Bauer & Springel, 2012; Hopkins, 2015).

A solution to kinetic energy build-up is to model the action of turbulent
eddies as a viscous process that diffuses momentum (and metals) and dissipates
kinetic energy via a diffusion equation and source term in the energy
equation, respectively. Usually the assumption is that the viscosity, or
diffusivity, depends on velocity fluctuations $v_{\mathrm{eddy}}$ near the
resolution scale $h$. The resulting diffusivity is $D\propto
hv_{\mathrm{eddy}}$, where $v_{\mathrm{eddy}}$ may be a characteristic
velocity, or velocity difference, within the neighbourhood of a fluid element
(Wadsley et al., 2008; Greif et al., 2009).

One particularly important choice of $v_{\mathrm{eddy}}$ is the Smagorinsky
model (Smagorinsky, 1963), which assumes that velocity shear fluctuations
drive dissipation and mixing through $v_{\mathrm{eddy}}\sim h|S^{*}|$, where
$|S^{*}|$ is the magnitude of the trace-free shear tensor. Wadsley et al.
(2008) showed, using SPH, that a Smagorinsky-like model for thermal energy
mixing improved resolved turbulent mixing in their hydrodynamical experiments
– a problem that plagued SPH for decades (Agertz et al., 2007). The
Smagorinsky model has been successfully used to treat metal and thermal energy
mixing in SPH (Wadsley et al., 2008; Shen et al., 2010; Shen et al., 2012;
Brook et al., 2014; Williamson et al., 2016; Tremmel et al., 2017; Wadsley et
al., 2017; Su et al., 2017) and has been extended to other Lagrangian
hydrodynamical methods, such as the MFM method (Rennehan et al., 2019). In
particular, the Smagorinsky model produces a reasonable description of
turbulence in isotropic scenarios (Colbrook et al., 2017), as well as
essential for understanding metal distribution functions of stellar
populations in dwarf galaxies (Escala et al., 2018), and understanding gas
flows in the CGM (Hafen et al., 2019, 2020). While the Smagorinsky model
improves mixing in Lagrangian simulations, it is important to consider two
assumptions in the model: (a) that shear fluctuations (i.e., changes in
$|S^{*}|$) always represent turbulence and (b) that the diffusive process acts
isotropically through the magnitude of the trace-free (i.e. ignoring
compression) shear tensor over the scale $h$. Are either of these assumptions
reasonable?

To address point (a), consider two incompressible fluids of different
densities in contact and in parallel opposing motion. Will one fluid encroach
the other? Will the flow become turbulent? Such a situation results in the
Kelvin-Helmholtz instability and, given certain conditions on the density and
velocities, we know that perturbations on the contact interface will grow
exponentially and mix the fluids, generating vorticity and turbulence
(Chandrasekhar, 1961). However, if the instability condition is not met then
there remains a large $|S^{*}|$ at the interace without turbulence. If we used
the simple assumption of a diffusivity $D\propto h|S^{*}|$ in a simulation we
would diffuse rapidly at the interface even though there are no turbulent
eddies! Piomelli & Liu (1995) recognized the problem and presented a solution:
a method of dynamically calculating the model coefficient at simulation time.
That method was employed in Rennehan et al. (2019) for the first time in
Lagrangian hydrodynamics and they found that the model significantly reduced
over-diffusion in non-turbulent shear flows, such as in rotating galactic
disks and the Kelvin-Helmholtz instability.

The second point (b) concerns the isotropy in $D$ and the discounting of
compression in $|S^{*}|$. Recently Hu & Chiang (2020) showed that a better
representation of sub-grid scale turbulence is obtained by treating the full
velocity tensor $\boldsymbol{\nabla}\otimes\boldsymbol{u}$ for the diffusivity
$\boldsymbol{D}$, now a tensor. The model is the gradient model (Clark et al.,
1979) and is a major change since diffusion now depends on the directionality
encoded in $\boldsymbol{\nabla}\otimes\boldsymbol{u}$ rather than equally in
every spatial direction. The trace of
$\boldsymbol{\nabla}\otimes\boldsymbol{u}$ is automatically included and,
therefore, compression is automatically handled111The trace of
$\boldsymbol{\nabla}\otimes\boldsymbol{u}$ is
$\boldsymbol{\nabla}\cdot\boldsymbol{u}$, the divergence of the flow. – an
important point for highly-compressible turbulence in cosmological flows.
However, Hu & Chiang (2020) used a common approach from the fluid mechanics
community to post-process driven turbulence simulations in order to check if
the model would have improved the results. Our goal is to implement the model
for the mesh-free finite mass method, a method growing in popularity, and
determine its feasibility at simulation time.

We introduce an implementation of the gradient model for Lagrangian
astrophysical simulations and additionally provide methods for computing the
model coefficient at simulation time. In Section 2 we provide a derivation of
the model, as well as a derivation of the dynamic procedure in Section 2.1
that allows calculation of the model coefficient at simulation time. Section 3
describes driven turbulence validation tests of the gradient model, included
at run time, for both eddy viscosity and metal mixing. As a first application,
we describe the qualitative impact of the eddy viscosity and metal mixing
model on cosmological gas phases in Section 4. We present our conclusions and
recommendations in Section 5.

##  2 The gradient model

In finite-mass Lagrangian hydrodynamics, such as the mesh-free finite mass
(MFM) method (Lanson & Vila, 2008a, b; Gaburov & Nitadori, 2011; Hopkins,
2015), the build-up of kinetic energy at the resolution scale demands an
additional dissipation mechanism. Additionally, metals follow the fluid mass
elements throughout the simulation volume and, therefore, the exchange of
metals between fluid elements due to sub-grid scale turbulent motion does not
occur. The crux of the issue is that discretisation of the fluid field leads
to damping out of the high-frequency turbulent fluctuations that should
continue down to the physical dissipation scale. It is useful to think of the
damping action as a high-pass filter acting on the fluid equations of motion.
By applying a general filter to the momentum conservation equation, it is
possible to derive the correct level of mixing that should occur between fluid
elements due to unresolved turbulence.

In general, the filtering action over a scalar field $f_{i}(\boldsymbol{r})$
can be represented as,

$\overline{f}_{i}(\boldsymbol{r})=\int_{D}f_{i}(\boldsymbol{r}^{\prime})G(|\boldsymbol{r}-\boldsymbol{r}^{\prime}|,h)d\boldsymbol{r}^{\prime},$
(1)

where $h$ is the smoothing scale over the domain $D$. We apply this to the
momentum equation to determine the correction terms due to unresolved
turbulence,

$\frac{\partial(\rho\boldsymbol{u})}{\partial
t}+\nabla\cdot(\rho\boldsymbol{u}\otimes\boldsymbol{u}+P\boldsymbol{I})=0,$
(2)

where $\rho$ is the gas density, $P$ is the pressure, and $\boldsymbol{u}$ is
the velocity vector. If we apply equation (1) to equation (2) it follows that,
assuming the filtering operation and derivatives commute,

$\frac{\partial(\overline{\rho\boldsymbol{u}})}{\partial
t}+\nabla\cdot(\overline{\rho\boldsymbol{u}\otimes\boldsymbol{u}}+\overline{P}\boldsymbol{I})=0.$
(3)

For simplicity, we switch to density-weighted variables such that
$\widetilde{\boldsymbol{u}}\equiv\overline{\rho\boldsymbol{u}}/\overline{\rho}$
and

$\frac{\partial(\overline{\rho}\tilde{\boldsymbol{u}})}{\partial
t}+\nabla\cdot(\overline{\rho}\widetilde{\boldsymbol{u}\otimes\boldsymbol{u}}+\overline{P}\boldsymbol{I})=0.$
(4)

The term $\widetilde{\boldsymbol{u}\otimes\boldsymbol{u}}$ is unknown at
simulation time because it relies on information below the resolution scale.
To put the equation in a more manageable form, we add
$\nabla\cdot(\overline{\rho}[\widetilde{\boldsymbol{u}}\otimes\widetilde{\boldsymbol{u}}-\widetilde{\boldsymbol{u}\otimes\boldsymbol{u}}])$
and rearrange,

$\frac{\partial(\overline{\rho}\widetilde{\boldsymbol{u}})}{\partial
t}+\nabla\cdot(\overline{\rho}\widetilde{\boldsymbol{u}}\otimes\widetilde{\boldsymbol{u}}+\overline{P}\boldsymbol{I})=\nabla\cdot(\overline{\rho}[\widetilde{\boldsymbol{u}}\otimes\widetilde{\boldsymbol{u}}-\widetilde{\boldsymbol{u}\otimes\boldsymbol{u}}]).$
(5)

Therefore, we define the sub-grid scalar flux $\boldsymbol{F}$,

$\boldsymbol{F}\equiv\overline{\rho}(\widetilde{\boldsymbol{u}\otimes\boldsymbol{u}}-\widetilde{\boldsymbol{u}}\otimes\widetilde{\boldsymbol{u}}),$
(6)

and retrieve a new equation,

$\frac{\partial(\overline{\rho}\widetilde{\boldsymbol{u}})}{\partial
t}+\nabla\cdot(\overline{\rho}\widetilde{\boldsymbol{u}}\otimes\widetilde{\boldsymbol{u}}+\overline{P}\boldsymbol{I})=-\nabla\cdot\boldsymbol{F}.$
(7)

The sub-grid scale momentum flux $\boldsymbol{F}$ is unknown at simulation
time and must be modelled yet is widely ignored in cosmological simulation
studies which usually focus only on the thermal energy and metal fluxes via
the Smagorinsky model.

There are a myriad of models in the literature for $\boldsymbol{F}$ (see
Garnier et al. 2009 for extensive lists) but here we take the direct approach
of using a Taylor series approximation following Clark et al. (1979) and Hu &
Chiang (2020). We expand222See Appendix A for a full derivation.
$\widetilde{\boldsymbol{u}\otimes\boldsymbol{u}}$ via a Taylor expansion as

$\widetilde{\boldsymbol{u}\otimes\boldsymbol{u}}\approx\widetilde{\boldsymbol{u}}\otimes\widetilde{\boldsymbol{u}}+\epsilon\nabla^{2}(\widetilde{\boldsymbol{u}}\otimes\widetilde{\boldsymbol{u}})$
(8)

and, therefore, the flux becomes

$\boldsymbol{F}=\overline{\rho}(\widetilde{\boldsymbol{u}\otimes\boldsymbol{u}}-\widetilde{\boldsymbol{u}}\otimes\widetilde{\boldsymbol{u}})\approx\overline{\rho}\epsilon\nabla^{2}(\widetilde{\boldsymbol{u}}\otimes\widetilde{\boldsymbol{u}}),$
(9)

where $\epsilon$ is a constant that depends on the kernel scale $h$ as
$\epsilon\propto h^{2}$ (Monaghan, 1989, 2002, 2011). Expanding the term on
the right hand side of equation (9) and keeping only the first derivative
terms we find

$\boldsymbol{F}=2\overline{\rho}Ch^{2}(\nabla\otimes\widetilde{\boldsymbol{u}})(\nabla\otimes\widetilde{\boldsymbol{u}})^{\mathrm{T}},$
(10)

where $C$ is our model parameter.

A similar result emerges when considering the mass flux of metals in a fluid,

$\frac{\partial(\rho Z)}{\partial
t}+\boldsymbol{\nabla}\cdot(\rho\boldsymbol{u}Z)=0.$ (11)

Applying the same approach as before gives an equation for the sub-grid flux
of metals,

$\boldsymbol{F}=2\overline{\rho}C_{\mathrm{Z}}h^{2}(\nabla\otimes\widetilde{\boldsymbol{u}})\cdot\mathbf{\nabla}Z$
(12)

The method for solving equation (7) is detailed in Hopkins (2017). From their
equation (2),

$\boldsymbol{F}=\boldsymbol{K}\cdot(\nabla\otimes\boldsymbol{q}),$ (13)

where $\boldsymbol{K}$ is the tensor describing the diffusive strength, and
$\boldsymbol{q}$ is the fluid field property. Therefore, we identify,

$\boldsymbol{K}\equiv
2\overline{\rho}Ch^{2}(\nabla\otimes\widetilde{\boldsymbol{u}}),$ (14)
$\boldsymbol{q}\equiv\widetilde{\boldsymbol{u}}.$ (15)

The model in equation (10) is known to lead to numerical instability due to
locally stretching flows. Balarac et al. (2013) showed that ignoring the
positive eigenvalues of $\boldsymbol{S}$ (the stretching effect) ensures the
model is always well behaved. Therefore, we follow Balarac et al. (2013) and
only keep the contribution due to the negative eigenvalues of the shear
tensor. However, we must first decompose
$\nabla\otimes\widetilde{\boldsymbol{u}}$ into the symmetric and anti-
symmetric parts,

$\nabla\otimes\widetilde{\boldsymbol{u}}=\frac{1}{2}(\nabla\otimes\widetilde{\boldsymbol{u}}+[\nabla\otimes\widetilde{\boldsymbol{u}}]^{\mathrm{T}})+\frac{1}{2}(\nabla\otimes\widetilde{\boldsymbol{u}}-[\nabla\otimes\widetilde{\boldsymbol{u}}]^{\mathrm{T}})\equiv\boldsymbol{S}+\boldsymbol{\Omega}.$
(16)

We further decompose $\boldsymbol{S}$ into the contribution from the positive
eigenvalues and negative eigenvalues as
$\boldsymbol{S}\equiv\boldsymbol{S}_{\oplus}+\boldsymbol{S}_{\ominus}$. It
then follows that

$\boldsymbol{S}_{\ominus}\equiv\sum_{k=1}^{3}\mathrm{min}(0,\lambda^{(k)})\boldsymbol{e}^{(k)}\otimes\boldsymbol{e}^{(k)},$
(17)

where $\lambda^{(k)}$ is the $k$th eigenvalue and $\boldsymbol{e}^{(k)}$ is
the corresponding eigenvector. Therefore, the new diffusion coefficient tensor
is

$\boldsymbol{K}=2\overline{\rho}Ch^{2}\boldsymbol{S}_{\ominus},$ (18)

where we have reduced the shear contribution into $\boldsymbol{S}_{\ominus}$.
Removing the negative eigenvalues appears to not be necessary for the metal
mixing case as we found no cases of numerical instability in all of our
hydrodynamical tests and cosmological simulations. However, we find that it is
absolutely necessary in cosmological simulations for the momentum flux.

It is important to note that our choice of $h$ differs from that in Hopkins et
al. (2018). We choose $h$ as the kernel radius of compact support as we found
better results in our tests in Section 3. Our value is approximately twice
that found in the FIRE studies (Su et al., 2017; Escala et al., 2018; Hafen et
al., 2019, 2020), but produces $4$ times as much dissipation and metal mixing
since the dependence is squared in equation (18). Our choice of $h$ is the
same as Wadsley et al. (2017) who use the radius of compact support for
turbulent mixing in the GASOLINE-2 code.

###  2.1 The dynamic gradient model

We apply the same procedure in Balarac et al. (2013) combined with the
density-weighted filtering procedure in Rennehan et al. (2019). Following
their notation, we replace $\tilde{f}$ with $\overline{f}$ since fluid
properties are inherently density-weighted in the mesh-free finite mass
method. Although we focus on the velocity fluctuations in the following
procedure, we note that it applies equally to the metal field. The resolved
fluctuations in the flow are,

$\boldsymbol{\mathcal{L}}=\widehat{\boldsymbol{\overline{u}}\otimes\overline{\boldsymbol{u}}}-\widehat{\boldsymbol{\overline{u}}}\otimes\widehat{\overline{\boldsymbol{u}}},$
(19)

where $\widehat{\overline{f}}$ represents the quantity filtered on twice the
resolution scale ($\sim 2\overline{h}$, see Section 2.4 of Rennehan et al.
2019). If we use only resolved quantities (i.e. doubly-filtered) then the
gradient model in equation (10) should reproduce $\boldsymbol{\mathcal{L}}$,

$\boldsymbol{\mathcal{L}}=\widehat{\boldsymbol{\overline{u}}\otimes\overline{\boldsymbol{u}}}-\widehat{\boldsymbol{\overline{u}}}\otimes\widehat{\overline{\boldsymbol{u}}}=2C\widehat{\overline{h}}^{2}\widehat{\overline{\boldsymbol{S}}}_{\ominus}(\boldsymbol{\nabla}\otimes\widehat{\overline{\boldsymbol{u}}})^{\mathrm{T}}.$
(20)

The above equation results in a solution for the one unknown parameter $C$,

$C=\frac{\boldsymbol{\mathcal{L}}\cdot\big{[}2\widehat{\overline{h}}^{2}\widehat{\overline{\boldsymbol{S}}}_{\ominus}(\boldsymbol{\nabla}\otimes\widehat{\overline{\boldsymbol{u}}})^{\mathrm{T}}\big{]}}{||2\widehat{\overline{h}}^{2}\widehat{\overline{\boldsymbol{S}}}_{\ominus}(\boldsymbol{\nabla}\otimes\widehat{\overline{\boldsymbol{u}}})^{\mathrm{T}}||^{2}}=\frac{1}{2}\frac{\boldsymbol{\mathcal{L}}\cdot\boldsymbol{\alpha}}{||\boldsymbol{\alpha}||^{2}},$
(21)

where we have defined
$\alpha\equiv\widehat{\overline{h}}^{2}\widehat{\overline{\boldsymbol{S}}}_{\ominus}(\boldsymbol{\nabla}\otimes\widehat{\overline{\boldsymbol{u}}})$.
Using this determined value of $C$, we use equation (18) as the diffusivity
tensor in the additional flux term from equation (10). Note that
$\overline{h}$ is the resolution scale and not the radius of compact support
of the kernel. The resolution scale is approximately half of the radius of
compact support, or the mean interparticle spacing.

Applying the same procedure to the metal field yields a separate equation for
$C_{\mathrm{Z}}$. The resolved fluctuations take a similar form to equation
(19),

$\boldsymbol{\mathcal{L}}_{\mathrm{Z}}=\widehat{\boldsymbol{\overline{u}}\,\overline{Z}}-\widehat{\boldsymbol{\overline{u}}}\,\widehat{\overline{Z}}.$
(22)

Note that $\boldsymbol{\mathcal{L}}$ is now a vector rather than a rank-2
tensor. Following the procedure we outline above results in an equation for
$C_{\mathrm{Z}}$

$C_{\mathrm{Z}}=\frac{\boldsymbol{\mathcal{L}}_{\mathrm{Z}}\cdot\big{[}2\widehat{\overline{h}}^{2}(\boldsymbol{\nabla}\otimes\widehat{\overline{\boldsymbol{u}}})\cdot\boldsymbol{\nabla}\widehat{\overline{Z}}\big{]}}{||2\widehat{\overline{h}}^{2}(\boldsymbol{\nabla}\otimes\widehat{\overline{\boldsymbol{u}}})\cdot\boldsymbol{\nabla}\widehat{\overline{Z}}||^{2}}=\frac{1}{2}\frac{\boldsymbol{\mathcal{L}}_{\mathrm{Z}}\cdot\boldsymbol{\beta}}{||\boldsymbol{\beta}||^{2}},$
(23)

where we have defined
$\boldsymbol{\beta}\equiv\widehat{\overline{h}}^{2}(\boldsymbol{\nabla}\otimes\widehat{\overline{\boldsymbol{u}}})\cdot\boldsymbol{\nabla}\widehat{\overline{Z}}$.

###  2.2 Comparison to the Smagorinsky model

The gradient model differs from the widely used Smagorinsky model in a subtle
yet important way. Returning to the definition of the sub-grid flux in
equation (10), the Smagorinsky model represents $\boldsymbol{F}$ as,

$\boldsymbol{F}=2\overline{\rho}(C_{\mathrm{s}}h)^{2}||\boldsymbol{S}^{*}||\boldsymbol{S}^{*},$
(24)

where $\boldsymbol{S}^{*}$ is the trace-free symmetric shear
tensor333$\boldsymbol{S}^{*}\equiv\boldsymbol{S}-\frac{1}{3}\mathrm{tr}(\boldsymbol{S})\cdot\boldsymbol{I}$.
The diffusivity tensor is isotropic,

$\boldsymbol{K}_{\mathrm{Smag}}\equiv
2\overline{\rho}(C_{\mathrm{s}}h)^{2}||\boldsymbol{S}^{*}||\,\boldsymbol{I}.$
(25)

and
$||\boldsymbol{K}_{\mathrm{Smag}}||=2\overline{\rho}(C_{\mathrm{s}}h)^{2}||\boldsymbol{S}^{*}||$
with $C_{\mathrm{s}}\sim 0.15$. The constant nature of $C_{\mathrm{s}}$
implies that the diffusivity scales with the magnitude of the symmetric shear
or, more directly, that the strength of turbulent fluctuations is only
determined by fluctuations in the shear strength. That is a good assumption in
purely turbulent flows but fails dramatically in laminar shear flows, where
the shear is not a good indicator of the presence of turbulence.

One solution to over-diffusion in the Smagorinsky model is to dynamically
calculate the coefficient $C_{\mathrm{s}}$ at simulation time based on the
local fluid properties. We showed in Rennehan et al. (2019) that the dynamic
procedure predicts much lower values of $C_{\mathrm{s}}$ in the majority of
our simple hydrodynamical tests. However, we did not consider the impact of
altering the isotropic nature of the diffusivity.

The gradient model has
$\boldsymbol{K}_{\mathrm{Grad}}\propto\boldsymbol{S}_{\ominus}$ with
$||\boldsymbol{K}_{\mathrm{Grad}}||=2\overline{\rho}Ch^{2}||\boldsymbol{S}_{\ominus}||$.
The constant $C$ is yet to be determined but the direction differs from the
dynamic and non-dynamic Smagorinsky models. The diffusivity itself no longer
acts in each direction equally but acts in the direction of the eigenvectors
of $\boldsymbol{S}_{\ominus}$. However, in a simple incompressible, low Mach
number turbulent flows we expect that
$||\boldsymbol{K}_{\mathrm{Smag}}||\sim||\boldsymbol{K}_{\mathrm{Grad}}||$
given that the velocity derivative tensor
$\boldsymbol{\nabla}\otimes\boldsymbol{\tilde{u}}$ is approximately isotropic
in that regime.

Our application of the dynamic method (Piomelli & Liu, 1995) to the gradient
model simultaneously allows $C=C(\boldsymbol{x},t)$ and the diffusivity to be
anisotropic, $\boldsymbol{K}_{\mathrm{Grad}}\propto\boldsymbol{S}_{\ominus}$.
This should be important for any complicated astrophysical flows such as those
we investigate in the following sections.

###  2.3 Models

Table 1: Turbulence models and parameters. The dynamic model calculates the model parameters at simulation and, therefore, the values listed below are the forced upper limit. Label | Dynamic | (Max.) Parameter Value | Anisotropic
---|---|---|---
None | N/A | N/A | N/A
Smag. | $\times$ | 0.15 | $\times$
Dyn. Smag. | ✓ | 0.20 | $\times$
FIRE | $\times$ | 0.05 | $\times$
Grad. | $\times$ | 0.22 | ✓
Dyn. Grad. | ✓ | 1.00 | ✓

Table 1 contains a compact description of our model set. In all cases where
there is a sub-grid turbulence model, we treat both metals and viscosity
simultaneously. There are three categories of models: no sub-grid model
(None), the Smagorinksy model, and the gradient model. The dynamic procedure
allows us to extend the Smagorinsky and gradient models with a model parameter
that depends on spatio-temporal coordinates. Additionally, we test the only
other calibration of the Smagorinsky model in the mesh-free finite mass method
from the FIRE collaboration (Escala et al., 2018).

For the Smagorinsky model, Smag., we use the theoretical value of
$C_{\mathrm{s}}$ $\sim 0.15$ and limit $C_{\mathrm{s}}$ to $0.20$ for the
dynamic Smagorinsky model (Dyn. Smag.) to avoid numerical instability. The
FIRE calibration of the Smagorinsky model (labelled FIRE) is $C_{\mathrm{s}}$
$\approx 0.046$ and we adopt $C_{\mathrm{s}}$ $=0.05$ for simplicity. For
these three versions of the Smagorinsky model the value of $C_{\mathrm{s}}$ is
the same for both metals and eddy-viscosity.

In the new gradient model we use fixed values of $C=0.22$ and
$C_{\mathrm{Z}}=0.22$ for the baseline comparison and label these as Grad. In
our other tests, we use the dynamic procedure outlined in Section 2.1 and
label these tests as Dyn. Grad. We derive our fixed values of $C$ and
$C_{\mathrm{Z}}$ from the approximate median value predicted by the dynamic
procedure in the driven turbulence tests in Section 3.

##  3 Homogeneous Turbulence

Turbulence is ubiquitous in astrophysical flows on a myriad of scales and Mach
numbers. Therefore, in this section, we investigate the impact of the gradient
model on the velocity statistics and metal distributions in homogeneous,
isotropic, driven turbulence at Mach numbers
$\mathcal{M}\in\\{0.3,0.7,2.1\\}$.

For each $\mathcal{M}$, our control (i.e., no sub-grid turbulence model)
simulation set comprises $5$ simulations with particle counts
$N\in\\{64^{3},128^{3},256^{3},512^{3},768^{3}\\}$ within a box of side length
$L=1$ with initial pressure, density, and specific internal energy of $P=1$,
$\rho=1$, $u=1000$, respectively444We use arbitrary code units in all of our
hydrodynamical tests.. Initially, we place equal mass particles on a uniform
grid and then subsequently mix the gas via the prescription of Bauer &
Springel (2012) and Hopkins (2015). We list our turbulent driving parameters
for GIZMO in Table 2 (cf. Table 1 in Bauer & Springel 2012). For more details
of our approach, see Section 3.1 of Rennehan et al. (2019).

Table 2: Parameters for our driven turbulence experiments. Units are arbitrary code units. $\sim\mathcal{M}$ | $\sigma$ | $\Delta t$ | $t_{\mathrm{s}}$ | $k_{\mathrm{min}}$ | $k_{\mathrm{max}}$ | $\tau_{\mathrm{mix}}$
---|---|---|---|---|---|---
0.3 | 0.014 | 0.005 | 1 | 6.27 | 12.57 | 3.33
0.7 | 0.045 | 0.005 | 1 | 6.27 | 12.57 | 1.67
2.1 | 0.21 | 0.005 | 0.5 | 6.27 | 12.57 | 0.50

Our interest lies in measuring the impact of (1) the eddy-viscosity model on
the velocity power spectra of these driven turbulence volumes and (2) the
convergence of metal distribution functions in these volumes. However, these
properties rely on driven turbulence volumes that are in statistical
equilibrium. To gauge whether our simulations are in equilibrium, we define a
mixing timescale $\tau_{\mathrm{mix}}\equiv L/\langle v\rangle$, where
$\langle v\rangle$ is the expected average velocity of the particles in each
volume. However, $\langle v\rangle=v_{\mathrm{s}}\mathcal{M}$, where $v_{s}=1$
is the isothermal sound speed of the gas and, additionally, $L=1$. Therefore,
$\tau_{\mathrm{mix}}=1/\mathcal{M}$. We evolved each control simulation for
several mixing timescales ($\sim 4\tau_{\mathrm{mix}}$) to ensure that the gas
is in steady-state statistical equilibrium and we confirmed the stability of
the Mach number. The mixing timescale and steady-state Mach numbers are listed
in Table 2. We measure the velocity power spectra following the same method as
Bauer & Springel (2012), which is available in the public version of GIZMO.

It may seem out of place to study low resolutions, such as $64^{3}$ and
$128^{3}$, in this study when it is clearly possible to study resolutions up
to $768^{3}$. The ultimate goal of our work is to apply the model to
cosmological simulations which have a huge dynamic range and, therefore, low
resolutions in individual galaxies. For example, IllustrisTNG (Pillepich et
al., 2018) and the RomulusC (Tremmel et al., 2019) simulations both have
particle mass resolutions of $\sim 10^{5}$ $\mathrm{M}_{\odot}$ for each gas
particle, at their best. Consider that in an $L^{*}$-galaxy, we expect perhaps
$\lesssim 10^{10}$ $\mathrm{M}_{\odot}$ of hot gas in the circumgalactic
medium (Anderson & Bregman, 2010). At the best resolutions we have today, that
gives $\sim 10^{5}$ particles per $L^{*}$-halo or $\sim 50^{3}$ particles.
Evidently, contemporary cosmological simulations that capture both hundreds of
$\mathrm{Mpc}$ on the large-scale as well as individual galaxies are far off
from $\gtrsim 256^{3}$ per galactic halo.

###  3.1 Velocity power spectra

A standard measure to determine if eddy viscosity models improve the accuracy
of hydrodynamical simulations is whether the velocity power spectra reproduces
the theoretically predicted Kolmogorov scaling, $E(k)\sim k^{-5/3}$. That
scaling holds for incompressible, low Mach number turbulence
($\mathcal{M}\lesssim 1$) but is shallower than the apparent scaling in
supersonic turbulence, $E(k)\sim k^{-2}$ (Federrath, 2013). In physical
turbulence, the dissipation scale is demarcated by a sharp decline from the
Kolmogorov slope on the smallest scales. In simulations, the resolution scale
forces dissipation to occur on a much larger scale than what would occur in
nature as the physical dissipation scale is unresolved. If the numerical
viscosity of a hydrodynamical method cannot rapidly dissipate that energy,
there will be a build-up of kinetic energy near the resolution scale that
causes an unphysical representation of turbulence.

Eddy-viscosity models introduce additional dissipation in the gas by
accounting for the unrepresented scales in the flow, or equivalently by
minimizing the error from the missing terms in the equations of motion. The
build-up of kinetic energy is usually observed as a "bump" in the velocity
power spectra where there is artificial correlation in the velocities on small
scales. However, before we discuss the impact of eddy viscosity models on the
power spectrum we must first test the convergence of the mesh-free finite mass
(MFM) method in simulations without eddy viscosity.

Figure 1: The velocity power spectra of our simulations with no sub-grid eddy
viscosity model at three Mach numbers from top to bottom, respectively:
$\mathcal{M}=0.3$, $0.7$, and $2.1$. We compensate the power spectra by $k$
for easy comparison with Bauer & Springel 2012. The coloured lines show
resolutions $64^{3}$, $128^{3}$, $256^{3}$, $512^{3}$, and $768^{3}$ from
lightest to darkest, respectively. In each panel, we show the predicted
scaling with a dotted line of arbitrary normalisation. Only in the highest
Mach number case do we see many of the scales in the inertial range
represented and clear convergence of the hydrodynamical method.

Fig. 1 shows the velocity power spectra for our set of simulations with
particle counts $64^{3}$, $128^{3}$, $256^{3}$, $512^{3}$, and $768^{3}$
coloured by lines from lightest to darkest, respectively. The power spectra
are compensated by $k$ for easy comparison to Bauer & Springel (2012). The
panels are ordered from lowest to highest Mach number from top to bottom —
$\mathcal{M}\sim 0.3$, $0.7$ and $2.1$, respectively. In each panel, the
dotted line represents the predicted scaling but at an arbitrary
normalisation.

From the top panel of Fig. 1, it is apparent that $64^{3}$ and $128^{3}$ do
not faithfully represent a turbulent gas as they are dominated by the bump.
More precisely, the $E(k)$ scaling is much too shallow compared to Kolmogorov
turbulence for a wide range of $k$. Our $256^{3}$ simulation shows an inkling
of the inertial range scaling but is slightly too steep below $k\lesssim 20$
and dominated by the bump at $k\gtrsim 30$. As we move up in resolution the
inertial range only grows slightly. At our highest resolution, the inertial
range spans $k\sim 40$ to $k\sim 60$ and the bump dominates the small scales.
We skip discussion of the middle panel as the results are qualitatively
equivalent between $\mathcal{M}\sim 0.3$ and $0.7$.

In the bottom panel of Fig. 1, we show the compensated power spectra for
supersonic turbulence at $\mathcal{M}\sim 2.1$. It is immediately evident that
the simulations converge much more rapidly to the proper scaling than in
subsonic turbulence. Still, at $64^{3}$ resolution the power spectra is
dominated by the build-up of kinetic energy near the resolution scale, with
the inertial range only beginning to appear at $128^{3}$ resolution. At our
highest resolution, $768^{3}$, there is arguably an entire order-of-magnitude
resolved in the inertial range before the build-up of kinetic energy dominates
at the smallest scales.

Figure 2: The velocity power spectrum of our turbulence volumes with various
eddy viscosity models. Rows are Mach numbers $0.3$, $0.7$, and $2.1$ from top
to bottom, respectively. The left column has no additional increase in
dissipation whereas the right column has an order-of-magnitude boost in
dissipation on particles with $\mathcal{M}<1$. The black curve shows the
$256^{3}$ simulation with no eddy viscosity model. The other solid lines show
the Smag. and Dyn. Smag. models, and the dashed and dotted lines show the
Grad. and Dyn. Grad. models, respectively. While all of the models improve the
inertial scaling, an additional boost factor ($\gamma\sim 10$) for subsonic
particles is required to reproduce the proper scaling at all Mach numbers.

Evidently, sub-grid eddy viscosity models are required for the mesh-free
finite mass (MFM) method at any resolutions that may be used in astrophysical
environments. We show the impact of our eddy viscosity models in Fig. 2. The
left column shows the compensated power spectra, $E(k)$, as a function of
wavenumber $k$. The rows represent the same three Mach numbers
$\mathcal{M}=0.3$, $0.7$, and $2.1$ from top to bottom, respectively. The
dotted black line shows the Kolmogorov scaling at each Mach number. All of the
simulations were run at $256^{3}$ resolution since that is the resolution
where, with no eddy viscosity model, we begin to see an extended inertial
range and distinguish the kinetic energy "bump". The coloured curves show the
eddy viscosity models, shaded from lightest to darkest: Dyn. Smag (solid
salmon), Smag. (solid magenta), Dyn. Grad. (dotted purple), and Grad. (dashed
purple). To explain the right column, we must first explain the results in the
left column.

The left column of Fig. 2 shows the velocity power spectra of the turbulent
gas in our simulated volumes with the same models presented in Section 2. The
black curve shows the control experiment at $256^{3}$ resolution, i.e. the
simulation with no eddy viscosity model and only numerical dissipation as in
Fig. 1. All Mach numbers show the same trend: the eddy viscosity models have
little impact on reducing the build-up of kinetic energy at small scales.
Especially important is that the sub-sonic ($\mathcal{M}\lesssim 1$)
simulations are much less improved than the supersonic case. However, the new
gradient model variants, Grad. and Dyn. Grad., dissipate slightly more rapidly
and allow for a steeper slope closer to the Kolmogorov scaling.

We did not expect that all of the eddy viscosity implementations would fail to
reduce the kinetic energy build-up a priori. The fact that there is not enough
dissipation suggests that some physical property in the diffusion tensor was
assigned incorrectly. As we outlined in Section 2.2, the diffusion strength
for the Smagorinsky and the gradient model classes should effectively scale
with each other
($||\boldsymbol{K}_{\mathrm{Smag}}||\sim||\boldsymbol{K}_{\mathrm{Grad}}||$)
in isotropic, homogeneous turbulence. Therefore, in both classes of models,
there are only two physically-motivated quantities that control the strength
of diffusion: the length-scale $h$ and the velocity tensor
$\boldsymbol{\nabla}\otimes\boldsymbol{u}$. The velocity tensor should not be
the issue since it has been verified through the hydrodynamical tests in
Hopkins (2015) and would cause the MFM method to fail drastically if the
velocity gradients were incorrect. That leaves $h$ as the issue, suggesting
that our estimate of the scale over which the eddy viscosity interactions
propagate is underestimated555Recall that our definition of $h$ is already
twice as large as Hopkins (2017), leading to $4$ times as much dissipation and
mixing.. Therefore, we introduce a boost factor $\gamma$ to the diffusion
tensor $||\boldsymbol{K}||$ (i.e.,
$||\boldsymbol{K}||\rightarrow\gamma||\boldsymbol{K}||$) in order to get a
more reasonable scaling in the inertial range.

We determined the boost factor by running a series of driven turbulence tests
with discrete $\gamma$ from $\gamma=1$ to $\gamma=100$. We did not perform a
quantitative fit to the Kolmogorov scaling as our interest is in the
approximate offset required to improve the inertial scaling and the exact
$\gamma$ is unimportant for the statistics of the flow. We additionally found
that we only need to correct the diffusion strength in particles that are
subsonic, $\mathcal{M}<1$, with the Mach number for each particle derived from
the current velocity of the particle divided by its thermal sound speed.

The right column of Fig. 2 shows the velocity power spectra of our simulations
with a dissipation boost factor of $\gamma\sim 10$ on only the subsonic
particles in each simulation. At all Mach numbers the additional dissipation
causes the kinetic energy to convert into thermal energy much more readily,
causing the disappearance of the additional power at small scales, $k\gtrsim
40$. The gradient models Grad. and Dyn. Grad. perform better than the
Smagorinksy variants, but only slightly. That is expected in isotropic
homogeneous turbulence since the dissipation strength is effectively the same,
and confirms that our implementation of the gradient model is a successful
eddy viscosity model.

It is very important that we emphasise the results we observe in Fig. 1 are
very similar to those found in Bauer & Springel (2012) and Hopkins (2015) for
the moving-mesh method (MM; as implemented in AREPO) and the mesh-free finite
mass method (MFM; as implemented in GIZMO), respectively. The only difference
is that we show results for much higher resolutions where we begin to see an
extended inertial range. Both the MFM and MM hydrodynamical methods produce
the same kinetic energy bump that exists in grid-based methods and, therefore,
we would expect an eddy viscosity model to also solve the problem in the MM
method, although we do not test that in this work. To reiterate, it is
necessary that all hydrodynamical simulations resolve the inertial range
combined with an immediate sharp drop-off in power, or they are not
reproducing what we physically observe as turbulence below a certain scale,
where the build-up of kinetic energy begins to dominate.

###  3.2 Metal mixing

After each simulation reached $\sim 4\tau_{\mathrm{mix}}$, we treated each
steady-state volume as new initial conditions for our metal mixing study. In
each volume, we gave the densest $50\%$ of particles a metal mass fraction of
$Z=1$ while keeping the rest $Z=0$. The metals in our simulations act as
passive scalars and have no impact on the flow properties. We will test the
model with more realistic metal distributions in Section 4.

First, we must determine if our simulations converge toward a solution for the
metal distribution as we increase resolution. We ran each of the metal
enriched volumes for an additional $4\tau_{\mathrm{mix}}$ to sample a wide
variety of metal distribution states. We expect a priori that by $\sim
2\tau_{\mathrm{mix}}$ the metal-enriched particles should be scattered
approximately homogeneously since a particle with the typical velocity
$\langle v\rangle$ should have crossed the volume twice in that time. Although
that is true for all resolutions, how can we compare each resolution on equal
footing after it has reached equilibrium?

The appropriate comparison involves smoothing the spatial distribution of
metals on the same scale in all of our simulations. The main assumption is
that our simulations with particle counts $\geqslant 128^{3}$ contain more
accurate information on the scales equivalent to our $64^{3}$ simulations.
Equivalently, if we degrade the resolution of the highest resolution
simulations to the lowest resolution, we should hope to obtain a result
similar to the lowest resolution simulation. To degrade the resolution for
each simulation, we first kernel-weight the particle data onto a grid with
resolution twice as fine as minimum smoothing length in the simulation,
$\Delta x_{\mathrm{sim,i}}$. Next, we smooth the grid data on a physical scale
equivalent to our $64^{3}$ simulation using a uniform top-hat
filter666Specifically, we use the uniform_filter function from the scipy
package in Python, with periodicity enabled. with width $w=\Delta
x_{\mathrm{low}}/\Delta x_{\mathrm{sim,i}}$, where $\Delta x_{\mathrm{low}}$
is always $\Delta x_{\mathrm{low}}\equiv 1/64$ since our box has length $L=1$.
See Appendix B for a comparison of the smoothed metal distributions.

Figure 3: The standard deviation of the Gaussian fit to the metal
distributions in our driven turbulence experiments as a function of
resolution, $N_{\mathrm{x}}$. Mach numbers are $\mathcal{M}=0.3$, $0.7$, and
$2.1$ from top to bottom, respectively. The stars show the simulations with no
sub-grid metal mixing model at resolutions given by the number associated with
each point. The simulations at $64^{3}$ resolution with sub-grid metal mixing
models are given by the figure labels. The dotted line shows an exponential
fit to the no metal mixing model simulations. Clearly, the gradient model
allows us, at $64^{3}$ resolution, to reproduce the spatial redistribution of
metals due to turbulence at a level better than $4-8$ times the resolution.

Fig. 3 shows the standard deviation ($\sigma_{\mathrm{Z}}$) of the smoothed
metal distribution as a function of $N_{\mathrm{x}}=\\{64,128,256,512,768\\}$
in our simulations at $t\sim 2\tau_{\mathrm{mix}}$, for Mach numbers $0.3$,
$0.7$, and $2.1$ from top to bottom, respectively. The stars correspond to the
simulations without a sub-grid metal mixing model at the resolution given by
their labels. The remaining symbols in the legend correspond to the
simulations at $64^{3}$ with a sub-grid metal mixing model (see Table 1 for a
description). The dotted line shows an exponential decay fit to test
convergence in the simulations without metal mixing. We expect decreasing
$\sigma_{\mathrm{Z}}$ with increasing resolution since hydrodynamical mixing
is more resolved and the metal value in each grid cell approaches the mean.

At $\mathcal{M}\sim 0.3$ in the top panel of Fig. 3, $\sigma_{\mathrm{Z}}$
follows an exponentially decreasing trend with resolution in the simulations
without metal mixing, as expected. However, there is not evidence for strong
convergence at our highest resolution — although it appears the curve is
beginning to flatten. The inverted triangle shows the results for the Dyn.
Smag. model and the left-pointing triangle shows the result for the Smag.
model. Both the dynamic and standard Smagorinsky models predict a more
reasonable $\sigma_{\mathrm{Z}}$, reproducing a metal distribution closer to a
resolution of $512^{3}$. The FIRE calibration of the Smagorinsky model, marked
as $\times$, shows little improvement in $\sigma_{\mathrm{Z}}$; the result is
effectively equivalent to having no model at all. The Dyn. Grad. and Grad.
models are marked by a plus sign and diamond, respectively. The gradient model
obviously improves $\sigma_{\mathrm{Z}}$ and can, at $64^{3}$ resolution, also
reproduce a metal distribution equivalent to a resolution $512^{3}$. The
trends are equivalent for $\mathcal{M}\sim 0.7$ turbulence, so we continue to
the next panel.

The bottom panel of Fig. 3 shows how the sub-grid metal mixing models impact
supersonic turbulence at $\mathcal{M}\sim 2.1$. There is much better
convergence of the metal distribution at $t\sim 2\tau_{\mathrm{mix}}$ than in
subsonic turbulence at this time. The qualitative trend remains the same for
the $64^{3}$ simulations with sub-grid metal mixing: both the gradient models
and Smagorinsky models perform equally as well. However, the width
$\sigma_{\mathrm{Z}}$ of the distributions are much wider in $\mathcal{M}\sim
2.1$ turbulence. With the exception of the FIRE calibration, the other metal
mixing models at $64^{3}$ resolution predict metal distributions similar to
$256^{3}$.

Our turbulence tests with metals demonstrate that sub-grid metal mixing models
are necessary in the mesh-free finite mass (MFM) method if one desires more
accurate metal distributions. The method we provided for calibrating the metal
mixing models is important and, we argue, must be investigated whenever one
implements a new metal mixing model into any hydrodynamics solver.
Particularly, one must calibrate the model parameter at the resolution they
desire for their mixing model if a dynamic procedure is not applied. That is
an important point: the dynamic procedure (Dyn. Smag. and Dyn. Grad.) allows
us to approximate the calibrated model parameter for the corresponding model
(Smag. and Grad.) without carrying out the calibration. However, the true
power of the dynamic procedure is in simulations with mixtures of non-
turbulent and turbulent gas at various Mach numbers. In those complex
environments the dynamic procedure automatically adjusts the model parameter
and, as we showed in Rennehan et al. (2019) for the Dyn. Smag. model,
drastically alters the resulting metal distributions.

As expected, in pure homogeneous, isotropic turbulence all of the sub-grid
metal mixing models improved the accuracy of the metal distributions when
using a proper calibration. That is expected in the homogeneous, isotropic
case since all of the models effectively act in the same way, on average. The
power of the dynamic procedure, and the new anisotropic model is in
cosmological simulations where many complex flows interact.

###  3.3 Applicability to other hydrodynamical methods

While our interest lies in the mesh-free finite mass (MFM) method, the models
in Section 2 and experiments in Section 3 are applicable to other
hydrodynamical solvers. Specifically, our results extrapolate with minor
modification to grid-based hydrodynamics. Only slight changes to the filtering
method must be implemented, as outlined in Schmidt (2015). More care must be
taken when applying the model to smoothed particle hydrodynamics (SPH).
However, as we outline below, there is much broader applicability to the
moving-mesh (MM) method.

First we consider eddy viscosity. In SPH, it is well-established that there is
a deficit in power near the resolution scale rather than a build-up of kinetic
energy (Bauer & Springel, 2012; Price, 2012; Hopkins, 2013, 2015). That fact
suggests that SPH reproduces turbulence better than the MFM or MM methods, but
produces results at a much lower effective resolution. The lack of power,
rather than the overabundance of power, implies that an eddy viscosity model
would only further degrade the resolution of SPH results and not improve the
inertial scaling. Therefore, we do not recommend eddy viscosity models for SPH
but rather the work of Di Mascio et al. (2017), who recently provided an SPH
equivalent.

For the MM method, the build-up of kinetic energy at the resolution scale is
present (Bauer & Springel, 2012) and of equivalent magnitude to our results in
the MFM method. Therefore, we recommend investigation into eddy viscosity
models for the MM method as they could improve the inertial scaling. In terms
of implementation, all of the derivations in Section 2 apply to the MM method.

Metal mixing using the Smagorinsky model has been studied in cosmological
simulations involving SPH but widely ignored in the MM method. However, no
calibration technique has been provided by the community for SPH and the
calibrations usually follow the theoretical value for the Smagorinsky model
(e.g., Shen et al. 2010; Williamson et al. 2016) or calibrations that depend
on sub-grid astrophysics models (e.g., Wadsley et al. 2017; Escala et al.
2018). The metal mixing calibration technique in Section 3 is completely
applicable to SPH since metals are treated equivalently to the MFM method, and
they are both constant mass methods. Additionally, our calibration technique
does not depend on the uncertainties within astrophysical sub-grid models —
only pure hydrodynamics.

For the MM method, all of Section 2 is applicable for metal mixing since the
MM method relies on transport equations such as equation (2) for advecting
metals throughout the fluid (Springel, 2010). In fact, Balarac et al. (2013)
find that the gradient model improves the inertial scaling in the power in the
metal field through identical transport equations. However, convergence tests
such as those in Fig. 3 are necessary in order to determine whether they are
truly required, or if numerical dissipation is adequate.

##  4 Cosmological Simulations

Understanding the evolution of galaxies is a complex enterprise involving
highly non-linear coupled physical processes. Not only do stellar feedback and
active galactic nuclei produce powerful outflows that drive turbulence locally
in the interstellar media of galaxies, but also in the gas reservoirs
surrounding galaxies. Turbulence also appears through the Kelvin-Helmholtz
instability in ram pressure stripping of galaxies moving through a hot medium,
and through the stellar winds from stars making their way out from the galaxy
into the circumgalactic medium.

The physical processes above occur on spatial scales much smaller than
currently possible to resolve in the average contemporary cosmological
simulation. For that reason, the majority of astrophysics in cosmological
simulations are encoded into parametrised sub-grid models that use the
information on the largest scales to predict what occurs below the resolution
of the simulation. The resulting calculations usually indicate how much mass
and energy should be injected (or removed) from the large-scale gas and
stellar components. However, there is no one correct way to approximate the
astrophysics on the sub-grid scale since it highly depends on the maximum
possible resolution, hydrodynamical method, as well as other complex numerical
effects. All of the issues with numerics and missing physics usually ends up
in one or more tunable free parameters in the model.

Assuming such a sub-grid astrophysical model is developed, how do we know that
it is correct? Or, at least approximating reality? Normally, one or more
trusted astronomical observation is used to test the validity of all of the
sub-grid astrophysics that may exist. A common example would be the galaxy
stellar mass function, or the $\mathrm{M}_{\mathrm{BH}}$ \-
$\sigma_{\mathrm{*}}$ relationship that links supermassive black hole masses
to the stellar masses of their host galaxies. However, two different
hydrodynamical methods may provide different parameter values for the same
sub-grid astrophysical models. Additionally, there may be two completely
different approaches to modelling the same physical phenomenon with no clear
mapping between free parameters! Calibrating sub-grid astrophysical models is
obviously a complicated endeavour and must be built on a strong hydrodynamics
base. How can we begin to trust that our understanding of the astrophysics of
small scale approximate reality if the hydrodynamics, as we showed in Section
3, does not reproduce reality?

Our goal is in determining whether the converged and separately calibrated
sub-grid turbulence models we presented in Section 2 have any significant
impact over the sub-grid astrophysical models that we use in large-volume
cosmological simulations. As a first step, we only investigate the broad,
qualitative impact on the gas properties in gaseous halos in a single set of
sub-grid astrophysical models. We stress that we do not intend to reproduce
the full galaxy population in a calibrated and predictive sense. Additionally,
we note that more testing is required across the all of sub-grid astrophysical
models that exist in the literature, as the turbulent mixing models may
interact in unexpected ways due to the non-linearity of the problem.

###  4.1 The simulations

For our comparison, we choose to use the SIMBA galaxy formation model. SIMBA
includes robust sub-grid models of star formation, cooling, stellar feedback,
chemical enrichment, active galactic nuclei feedback, and dust evolution — all
evolved in concert with the mesh-free mass method (MFM) (Davé et al., 2016,
2019). For this study, we implemented the SIMBA models into the public version
of GIZMO as described in Davé et al. (2019) and we point the interested reader
to that study for the details of the sub-grid models. We follow the approach
of Schaye et al. (2015) and calibrate our implementation of the SIMBA model
only to the galaxy stellar mass function and the $M_{\mathrm{BH}}$-$M_{*}$
relationship at $z=0$ for the purposes of this study.

Table 3: Cosmological parameters and simulation information for our simulation set. Our parameters follow Davé et al. 2019 except begin at a higher redshift. Cosmological Parameters |
---|---
$\Omega_{\mathrm{m,0}}$ | $0.3$
$\Omega_{\mathrm{\Lambda,0}}$ | $0.7$
$\Omega_{\mathrm{b,0}}$ | $0.048$
$h$ | $0.68$
$\sigma_{\mathrm{8}}$ | $0.82$
$n_{\mathrm{spec}}$ | $0.97$
Simulation Information |
$z_{\mathrm{begin}}$ | 249
$N_{\mathrm{part}}$ | $2\times 256^{3}$
$L_{\mathrm{side}}$ |  $25$ $\mathrm{cMpc}\,h^{-1}$
$m_{\mathrm{part,gas}}$ |  $1.26\times 10^{7}$ $\mathrm{M}_{\odot}h^{-1}$
$m_{\mathrm{part,dark}}$ |  $6.88\times 10^{7}$ $\mathrm{M}_{\odot}h^{-1}$
$\epsilon_{\mathrm{soft,min}}$ |  $0.5$ $\mathrm{kpc}\,h^{-1}$

We run 6 cosmological-scale volumes of side-length $L=25$
$\mathrm{cMpc}\,h^{-1}$ ($\sim 37$ $\mathrm{cMpc}$) in order to compare our
various mixing models. The simulations begin from initial conditions generated
with MUSIC (Hahn & Abel, 2011) at a redshift of $z=249$, with a standard
$\Lambda$ Cold Dark Matter cosmology (see Table 3 for values). The mass
resolution in gas and dark matter follows the SIMBA simulations with
$m_{\mathrm{part,gas}}=1.26\times 10^{7}$ $\mathrm{M}_{\odot}h^{-1}$ and
$m_{\mathrm{part,dark}}=6.88\times 10^{7}$ $\mathrm{M}_{\odot}h^{-1}$,
respectively. We use adaptive gravitational softening (Hopkins, 2015; Hopkins
et al., 2018) to compute the softening lengths of all of our particles, and
enforce a minimum softening length of $\epsilon_{\mathrm{soft,min}}=0.5$
$\mathrm{kpc}\,h^{-1}$.

###  4.2 Global metal mixing

While sub-grid scale turbulence models maximally impact the smallest scales in
a cosmological simulation, their integrated effect impact the properties of
the largest scales, such as global metal distribution functions (Shen et al.,
2010; Escala et al., 2018; Rennehan et al., 2019). Therefore, in this Section,
we examine the impact of the gradient model on the metal distribution
functions (MDFs) in the circumgalactic medium (CGM) and warm-hot intergalactic
medium (WHIM) — both known to be turbulent environments (Iapichino et al.,
2011; Iapichino et al., 2013; Tumlinson et al., 2017).

Our definition of CGM and WHIM depends on separating gas that is bound to
halos from that which is unbound, at a given epoch. A good estimation comes
from Davé et al. (2010),

$\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,$
(26)

where $\Omega_{\mathrm{b}}(z)$ is the baryon fraction as a function of
redshift, $\Omega_{\mathrm{m}}(z)$ the matter fraction,
$\rho_{c}(z)=3(H(z))^{2}/(8\pi G)$, and $H(z)$ the redshift-dependent Hubble
function. All gas above $\rho_{\mathrm{bound}}(z)$ we consider bound to halos,
and have confirmed that the approximation holds well.

We define the CGM to be all gas in the volume that is above
$\rho_{\mathrm{bound}}(z)$ in equation (26) and below the star formation
density threshold, $\rho_{\mathrm{*,crit}}=4.4\times 10^{-25}$ g cm-3, at any
temperature. That includes gas in the intragroup medium of our most massive
halos in the $(25\,\mathrm{cMpc}\,h^{-1})^{3}$ volumes. The WHIM is all gas
that is below $\rho_{\mathrm{bound}}(z)$ and above a temperature of $T=10^{5}$
$\mathrm{K}$.

Fig. 4 shows the metal distribution functions (MDFs) for our two gas phases in
columns: WHIM (left) and CGM (right), and at $z=6$, $4$, $2$, and $0$ in rows
from top to bottom, respectively. These are probability density functions, and
were constructed by binning the particle metallicities in the range
$10^{-6}<\log(Z/Z_{\odot})<10^{1}$, where $Z_{\odot}=0.0134$ (Asplund et al.,
2009). The black curves show the control simulation, None, with no sub-grid
metal mixing. The coloured curves show the simulations with sub-grid metal
mixing and are, from lightest to darkest: Dyn. Smag (solid salmon), Smag.
(solid magenta), FIRE (dotted magenta), Dyn. Grad. (dotted purple), and Grad.
(dashed purple). See Table 1 for more details.

Figure 4: The metal distribution function across global phases: the warm-hot
intergalactic medium (WHIM; left column) and the circumgalactic medium (CGM;
right column). Redshifts are given by the rows from top to bottom: $z=6$, $4$,
$2$, and $0$, respectively. The gradient model variants mix metals much more
rapidly in the early stages of the simulation, suggesting that the choice of
model will impact studies that focus on enrichment timing of halo gas.

First we focus on the WHIM. At $z=6$, there are two distinct components across
all of our model variants. The peak at $Z\sim 10^{-1}$
$\mathrm{Z}_{\mathrm{\odot}}$ is the highly enriched interstellar medium (ISM)
gas that recently joined the WHIM via stellar winds from the integrated star
formation in the early universe. The lower distribution is gas that has mixed
into the WHIM but did not recycle through the ISM, missing the opportunity for
further enrichment via supernova feedback. It is important to note that since
there is no mixing in the None case, when a particle leaves the ISM it cannot
change its metallicity. The models that include sub-grid mixing show varying
spread in the MDFs, with the FIRE and Dyn. Smag. showing the widest spread.
The Dyn. Grad. and Grad. models show the tightest distributions, with the two
peaks in the distributions seemingly merging at $Z\sim 10^{-1.5}$
$\mathrm{Z}_{\mathrm{\odot}}$. The Smag. model matches the Dyn. Smag. model at
$Z\gtrsim 10^{-3}$ $\mathrm{Z}_{\mathrm{\odot}}$, but is biased toward higher
metallicities below that threshold.

The next 3 panels in the left column of Fig. 4 show the MDFs in the WHIM for
$z=4$, $2$, and $0$, from top to bottom, respectively. The trend for all
models is to approach a singly-peaked distribution as the simulation evolves
through cosmic time. Most of the evolution in the MDFs occurs from $z=6$ to
$z=2$ after which the distributions are mostly stationary. The transition from
$z=6$ to $z=4$ demonstrates how rapidly the WHIM evolves at high redshift, and
how each sub-grid metal mixing model impacts the MDFs with different mixing
rates. Specifically, the Dyn. Grad. and Grad. models predict similar
distributions at $z=4$, and produce the tightest MDFs compared to all of the
other models. In fact, there is the same trend at $z=4$ as at $z=6$ — the
gradient model variants (Dyn. Grad. and Grad.) predict tighter MDFs, followed
by wider distributions in the Smagorinsky models (Smag., Dyn. Smag., and FIRE,
respectively).

There are similar trends in the CGM, as we see in the right column of Fig. 4.
To reiterate, the panels show redshifts $z=6$, $4$, $2$, and $0$ from top to
bottom, respectively.

At $z=6$, there is a clear distinction between the distribution at $Z\lesssim
10^{-2}$ $\mathrm{Z}_{\mathrm{\odot}}$ and $Z\gtrsim 10^{-2}$
$\mathrm{Z}_{\mathrm{\odot}}$ in the None case, and the Smagorinsky variants
(Smag., Dyn. Smag., and FIRE). Stellar feedback drives the peak at $Z\sim
10^{-1}$ $\mathrm{Z}_{\mathrm{\odot}}$ similarly to the WHIM at this redshift,
whereas the distribution at $Z\lesssim 10^{-2}$ $\mathrm{Z}_{\mathrm{\odot}}$
is from the very first generations of stars. By this time the Dyn. Grad. and
Grad. models have mixed the most rapidly to create a single broad distribution
in their MDFs. All of the models with sub-grid metal mixing have much more gas
mass enriched above $Z>10^{-6}$ $\mathrm{Z}_{\mathrm{\odot}}$ than the None
case, especially compared to the deficit at $Z\sim 10^{-1.5}$
$\mathrm{Z}_{\mathrm{\odot}}$ in the None case. The Smagorinsky models vary in
mixing rate as Smag., Dyn. Smag., and FIRE, from fastest to slowest,
respectively. The MDFs in the CGM at redshifts $z=4$ to $z=0$ demonstrate the
same trends as in the WHIM phase at the same redshifts: the gradient models
mix much more rapidly at early stages than the Smagorinsky models. At $z=0$
the Dyn. Grad., Grad., and Smag. models predict the same distribution in the
global CGM phase, whereas the Dyn. Smag. and FIRE models predict slightly less
enriched gas.

The MDFs in the turbulent WHIM and CGM show the importance of sub-grid metal
mixing models in cosmological simulations as well as the importance of model
choice. In all cases we include metal mixing, the MDFs are significantly
tighter at all redshifts we measure and significantly tighter for the CGM at
$z=0$. This is contrary to the study in Su et al. (2017) that found metal
mixing to be relatively unimportant on cosmological scales. However, we find
that the FIRE calibration is much too low to reproduce the correct converged
hydrodynamical mixing of metals (see Section 3.2). With our new calibrations
of the Smagorinsky model, Smag., and the new gradient models, Dyn. Grad. and
Grad., we see significant differences at all redshifts.

A common theme in theoretical galaxy evolution is that equivalent results
between different models at $z=0$ does not necessarily imply a similar
integrated history. The evolutionary paths for each sub-grid metal mixing
model all evolve at slightly different rates as we would expect based on their
diffusivities from Section 2. The lesson from our results is that the metal
mixing model choice impacts the early development phases of galaxies rather
than the long term equilibrium stages. At higher redshifts, $z\gtrsim 2$, gas
is collapsing to form galaxies, while stellar feedback and supermassive black
holes are driving outflows out of the potential wells and forcing turbulence.
Our inclusion of the full diffusion tensor in the Dyn. Grad. and Grad. models
allows the gas that is compressing from feedback and infall to further mix its
metal mass with nearby neighbours, tightening the MDFs.

The Smagorinsky variants also improve the results and allow metal mixing
between gas particles but produce broader distributions, notwithstanding the
Dyn. Smag. and Smag. models showing good matches in the convergence tests of
Section 3.2. This is an important point: the simple turbulence tests in
Section 3.2 showed agreement between the gradient and Smagorinsky models
(except for the lower FIRE calibration) but now we see disagreement in complex
cosmological environments. Evidently, ignoring the trace of the velocity
tensor is not the correct approach for cosmological contexts and we recommend
either the Grad. or Dyn. Grad. models as we see no difference with the dynamic
procedure applied to the constant coefficient gradient case.

###  4.3 The impact of eddy viscosity

While the velocity power spectra results in Section 3.1 show that eddy
viscosity is required in Lagrangian finite mass methods, we find no
significant impact of eddy viscosity on the average gas and galaxy properties
in our cosmological simulations. However, we should expect these results
considering that the eddy viscosity models we tested in Section 3.1 had no
impact on the largest scales of the simulation (as they should not). In terms
of global gas properties, we investigated the vorticity, temperature, and
density distributions of the warm-hot intergalactic medium (WHIM) and
circumgalactic medium (CGM) phases777Definitions in Section 4.2. and found
only minor differences in data binned by galaxy stellar mass. Additionally,
there were little differences between stellar mass distributions in our small
volume, low resolution tests. To compare the impact of eddy viscosity between
the models we need to investigate the small scales of the simulations in a
controlled manner. For that reason, we will use the most massive galaxy in our
cosmological simulations as a qualitative case study of the impact of the eddy
viscosity model on the halo gas since this galaxy represents the same system
in each simulation we ran.

First, we investigate the temperature projections of the most massive (in
stellar mass) galaxy at redshifts $z=2$, $1$, and $0$. We confirmed that the
most massive galaxy at $z=2$ ends up as part of the most massive galaxy at
$z=1$ and $z=0$. We restrict our analysis to a radius of $500$ $\mathrm{kpc}$
(physical) for the purposes of this introductory study. Additionally, given
that we see the Smag., Grad., and Dyn. Grad. models as the best choices from
the results in Section 4.2, we restrict our analysis to the None, Smag., and
Grad. models. The Grad. model is much less computationally expensive than the
Dyn. Grad. model and is, therefore, the better choice for cosmological
simulations888Further study is required for tests that rely on getting the
small scale structure as correct as possible, such as cold clouds interacting
with a hot medium..

Figure 5: Temperature projections of the most massive halo in three of our
cosmological simulations at redshifts $z=2$, $1$, and $0$ in rows from top to
bottom, respectively. The columns show our None, Smag., and Grad. models from
left to right, respectively. Each of the panels represents a $1$
$\mathrm{Mpc}$ by $1$ $\mathrm{Mpc}$ (physical) region centred on the most
massive galaxy at each redshift. The Smag. simulation shows the most small-
scale structure at all redshifts, and a smoother distribution of temperature
at high-redshift compared to the None case. The Grad. model produces less
small-scale structure than any other model, and much more hot gas at high
redshift. Surprisingly, the Grad. model also produces more extended tails from
sub-structure moving through the hot halo at lower redshift, suggesting it may
impact future studies of jellyfish galaxies. Figure 6: Velocity magnitude
projections of the most massive halo at $z=0$ in our cosmological simulations.
From left to right we show the None, Smag., and Grad. simulations,
respectively. The images are centred on the most massive galaxy and each show
a $1$ $\mathrm{Mpc}$ by $1$ $\mathrm{Mpc}$ (physical) region. The Smag. model
shows more substructure than the None and Grad. models. The Grad. model shows
many extended tails from substructure moving through the hot halo compared to
the control and Smag. models, in addition to an overall smoother distribution
in velocity space.

Fig. 5 shows the density-weighted temperature projections of the most massive
galaxy at $z=2$, $1$, and $0$ in rows from top to bottom, respectively. The
columns show the None, Smag., and Grad. models from left to right,
respectively.

First, we compare the results at $z=2$ across mixing model variants. The dark
clumps in the None case are individual cold galaxies that are being fed into
the main structure via filaments. Ongoing stellar and AGN feedback leads to
the temperature increase in the central region, while the extended $T\sim
10^{6}$ $\mathrm{K}$ halo is a mixture of gravitational heated gas and gas
heated by previous generations of stellar feedback and AGN. There are many
filamentary structures visible feeding the main galaxy (centred), and there is
some small scale structure visible surrounding the galaxies. The Smag. case
resembles the None case, except there less cold gas in the upper in-falling
structure. Additionally, the satellite galaxies have more small-scale
structure in the cold gas surrounding them. Importantly, the distribution of
hotter gas appears mores smoothly distributed throughout the volume due to an
increase in conversion of kinetic energy in thermal energy via the eddy
viscosity.

The temperature projection of the Grad. model simulation is strikingly
different at $z=2$ than both the Smag. and None cases. While in the None and
Smag. cases there appears to be large-scale gas at $T\lesssim 10^{5}$
$\mathrm{K}$ at the boundary of the $500$ $\mathrm{kpc}$ radius, there is no
such gas in the Grad. case. There is much less small-scale structure in the
Grad. case, and the filaments are much smoothly distributed in space. We
conclude that the stellar feedback, gravitational heating, and AGN feedback
are much more effective at heating the gas on small scales in the Grad.
simulation, since the eddy viscosity is the strongest.

At $z=1$, in the middle row of Fig. 5, a similar picture emerges. The Smag.
case has much more small scale structure than the None case as is evident in
the central region and to the leftmost satellite galaxy structure. There is
still a much smoother distribution of cold gas that extends outwards from the
satellite galaxies. There are only a few centrally-located cold gas clumps in
the Grad. case compared with the Smag. case and there is effectively no
fragmentation identifiable in the satellite galaxy structure to the left of
panel.

The last row of Fig. 5 shows $z=0$ across the three model comparisons. By this
redshift, there is little structure remaining in the group-sized halo and the
temperature distribution of the intragroup gas is very smooth. Similarly to
the other redshifts, the Smag. shows the most fragmentation of cold gas,
followed by the None case, and then the Grad. case. All three show a satellite
being stripped of gas in the upper right of the panels but the Smag. and Grad.
models each differ in an important unique way compared to the None case.

It is much easier to see the differences in ram pressure stripping in velocity
space rather than temperature space. Fig. 6 shows the density-weighted
velocity magnitude projections of the most massive galaxy at $z=0$ for the
three eddy viscosity models None, Smag., and Grad. from left to right,
respectively. It is apparent across all three cases that the substructure is
moving at least a factor of $\approx 5$ faster than the background gas yet
only the Grad. case shows the most clearly defined long stripping tails from
the cold gas in the galaxies. The Smag. case does not have the cold gas in the
satellite galaxy marked by the arrows in the None and Grad. cases, as the gas
has been completely stripped away. The Grad. case produces the cold gas in
that satellite galaxy, and the tail is much more extended than in the None
case. In fact, it is evident upon close inspection of the cold gas structures
in the halo that the Grad. case produces tails from the cold gas in the
galaxies more readily than the None or Smag. cases. That has implications for
the study of ram pressure stripping in general, and should be further
investigated in the future.

##  5 Conclusions

Turbulence is a key physical process in the study galaxy evolution and one of
many highly complex non-linear interactions that must be understood to advance
our knowledge of the Universe. The complexity demands the use of simulations
that combine astrophysical sub-grid models with hydrodynamics and gravity in
an expanding universe. All hydrodynamical simulations are known to require
additional sub-grid models to accurately treat the impact sub-grid turbulence,
yet these models have been widely ignored in the astrophysical community. We
have, for the first time in Lagrangian hydrodynamics, implemented and studied
the gradient model (Clark et al., 1979) – an anisotropic sub-grid turbulence
model for viscosity and metal mixing (Hu & Chiang, 2020). The model is based
on directly modelling the error terms that arise from discretisation of the
fluid field via a Taylor series expansion, including the compression terms
that are missing from the standard Smagorinsky model. We additionally
implemented a dynamic procedure that computes the model parameter on-the-fly
for the gradient model following the approach of Rennehan et al. (2019). We
used the mesh-free finite mass method in the GIZMO code (Hopkins, 2015) as our
numerical hydrodynamics solver for all of our experiments.

We ran driven turbulence simulations at Mach numbers $\mathcal{M}=0.3$, $0.7$,
and $2.1$ to validate the gradient model and compare with the popular
Smagorinsky model. Hu & Chiang (2020) recently showed, by post-processing
driven their turbulence simulations, that the gradient model should be able to
reduce the build-up of kinetic energy near the resolution scale in isotropic,
homogeneous turbulence and better reproduce the sub-grid metal flux. We
confirmed these results in Section 3 by using the gradient model at simulation
time.

Our analysis of the velocity power spectra in driven turbulence produced
unexpected results. We found that the gradient and Smagorinsky models, and
their dynamic variants, predicted insufficient dissipation to reduce the
artificial build-up of kinetic energy near the grid-scale in our $256^{3}$
simulations. For that reason we introduced a boost factor $\gamma$ for the
dissipation strength (i.e.
$||\boldsymbol{K}||\rightarrow\gamma||\boldsymbol{K}||$) and experimented with
various values in the range $\gamma\in[1,100]$. We found that a factor of
$\gamma\sim 10$ is sufficient to reduce the build-up of kinetic energy and is
required for all of the Smagorinsky and gradient model variants. Additionally,
we found that the boost factor only needs to be applied to the subsonic
($\mathcal{M}<1$) particles in our simulations to produce the correct
statistics in supersonic turbulence (our $\mathcal{M}\sim 2.1$ test). The true
boost factor is higher since we used the maximum interaction distance between
neighbouring gas particles to calculate the diffusion tensor, leading to a
$\sim 4$ times additional boost (total $\sim 40$) over other the default
implementation in GIZMO (Hopkins et al., 2018).

Our converged metal mixing simulations in Section 3.2 show that when we use
the gradient and Smagorinsky models at lower resolution ($64^{3}$ particles)
we are able to produce MDFs that are equivalent to $4$ to $12$ times the
resolution. That is true for both the constant-coefficient and dynamic
variants of the gradient and Smagorinsky models with standard parameter
values, with the additional factor of $\sim 4$ boost from using the maximum
interaction distance in the kernel. However, lower calibrations such as those
from the FIRE simulations do not produce the correct rate of mixing as they
are at least a factor of $\sim 20$ too low. We posit that this is because of
the common calibration approach in cosmological simulations: calibration in
tandem with the full suite of astrophysical sub-grid models.

We argue that calibration of fundamental hydrodynamics models, such as the
metal mixing model here, _must be done in the absence of sub-grid
astrophysical models_. There must be a strong hydrodynamics base before the
complexity of astrophysics is built on top. Note that our dynamic Smagorinsky
and dynamic gradient models _do not require calibration_ and produce accurate
predictions for the rate of metal mixing in isotropic, homogeneous turbulence.

As a first application of the new gradient model, we investigated a set of
cosmological simulations to determine if there is any dominant impact on the
galaxy and gas properties. We investigated the metal mixing and eddy viscosity
separately in Sections 4.2 & 4.3, respectively.

We found that the choice of sub-grid metal mixing model strongly impacts the
MDF evolution in the warm-hot intergalactic medium (WHIM) and circumgalactic
medium (CGM). We found that the gradient and dynamic gradient models mix
metals much more rapidly than the Smagorinsky variants and produce tighter
MDFs up until $z\sim 1$ when they approach a similar distribution down to
$z=0$. In our simulation without sub-grid metal mixing, the MDFs in the WHIM
and CGM are significantly broader than any of the simulations with sub-grid
metal mixing demonstrating that, at the very minimum, including any sub-grid
metal mixing model is an improvement. The most important result we discovered
is that the metal mixing models are most impactful during the tempestuous
early stages of galaxy evolution. On very long timescales, the equilibrium
distributions match quite closely across the models.

Including eddy viscosity in our cosmological simulations did not significantly
impact the galaxy properties we investigated when averaging in bins of stellar
mass after $z=2$. The galaxy stellar mass function was relatively unchanged,
along with only slight variations in the stellar mass to halo mass function.
We also found that including eddy viscosity did not significantly impact the
averaged gas distributions of vorticity, temperature, and density across
galaxies of similar stellar mass. We expected a priori that the large-scale
properties of galaxies would be unaffected as the sub-grid eddy viscosity
mainly impacts the small-scale. For that reason, we investigated a single
system that could be linked across all of our cosmological simulations to gain
a qualitative view of the impact.

In Section 4.3 we showed the temperature projections of the most massive halo
traced from $z=2$, $1$ and $0$ in our cosmological simulations for three eddy
viscosity models: no model, the standard Smagorinsky model, and the new
constant-coefficient gradient model. We found that the Smagorinsky model
produced much more fragmentation in the halo gas of the most massive galaxy on
the small scale compared to having no eddy viscosity, at all redshifts. We
also found that the spatial temperature distribution was much smoother at
$z=2$ when stellar and active galactic nuclei feedback was much stronger,
showing that the small-scale kinetic energy was being efficiently converted
into thermal energy. Although the constant-coefficient gradient model
seemingly dissipates faster based on our results in Section 3.1, we observed
that its inherent anisotropy does not lead to the same fragmentation we saw in
the Smagorinsky model. At high redshift, $z=2$ the gradient model produced a
much more widespread hot gas. The filamentary structure at all redshifts was
much smoother and, after $z=1$, the satellite galaxies in the halo had many
more clearly defined tails due to an improved treatment of ram-pressure
stripping.

Sub-grid metal mixing and eddy viscosity models have a strong impact on galaxy
evolution simulations. In this work, we showed in the simplest case of
isotropic, homogeneous turbulence that the all of the models tested here
improved the accuracy of metal mixing and turbulent kinetic energy dissipation
in the mesh-free finite mass method. The most significant differences between
model choice appeared at high redshift in the early stages of galaxy
evolution, before any equilibrium is reached. Given that contemporary
cosmological simulations have resolutions of less than $\sim 50^{3}$ particles
in a typical galaxy, we recommend that future studies must at least use the
constant-coefficient gradient model as it is (a) computationally inexpensive
compared to the dynamic version, while producing similar results and (b)
includes the full velocity tensor in the diffusion tensor to give the most
accurate solution for sub-grid turbulence. Our recommendation is especially
pertinent given the recent push to study higher redshift systems driven by the
upcoming launch of the James Webb Space Telescope, as our theoretical
understanding of enrichment and thermodynamic histories will depend directly
on sub-grid turbulence model choice.

## Acknowledgements

This research was enabled in part by support provided by WestGrid and
Compute/Calcul Canada. The simulations in this research were made possible by
SciNet and the Niagara supercomputing cluster. DR acknowledges the support of
the Natural Sciences and Engineering Research Council of Canada (NSERC),
[funding reference number 534263] and through the Discovery Grant program. DR
thanks Arif Babul, Belaid Moa, Ondrea Clarkson, Austin Davis, Drummond
Fielding, Phil Hopkins, and Wolfram Schmidt for useful advice during the
course of this research. DR gives special acknowledgement to Arif Babul and
Belaid Moa for providing invaluable support to this research, without which it
would not have been possible.

Our analysis was performed using the Python programming language (Python
Software Foundation, https://www.python.org). The following packages were used
throughout the analysis: h5py (Collette, 2013), numpy (Harris et al., 2020),
scipy (Virtanen et al., 2020), yt (Turk et al., 2011), and matplotlib (Hunter,
2007). Prototyping of the analysis scripts was performed in the IPython
environment (Perez & Granger, 2007).

## Data availability

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

## References

  * Agertz et al. (2007) Agertz O., et al., 2007, Monthly Notices of the Royal Astronomical Society, 380, 963
  * Anderson & Bregman (2010) Anderson M. E., Bregman J. N., 2010, The Astrophysical Journal, 714, 320
  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, Annual Review of Astronomy and Astrophysics, 47, 481
  * Balarac et al. (2013) Balarac G., Le Sommer J., Meunier X., Vollant A., 2013, Physics of Fluids, 25
  * Bauer & Springel (2012) Bauer A., Springel V., 2012, Monthly Notices of the Royal Astronomical Society, 423, 2558
  * Bennett & Sijacki (2020) Bennett J. S., Sijacki D., 2020, Monthly Notices of the Royal Astronomical Society, 499, 597
  * Bourne & Sijacki (2017) Bourne M. A., Sijacki D., 2017, Monthly Notices of the Royal Astronomical Society, 472, 4707
  * Brook et al. (2014) Brook C. B., Stinson G., Gibson B. K., Shen S., Macciò A. V., Obreja A., Wadsley J., Quinn T., 2014, Monthly Notices of the Royal Astronomical Society, 443, 3809
  * Buie et al. (2020) Buie E., Fumagalli M., Scannapieco E., 2020, The Astrophysical Journal, 890, 33
  * Chandrasekhar (1961) Chandrasekhar S., 1961, Hydrodynamic and hydromagnetic stability. Oxford University Press, London and New York
  * Clark et al. (1979) Clark R. A., Ferziger J. H., Reynolds W. C., 1979, Journal of Fluid Mechanics, 91, 1
  * Colbrook et al. (2017) Colbrook M. J., Ma X., Hopkins P. F., Squire J., 2017, Monthly Notices of the Royal Astronomical Society, 467, 2421
  * Collette (2013) Collette A., 2013, Python and HDF5. O’Reilly
  * Davé et al. (2010) Davé R., Oppenheimer B. D., Katz N., Kollmeier J. A., Weinberg D. H., 2010, Monthly Notices of the Royal Astronomical Society, 408, 2051
  * Davé et al. (2016) Davé R., Thompson R., Hopkins P. F., 2016, Monthly Notices of the Royal Astronomical Society, 462, 3265
  * Davé et al. (2017) Davé R., Rafieferantsoa M. H., Thompson R. J., Hopkins P. F., 2017, Monthly Notices of the Royal Astronomical Society, 467, stx108
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, Monthly Notices of the Royal Astronomical Society, 486, 2827
  * Dekel et al. (2009) Dekel A., Sari R., Ceverino D., 2009, The Astrophysical Journal, 703, 785
  * Di Mascio et al. (2017) Di Mascio A., Antuono M., Colagrossi A., Marrone S., 2017, Physics of Fluids, 29, 035102
  * Elmegreen & Scalo (2004) Elmegreen B. G., Scalo J., 2004, Annu. Rev. Astron. Astrophys, 42, 211
  * Escala et al. (2018) Escala I., et al., 2018, Monthly Notices of the Royal Astronomical Society, 474, 2194
  * Federrath (2013) Federrath C., 2013, Monthly Notices of the Royal Astronomical Society, 436, 1245
  * Fielding et al. (2017) Fielding D., Quataert E., McCourt M., Thompson T. A., 2017, Monthly Notices of the Royal Astronomical Society, 466, 3810
  * Fielding et al. (2018) Fielding D., Quataert E., Martizzi D., 2018, Monthly Notices of the Royal Astronomical Society, 481, 3325
  * Gaburov & Nitadori (2011) Gaburov E., Nitadori K., 2011, Monthly Notices of the Royal Astronomical Society, 414, 129
  * Garnier et al. (2009) Garnier E., Adams N., Sagaut P., 2009, Large Eddy Simulation for Compressible Flows. Scientific Computation, Springer Netherlands, Dordrecht, doi:10.1007/978-90-481-2819-8
  * Genel et al. (2014) Genel S., et al., 2014, Monthly Notices of the Royal Astronomical Society, 445, 175
  * Gingold & Monaghan (1977) Gingold R. A., Monaghan J. J., 1977, Monthly Notices of the Royal Astronomical Society, 181, 375
  * Greif et al. (2009) Greif T. H., Glover S. C. O., Bromm V., Klessen R. S., 2009, Monthly Notices of the Royal Astronomical Society, 392, 1381
  * Guedes et al. (2011) Guedes J., Callegari S., Madau P., Mayer L., 2011, The Astrophysical Journal, 742, 76
  * Hafen et al. (2019) Hafen Z., et al., 2019, Monthly Notices of the Royal Astronomical Society, 488, 1248
  * Hafen et al. (2020) Hafen Z., et al., 2020, Monthly Notices of the Royal Astronomical Society, 494, 3581
  * Hahn & Abel (2011) Hahn O., Abel T., 2011, Monthly Notices of the Royal Astronomical Society, 415, 2101
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hernquist & Katz (1989) Hernquist L., Katz N., 1989, The Astrophysical Journal Supplement Series, 70, 419
  * Hopkins (2013) Hopkins P. F., 2013, Monthly Notices of the Royal Astronomical Society, 428, 2840
  * Hopkins (2015) Hopkins P. F., 2015, Monthly Notices of the Royal Astronomical Society, 450, 53
  * Hopkins (2017) Hopkins P. F., 2017, Monthly Notices of the Royal Astronomical Society, 466, 3387
  * Hopkins et al. (2014) Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, Monthly Notices of the Royal Astronomical Society, 445, 581
  * Hopkins et al. (2018) Hopkins P. F., et al., 2018, Monthly Notices of the Royal Astronomical Society, 480, 800
  * Hu & Chiang (2020) Hu C.-Y., Chiang C.-T., 2020, The Astrophysical Journal, 900, 29
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Iapichino et al. (2011) Iapichino L., Schmidt W., Niemeyer J. C., Merklein J., 2011, Monthly Notices of the Royal Astronomical Society, 414, 2297
  * Iapichino et al. (2013) Iapichino L., Viel M., Borgani S., 2013, Monthly Notices of the Royal Astronomical Society, 432, 2529
  * Karen Yang & Reynolds (2016) Karen Yang H.-Y., Reynolds C. S., 2016, The Astrophysical Journal, 829, 90
  * Khatri & Gaspari (2016) Khatri R., Gaspari M., 2016, Monthly Notices of the Royal Astronomical Society, 463, 655
  * Landau & Lifshitz (1987) Landau L., Lifshitz E., 1987, Course of theoretical physics. Vol. 6: Fluid Mechanics. Pergamon Press
  * Lanson & Vila (2008a) Lanson N., Vila J.-P., 2008a, SIAM Journal on Numerical Analysis, 46, 1912
  * Lanson & Vila (2008b) Lanson N., Vila J.-P., 2008b, SIAM Journal on Numerical Analysis, 46, 1935
  * Li et al. (2020) Li M., Li Y., Bryan G. L., Ostriker E. C., Quataert E., 2020, The Astrophysical Journal, 898, 23
  * Lochhaas et al. (2020) Lochhaas C., Bryan G. L., Li Y., Li M., Fielding D., 2020, Monthly Notices of the Royal Astronomical Society, 493, 1461
  * Lucy (1977) Lucy L. B., 1977, The Astronomical Journal, 82, 1013
  * Monaghan (1989) Monaghan J. J., 1989, Journal of Computational Physics, 82, 1
  * Monaghan (2002) Monaghan J. J., 2002, Monthly Notices of the Royal Astronomical Society, 335, 843
  * Monaghan (2011) Monaghan J. J., 2011, European Journal of Mechanics, B/Fluids, 30, 360
  * Naab & Ostriker (2017) Naab T., Ostriker J. P., 2017, Annual Review of Astronomy and Astrophysics, 55, 59
  * Pan et al. (2013) Pan L., Scannapieco E., Scalo J., 2013, The Astrophysical Journal, 775, 111
  * Perez & Granger (2007) Perez F., Granger B. E., 2007, Computing in Science & Engineering, 9, 21
  * Pillepich et al. (2018) Pillepich A., et al., 2018, Monthly Notices of the Royal Astronomical Society, 475, 648
  * Pinto et al. (2015) Pinto C., et al., 2015, Astronomy & Astrophysics, 575, A38
  * Piomelli & Liu (1995) Piomelli U., Liu J., 1995, Physics of Fluids, 7, 839
  * Poole et al. (2006) Poole G. B., Fardal M. A., Babul A., McCarthy I. G., Quinn T., Wadsley J., 2006, Monthly Notices of the Royal Astronomical Society, 373, 881
  * Pope (2000) Pope S. B., 2000, Turbulent Flows, doi:10.1017/CBO9780511840531
  * Prasad et al. (2015) Prasad D., Sharma P., Babul A., 2015, The Astrophysical Journal, 811, 108
  * Prasad et al. (2018) Prasad D., Sharma P., Babul A., 2018, The Astrophysical Journal, 863, 62
  * Price (2012) Price D. J., 2012, Monthly Notices of the Royal Astronomical Society: Letters, 420, 33
  * Rennehan et al. (2019) Rennehan D., Babul A., Hopkins P. F., Davé R., Moa B., 2019, Monthly Notices of the Royal Astronomical Society, 483, 3810
  * Ruggiero & Lima Neto (2017) Ruggiero R., Lima Neto G. B., 2017, Monthly Notices of the Royal Astronomical Society, 468, 4107
  * Sagaut (2006) Sagaut P., 2006, Large Eddy Simulation for Incompressible Flows. Scientific Computation, Springer-Verlag, Berlin/Heidelberg, doi:10.1007/b137536
  * Scalo & Elmegreen (2004) Scalo J., Elmegreen B. G., 2004, Annual Review of Astronomy and Astrophysics, 42, 275
  * Scannapieco & Brüggen (2008) Scannapieco E., Brüggen M., 2008, The Astrophysical Journal, 686, 927
  * Schaye et al. (2015) Schaye J., et al., 2015, Monthly Notices of the Royal Astronomical Society, 446, 521
  * Schmidt (2015) Schmidt W., 2015, Living Reviews in Computational Astrophysics, 1, 64
  * Schmidt et al. (2014) Schmidt W., et al., 2014, Monthly Notices of the Royal Astronomical Society, 440, 3051
  * Semenov et al. (2016) Semenov V. A., Kravtsov A. V., Gnedin N. Y., 2016, The Astrophysical Journal, 826, 200
  * Shen et al. (2010) Shen S., Wadsley J., Stinson G., 2010, Monthly Notices of the Royal Astronomical Society, 407, 1581
  * Shen et al. (2012) Shen S., Madau P., Aguirre A., Guedes J., Mayer L., Wadsley J., 2012, The Astrophysical Journal, 760, 50
  * Shen et al. (2013) Shen S., Madau P., Guedes J., Mayer L., Prochaska J. X., Wadsley J., 2013, The Astrophysical Journal, 765, 89
  * Simons et al. (2020) Simons R. C., et al., 2020, The Astrophysical Journal, 905, 167
  * Smagorinsky (1963) Smagorinsky J., 1963, Monthly Weather Review, 91, 99
  * Sokołowska et al. (2018) Sokołowska A., Babul A., Mayer L., Shen S., Madau P., 2018, The Astrophysical Journal, 867, 73
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, Annu. Rev. Astron. Astrophys, 53, 51
  * Springel (2010) Springel V., 2010, Monthly Notices of the Royal Astronomical Society, 401, 791
  * Su et al. (2017) Su K.-Y., Hopkins P. F., Hayward C. C., Faucher-Giguère C.-A., Kereš D., Ma X., Robles V. H., 2017, Monthly Notices of the Royal Astronomical Society, 471, 144
  * Tremmel et al. (2017) Tremmel M., Karcher M., Governato F., Volonteri M., Quinn T. R., Pontzen A., Anderson L., Bellovary J., 2017, Monthly Notices of the Royal Astronomical Society, 470, 1121
  * Tremmel et al. (2019) Tremmel M., et al., 2019, Monthly Notices of the Royal Astronomical Society, 483, 3336
  * Tumlinson et al. (2013) Tumlinson J., et al., 2013, The Astrophysical Journal, 777, 59
  * Tumlinson et al. (2017) Tumlinson J., Peeples M. S., Werk J. K., 2017, Annual Review of Astronomy and Astrophysics, 55, 389
  * Turk et al. (2011) Turk M. J., Smith B. D., Oishi J. S., Skory S., Skillman S. W., Abel T., Norman M. L., 2011, The Astrophysical Journal Supplement Series, 192, 9
  * Vazza et al. (2010) Vazza F., Brunetti G., Gheller C., Brunino R., 2010, New Astronomy, 15, 695
  * Vazza et al. (2012) Vazza F., Roediger E., Brüggen M., 2012, Astronomy & Astrophysics, 544, A103
  * Vazza et al. (2017) Vazza F., Jones T. W., Brüggen M., Brunetti G., Gheller C., Porter D., Ryu D., 2017, Monthly Notices of the Royal Astronomical Society, 464, 210
  * Vazza et al. (2018) Vazza F., Angelinelli M., Jones T. W., Eckert D., Brüggen M., Brunetti G., Gheller C., 2018, Monthly Notices of the Royal Astronomical Society: Letters, 481, L120
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, Monthly Notices of the Royal Astronomical Society, 444, 1518
  * Wadsley et al. (2008) Wadsley J. W., Veeravalli G., Couchman H. M. P., 2008, Monthly Notices of the Royal Astronomical Society, 387, 427
  * Wadsley et al. (2017) Wadsley J. W., Keller B. W., Quinn T. R., 2017, Monthly Notices of the Royal Astronomical Society, 471, 2357
  * Wang et al. (2020) Wang C., Ruszkowski M., Pfrommer C., Oh S. P., Yang H. Y. K., 2020, arXiv preprint arXiv:2012.11085
  * Werk et al. (2016) Werk J. K., et al., 2016, The Astrophysical Journal, 833, 1
  * Williamson et al. (2016) Williamson D., Martel H., Kawata D., 2016, The Astrophysical Journal, 822, 91
  * Wittor et al. (2017) Wittor D., Jones T., Vazza F., Brüggen M., 2017, Monthly Notices of the Royal Astronomical Society, 471, 3212
  * Zhuravleva et al. (2014) Zhuravleva I., et al., 2014, Nature, 515, 85
  * Zhuravleva et al. (2015) Zhuravleva I., et al., 2015, Monthly Notices of the Royal Astronomical Society, 450, 4184
  * Zhuravleva et al. (2018) Zhuravleva I., Allen S. W., Mantz A., Werner N., 2018, The Astrophysical Journal, 865, 53
  * ZuHone et al. (2013) ZuHone J. A., Markevitch M., Ruszkowski M., Lee D., 2013, The Astrophysical Journal, 762, 69

##  Appendix A FILTERING APPROXIMATION

Consider a scalar field $f_{i}(\boldsymbol{r})$ where $\boldsymbol{r}=(x,y,z)$
in a general coordinate system. We expand $f_{i}(\boldsymbol{r}^{\prime})$
about $\boldsymbol{r}$ in equation (1) via a Taylor expansion,

$\begin{split}f_{i}(\boldsymbol{r}^{\prime})=f_{i}(\boldsymbol{r})+[(\Delta\boldsymbol{R})\cdot\boldsymbol{\nabla}f(\boldsymbol{r})]+\frac{1}{2}[(\Delta\boldsymbol{R})\otimes(H(\boldsymbol{r})\cdot(\Delta\boldsymbol{R}))]\\\
+\mathcal{O}(|\Delta\boldsymbol{R}|^{3}),\end{split}$ (27)

where $\Delta\boldsymbol{R}=\boldsymbol{r}^{\prime}-\boldsymbol{r}$ and
$H(\boldsymbol{r})$ is the Hessian matrix. Putting this into equation 1 gives,

$\begin{split}\overline{f}_{i}(\boldsymbol{r})=f_{i}(\boldsymbol{r})+\frac{1}{2}H(\boldsymbol{r})\int_{\mathrm{kern}}((\Delta\boldsymbol{R})\otimes(\Delta\boldsymbol{R}))G(|\Delta\boldsymbol{R}|,h)d\boldsymbol{r}^{\prime}\\\
+\mathcal{O}(|\Delta\boldsymbol{R}|^{4}).\end{split}$ (28)

using the fact that the kernel function is normalized and the integral of an
odd function over the domain is zero. The integral in equation 28 can be
tabulated for a specific function $G=G(|\Delta\boldsymbol{R}|,h)$, so let’s
define,

$\epsilon(\boldsymbol{r})\equiv\frac{1}{2}\int_{\mathrm{kern}}((\Delta\boldsymbol{R})\otimes(\Delta\boldsymbol{R}))G(|\Delta\boldsymbol{R}|,h)d\boldsymbol{r}^{\prime}.$
(29)

We then have,

$\overline{f}_{i}(\boldsymbol{r})\approx
f_{i}(\boldsymbol{r})+\epsilon(\boldsymbol{r})\cdot H(\boldsymbol{r}).$ (30)

However, $\epsilon(\boldsymbol{r})$ is isotropic
$\epsilon(\boldsymbol{r})=\epsilon\,\boldsymbol{I}$, so we can write

$\overline{f}_{i}(\boldsymbol{r})\approx
f_{i}(\boldsymbol{r})+\epsilon\boldsymbol{\nabla}^{2}f_{i}(\boldsymbol{r}).$
(31)

##  Appendix B METAL DISTRIBUTIONS

Fig. 7 shows the normalised histograms of the filtered metal field from our
driven turbulence simulations in Section 3.2. The panels are Mach numbers
$\mathcal{M}=0.3$, $0.7$, and $2.1$ from top to bottom, respectively. The
black curves with markers show the convergence of the filtered metal field for
resolutions $64^{3}$, $128^{3}$, $256^{3}$, $512^{3}$, and $768^{3}$. The
coloured lines show, from lightest to darkest: Dyn. Smag (solid salmon), Smag.
(solid magenta), FIRE (dotted magenta), Dyn. Grad. (dashed purple), and Grad.
(solid purple) at $64^{3}$ resolution. We obtained all of the information for
this figure after approximately two mixing timescales, $t\sim
2\,\tau_{\mathrm{mix}}$, where $\tau_{\mathrm{mix}}=1/\mathcal{M}$ (see
Section 3.2 for more details). Fig. 3 we used the Gaussian fits to these
distributions and found the widths of the distribution, $\sigma_{\mathrm{Z}}$.

Figure 7: The normalised histograms of the filtered metal mass fraction field
from the simulations in Section 3.2. The panels show Mach numbers
$\mathcal{M}=0.3$, $0.7$, and $2.1$ from top to bottom, respectively after two
mixing timescales. The black curves show the simulations with no sub-grid
metal mixing. The coloured curves show the simulations at $64^{3}$ resolution
with a sub-grid metal mixing model given by the label. All of the sub-grid
metal mixing models show improvement except for FIRE which lags due to the
lower calibration coefficient.