Mean fields and fluctuations in compressible magnetohydrodynamic flows

We apply Gaussian smoothing to obtain mean magnetic field, density, velocity, and magnetic and kinetic energy densities from our numerical model of the interstellar medium, based on three-dimensional magnetohydrodynamic equations in a shearing box in size. The interstellar medium is highly compressible, as the turbulence is transonic or supersonic; it is thus an excellent context in which to explore the use of smoothing to represent physical variables in a compressible medium in terms of their mean and fluctuating parts. Unlike alternative averaging procedures, such as horizontal averaging, Gaussian smoothing retains the three-dimensional structure of the mean fields. Although Gaussian smoothing does not obey the Reynolds rules of averaging, physically meaningful and mathematically sound central statistical moments are defined as suggested by Germano [Turbulence – The filtering approach. J. Fluid Mech. 1992, 238, 325–336]. We discuss methods to identify an optimal smoothing scale ℓ and the effects of this choice on the results. From spectral analysis of the magnetic, density and velocity fields, we find a suitable smoothing length for all three fields, of . Such a smoothing scale is likely to be sensitive to the choice of simulation parameters; this may be considered in future work, but here we just explore the methodology. We discuss the properties of third-order statistical moments in fluctuations of kinetic energy density in compressible flows, and suggest their physical interpretation. The mean magnetic field, amplified by a mean-field dynamo, significantly alters the distribution of kinetic energy in space and between scales, reducing the magnitude of kinetic energy at intermediate scales. This intermediate-scale kinetic energy is a useful diagnostic of the importance of SN-driven outflows.


Introduction
The injection of thermal and kinetic energy by stellar winds and supernova (SN) explosions drives transonic turbulence in the interstellar medium (ISM) and produces an inhomogeneous, multiphase system (Elmegreen and Scalo 2004, Scalo and Elmegreen 2004, Mac Low and Klessen 2004).The outer scale of turbulent motions in the ISM consistently suggested by observations, theory and simulations is of order 10-100 pc, and the turbulent scales extend to a fraction of a parsec (Armstrong et al. 1995).
Understanding the properties and nature of a turbulent flow requires the separation of mean and fluctuating quantities.Such a separation is well understood for statistically homogeneous random flows where a number of averaging procedures are available.Volume or area averaging are most important in astronomy, while numerical simulations provide a further opportunity to average over time.Under favourable conditions (defined by ergodic theorems and hypotheses), the resulting averages are equivalent to the statistical ensemble averages employed in theory (e.g., Monin and Yaglom 2007a, Panchev 1971, Tennekes and Lumley 1972).The ensemble averages are rarely accessible in applications, as their calculation requires the availability of a large number of statistically independent realizations of the random processes.
Space and time averaging procedures are consistent with ensemble averaging f + g = f + g , f g = f g and f = f , where f and g are random functions and angular brackets denote averaging (e.g., Sect.3.1 in Monin and Yaglom 2007a).Volume and time averaging only satisfy the Reynolds rules in an approximate manner when the scales of variations of the mean quantities and the fluctuations differ significantly (the requirement of scale separation between the averaged quantities and the fluctuations) and the averaging scale is large in comparison with the scale of the fluctuations and small in comparison with that of the mean quantities.In practice, the mean quantities need to be homogeneous or time-independent for the ensemble and volume (or time) averages to be consistent with each other.
The outer scale of the interstellar turbulence is comparable to the scale height of the gas density distribution in spiral galaxies (about 0.1 kpc and 0.5 kpc for the cold and warm diffuse HI, respectively).Therefore, the interstellar turbulent flow cannot be considered statistically homogeneous apart from along the horizontal directions.However, numerical simulations of the supernova-driven, multi-phase ISM have relatively small horizontal domains of order 1 kpc × 1 kpc or less (e.g., Korpi et al. 1999, Joung andMac Low 2006, de Avillez and Breitschwerdt 2007, 2012a,b, Gressel et al. 2008, Federrath et al. 2010, Hill et al. 2012, Gent et al. 2013b,a, Gressel et al. 2013, Bendre et al. 2015, Walch et al. 2015, Girichidis et al. 2016b,a, 2018).Meanwhile, the ISM has a wide range of density and velocity structures (e.g., those related to gas clouds, galactic outflows and spiral patterns) that cover continuously the range of scales from 1 pc to 10 kpc.Therefore, scale separation between the random and large-scale ISM flows is questionable at best.This poses difficulties for the interpretation of numerical simulations.Similar difficulties arise in the interpretation of observations, but numerical simulations have exposed the problems especially clearly.
MHD simulations of galaxies employ domains of scale 10 kpc and larger, and contain kiloparsecscale structures that would be associated with mean-field dynamics (e.g., Slyz et al. 2003, Dobbs and Price 2008, Hanasz et al. 2009, Kulesza-Żydzik et al. 2009, Dobbs and Pringle 2010, Pakmor et al. 2016, 2017, Rieder and Teyssier 2016, 2017).Such simulations can also support a small-scale magnetic field driven by dynamo action (Pakmor et al. 2016, 2017, Rieder and Teyssier 2016, 2017).However, since computational limits restrict the resolution of such simulations to order 10 pc, these simulations are unable to fully model physics such as the injection of SN energy into the ISM and must rely on a simplified prescription.Thus, whilst scale separation may be obtained and lead to useful results for the magnetic field, such results for the velocity field may not be reliable.
Simulations on the scale of interstellar clouds have provided powerful insight into the parsec and sub-parsec physics in the ISM (e.g., Klessen et al. 2000, Heitsch et al. 2001, Brandenburg et al. 2007).Such simulations drive turbulence at specific wavenumbers to simulate ISM turbulence.Whilst scale separation may be possible, it would be a separation between the parameterised injected flow and any resultant smaller-scale flows.Thus, the results from these scale separations may not be fully informative.
The division of the Navier-Stokes and magnetohydrodynamic (MHD) equations into evolution equations for the mean flow and the fluctuations has been explored for both ensemble averaging and filtering of the fluctuations (also known as coarse-graining); i.e., volume averaging via convolution with a compact kernel.The Reynolds rules of averaging are not satisfied for this procedure but this is not an obstacle to developing a mathematically sound formalism that leads to evolution equations for averaged quantities and their moments (Germano 1992).The most widely known application of this technique is to subgrid models for large eddy simulations of turbulent flows (Meneveau 2012).Eyink (2018) and Aluie (2017) provide details and a review of this approach to hydrodynamic and MHD turbulence (see also Eyink 1995Eyink , 2015)).
An important advantage of the filtering approach is that, together with ensemble averaging, it does not require scale separation between the mean fields and their fluctuations (e.g., Aluie 2017).The separation of the mean and fluctuating quantities in a random flow is of crucial significance in the theory of mean-field turbulent dynamos.Mean-field dynamo theory is based on ensemble averaging, but numerical simulations rely on various volume and time averaging procedures.(We note that the algebraic form of the mean-field equations is the same for a wide range of averaging methods: Germano 1992, Aluie 2017.)For example, the separation of the magnetic field into mean and fluctuating components often involves averaging over the whole computational volume or, in systems stratified along the z-direction due to gravity, averaging in the (x, y)-planes (horizontal averaging; see Brandenburg and Subramanian 2005).The resulting mean magnetic field is either perfectly uniform or only dependent on z.These constraints on the form of the mean magnetic field are often artificial and unphysical.An inhomogeneous system, such as the ISM, is expected to produce a spatially complex mean field, which is ignored by these simple volume or horizontal averaging techniques.A further complication with horizontal averaging, when periodic boundary conditions are used in x and y, is that B z must vanish to guarantee the solenoidality of the mean magnetic field (e.g., Gent et al. 2013a).The main advantage of horizontal averaging is to obey the Reynolds rules, achieved often at the expense of physical validity.Another, Reynolds rules compliant, option is to use azimuthal averaging in global simulations of dynamo action in a rotating spherical object to obtain an axially symmetric mean magnetic field (see Simard et al. 2016, for a review).This approach is easier to justify but still it excludes physically admissible azimuthal variations of the mean field.Furthermore, the kinematic mean-field dynamo action, with homogeneous transport coefficients α and β (representing the αeffect and turbulent diffusion), in infinite space produces an inhomogeneous mean magnetic field that varies at all wavenumbers below α/β, with the dominant mode having the wavenumber α/(2β) (e.g., Sokoloff et al. 1983).The spatial structure of any mean field is controlled by the physical properties of the system, rather than by the size of the computational domain.
We discuss an alternative approach to averaging based on Gaussian smoothing as suggested by Germano (1992), and employ it to obtain the mean fields in simulations of the multi-phase, supernovadriven ISM.Averaging with a Gaussian (or another) kernel is inherent in astronomical observations, where such smoothing is applied either during data reduction or stems from the finite width of a telescope beam.This approach has been applied by Gent et al. (2013a) to the simulated magnetic field; here we extend it to the velocity and density fields and, importantly, energy densities, which represent higher-order statistical moments.In particular, kinetic energy density in a compressible flow represents a third-order statistical moment and requires special attention.
A detailed exploration of the sensitivity of our results to the parameter space is beyond the scope of this paper.Rather, we seek to demonstrate that the extension of this approach to also concern the density and velocity, and thus also kinetic energy, does indeed lead to mathematically and physically meaningful results.A summary of our numerical model of the ISM is presented in Section 2, and Section 3 introduces averaging based on Gaussian smoothing.Various approaches to the selection of the smoothing length are discussed in Section 4. Section 5 analyses the behaviour of magnetic and kinetic energy densities.Section 6 details the effects of the amplified mean field on the magnetic and kinetic energies.Section 7 compares Gaussian smoothing with horizontal averaging, to show the advantages of the former in the current context.Finally, section 8 summarises and concludes our discussion.

A numerical model of the multiphase ISM
We use our earlier numerical model of the ISM, described in detail by Gent et al. (2013a,b).The model involves solving, with the Pencil Code (Pencil Code Collaboration et al. 2021), the full, compressible, non-ideal MHD equations with parameters generally typical of the Solar neighbourhood in a threedimensional Cartesian, shearing box with radial (x) and azimuthal (y) extent L x = L y = 1.024 kpc and vertical (z) extent L z = 1.086 kpc on either side of the midplane at z = 0. Our numerical resolution is ∆ = ∆x = ∆y = ∆z = 4 pc, using 256 grid points in x and y and 544 in z.Gent et al. (2013a) and, in greater detail, Gent et al. (2020) demonstrate that this resolution is sufficient to reproduce the known solutions for expanding SN remnants in the Sedov-Taylor and momentum-conserving phases.Details of the numerical implementation and its comparison with some other similar simulations can be found in Appendix A.
The mass conservation, Navier-Stokes, heat and induction equations are solved for mass density ρ, velocity u, specific entropy s, and magnetic vector potential A (such that B = ∇ × A).The Navier-Stokes equation includes a fixed vertical gravity force with contributions from the stellar disc and dark halo.The initial state is an approximate hydrostatic equilibrium.The Galactic differential rotation is modelled by a background shear flow U = (0, −qΩx, 0), where q is the shear parameter and Ω is the Galactic angular velocity.Here we use q = 1, as in a flat rotation curve (i.e., with the rotation speed independent of the cylindrical radius), and Ω = 50 km s −1 kpc −1 , twice that of the Solar neighbourhood in order to enhance the mean-field dynamo action and thus reduce the computational time.The velocity u is the perturbation velocity in the rotating frame, that remains after the subtraction of the background shear flow from the total velocity.However, it still contains a large-scale vertical component due to an outflow driven by the SN activity.
Both Type II and Type I supernova explosions (SNe) are included in the simulations.These differ only in their vertical distribution and frequency.The frequencies used correspond to those in the Solar neighbourhood.We introduce Type II SNe at a mean rate per surface area of ν II = 25 kpc −2 Myr −1 .Type I SNe have a mean rate per surface area of ν I = 4 kpc −2 Myr −1 .The SN sites have uniform random distribution in the horizontal plane.Their vertical positions have Gaussian distributions with scale heights h II = 0.09 kpc and h I = 0.325 kpc for SN II and SN I, respectively.
No spatial clustering of the SNe is included since the size of superbubbles produced by SNe clustering are comparable to the horizontal size of the computational domain.Simulations in a domain of significantly larger size are required to capture the effects of the SN clustering.de Avillez and Breitschwerdt (2007) include SN clustering in their simulations and obtain the correlation scale of the random flows of 75 pc, comparable to those obtained from the correlation analysis of this model (see Hollins et al. 2017).Each SN is initialised as an injection of 0.5 × 10 51 erg of thermal energy and a variable amount of kinetic energy that depends on the local gas density fluctuations and the turbulent ambient velocities, and has the mean net value 0.4 × 10 51 erg.
We include optically thin radiative cooling with a parameterised cooling function.For T < 10 5 K, we adopt a power-law fit to the 'standard equilibrium' pressure-density curve of Wolfire et al. (1995), as given in Sánchez-Salcedo et al. (2002).For T > 10 5 K, we use the cooling function of Sarazin and White (1987).Photoelectric heating is also included as in Wolfire et al. (1995); it decreases with |z| on a length scale comparable to the gas scale height near the Sun.The system exhibits distinct hot, warm and cold gas phases identifiable as peaks in the joint probability distribution of the gas density and temperature.
Shock-capturing kinetic, thermal and magnetic diffusivities (in addition to background diffusivities) are included to resolve shock discontinuities and maintain numerical stability in the Navier-Stokes, heat and induction equations.Periodic boundary conditions are used in y, and sheared-periodic boundary conditions in x.Open boundary conditions, permitting outflow and inflow, are used at the vertical boundaries at z = ±L z .Gent et al. (2013a,b) provide further details on the boundary conditions used and on the other implementations briefly described above.
Starting with a weak azimuthal magnetic field at the midplane, the system is susceptible to the dynamo instability.Dynamo action can be identified with exponential growth of magnetic field saturating after about 1.4 Gyr at a level of 2.5 µG, comparable to observational estimates for the solar neighbourhood (Gent et al. 2013a).The magnetic field has energy at a scale comparable to the size of the computational domain, suggesting a mean-field dynamo action (Gent et al. 2013b).
We analyse snapshots in the range 0.8 ≤ t ≤ 1.725 Gyr.Three distinct temporal stages can be identified in the dynamo action and the magnetic field.With magnetic energy low, compared to the thermal and kinetic energies, 0.8 ≤ t < 1.1 Gyr hosts the kinematic phase of the mean-field dynamo.The dynamo adjusts itself to a non-linear stage at 1.1 ≤ t < 1.45 Gyr as the magnetic energy reaches approximate equipartition with kinetic energy of the random flow.Finally, at 1.45 ≤ t ≤ 1.725 Gyr, the mean-field dynamo saturates and the magnetic energy slightly exceeds the kinetic energy (see Gent et al. 2013b).Since the evolution of the magnetic field is expected to significantly affect the structure of the gas density and velocity, each stage is considered separately.The results are illustrated in the figures shown below using the snapshot at t = 1.6 Gyr.

Mean fields and fluctuations in a compressible random flow
Averaging procedures can be used to represent a physical variable f as a superposition of its mean f and fluctuations f : f = f + f .Ensemble averaging is used in most theoretical contexts.Ensembleaveraged quantities do not need to be independent of any spatial or temporal variable.However, volume and time averaging are often the only options available in simulations and observations, and those averages clearly are independent of spatial and time variables, respectively.Consider for example, the average over a volume V, It satisfies the Reynolds rules of averaging, including leading to for the random variables f and g.This allows evolutionary equations for the central moments f g V , f g h V (with another random variable h), etc., to be derived by averaging the governing equations using relations, such as for the velocity field u (e.g., Monin and Yaglom 2007a), In numerical simulations, V is often the whole computational domain, or some significant part of it, or a (thin) slice parallel to one of the coordinate planes, as in averages over a horizontal plane (x, y), or azimuthal averaging.Such averages are constrained to be partially or fully independent of position, in all three directions in the case of volume averages, in two dimensions for horizontal averages and in the azimuth for axial averages.As we discuss in Section 1, these constraints may be, and often are, unreasonably restrictive.Moreover, any observational data obtained with a finite resolution represent a convolution of the quantity observed with the telescope beam, and are free to vary with position.It is therefore desirable to apply to numerical results an averaging procedure, compatible with the observational procedures, in a manner that does not impose unjustifiable restrictions on the averaged quantities.This is the goal of this paper.
A local mean part of a random field f (x), denoted f , is obtained by spatial smoothing (filtering) of its fluctuations at scales l < , with a certain smoothing length , using a smoothing kernel G : where integration extends over the whole volume where f (x) is defined.The filtering kernel is normalized, and assumed to be symmetric, To ensure that fluctuations in kinetic energy density are positive definite, the kernel must be positive for all x (Aluie 2017, and references therein).The fluctuation field is obtained as (with the link between the prime and the scale being understood).This procedure retains the spatial structure of both the mean field and the fluctuations.We discuss below physically motivated choices for the smoothing length .Thus defined, the averaging procedure does not satisfy the Reynolds rules outlined in equations ( 2) and (3).In particular, the mean of the fluctuations does not vanish, repeated averaging affects the mean field f (x , and the mean and fluctuating fields are not uncorrelated: As a consequence, the standard relations between statistical moments of total fields and their fluctuations, shown in equation ( 4), are no longer valid.
To address these complications, Germano (1992) introduced generalised statistical moments µ( f, g) µ( f, g, h) . . ., of random fields f (x), g(x) and h(x) to ensure that the mathematical soundness and simplicity of the averaged governing equations is regained for both the mean fields and their statistical moments.In fact, relations between the statistical moments are quite similar to the standard ones of equation (4).For example, the generalised statistical moments of the velocity field u(x) are defined as Statistical moments of the fluctuations are obtained from the moments of the total fields and their averages as, for example, As in equation ( 8), we have u i u j 0 and u i u j 0. In addition, u i u j u i u j since u i u i .As a consequence, u i u j u i u j − u i u j = µ(u i , u j ).Replacing statistical moments of the fluctuations such as u i u j wherever they appear with generalised central moments such as µ(u i , u j ), leads to governing equations for the fluctuations in a mathematically simple form practically identical to that obtained under ensemble averaging (see Aluie 2017, for the case of MHD equations).The algebraic structure of the closure is the same, regardless of the choice of the filter G.Such a property is called the averaging invariance of the turbulent equations (see Germano 1992).
In application to the ISM simulations, we consider the decomposition of the physical fields into mean and fluctuating components with the mean fields obtained via filtering with a Gaussian kernel, where is the smoothing length.We perform this analysis for magnetic field B, gas density ρ and velocity u.All averages are denoted with the subscript and fluctuations with the prime, with the exception of magnetic field fluctuations denoted b: 4. The smoothing scale and Fourier spectra The challenge in applying the filtering approach in our context is to determine an appropriate smoothing length or its admissible range.We note that the mean and fluctuating parts of different variables, Table 1.Notation for the total (T), mean (M) and fluctuating (F) fields and their respective Fourier spectra, arising from the filtering decomposition.See section 4 for definitions.Scale separation between the mean and fluctuation fields is required neither by theory based on ensemble averages nor by the filtering technique.Nevertheless, it is natural to expect some difference in scales between the two.For example, the scale of the mean field in a turbulent dynamo is controlled by deviations of the random flow from mirror symmetry and mean velocity shear, whereas turbulent scales depend on the nature of the driving forces.Given the fundamental difference between the two groups of physical effects, it is unlikely that the two parts of magnetic field have similar scales.Since deviations from mirror symmetry are usually weak, the scale of the mean field is expected to be correspondingly large and to exceed the turbulent scale.Arguments of this kind are used to justify the two-scale approach in mean-field magnetohydrodynamics (Moffatt 1978, Krause and Rädler 1980, Zeldovich et al. 1983).However, numerical simulations of dynamo systems (including those discussed here) are performed in domains that are only moderately larger than the integral scale of the simulated random flow (Brandenburg and Subramanian 2005, and references therein) which precludes any strong scale separation between the simulated mean and fluctuating fields.Nevertheless, evidence for such separation is usually sought, in the form of a pronounced minimum in the Fourier spectra at a scale exceeding the presumed integral scale of the fluctuations (often, the scale at which the random flow is driven by an explicit force) and the domain size.In application to the magnetic field, Gent et al. (2013b) demonstrate that the situation can be more subtle and, despite a pronounced difference of the two scales (by a factor of two), the Fourier spectrum of the total magnetic field may not have a noticeable minimum between them.

Spectrum Energy density Energy
The Fourier spectrum of the total magnetic field B is given by where We also consider the integral scale of each field (Sect.12.1 in Monin and Yaglom 2007b), calculated using the appropriate power spectrum S (k), where ∆ is the grid spacing and D the size of the computational domain.Since both the mean field and the fluctuations are inhomogeneous, equation ( 15) can be used to derive the characteristic scales of both the mean and fluctuating fields: e.g., L B for the mean magnetic field, such that The spectra and lengths scales for ρ, u and their respective mean and fluctuations are defined in a similar manner and denoted S ρ (k), S ρ (k), S ρ (k), S u (k), S u (k) and S u (k), with the corresponding length scales L ρ , L ρ , L ρ , L u , L u and L u .The notation is summarized in Table 1.
As discussed in Sections 4.1, 4.2 and 4.3, none of the Fourier spectra of B, ρ and u, have a local minimum.Nonetheless, each variable has distinct, well separated length scales for the mean and fluctuating fields.The optimal smoothing scale for each variable is obtained under the requirements that: (i) the major maxima in the Fourier spectra of the mean fields and fluctuations in each variable occur on different sides, along the wavenumber axis, of the wavenumber where they intersect; and (ii) that the ratio of the integral scales of the mean fields and the fluctuations is (approximately) maximized.
The power spectrum S B (k) is equivalent, up to a constant factor of 1/(8π), to the magnetic energy spectrum M(k) = S B (k)/(8π), and the total magnetic energy can be obtained as an integral over the relevant wavenumber range, E B = k M(k) dk.However, unlike the case of incompressible flows, the power spectrum of the velocity field cannot be directly equated to the kinetic energy density because of the contribution from the gas density fluctuations.
Calculation of energy densities due to the mean fields and fluctuations should be done with care.To illustrate the general approach, consider magnetic energy.In the filtering approach, the energy densities sum as required from equation ( 9), with the following definitions: (with summation over repeated indices understood in the final definition).Note that e b b 2 /(8π); this is discussed further below.We introduce a distinct notation for the volume integrals of these energy densities, to allow the meaningful summation of the energies: where dV = d 3 x.These energy densities and their volume integrals are summarised in Table 1, and discussed further in section 5.
It is important to appreciate the distinction, within the filtering approach, between the quantities introduced above (which allow a meaningful decomposition of energy), and the 'naive' energies obtained directly from the decomposed parts of the field (which differ in some cases, and which do not meaningfully sum).Although energies can be obtained from the latter -via volume integrals of the squared quantities, or via wavenumber integrals of the power spectra -these quantities do not sum to a valid decomposition, and so are not here identified with the energies of the mean and fluctuating parts, or of their sum.
For clarity, we do briefly discuss these alternative definitions of energies here, but we do not include them in Table 1, or use these definitions in the following.The total magnetic energy E B satisfies As noted above, and discussed in more detail in section 5, e b must be defined in terms of the generalised second moment with i = j (from equation ( 9)), e b = µ(B i , B i )/(8π), so that e b b 2 /(8π).
To summarise the comparison of energies between the two approaches: and E b E b ; and as noted before, while As a result, we focus in the following on the filtered energies (E B , etc.), and do not further refer to the naive energies (E B , etc.).We do analyse the Fourier spectra of the basic physical variables in Sections 4.1-4.3 to identify appropriate smoothing lengths , which can be different for different variables; but then use the filtering approach to derive and discuss the corresponding energy densities in Section 5. in the three stages of magnetic field evolution, kinematic 0.8 ≤ t < 1.1 Gyr (solid, blue), transitional 1.1 ≤ t < 1.45 Gyr (dash-dotted, green) and non-linear 1.45 ≤ t ≤ 1.725 Gyr (dashed, red).

Magnetic field
In figure 1 we show the effect of varying the smoothing length on the power spectra of the mean magnetic field and its fluctuations.A short = 20 pc (figure 1b) assigns the majority of the energy to the mean field including a large proportion at small scales, while too long = 140 pc (figure 1c) assigns most of the energy to the fluctuations even at very large scales.At wavelength λ the mean and fluctuation power spectra intersect, S B (λ) = S b (λ), with k > λ mainly characteristic of fluctuations and k < λ of mean field.We also compute the integral scale L for the decomposed fields, and would expect L B > λ > L b .Instead for = 20 pc, λ = 0.09 kpc < L b = 0.17 kpc, and for = 140 pc L B = 0.92 kpc < λ = 1.09 kpc, which are both physically inconsistent.
A more satisfactory picture emerges when = 50 pc, shown in figure 1a.L B = 0.65 kpc, λ = 0.3 kpc and L b = 0.27 kpc, such that L b < λ < L B .Thus, = 50 pc could be adopted as an appropriate smoothing length for the magnetic field: then the mean field dominates at scales around L B whereas the fluctuations contribute most of the power at scales around L b .
The ratio of L B and l b as a function of is shown in figure 1d for the three stages of the magnetic field evolution.When magnetic field is still weak, there is a pronounced maximum at = 65 pc which becomes less prominent as the magnetic field growth saturates.Thus, the requirement that L b < λ < L B is compatible with the maximum scale separation between the mean field and the fluctuations.The ratio reaches an asymptotic value in the range 0.3-0.4 at ≈ 90 pc.

Gas density
Using the same arguments as for magnetic field, we conclude that = 50 pc is a suitable smoothing length for the density distribution, as also shown in figure 2. Indeed, when = 50 pc (figure 2a), we obtain L ρ = 0.62 kpc > λ = 0.31 kpc > L ρ = 0.27 kpc.In contrast for = 20 pc (figure 2b), λ = 0.11 kpc < L ρ = 0.17 kpc and, for = 140 pc (figure 2c), L ρ = 0.91 kpc < λ = 0.95 kpc.The ratio of L ρ and L ρ as a function of is shown in figure 2d.Its maximum is reached at values of increasing from 65 pc to 75 pc as the magnetic field saturates, suggesting a suitable smoothing length of approximately 70 pc.

Gas velocity
Figures 3a-d illustrate similar arguments for the velocity field u (we recall that u represents deviations from the overall shearing flow and contains a systematic vertical outflow velocity).When = 50 pc, L u = 0.66 kpc > λ = 0.3 kpc > L u = 0.27 kpc.However, for = 20 pc we have λ = 0.12 kpc < L u = 0.16 kpc, whilst for = 140 pc we have L u = 0.92 kpc < λ = 1.12 kpc, each in conflict with the demarcations of mean and fluctuation field.Unlike for the magnetic field and gas density spectra, the ratio of length scales in figure 3d does not have any pronounced maxima, as it increases monotonically with for t < 1.45 Gyr, and has a very broad maximum at = 90-100 pc for t > 1.45 Gyr.
It is clear from each of figures 1, 2 and 3, that the spectral properties of each of these fields are distinct.In addition, the properties of each field vary in time.The simulation times considered here, 0.8 ≤ t ≤ 1.725 Gyr, are all well after the SN-driven hydrodynamics has reached a statistically-steady state, which occurs at about 400 Myr.Thus, we are confident that any changes in time result from the evolution of the mean-field dynamo, which evolves over a time-scale of order Gyr.
It would therefore seem most appropriate to select different smoothing lengths to obtain the fluctuations, depending on both the variable considered and the simulation time.However, complications would then arise with the interpretation of results obtained from such choices.The sensitivity of the results to any change in smoothing length would have to be considered.Theories based on a filtering approach to the MHD equations requires a consistent filter as the averaging operator.Hence, applying different smoothing lengths for each variable would introduce new difficulties when trying to interpret the mean fields and moments of the fluctuating fields as solutions of the filtered equations.In addition, complications could arise when selecting a smoothing scale for moments computed from multiple basic variables, such as the kinetic energy ρu 2 .A time dependent smoothing length could be used, interpreted as a change in the grid scale of such a simulation.
We suggest it most useful to identify an appropriate value of that can be used as a smoothing length for all three variables throughout the times considered.We adopt = 75 pc as the smoothing length for magnetic field, gas density and gas velocity, since for magnetic field and gas density the local maxima in the ratios of the mean and fluctuating length scales occur close to 75 pc.For the gas velocity, the value of this ratio at 75 pc is above 90% of the asymptotic value in each stage, whilst the value at 75 pc in the saturated stage is very similar to the value at the broad local maximum.
In the context of the numerical model described in Section 2, this smoothing length is comparable to the correlation scales for the density fluctuations, random velocity and magnetic fields as obtained by Hollins et al. (2017).Such a result is sensitive to the choice of simulation parameters.For example, an increased supernova rate results in increased gas compressions, leading to stronger local density gradients and the formation of more filamentary structures.Additionally, stronger local changes in velocity would be observed.Thus, the length scales of the small-scale density and velocity would likely decrease, and so the optimal choice of would be reduced.Similar effects are likely to be observed with the reduction of the kinematic viscosity and the thermal and magnetic diffusivities.
Such a detailed exploration of the parameter space, involving multiple simulations with different parameter choices, is beyond the scope of this paper.However, in this section, we have successfully demonstrated that an optimal choice of the smoothing length for each of the fundamental physical fields can be reasonably obtained within a detailed simulation of the ISM.

Energy densities
Magnetic and kinetic energy densities have to be derived using the generalized central moments, as discussed in Section 3. The required moments are derived in Appendix B. Since the mean and fluctuating fields are sensitive to the choice of smoothing length, the resultant energies will also depend on .The maximum admissible value of is half the horizontal extent of the simulation domain.We derive the energy densities obtained with various smoothing lengths in the range 0 < < 0.5 kpc Figure 4. a) Volume averages of the mean magnetic energy density e B V at times 0.8 ≤ t < 1.1 Gyr (green, dash-dotted), 1.1 ≤ t < 1.45 Gyr (black, solid) and t ≥ 1.45 Gyr (red, dotted); also the fluctuating magnetic energy density e b V at 1.1 ≤ t < 1.45 Gyr (blue, dashed), as functions of the smoothing length .These are normalised by the volume average of the smoothed magnetic energy density, e B V , with the volume averaging over the region |z| < 0.5.b) Derivatives of e b V , normalised by e B V , with respect to at 0.8 ≤ t < 1.1 Gyr (green,dash-dotted), 1.1 ≤ t < 1.45 Gyr (blue, dashed) and t ≥ 1.45 Gyr (red, dotted).4a but for the volume average of the mean kinetic energy density e s V at 0.8 ≤ t < 1.1 Gyr (green, dash-dotted), 1.1 ≤ t < 1.45 Gyr (black, solid) and t ≥ 1.45 Gyr (red, dotted); with the volume average of the fluctuating kinetic energy density e st + e t V at 1.1 ≤ t < 1.45 Gyr (blue, dashed).These are normalised by the volume average of the smoothed kinetic energy e k V .(b) As figure 4b but for the derivative of e st + e t V , with respect to (normalised by e k V ); at 0.8 ≤ t < 1.1 Gyr (green, dash-dotted), 1.1 ≤ t < 1.45 Gyr (blue, dashed) and t ≥ 1.45 Gyr (red, dotted).and discuss the results in this section.As previously, we consider the three stages of the mean-field dynamo separately, and present results averaged over the snapshots within each stage.

Magnetic energy
The total magnetic energy density is given by with the energy density of the fluctuating magnetic field obtained as This ensures the energies of the mean and fluctuating magnetic fields sum to the energy of the (filtered) total magnetic energy, i.e.where e B = |B | 2 /(8π) is the energy density of the mean magnetic field.We note that e b |b| 2 /(8π), but it can be shown, by expanding B(x ) in a Taylor series around x, that e b = |b| 2 /(8π) + O( 2 /L 2 B ). Thus, the difference between the volume and filtering averages decreases as /L B → 0. This fact, also true for any other variable, suggests one consideration for the choice of might be to maximise the ratio for L B / .In practice, however, this would simply lead to → 0; i.e. all the signal in the mean field, and effectively no decomposition.
The larger is , the smaller part of the total field is deemed to be a mean field, and e B V monotonically decreases with whilst e b monotonically increases, as shown in figure 4a.The rate of variation of e b V / e B V with , shown in figure 4b -and also of e B V / e B V , not shown -becomes relatively small when > 50 pc.This confirms that the appropriate choice for the smoothing length is > 50 pc.(The difference between figure 4a and figure 2a of Gent et al. (2013a) is caused by a downsampling to a grid ∆x = 8 pc used in the Fourier transform for that calculation in Gent et al. (2013a).) The mean magnetic energy grows with time due to dynamo action, and the value of for which the two energies are equal to each other increases.At late times, the mean magnetic field is energetically dominant over the fluctuating magnetic field for all .

Kinetic Energy
In a compressible flow, the mean kinetic energy density is represented by a third-order moment involving the density and velocity fields.Under ensemble (or volume) averaging, the mean kinetic energy density is conveniently -and physically meaningfully -represented (see Section 6.4 in Monin and Yaglom 2007a) as where e s is the energy density of the mean flow, e t is the energy density of the fluctuations and e st represents the transport of momentum ρ u i by the mean flow (summation over repeated indices is understood here and below).An equivalent decomposition is appropriate under the filtering approach as well: where the moments involved are derived in Appendix B in explicit integral forms: where ∆ρ (x, x ) = ρ(x ) − ρ (x) and ∆u (x, x ) = u(x ) − u (x).
Figure 5 shows how various parts of the kinetic energy density depend on the smoothing length .The behaviour of the volume averages of these contributions to the kinetic energies is much less straightforward than for magnetic energy, except for t > 1.45 Gyr where similar monotonic dependence on is observed.Additionally for both 0.8 ≤ t < 1.1 Gyr and 1.1 ≤ t < 1.45 Gyr, we observe that the fluctuating kinetic energy e st + e t V is equal to zero within errors for 50 ≤ ≤ 100 pc.This results from cancellation between e st V and e t V , with e st V significantly negative, as confirmed by figure 7. The quantity e st = u i µ(ρ, u i ) is dominated by the contribution of the z-component of the velocity field (i = 3) since u z is much larger than the x-and ycomponents because of a systematic gas outflow from the midplane.
The gas involved in the outflow is hotter and less dense than on average, leading to large negative values of − ρ u z for z > 0 and, hence, of u z µ(ρ, u z ) = u z ( ρu z − ρ u z ) (the dominant component of e st ).
For z < 0 kpc, the mean vertical velocity driven by the supernovae u z is large and negative, resulting in large, positive values of µ(ρ, u z ) = ρu z − ρ u z .Thus, the opposite signs of u z and µ(ρ, u z ) result in large, negative values of e st for negative z.These large, negative values for e st appear to dominate the kinetic energy statistics at earlier times.This is discussed in more detail below.
The variation with of the fluctuating kinetic energy produces a more complicated pattern than for fluctuating magnetic energy, see figure 5.The values of for which the variation is weak are > 300 pc.Such a smoothing length is much larger than any estimate of the correlation scale of the random motions, and the optimal smoothing lengths of both ρ or u.As a result, the criterion that the variation of the fluctuating kinetic energy must be weak is not an appropriate method for choosing suitable smoothing lengths for either ρ or u.

Influence of the mean-field dynamo
Figures 4 and 5 both suggest that the structures of magnetic and kinetic energy distributions vary with the state of the mean-field dynamo.We first examine the vertical structure of both energies, comparing the three temporal stages discussed previously, to demonstrate the changes in structure caused by the dynamo.
At early times, when the fluctuating magnetic field dominates the mean field, the magnetic field is strongest at z = 0 kpc, but significant local maxima (in each of e B , e B and e b ) develop at |z| = 0.3 kpc.This is in line with the kinetic energy, where e k is maximal at the midplane, but e s and e st have extrema at |z| = 0.3 kpc; see figures 6a and 7a.
As the mean field dynamo saturates, the mean magnetic field dominates compared to the fluctuating field.The vertical profile of the smoothed total magnetic energy is increasingly dominated by the mean magnetic energy, with the peaks at |z| = 0.3 kpc now dominating; see figures 6b,c.The increasing mean magnetic field significantly alters the vertical profile of the kinetic energy, as shown in figures 7b,c.All the components in the division of kinetic energy are ultimately concentrated near the midplane (so that the peaks at z = 0.3 kpc have been suppressed), and the maximum value of e k decreases.
Strong mean magnetic fields generated via dynamo action in the same ISM simulations have been shown to suppress outflows of the hot gas (see Evirgen et al. 2017), which are associated with high values of kinetic energy.This would lead to a vertical profile of kinetic energy with the characteristics present in figure 7c.The most dramatic change is the effect on the 'intermediate scale' component of the kinetic energy, e st .As the magnetic field strength increases, the magnitude of the horizontal average of e st decreases significantly, becoming small except near to the midplane.As a result, the kinetic energy is approximately split between the mean and small-scale energies e s and e t .
As this change appears to be the most significant, we focus on horizontal planes from the snapshots t = 0.8 Gyr and t = 1.6 Gyr at which the vertical profiles of e st shown in figure 7 show the most profound differences.
In the kinematic stage of the mean-field dynamo, there are regions in which e st differs from zero significantly, whilst e k is uniform by comparison (see figure 8).The mean and turbulent kinetic energies, e s and e t respectively (not shown here), are also significant in the same regions as e st .The contribution to e st from the z-component of u, u z µ(ρ, u z ), comprises a large fraction of the total kinetic energy (about 80%) and so the vertical flow is dominant in e st at this stage.The values for which e st is highest strongly coincide with regions of large positive u z , which are the regions of hot gas outflows.Thus, at the kinematic stage of the mean-field dynamo, e st is strongly correlated with the outflows of hot gas.In this model, the mean magnetic field is absent from the regions of hot gas, as demonstrated by Evirgen et al. (2017).Thus, the mean magnetic field also avoids regions in which e st is significant in magnitude.The effect of the amplified mean magnetic field on the kinetic energies is demonstrated in figure 9.The values of e st are reduced significantly and e st appears more uniform.By contrast, e k is now more significant and the non-uniform structure of e k is more pronounced.The contribution to e st from the vertical flow is also dramatically reduced, and no longer dominates.The mean vertical velocity is reduced both in maximal value and in the size of regions in which u z is significant, indicative of the suppression of the hot gas outflow.Thus, the partial suppression of the hot gas outflow by the mean magnetic field has both significantly reduced the value of e st and resulted in the dynamics of the overall kinetic energy becoming independent of the behaviour of e st .

Comparison with horizontal averaging
To emphasise the physical relevance of the filtered quantities as representations of the mean and fluctuating fields, as compared with the corresponding quantities from horizontal averaging, in this section we perform a direct comparison.
Figure 10 shows a representative comparison of the two decompositions.In both panels, the solid blue curve gives the total vertical velocity at fixed y and z; panel (a) also shows the fluctuating parts obtained via both averaging methods, and panel (b) shows the corresponding mean parts.We notice that here the velocity is small in magnitude for x < 0. However, there is a region of higher vertical velocity for x > 0, which corresponds to gas flowing away from the midplane.This region has the typical characteristics of a hot gas structure.
Since horizontal averaging takes an average value for the entire horizontal plane, panel (b) shows the mean vertical gas velocity as constant with respect to x, with U h,z = u z xy ≈ 20 km s −1 .The systematic flow seen in x > 0 kpc has approximate velocity 60-80 km s −1 .Consequently, the fluctuating field, derived using horizontal averaging, retains a large proportion of this systematic flow.This is apparent in the similarity of the curves and magnitudes of the total vertical velocity and fluctuating vertical velocity defined using horizontal averaging.In general, where there are significant local deviations within a horizontal layer (as with the localised hot gas structures with large vertical velocity here), horizontal averaging will not capture the systematic features well.The horizontal mean calculated may be unrepresentative of the flow in most of the layer (as in panel (b)), and a significant proportion of the fluctuating part may simply be correcting for this locally inappropriate mean (as in panel (a)), rather than accurately reflecting local fluctuations.Here this leads to downwards fluctuating flows of up to 20 km s −1 ; in addition to inflating the magnitude of the fluctuating velocity, this By contrast, the mean vertical velocity obtained with Gaussian smoothing, shown by the dashed magenta curve of panel (b), extracts the systematic flow well, retaining information on the region of hot gas moving away from the plane.In addition, this mean retains features of the systematic vertical velocity more accurately in regions where that velocity is smaller (close to zero).The fluctuating vertical velocity derived using Gaussian smoothing retains a smaller proportion of the systematic background trend, so its magnitude remains commensurately small.Figure 10.Vertical component of gas velocity as a function of x (galactocentric radius) at fixed y and z, to compare horizontal averaging and Gaussian smoothing.In Panel (a), we present the total vertical gas velocity U z = u z (blue solid); fluctuating vertical velocity defined using horizontal averaging (green dash dotted), u h,z = u z − u z xy ; fluctuating vertical velocity defined using Gaussian smoothing (magenta dashed), u l,z = u z = u z − u z .In Panel (b), we present the total vertical gas velocity (blue solid); mean vertical velocity defined using horizontal averaging (green dash dotted), U h,z = u z xy ; mean vertical velocity defined using Gaussian smoothing (magenta dashed), U l,z = u z .
Given that the system discussed here, like many other astrophysical systems, contains systematic features such as those discussed above -driven by the bubbles of hot gas associated with the SN activity in the case of the ISM, which are well-defined local 3D features that cannot be well represented by a horizontal average across the whole layer -it is perhaps not surprising that Gaussian smoothing extracts such systematic features more successfully than horizontal averaging.This will also be true in many other instances with pronounced local 3D structures.(It should be noted that the effect will not be so pronounced where the horizontal layer is large enough to contain multiple such structures -involving both outflows and inflows -which might then average out over the full layer.But even in such cases, the horizontally-averaged mean cannot capture these systematic features, which should sensibly be characterised as part of the large-scale field, rather than as fluctuations.) Gaussian smoothing (or filtering with any other kernel) is not technically or conceptually more difficult to implement and interpret than any other averaging method, but it requires an appropriate filtering scale to be identified.The Gaussian kernel used here is isotropic but the filtering approach also offers the opportunity to use anisotropic kernels wherever appropriate.This can be specially important in the case of the velocity field; for example, SN remnants expand more strongly along the density gradient in galactic discs, and the hot gas is buoyant.It is reasonable to expect that the magnetic field reflects such features of the gas flow and application of an anisotropic filtering kernel might lead to a clearer physical picture, although this would introduce an additional parameter, the degree of the kernel anisotropy.This situation is quite similar to that with the anisotropic wavelet decomposition (e.g., Patrikeev et al. 2006).

Discussion
We have applied Gaussian smoothing to obtain mean fields for magnetic field, density and velocity in a simulation of the multiphase interstellar gas: a complex, partially ordered magnetohydrodynamic system that supports the mean-field dynamo action.The optimal smoothing lengths were obtained by spectral analysis of each field independently.We find = 75 pc, approximately 19 grid cells, as an appropriate smoothing length to use for each of these fields.Such a result is likely to be sensitive to the choice of simulation parameters, and should be investigated with a subsequent detailed exploration of the parameter space.However, we have successfully demonstrated that an optimal value of for each of the three fundamental physical fields is indeed possible in our simulation of the ISM.We have also shown that Gaussian smoothing, unlike horizontal averaging, retains large-scale 3D features of the mean fields.The filtering approach allows for a more physically-meaningful decomposition into mean and fluctuating parts for each variable and for their higher statistical moments, such as the magnetic and kinetic energy densities.
It is natural to expect that the magnetic, density and velocity fields can have distinct spatial structures since they are controlled by different physical processes, even though they do not evolve independently.The Gaussian smoothing successfully reveals such differences (e.g., the different quantities have different correlations lengths), as shown in Hollins et al. (2017) for a filtering length of 50 pc.(The change to 75 pc suggested here does not affect this significantly.)Gressel and Elstner (2020) apply kernel filtering to acquire the mean magnetic field from simulations of SN driven turbulence and explore the scale separation between mean and random fields in their simulations.They show that the diamagnetic transport of the mean field is stronger when the filtering scale is smaller.
We examine the mean and fluctuating magnetic and kinetic energies, using the generalised central moments to define the energy density of the fluctuations.We examine the dependencies of the energies on , and on the magnetic dynamo saturation.This allows us to identify the key physical processes affecting the mean and fluctuating fields.Amplification of the mean magnetic field by dynamo action has a significant impact on the contributions to the magnetic and kinetic energies due to the mean fields, fluctuations and, in the case of the kinetic energy, advection of the fluctuations by the mean flow.The growing mean magnetic field shifts the maximum of the vertical profile of kinetic energy towards the disc midplane.part of the kinetic energy density associated with the advection of the velocity fluctuations by the mean flow, e st , is closely correlated with a systematic gas outflow and is partly suppressed by the growing mean magnetic field.This results in a dramatic reduction in e st at late times in the simulation, when the kinetic energy is mostly associated with the large-scale flow and the velocity fluctuations.
is the Fourier transform of B and • k denotes the average value within a spherical shell of thickness δk with radius k = |k|.The power spectra for the mean and random fields, S B (k) and S b (k), are similarly defined in terms of B (k) and b(k), the Fourier transforms of B and b: B dV (where e B = B 2 /(8π), as above); but E B E B .The integral of the mean field satisfies E B = 1/(8π) S B (k) dk = V e B dV (where e B = B 2 /(8π), as above), and does equate to the mean energy introduced above, E B = E B .But the corresponding quantities for the fluctuating field, E b = 1/(8π) S b (k) dk = 1/(8π) V b 2 dV, and E b = V e b dV, are not equal: E b E b .

Figure 1 .
Figure 1.Fourier spectra of the total magnetic field S B (k) (solid, blue), its mean part S B (k) (dash-dotted, green) and the fluctuations S b (k) (dashed, red) for = 50 pc at t = 1.6 Gyr for various values of the smoothing length : (a) = 50 pc, (b) = 20 pc and (c) = 140 pc.The vertical dotted lines indicate (from left to right) the wavenumbers corresponding to the scale of the mean field L B , its fluctuations L b , the smoothing length and the resolution of the simulations ∆. (d): Ratio of the integral scales L b and L B as a function of the smoothing lengthin the three stages of magnetic field evolution, kinematic 0.8 ≤ t < 1.1 Gyr (solid, blue), transitional 1.1 ≤ t < 1.45 Gyr (dash-dotted, green) and non-linear 1.45 ≤ t ≤ 1.725 Gyr (dashed, red).

Figure 5 .
Figure 5. (a) As figure4abut for the volume average of the mean kinetic energy density e s V at 0.8 ≤ t < 1.1 Gyr (green, dash-dotted), 1.1 ≤ t < 1.45 Gyr (black, solid) and t ≥ 1.45 Gyr (red, dotted); with the volume average of the fluctuating kinetic energy density e st + e t V at 1.1 ≤ t < 1.45 Gyr (blue, dashed).These are normalised by the volume average of the smoothed kinetic energy e k V .(b) As figure4bbut for the derivative of e st + e t V , with respect to (normalised by e k V ); at 0.8 ≤ t < 1.1 Gyr (green, dash-dotted), 1.1 ≤ t < 1.45 Gyr (blue, dashed) and t ≥ 1.45 Gyr (red, dotted).
e B = e B + e b ,

Figure 6 .
Figure 6.Vertical profiles of the horizontal averages of the smoothed total magnetic energy, e B xy (blue, solid), mean magnetic energy, e B xy (green, dashed), and fluctuating magnetic energy, e b xy (red, dash-dotted); at times (a) 0.8 ≤ t < 1.1 Gyr, (b) 1.1 ≤ t < 1.45 Gyr and (c) t ≥ 1.45 Gyr.The smoothing length applied for each snapshot is = 75 pc.

Figure 7 .
Figure 7.As in figure 6 but for the smoothed total kinetic energy density e k xy (blue, solid), mean kinetic energy density e s xy (green, dashed), 'intermediate scale' kinetic energy density e st xy (cyan, dash-dot-dotted), fluctuating kinetic energy density e t xy (red, dot-dashed), and the sum, e st + e t xy (purple, dotted); at times (a) 0.8 ≤ t < 1.1 Gyr, (b) 1.1 ≤ t < 1.45 Gyr and (c) t ≥ 1.45 Gyr.As in figure 6, the smoothing length applied is = 75 pc.

Figure 8 .Figure 9 .
Figure 8. Horizontal slices of of the smoothed total kinetic energy density e k (top-left panel), the 'intermediate-scale' kinetic energy density e st (top-right panel), u z µ(ρ, u z ) the vertical contribution to e st (bottom-left panel), all in units of 10 −13 erg cm −3 , and the mean vertical velocity u z (bottom-right panel) in km s −1 ; at z = 290 pc from the snapshot t = 0.8 Gyr.The smoothing length used is = 75 pc.
e st , e t E k E s E st , E t e.g., B, ρ and u, can have different spatial properties and, hence, different smoothing lengths may be required to separate the fluctuations in different variables.For example, Hollins et al. (2017) find that the correlation lengths of the three variables are different in the simulations discussed here.Unlike applications to subgrid turbulence models, where is identified with the spatial resolution of a simulation, the choice of in the present context is motivated by physical considerations.Following Gent et al. (2013b), we select using the spectral structure of each variable as discussed below.