Parameters of the Supernova-Driven Interstellar Turbulence

Galactic dynamo models take as input certain parameters of the interstellar turbulence, most essentially the correlation time $\tau$, root-mean-square turbulent speed $u$, and correlation scale $l$. However, these quantities are difficult, or, in the case of $\tau$, impossible, to directly observe, and theorists have mostly relied on order of magnitude estimates. Here we present an analytic model to derive these quantities in terms of a small set of more accessible parameters. In our model, turbulence is assumed to be driven concurrently by isolated supernovae (SNe) and superbubbles (SBs), but clustering of SNe to form SBs can be turned off if desired, which reduces the number of model parameters by about half. In general, we find that isolated SNe and SBs can inject comparable amounts of turbulent energy into the interstellar medium, but SBs do so less efficiently. This results in rather low overall conversion rates of SN energy into turbulent energy of $\sim1$-$3\%$. The results obtained for $l$, $u$ and $\tau$ for model parameter values representative of the Solar neighbourhood are consistent with those determined from direct numerical simulations. Our analytic model can be combined with existing dynamo models to predict more directly the magnetic field properties for nearby galaxies or for statistical populations of galaxies in cosmological models.


Introduction
Turbulence affects a wide range of physical processes in the interstellar medium (ISM) of spiral galaxies, including the turbulent dynamo. Its spectrum extends over a wide range of scales [1,2], apparently maintained by diverse physical effects. Turbulent flows are characterized by certain physical parameters, which are defined through some sort of averaging. These include the velocity correlation scale (or integral scale) l, which is similar to, but smaller than, the scale of the force driving the turbulence, the root-mean-square (RMS) turbulent speed u, and the turbulent correlation time τ.
In this work, our aim is to obtain scaling relations for the turbulence parameters l, u and τ in terms of other quantities like gas density, sound speed, and supernova rate density. The latter parameters are often readily computed from observations or models, so such relations can provide a missing link. We are motivated by one application in particular: the galactic dynamo, which is responsible for amplifying a galaxy's magnetic field up to an energy comparable with that of the turbulence. The properties of the magnetic field as it evolves and saturates are predicted to depend on the values of the parameters u, τ and l [3][4][5][6][7]. Table 1. The underlying ISM parameters with adopted values. 'Usage' designates whether the parameter is needed to model isolated SN driving or SB driving. (If either SBs or isolated SNe are neglected, then f SB = 0 or f SB = 1, respectively.) The range suggests plausible values for Milky Way-like galaxies, but the model is not restricted to this range. The fiducial values refer to typical estimates for the Solar neighbourhood of the Milky Way, averaged across the Galactic disc [40]. thin cloud layer in Section 5. In Section 6, we discuss the limitations of our model and opportunities for extending it. Finally, we summarize and conclude in Section 7.

Estimation of Interstellar Turbulence Parameters
Our approach for calculating the turbulence parameters l, u, and τ is summarized as follows. Using a standard scaling relation for the outer radius and speed of an SNR, and assuming that the SNR injects its energy into the ISM once it slows to the ambient sound speed, we derive a scaling relation for the injection scale l SN of turbulence driven by isolated SNe. We then repeat this procedure for SBs, which have their own standard scaling relation, to obtain an injection scale l SB ; however, we in addition include the possibility that the SB blows out of the disc before it can slow to the ambient sound speed. Combining driving by SNe and SBs additively allows us to write down the energy injection rate density. We then assume a simple form for the spectrum which allows us to compute analytically the integral scale l. Given l, we can write the spectral energy transfer rate density, equal to the turbulent energy dissipation rate density, in terms of the turbulent velocity u. We next balance energy injection and dissipation rates, and solve for u to obtain a scaling relation for the latter. The next step is to estimate the correlation time τ, which we take to be equal to the smaller of the eddy turnover time τ e and the average time for the flow at a given position to renovate due to the passage of an SN or SB blast wave, τ r . Table 1 summarizes the independent parameters used in our scaling relations. These are free to take on any values, but we have included a rough range. If desired, our model can be further simplified by including only one of the driving channels, either isolated SNe or SBs.

SNRs
Toward the end of its life, an SNR experiences a momentum-conserving snowplough (MCS) phase, such that its shell expands according to [41,42] R SN At a andṘ SN aAt a−1 , where A = R 1−a 0 (Ṙ 0 /a) a , a = 1/4, R 0 ≡ R(t 0 ),Ṙ 0 ≡Ṙ(t 0 ), the subscript '0' denotes the onset of the MCS, when the interior has cooled, and t t 0 has been assumed. Using the Sedov-Taylor similarity solution (R SN ∝ t 2/5 ) and a prescription for radiative cooling, Woltjer [41] (see also Reference [42] Cioffi et al. [44] argue that while the MCS phase is approached asymptotically at late times, the SNR will typically merge with the ISM before entering a full-fledged MCS phase. They derive a slightly different solution, valid for t > t 0 = 0. which are very similar to equations (2). Below we choose to adopt equations (2), but it is trivial to replace these with equations (3), with a = 3/10. This leads to only minor differences in the results. The ISM contains magnetic fields and is multi-phase, hence inhomogeneous; these details are neglected in our model. Kim and Ostriker [45] simulated the expansion of a radiative SNR into a uniform medium or a two-phase ambient medium containing cold clouds embedded in a warm neutral medium. For the former, they also considered the case where the environment was filled with an initially uniform magnetic field. They found that magnetic fields do not affect the evolution of the SNR unless the field is very strong (plasma β ∼ 0.1). Further, they found that the radial momentum injection by the SNR into the environment is only 5% smaller for a two-phase medium than for a uniform medium with the same mean density. Meanwhile, Martizzi et al. [46] found that the asymptotic radial momentum of an SNR is typically smaller by about 30% in an inhomogeneous medium compared to a homogeneous medium of the same mean density. This decrement is caused by extra cooling due to the inter-mixing of cold clouds and hot shocked gas. As these differences are relatively small, we model the ambient medium as being of uniform density. These studies do suggest, however, that it is important to include the different ISM phases in estimates of n, and in our model, n could be thought of as a volume-averaged value in the galactic disc.

SBs
Following Mac Low and McCray [47], we use the similarity solution of Weaver et al. [48] and Weaver et al. [49] to model the evolution of an SB in a homogeneous medium: where A = (125L SB /154πρ) 1/5 , α = 3/5, and ρ = (14/11)m H n has been assumed to allow for helium. Further, t 10 = t/(10 Myr) and L 38 = L SB /(10 38 erg s −1 ), where L SB is the equivalent mechanical luminosity of the SNe in the OB association, and where N SB is the number of SNe contributing to the SB over its lifespan t SB and η is the fraction of the injected energy converted into bulk kinetic energy of the SB. A detailed discussion of the features and applicability of the above similarity solution is given by Mac Low and McCray [47] and Breitschwerdt et al. [39]. In particular, radiative cooling of the SB interior can be neglected for the Milky Way parameter values but may be important in denser media, thicker discs, or smaller OB associations. We have included an efficiency factor η ≤ 1 to account for the fact that not all of the energy from SNe ends up as mechanical energy of the SB. The value of η is currently not well-constrained. Yadav et al. [50] compare 3D and 1D simulations of varying resolution. They conclude that η (which in their definition is the ratio of the total combined thermal and bulk kinetic energy to energy injected by SNe) increases with resolution and is not converged at the highest resolutions. They attribute this to the exclusion in the simulations of explicit diffusion processes needed to obtain radiative layers thick enough to be resolved. Nevertheless, their work suggests that η is a decreasing function of the ambient density n and time since the onset of SB expansion. These results are generally consistent with those of El-Badry et al. [12], who use 1D simulations and analytical solutions to model the evolution of an SB. To model mixing of hot and cold gas due to 3D hydrodynamic instabilities, which would affect cooling in the shell-bubble interface, they include an explicit diffusion term with adjustable diffusivity. 2 We adopt η = 0.1 as a fiducial value for an ambient density of n = 0.1 cm −3 , and η = 0.05 for n = 1 cm −3 . These values are generally consistent with simulations, and are also comparable to mechanical efficiencies estimated for isolated SNe (e.g., Ch. 7 of Dyson and Williams 43, and Section 4 below).
The solutions presented above have been obtained for a homogeneous ambient medium, and, strictly speaking, they do not apply when R SB becomes comparable to the pressure scale height. We address this limitation below by considering the case where the SB blows out of the disc. In 3D MHD numerical simulations with a realistic stratification of the ambient medium and an imposed magnetic field, Stil et al. [52] find that the above self-similarity solution approximates the SB evolution rather accurately. On the other hand, Ntormousi et al. [53] simulate two SBs colliding within a turbulent ISM. While the expected R SB ∝ t 3/5 expansion is reproduced without magnetic fields or for a mean magnetic field that is parallel to the collision axis, they find that the presence of a mean magnetic field component perpendicular to the collision axis can result in a flatter SB expansion with time. Furthermore, Breitschwerdt et al. [39] argue, using MHD simulations of the local ISM, that analytical solutions are rather inadequate for modeling SB evolution in a multi-phase ISM. In this work we choose, for simplicity and tractability, to treat the ISM as a uniform medium and to ignore the effects of magnetic fields.

SNRs
An SNR merges with the ISM at t ≈ t SN s , where t SN s is the time at which the expansion velocity of the SNR shell equals the ambient sound speed c s . 3 Then we obtain  2 Fierlinger et al. [51] performed 1D simulations of SN explosions of massive single stars in a dense ambient medium ∼ 100 cm −3 , and found that the preceding wind-blown bubble phase leads to a strong overall reduction in radiative cooling losses. Consistent with the works mentioned above, they find that the importance of cooling depends on resolution, and argue that numerical diffusive mixing in their simulations approximates mixing by turbulent diffusion in nature.
where n 0.1 = n/0.1 cm −3 and c 10 = c s /(10 km s −1 ). The turbulent driving scale can be taken to be of order the SNR radius at that time,

SBs
Similarly, the time t SB s when the expansion speed of an SB reduces to the ambient sound speed, and hence it comes into pressure balance with the surrounding medium, follows from the second of equations (4), assumingṘ SB = c s . Unlike a typical SNR, an SB may be able to blow out of the disc at some time t SB b < t SB s . Then the SB time scale appropriate to the turbulent flow in the disc is given by t SB b rather than t SB s . For the mechanical luminosity L SB , this leads to The self-similar solution presented above yields, fromṘ SB (t SB s ) = c s , where η 0.1 = η/0.1 and N 100 = N SB /100. This corresponds to The evolution of an SB in a stratified ISM [54] is discussed by Mac Low and McCray [47] (see also Refs. [55] and [56]). These authors identify the time when an SB blows out of the disc with the moment when its shell begins to accelerate away from the midplane. They suggest that this occurs when the SB radius at the midplane is about equal to the gas density scale height H, and the SB extends to a height z ≈ 3H above the galactic midplane. Mac Low and McCray [47] adopt H = 500 pc, their estimate of the scale height of the Lockman layer of neutral hydrogen. They argue convincingly that it is this diffuse layer, rather than the thinner cloud layer with scale height ∼ 100 pc, that has the more important influence on the dynamics of the SB. More recent estimates of the local scale height of the Lockman layer suggest 300-400 pc [40], and we adopt H = 400 pc as our fiducial value.
As with turbulence driven by isolated SNe, we identify the injection scale of turbulence with the final radius reached by an SB near the midplane, where the second case corresponds to blowout, and ξ is a parameter of order unity. The scale ξ H is the turbulent driving scale in the case that the SB blows out of the disc. Note that the vertical size of the SB when it blows out is not required, though it is expected to be equal to a few × H [47]. Below we set ξ = 1 for simplicity and because this is consistent with the horizontal extents of SBs at blowout in the simulations of Mac Low and McCray [47], but we return to consider the case ξ = 1/3 in Section 5.
Hence, l SB = R SB (t SB s ) is replaced by l SB = ξ H if the time to attain blowout t SB b is less than the time to slow to the ambient sound speed t SB s . For a flared disc, we might have H ∼ 0.2 kpc near the disc centre, and H ∼ 1 kpc near the outskirts of the star-forming disc. Since R SB (t SB b ) = ξ H, we then have where H 400 = H/(0.4 kpc). For the equivalent mechanical luminosity L SB , we then obtain which is close to the values typically assumed. We assume that most of the energy is injected into the ISM once the SNR or SB reaches its maximum radius, respectively l SN or l SB , and fragments. The mass of the outer shell is equal to that swept up from the ambient medium. A fraction f SB of SNe are assumed to contribute to SBs, so that where ν is the rate per unit volume for all SNe in the galaxy. The rate per unit volume of SBs is then Higdon and Lingenfelter [57] estimate that the fraction of SNe occurring in OB associations is ∼ 3/4 for the Milky Way, and we adopt this as our fiducial value for f SB .
If the SB blows out, we assume that it immediately loses pressure support and slows to the ambient sound speed in the disc. The energy injected is assumed to be equal to the bulk kinetic energy of swept up ambient matter ∼ 2 3 πρξ 3 H 3 c 2 s . Likewise, the SB would have injected ∼ 2 3 πρR 3 SB (t SB s )c 2 s had it been able to decelerate to the sound speed and break up without blowing out. Therefore, the efficiency of conversion of SN energy to turbulent energy in the disc by SBs is reduced by the factor ξ 3 H 3 /R 3 SB (t SB s ). This implies that the value of the scale height can be critical in determining the efficiency of energy conversion.

Energy Conversion Efficiency
The efficiency of turbulent energy injection is given by whereĖ i is the rate of turbulent energy injection per unit volume and νE 0 is the rate of total energy per unit volume released by SNe. The injection rateĖ i is given by the sum of the contributions from isolated SNe and SBs, where the swept up material is assumed to have the same average density ρ for both isolated SNe and SBs. While SBs expand to larger heights, where the ambient density is smaller, it is the lateral expansion within the disc that is most relevant for our calculation, so we regard this assumption as reasonable. For fiducial parameter values this givesĖ i ∼ 2.5 × 10 37 erg s −1 kpc −3 , which leads to an overall efficiency = 0.016, or 3.4 × 10 37 erg s −1 kpc −3 and 0.022 if we instead adopt n 0.1 = 10. This value of is still a few times lower than the value ≈ 0.1 obtained by Thornton et al. [58]. The latter could be an overestimate because their simulations end long before the shell velocity has reduced to the ambient sound speed [51]. However, the values obtained in our model agree closely with those obtained by Bacchini et al. [10], who quote a median SN efficiency, for the galaxies modeled, of ∼ 0.015 in the warm atomic gas, and a range of ∼ 0.01-0.03.

Energy input into an outflow
When an SB blows out, the part of its mechanical energy that is not deposited into the ISM is, in principle, available for driving a galactic outflow. This may take the form of a galactic wind if the gas escapes, or fountain flow if it returns to the disc. Thus, we estimate the power per unit volume of the ISM made available to drive outflows aṡ , right) is the SB expansion speed at the time at which blowout occurs, t SB b (Equation 12), and where we have made use of Equation (11) for l SB . For our fiducial values, we obtaiṅ which is equal to about half the value deposited into the ISM by isolated SNe and SBs, and about 1.25 times the value deposited by SBs alone. The estimate is sensitive to the degree of SN clustering, disc thickness (and, hence, galactocentric distance) among other parameters. We leave further exploration of the outflow properties for future work.

Relative Contribution from Isolated SNe and SBs
The ratio of the rates of energy per unit volume injected by the two mechanisms iṡ Thus, for typical parameter values isolated SNe and SBs inject comparable amounts of energy, but the ratio can vary rather strongly between and within galaxies. Likewise, Norman and Ferrara [21] find that SBs and isolated SNe contribute about equal energies (see their Figure 1 and discussion following their Equation (3.13)). Their model, however, is for highly supersonic turbulence, obtained for the case where heating of the ISM by turbulent dissipation is negligible. Turbulence in the warm phase of the ISM is known to be transonic [e.g. 59-61].

Turbulent Correlation Scale
We are now in a position to estimate the correlation or integral scale l. We adopt a simplified spectral model with constant spectral index −γ. The energy per unit volume injected at scale l 0 with wavenumber k 0 = 2π/l 0 is given by where the spectral energy distribution E(k) = E 0 (k/k 0 ) −γ for k ≥ k 0 and E(k) = 0 for k < k 0 , and E 0 is a constant. We have assumed that E(k) → 0 abruptly for k < k 0 and we also assume that γ > 0, so the integral converges. Defining k ≡ k/k 0 we then have Applying these expressions individually to the spectrum of turbulence driven by SBs or isolated SNe, and assuming each spectrum to have the same spectral index −γ, gives and Here, E SB and E SN are constants, k SB = 2π/l SB , k SN = 2π/l SN , and the injection scales l SB and l SN , with l SN < l SB , can be estimated from equations (11) and (7), respectively. Since the total energy density is simply equal to the sum of the two contributions, it follows from Equation (20) and the first equalities of equations (22) and (23) that The correlation scale for the overall spectrum can be approximated as with C a dimensionless parameter of order unity. For solenoidal (zero divergence) turbulence, Monin and Yaglom [62] (Chapter 12) obtain C = 3/8 and 3/4 for longitudinal and transverse correlations, respectively. For potential (zero curl) turbulence, they obtain C = 0 and 3/8 for longitudinal and transverse correlations, respectively. Interstellar turbulence is in reality expected to contain potential (or compressible) modes as well as solenoidal modes, with the latter modes containing about twice as much energy as the former [63], and they are coupled to one another [see, e.g. Reference 15, for a discussion]. We choose C = 3/4 as our fiducial value but also provide some examples which use C = 3/8. To make progress, we need an expression that relates E SN and E SB . It follows from comparison of equations (22) and (23) that In a steady state, the energy injection and dissipation rates balance, and can be assumed to be equal to the spectral energy transfer rate. If we assume that this holds separately for energy injected by isolated SNe and SBs, we can writeĖ i SN ∼ E SN /τ e andĖ i SB ∼ E SB /τ e . Thus, we assume that the ratio of the turbulent energy densities supplied by isolated SNe and SBs is equal to the ratio of the energy injection rate densities, i.e. E SN /E SB =Ė i SN /Ė i SB . Then from Equation (26) we obtain where the ratioĖ i SN /Ė i SB is obtained from Equation (19). Below we assume γ > 1. Evaluating the integrals in Equation (25) assuming the spectral energy distribution (24) we obtain where we have made use of Equation (27). Substituting the evaluated integrals into Equation (25), we obtain Henceforth we choose γ = 5/3, appropriate for Kolmogorov turbulence. 4 With the choice γ = 5/3, we have (γ − 1)/γ = 2/5, l → 2 5 Cl SB asĖ i SN → 0 or l SN → l SB , and l → 2 5 Cl SN asĖ i SN → ∞. For our fiducial values, Equation (29) with γ = 5/3 and C = 3/4 gives l = 74 pc.
In Figure 1, we plot an example of the energy spectrum using the fiducial parameter values of Table 1 and γ = 5/3. SB and isolated SNe spectra are shown by dashed curves, with the combined spectrum shown by a solid curve. The top panel shows the simple energy spectrum (24). By construction, the ratio of the areas under the isolated SNe and SB spectra is equal toĖ i SN /Ė i SB . On the bottom we illustrate the more realistic modified von Karman spectrum employed by Wilkin et al. [66], where the dissipative wave number k d is chosen, arbitrarily, as 10 4 kpc −1 , and γ = 5/3. The form of the spectrum at k → 0, E ∝ k 4 , is characteristic of both solenoidal and potential velocity fields [p. 52 in Reference 62]; this part of the spectrum has negligible effect on the results. (For the von Karman spectrum, the ratio of the areas differs only slightly fromĖ i SN /Ė i SB in the fiducial case.) In each panel of Figure 1, the integral scale of turbulence l is represented by a vertical solid line. 4 The arguments can be adapted to other spectra. One possibility is to model the spectrum as a broken power law with a different spectral index, say γ = 2, for scales larger than the sonic scale l s , and γ = 5/3 for scales smaller than l s [64,65]. However, the effect on the estimates provided would be rather negligible, particularly considering the myriad of other uncertainties and assumptions on which our model relies.  Table 1 have been used, along with γ = 5/3. Dotted horizontal lines reference the peaks of the isolated SN and SB spectra, and vertical dashed lines their respective injection scales l SN and l SB . The wavelength corresponding to the correlation scale of turbulence, 2π/l, with l given by Equation (29) with C = 3/4, is marked by a vertical solid line.

RMS Turbulent Velocity
We now estimate the RMS turbulent velocity. The energy dissipation rate per unit volume is equal to the spectral energy transfer rate, 5Ė where the eddy turnover time (comparable to the lifetime of the largest eddies) is estimated as the ratio of the integral scale and RMS turbulent speed, Assuming a statistical steady state, equating the injection rate per unit volume from Equation (17) and dissipation rate per unit volume from Equation (31), and solving for u, we obtain where the right-hand side is independent of the density. This result is consistent with that obtained by, for example, Schober et al. [28], who also obtain the scaling u ∝ ν 1/3 . Krumholz et al. [9] compute the 1D turbulent velocity dispersion that can be sustained by star formation feedback (SNe) alone. Multiplying the approximate range they obtain by √ 3, assuming isotropy, gives u ∼ 10-17 km s −1 . Their no-transport models assume that turbulence is driven by star formation feedback alone, and that this is dominated by SN feedback, as assumed in our model. Assuming the star formation rate surface density to be proportional to ν, 6 their fixed star formation efficiency per free-fall time no-transport model predicts u to be independent of ν, while their fixed Toomre Q no-transport model predicts u ∝ ν 1/2 . The power law index of our model (1/3) thus lies in between that of those models (0 and 1/2). We caution that in our model u ∝ ν 1/3 only assuming that certain other parameters, like c s , are independent of ν.
For the fiducial parameter values Equation (33)

Correlation Time
We estimate the correlation time τ as the minimum of the turnover time τ e of energy-carrying eddies, and the time τ r for the flow to renovate due to the passage of an SN or SB blast wave, This prescription guarantees that the Strouhal number, which we define here as St = τ/τ e , will be less than or equal to 1 in our model.
The renovation rate is equal to the sum of the rates from isolated SNe and SBs, We take τ r SN to be equal to the average time between successive SNR shells (blast waves) passing through a given point x. In a statistical steady state, the rate at which blast waves cross x is equal to the rate of isolated SN explosions within a sphere of radius l SN = A(t SN s ) a , since more distant SNRs will break up before reaching x. We then obtain 4 3 πl 3 SN ν SN τ r SN = 1, or, solving for τ r SN , where ν 50 = ν/(50 kpc −3 Myr −1 ) and we have made use of equations (7) and (14) in the last equality. The same result can be obtained by computing how long it takes for SNRs to fill the volume, neglecting overlapping. Using the same approach to estimate the renovation time due to SBs, we find where we have made use of equations (11) and (15). For the fiducial parameter values, t SB b < t SB s , so τ r SB = 9.9 Myr, and τ r = [1/(6.8 Myr) + 1/(9.9 Myr)] −1 = 4.0 Myr. On the other hand, τ e = l/u = 5.9 Myr for fiducial parameter values. In this case τ = τ r since τ r < τ e .

Graphical Example
We illustrate the fiducial case (Table 1) in Figure 2, where we plot solutions as functions of H. Solid curves show the full solutions, while dashed and dotted curves show the contributions from SBs and isolated SNe, respectively. The breaks in the curves at H = 0.53 kpc are caused by the transition from SBs blowing out of the disc when the disc is thin, to being unable to blow out when the disc is thick. In Figure 2f, there is a second break at H = 0.30 kpc where τ = min(τ r , τ e ) transitions from τ = τ e for a thinner disc to τ = τ r for a thicker disc.

Exploration of the Parameter Space
We plot several example solutions in Figure 3, and refer to rows and columns by number starting from top left. Panels show a given turbulence parameter as a function of disc scale height H (see the figure caption for a detailed description). The fiducial case is illustrated by the solid black curves of the left-hand column, at H = 0.4 kpc.
The correlation scale l (top row) lies in the range 20 pc < l < 175 pc. It increases with scale height if H is small enough that SBs blow out, and is otherwise independent of H. The dependence is stronger than linear because of the termĖ i SN /Ė i SB ∝ H −3 in the denominator of Equation (29) (the same ratio is suppressed by the factor l SN /l SB in the numerator, so is less consequential there). Since l is independent of the SN rate density ν, dashed-dotted red and solid black curves overlap, as do short-dashed blue and long-dashed orange curves. Since l ∝ C, its value is somewhat sensitive to the choice of C (compare first and fourth columns).
The turbulent velocity u (second row) is found to occupy the range 3 km s −1 ≤ u ≤ 23 km s −1 , but u ≥ 7 km s −1 for n = 0.1 cm −3 . These values are in good agreement with observations and simulation results [e.g. Reference 36]. If other sources of turbulence, in addition to SNe, were important, then there would be additional contributions toĖ i , soĖ i would increase. If, further, such contributions injected energy on similar or larger scales compared to SNe, so that l remained roughly the same or larger, then, from the first equality in Equation (33), we would expect u to increase. Therefore, if additional sources of turbulence were important, our model would predict a larger value of u. Gravitationally-driven inward radial transport is one such source that is likely to be important, at least at high cosmological redshift [9]. The lower end of the range of u corresponds to dense gas (compare third and fifth columns), mainly owing to smaller l SN and l SB . The upper end of the range corresponds to large H, which prevents or delays blowout, allowing more power to be deposited into the ISM. Also, u increases with ν since more power is deposited (u ∝ ν 1/3 ). Note that doubling the sound speed from 10 km s −1 to 20 km s −1 can actually lead to a slight reduction in u because expanding SBs slow to the ambient sound speed earlier (compare first and second columns). If all SNe were to reside in SBs, and the SBs did not blow out, then l ∝ l SB and we would have u ∝ c −2/9 s . If all SNe were isolated, then we would have l ∝ l SN and u ∝ c 2/9 s . On the other hand, at small H, SBs blow out and the injection scale from SBs is independent of c s , while the injected energy from SBs scales as c 2 s , so u always increases with c s in this case. For the correlation time (third row), we find 1 Myr ≤ τ ≤ 13 Myr, and τ ≤ 6 Myr for n = 0.1 cm −3 , which agrees with order-of-magnitude estimates [6] and simulations [37]. The smallest values of τ are obtained for large values of ν, and the largest values of τ are obtained for large values of n. In the fourth row, we plot St = τ/τ e , which is ≤ 1 by design in our model. This quantity is found to lie in the range 0.2-1, which implies that τ r and τ e are generally comparable to one another. At small scale height τ e < τ r since the expansion of SBs, and hence the renovation rate, is reduced by blowout. For large values of ν, τ r < τ e is more likely. As n increases, τ r /τ e also increases, leading to τ = τ e for the examples with n = 1 cm −3 . Our results suggest that the common assumption that the correlation time is equal to the eddy turnover time is generally reasonable to within a factor of 2-4.
It is worth pointing out that separate relations would exist between the underlying parameters of our model. For example, ν and H are likely to be inversely related. Firstly, for a flared disc H is larger in the disc outer region, where the SN rate surface density is low. Secondly, the SN rate volume density decreases with H for a given SN rate surface density. Likewise H and n are expected to be inversely related, and ν would be expected to increase with n.
We next turn to the ratio of energies injected by SBs and isolated SNe,Ė i SB /Ė i SN , plotted in the fifth row. As this ratio can be less than or greater than unity, both types of driving can be important. The ratio is independent of the overall supernova rate ν, but depends on the fraction f SB of SNe contributing to SBs, from Equation (19). In the case where H is small enough that SBs blow out,Ė i SN /Ė i SB varies linearly with N SB and approximately linearly with E 0 . This is because the more energetic the SB, the fraction of its energy that is "wasted" when the supernova blows out of the disc increases. In particular, for N SB = 1000, n = 0.1 cm −3 , and H < 0.5 kpc, isolated SNe strongly dominate the energy injection. This becomes clearer by computing the efficiency of conversion of SN energy into turbulent energy, shown in the bottom row. For the examples plotted, the efficiency ranges between about 1% and 4%, and is smaller when SBs experience blowout.

Relative Importance of Isolated SNe and SBs and Effect of SN Clustering
It is interesting to consider how the results change if all SNe are assumed to be isolated, or conversely, if all SNe are assumed to reside in SBs. Thus, in Figure 4 we adopt the same underlying parameter values as for Figure 3, but now we plot the cases f SB = 0 (thin lines) and f SB = 1 (thick lines).
The correlation scale l is considerably smaller for the pure isolated SN case (top row), generally in the range 20-50 pc (in this case the thin lines coincide in each panel as l does not depend on N SB or ν). For the pure SB case, l increases linearly with H (Equation 11) at small enough H such that blowout occurs, and can reach values as large as 300 pc.
The value of u (second row) is similar for the pure SB and pure SNR cases, u ∼ 10-20 km s −1 for n = 0.1 cm −3 and u ∼ 5-7 km s −1 for n = 1 cm −3 , if blowout is avoided in the SB case. However if blowout does take place, values of u can be much smaller in the pure SB case (as low as 2 km s −1 for H = 200 pc and N SB = 1000), due to inefficient energy conversion from SNe to turbulence. The pure SB low H case is not very realistic though, because if even a small fraction of SNe were isolated, that channel could dominate the energy injection, raising u to values closer to the thin lines, as can be seen in Figure 3 where 1/4 of SNe are isolated.
The correlation time, like the correlation scale, is much smaller in the pure isolated SN case as compared to the pure SB case, as seen in the third row. In the fourth row, we see that τ is sometimes equal to τ e , and sometimes to τ r , but for large c s or large n, usually τ = τ e . The correlation time increases with stronger clustering of SNe (compare solid black to long-dashed orange and dash-dotted red to short-dashed blue), and with reduced SN rate ν (compare red with black and blue with orange). 7 The ratio of injected energies for the pure isolated SN and pure SB cases, shown in the fifth row, are similar to those shown in Figure 3, except lower by a factor of three. This is because in the model plotted in Figure 3, f SB /(1 − f SB ) = 3 times more SNe reside in SBs compared to those that are isolated.
In the bottom row we see that the conversion efficiency of SN energy to turbulent energy is very low for the case f SB = 1 if H is small. The efficiency for f SB = 1 becomes larger with increasing H, but even neglecting blowout, SBs are less efficient than isolated SNe in transferring SN energy to turbulence. Hence, clustering of SNe leads to less efficient turbulence driving in our model. The efficiencies obtained for the pure isolated SN case ( f SB = 0) are generally consistent with the estimate of 4% made by Dyson and Williams [43,Ch. 7].
If we ignore clustering of SNe ( f SB = 0) so that all SNe are isolated, we obtain, for fiducial parameter values, l = 42 pc, u = 13 km s −1 , τ = τ r = 2 Myr and = 0.037, whereas if all SNe are assumed to reside in SBs ( f SB = 1) we obtain l = 120 pc, u = 12 km s −1 , τ = τ r = 7 Myr and = 0.008. If instead we retain f SB = 3/4 but only the isolated SNe are assumed to drive turbulence and those residing in SBs are neglected, l and are unchanged from the f SB = 0 case but u = 8 km s −1 and τ = τ e = 5 Myr, whereas if isolated SNe are neglected and only the SBs are assumed to drive turbulence, l and are unchanged from the f SB = 1 case but u = 11 km s −1 and τ = τ r = 10 Myr.

Varying the Parameter ξ
In our estimates above we parameterized the driving scale of an SB that blows out as ξ H, with H the gas scale height. We then made the choice ξ = 1, since an SB which blows out has expanded to a horizontal scale ∼ H, according to simulations of Mac Low and McCray [47] (c.f. their figure 8). Those authors made use of a two-component density model from Lockman et al. [54] which consists of a diffuse exponential component of scale height H (the Lockman layer) and a thin Gaussian cloud layer with scale height a few times smaller than H. They found in their simulations that while the vertical extent of an SB at blowout is a few times larger than H, it experiences a pinch due to the greater density at the midplane so that its horizontal expanse is roughly equal to H. In this section we consider adopting a smaller value of ξ, corresponding to a more severe pinch near the midplane. In Equation (11) it is implied that an SB cannot expand to a horizontal scale > ξ H because once it reaches that scale it blows out. Before blowout, SB expansion is independent of ξ, so a smaller ξ causes SBs to blow out earlier.

Motivation for Reducing ξ
It appears to violate self-consistency that turbulence can be driven in our model at a scale a few times larger than the cloud layer, since it would seem that the cloud layer would then be disrupted by such turbulence. However, the current understanding of interstellar clouds is still fragmentary, and it is possible that they are transient and formed by compressions in the SN-shocked transonic gas [e.g 14,18]. In any case, the entrainment of clouds by SBs can lead to enhanced cooling in the SB interior. Cloud evaporation in the SB interior can occur through thermal conduction, while mixing between clouds and hot gas can lead to radiative cooling. These processes cause a reduction in the interior pressure of the SB [47,67,68]. Since SBs interact energetically with the thin cloud layer, which has a scale height h ∼ 130 pc in the Solar neighbourhood [40,69], the driving scale could plausibly be capped at h ∼ 1 3 H by setting ξ to be ∼ 1/3 instead of 1. In this scenario an SB would experience a more severe pinch than predicted by Mac Low and McCray [47], presumably leading to narrower chimneys, as seen in simulations by de Avillez and Breitschwerdt [70], where the horizontal size of the chimneys of hot gas produced by clusters of SNe in the local ISM is about 150-200 pc.

Results
With this in mind, we present results of a model that is the same as that presented in Sections 2 and 3 except that ξ = 1/3 instead of 1. Results of this alternative small-ξ model are shown in Figure 5, which is otherwise similar to Figure 3. The efficiency is lower in the ξ = 1/3 model (bottom row) as compared to the ξ = 1 model because SBs blow out early in their expansions and thus transfer very little energy to the ISM. In fact, the results are very similar to those obtained if turbulence driving by SBs is completely neglected, but f SB = 3/4 is still assumed. This case has not been plotted, but it is just the H → 0 limit of the plots. When blowout occurs, the differences between the ξ = 1/3 and ξ = 1 cases are more significant for larger values of H, where energy injected by SBs is a larger fraction of the total injected energy. Only in the last column of Figure 5, where n = 1 cm −3 , can blowout be suppressed at large H due to the high ambient density.
Adopting ξ = 1/3 generally leads to smaller l due to the reduction in large-scale driving, and smaller u due to the reduction in energy input and in l, compared to ξ = 1. The value of τ is not greatly affected. For ξ = 1/3 these quantities depend only weakly on H and N SB , since these parameters only affect SBs. The values obtained when ξ = 1/3 and other parameters fiducial are l ≈ 42 pc, u ≈ 9 km s −1 , τ = τ e ≈ 5 Myr and ≈ 0.009 (first column solid black at H = 0.4 kpc), as compared with l ≈ 74 pc, u ≈ 12 km s −1 , τ = τ r ≈ 4 Myr and ≈ 0.016 for ξ = 1 (Figure 3).

Two-layer Model
A further possibility would be to generalize the model to consist of two distinct gas layers, namely the cloud and Lockman layers, each with a different set of properties. This treatment would allow for the possibility that average turbulence parameters vary with distance from the midplane |z|. Naturally l (and ξ) would increase with |z| because SBs would bulge out into the more tenuous outer layer. However, in a two-layer model one would need to consider how SN energy gets divided between the two layers. One would also need to consider the possibility of individual SNe blowing out of the thin layer and into the thick layer. Another possibility is to include explicitly the stratification of the ISM.

Limitations and Opportunities for Extending the Model
On account of keeping the model reasonably simple, certain details have necessarily been omitted. Expanding SN shocks may reflect off of interstellar clouds to produce secondary shocks [71]. The distribution function P(µ) of secondary shocks with Mach number ≥ µ was calculated by Bykov and Toptygin [72], who assumed secondary shocks to be weak. The distribution function derived diverges unphysically for µ → 1, which leads to a vanishing renovation time (Equation 36). Clearly τ r would decrease if secondary shocks were included, but including them would require a new model that remains valid as µ → 1, which is beyond the scope of our simple treatment here. On the other hand, we have also neglected mutual isolated SNe-SNe, isolated SNe-SB or SB-SB interactions, as we have assumed that every spherical shock propagates independently of neighbouring shocks. Such interactions probably make shock expansion less efficient and would thus lead to an increase in τ r , which could help to offset the omission of secondary shocks.
We have also chosen to omit the possible influence of galactic shear on the correlation time. In their galactic dynamo model, Zhou and Blackman [73] adopt τ −1 = (τ e ) −1 + qΩ, where Ω(r) is the large-scale angular velocity of gas at radius r, and q(r) = −d ln Ω/d ln r. Here q = 1 corresponds to a locally flat rotation curve. In the vicinity of the Solar neighbourhood, Ω ∼ 20-30 km s −1 kpc −1 , which translates to a shearing timescale of ∼ 30-50 Myr. This is large compared to our estimates for τ, so shear would not cause an important reduction in τ. However, this effect might be important for parts of the phenomenologically relevant parameter space. A more detailed model would also allow for a distribution of N SB , rather than considering only the extreme case of isolated SNe and SBs each containing the same number of SNe, as we have done here [74][75][76]. Whereas we have treated isolated SNe and SBs as different kinds of object with different similarity solutions governing their evolutions, there would in reality be a transition from one type to the other as N SB is increased from 1 to larger values.
Another possibility is to include effects of the magnetic field on the expansion of SBs [e.g. 77], which would provide nonlinear feedback in a dynamo model.
We have assumed that SB blowout occurs instantaneously when R SB = ξ H, with ξ = 1 for our fiducial model. At this time, the SB merges with the ISM, driving turbulence. A more sophisticated treatment would model the blowout phase from the time the SB begins to lose pressure support to the time it breaks up near the midplane. Using the similarity solution for the SB temperature as a function of radius [47], we compute the sound crossing time for the SB at blowout to be ∼ 4 Myr for our fiducial parameter values, as compared with t SB b = 15 Myr, so modeling this phase could lead to important differences. Our model assumes a uniform ISM of a given disc semi-thickness H, so does not differentiate between the various ISM phases. Thus, it effectively averages over the individual locations of isolated SNe, both in terms of distance from the midplane and occurrence inside or outside molecular clouds. Furthermore, it does not distinguish between the average properties of the ambient medium encountered by SBs as opposed to that encountered by isolated SNe, or between the average properties of the medium in which SNRs and SBs expand as opposed to that in which the turbulence is primarily driven. Refinements to the model to address theses shortcomings would entail introducing more parameters, and it is not clear whether this would be warranted given the various uncertainties.
More fundamentally, our model could also be extended to include turbulence driving by mechanisms other than SN feedback.

Summary and Conclusions
Various astrophysical phenomena, including galactic dynamos, are sensitive to the parameters of interstellar turbulence. However, determining the values of these parameters, as well as their dependencies on other galactic or interstellar medium parameters, has remained an elusive goal. We have applied standard similarity solutions for SNRs and SBs to a model of SN-driven interstellar turbulence, and obtained simple analytic expressions for the velocity correlation scale l, RMS turbulent speed u, velocity correlation time τ, and SN to turbulent energy transfer efficiency .
Our main motivation for this work was to extend dynamo models such that they can be parameterized by quantities that are more accessible than l, u and τ, namely, the underlying model parameters summarized in Table 1. In this way, l, u and τ would become intermediate quantities computed within the dynamo model, rather than input parameters. As it is beyond the scope of the present work, we leave the implementation of this idea for the future. This could be done using already available dynamo models [e.g. [78][79][80][81], but the underlying ISM parameters of Table 1, including their spatial and temporal dependencies, would need to be modelled independently and/or constrained using observations. The relations between turbulence parameters and underlying ISM parameters that we have derived are testable and may have many applications besides dynamos. However, care must be taken to account for separate relations between the underlying parameters themselves. No attempt is made to include these additional dependencies in our model.
In Figure 3, we present solutions over a large region of the underlying parameter space. We obtain the following approximate ranges for the quantities computed: l ∼ 20-175 pc, u ∼ 3-23 km s −1 , τ ∼ 1-13 Myr and ∼ 0.01-0.04. For our fiducial set of underlying parameter values, applicable to the Solar neighbourhood, we obtain l ≈ 74 pc, u ≈ 12 km s −1 , τ ≈ 4 Myr and ≈ 0.016, which can be compared with "canonical" order-of-magnitude estimates of l ∼ 100 pc, u ∼ 10 km s −1 , τ ∼ 10 Myr and ∼ 0.04, and recent local ISM simulations of Hollins et al. [37] which find l ≈ 74 pc, u ≈ 13 km s −1 and τ ≈ 5 Myr. The close agreement with the latter results is partly just a coincidence, but our results are in broad agreement with other simulations as well [e.g. 34,35,82].
Our model includes both isolated SNe and SBs, with a fraction f SB = 3/4 of SNe residing in SBs. Turbulent energy is deposited in the ISM once the expansion velocity of an SNR or SB reaches the ambient sound speed. We find that isolated SNe and SBs can both contribute significantly to the turbulent energy injection in the ISM. This suggests that both of these sources should be accounted for in models of interstellar turbulence. 8 However, one is free to choose to include only one or the other of these components. In particular, the model is simpler and contains many fewer parameters if clustering of SNe to form SBs is neglected ( f SB = 0).
The evolution of SBs depends on the density scale height H of the warm diffuse ambient gas (the Lockman layer in the Milky Way), and SB blowout happens in our model when the SB radius in the midplane is of order the ambient scale height [47]. For H 0.5 kpc, SBs tend to blow out of the disc, and consequently, a smaller fraction of their energy ends up in interstellar turbulence. Blowout is assumed to result in a rapid loss of pressure support near the midplane, and a sudden reduction of the expansion speed to the ambient sound speed, at which point the SB merges with the ISM. If blowout happens early, i.e. when the SB expansion speedṘ SB c s , then the energy deposited into the ISM is only a small fraction ∼ H 3 /R 3 SB (t SB s ) of what would have been injected had the SB been able to expand up to the time t SB s wheṅ R SB (t SB s ) = c s . In this case, isolated SNe dominate the energy injection into turbulence in our model. We also computed the fraction of SN energy that ends up in turbulence, and found this efficiency factor to typically lie in the range = 0.01-0.04. This agrees closely with the range found by Bacchini et al. [10] for the warm atomic gas. The efficiency of conversion of SN energy to turbulent energy is found to be lower when SBs are included, compared to a model for which all SNe are isolated. This is because SBs can blow out and also because SBs drive turbulence less efficiently than isolated SNe even when they do not blow out, as seen in the bottom row of Figure 4. Since most SNe can be contained within SBs, the reduction in the overall efficiency can be significant. When SN clustering is absent ( f SB = 0) in our model, all SNe are isolated. This results in values of l and τ that are smaller than the f SB = 3/4 case, and values of that are larger, while the opposite is true if f SB = 1. It could be argued that large driving scales comparable to H cannot be present because the resulting turbulence would then disrupt the thin HI cloud layer, which has a scale height of ∼ 130 pc in the Solar neighbourhood. This motivated us to consider a variation on the model in which an SB blows out earlier, when its radius in the midplane is equal to ξ H with ξ = 1/3 instead of 1. As ξ is reduced, SBs end up being less important in transferring energy to turbulence in the disc. For ξ = 1/3, we obtain values of l up to a few times lower than in the ξ = 1 case, and values of u up to two times smaller, but τ remains about the same.   8 Yoo and Cho [83] study MHD simulations with forcing on two scales and find that even a relatively small amount of energy injection on the larger scale can have important effects on the properties of the turbulence.