Abstract
We study supernovadriven galactic outflows as a mechanism for injecting turbulence in the intergalactic medium (IGM) far from galaxies. To this aim, we follow the evolution of a 10^{13}M_{⊙} galaxy along its merger tree, with carefully calibrated prescriptions for star formation and wind efficiencies. At z ≈ 3, the majority of the bubbles around galaxies are old (ages > 1 Gyr), that is, they contain metals expelled by their progenitors at earlier times; their filling factor increases with time, reaching about 10 per cent at z < 2. The energy deposited by these expanding shocks in the IGM is predominantly in kinetic form (mean energy density of 1 μeV cm^{−3}, about two to three times the thermal one), which is rapidly converted in disordered motions by instabilities, finally resulting in a fully developed turbulent spectrum whose evolution is followed through a spectral transfer function approach. The derived mean IGM turbulent Doppler parameter, b_{t}, peaks at z ≈ 1 at about 1.5 km s^{−1} with the maximum b_{t}= 25 km s^{−1}. The shape of the b_{t} distribution does not significantly evolve with redshift but undergoes a continuous shift towards lower b_{t} values with time, as a result of bubble ageing. We also find a clear trend of decreasing b_{t} with and a more complex dependence on R_{s} resulting from the age spread of the bubbles. We have attempted a preliminary comparison with the data, hampered by the scarcity of the latter and by the challenge provided by the subtraction of peculiar and thermal motions. Finally, we comment on the implications of turbulence for various cosmological studies.
1 INTRODUCTION
The intergalactic medium (IGM) is a pervasive, diffuse cosmic baryonic component out of which galaxies form by accretion. In turn, the IGM is replenished with gas and heavy elements carried by galactic outflows. Such dynamical processes hence encode a record of the complex galaxy–IGM interplay.
Observations of metal absorption lines (e.g. C iii, C iv, Si iii, S iv and O vi) in quasar spectra are a direct probe of this evolution. They show that regions of enhanced IGM density far from large galaxies are polluted with nonnegligible amounts of heavy elements (e.g. Songaila & Cowie 1996; Davé et al. 1998; Schaye et al. 2000a; Pettini et al. 2003; Schaye et al. 2003; Aguirre et al. 2004; Aracil et al. 2004; Simcoe, Sargent & Rauch 2004). In addition, Lyα forest clouds with the neutral hydrogen column density are metal enriched to Z ≈ 10^{−3}–10^{−2} Z_{⊙} (Cowie et al. 1995; Songaila & Cowie 1996; Ellison et al. 2000); finally, the carbon and silicon abundances of such systems remain roughly constant throughout the redshift range 1.5 < z < 5.5 (Songaila 2001; Pettini et al. 2003), although the first signs of a slight decrease above z = 6 have been now tentatively reported (RyanWeber et al. 2009). More recently, at variance with the constant behaviour picture, D’Odorico et al. (2010) reported a significant increase in the C iv abundance towards lower redshifts z ≲ 2. This result has been also confirmed by Cooksey et al. (2010) for z ≲ 1.
Since metals are created in stars inside galaxies, their transport into these distant regions must rely on some yet unknown mechanisms. The most likely one is supernovadriven (SNdriven) galactic outflows (Mac Low & Ferrara 1999; Ferrara, Pettini & Shchekinov 2000, hereinafter FPS; Aguirre et al. 2001b; Oppenheimer & Davé 2006; Cen & Chisari 2010), although different alternatives have been proposed, that is, dust sputtering (Aguirre et al. 2001a; Bianchi & Ferrara 2005) or merging (Gnedin 1998). This hypothesis is supported by a series of observations of galactic winds (e.g. Heckman 2001; Pettini et al. 2001; Adelberger et al. 2003; Shapley et al. 2003; Erb et al. 2006) from which one can conclude that (on some spatial scale) galactic wind velocities range from hundreds to thousands of km s^{−1} and the massloss rates are comparable to star formation rates (SFRs).
A complementary piece of information comes from the wellestablished fact that the interstellar medium (ISM) in galactic discs is turbulent (Larson 1981; Scalo 1987; Dickey & Lockman 1990; Elmegreen & Scalo 2004). In the Milky Way (MW), the observed vertical distribution of gas, and especially the cold neutral H i component, cannot be supported by thermal pressure alone (Lockman & Gehman 1991; Ferrara 1993) and an additional pressure contribution must be advocated. Furthermore, pulsar scintillation experiments have convincingly demonstrated that electron density fluctuations in the ionized ISM are characterized by a Kolmogorov spectrum on scales ranging from au to kpc (Armstrong, Rickett & Spangler 1995; Han, Ferriere & Manchester 2004). The Kolmogorov spectrum is the spectrum of the density fluctuations which arise from turbulent motions in a compressible gas with constant dynamical energy exchange between irregularities on different scales and under very general conditions, it assumes the form of a power law (∼k^{−11/3}, where k is the wavenumber associated to a given scale L) within the inertial range. In other spiral galaxies, the vertical velocity dispersions (derived essentially from H i observations) are observed to vary radially from 12–15 km s^{−1} in the central parts to 4–6 km s^{−1} in the outer regions and, in most cases, these values exceed those expected by the thermal broadening of the H i emission line (Dib, Bell & Burkert 2006).
Often invoked to provide the required additional pressure are turbulent motions of the interstellar gas, magnetic fields and cosmic rays (Boulares & Cox 1990; Lockman & Gehman 1991; Dib, Bell & Burkert 2006; Joung, Mac Low & Bryan 2009). Turbulence is likely to be the mostrelevant pressure contribution; its driving may be due to different sources (star formation, stellar outflows, instabilities) acting on specific characteristic scales and injecting different amounts of kinetic energy into the medium; among these, SN explosions are thought to be the dominant agent, at least in the starforming regions of galaxies (McKee & Ostriker 1977; Norman & Ferrara 1996; MacLow 2004; de Avillez & Breitschwerdt 2007). If both turbulence and metal enrichment are predominantly driven by the same physical phenomenon, that is, kinetic energy deposition by SN shocks, then a natural expectation is that the IGM gas in enriched regions should show signatures of a turbulent regime.
A direct way to measure the turbulence in the IGM is to look for velocity differences on the smallest spatial scales accessible to observations; to this scope, Rauch, Sargent & Barlow (2001) in a pioneering experiment used lensed quasars. They observed adjacent C iv profiles in paired lines of sight and found velocity differences of δv ≈ 5 km s^{−1} on scales of 300 pc at redshift ∼2.7. They apply the Kolmogorov steadystate assumption to calculate an energy transfer rate ε≈ 10^{−3} cm^{2} s^{−3} and a turbulent dissipation timescale as short as ≈10^{8} yr, thus implying a recurrent stirring of the gas by some sort of galacticscale feedback.
Another hint of a significant turbulent contribution to the IGM kinetic budget comes from the fact that the median Doppler parameters measured in the Lyα forest are significantly larger than those predicted by cosmological simulations. Again, this implies that some energy in nonthermal form must be injected in the gas to explain the observed line broadening (Meiksin, Bryan & Machacek 2001; Meiksin 2009).
Quite surprisingly, the properties of intergalactic turbulence have received so far relatively little attention in spite of the insights that its detailed understanding might provide into the galaxy–IGM interplay (FPS; Oppenheimer & Davé 2009; Scannapieco & Brüggen 2010). Oppenheimer & Davé (2009) trying to model lowz O vi absorbers concluded that their properties fit observables if and only if subresolution turbulence (in practice, a densitydependent turbulent b parameter) is added at a level which increases with the O vi absorber strength (and hence presumably closer to highdensity regions hosting galaxies). In addition, stronger absorbers arise from more recent outflows. It is important to note that turbulence, if present, can also profoundly modify the enrichment patterns through diffusion processes, as pointed out already by FPS. Additional important information might come from investigations of structure formation, during which the released gravitational energy is transferred to the IGM in different forms (e.g. gas entropy, instabilities, cosmicray acceleration); in particular, turbulence is induced by the vorticity cascade originating at cosmological shocks. This has been investigated throughout cosmological simulations by Ryu et al. (2008) and Zhu, Feng & Fang (2010) who derived the average magnetic field strength and turbulent pressure, respectively, in the overdense IGM regions outside clusters/groups.
In this work, we try to fill the aforementioned gap by investigating in detail the SNdriven outflow scenario as a mechanism to pump and sustain turbulence in the IGM surrounding highredshift galaxies. We investigate the evolution of winds via a semianalytic approach following both the galaxy evolution along its hierarchical growth and the expansion of SNdriven superbubbles as they escape from the halo potential well. With such results in hand, we model the turbulence evolution within wind shells using a novel approach based on the spectral transfer equation. This approach allows to derive the turbulent energy density deposited by galactic winds in the IGM, along with its spectral features and dissipation timescales, and to predict the corresponding thermal/kinetic evolution of the IGM.
2 THE MODEL
In this section, we describe how we determine the evolution of the turbulent energy deposited by galactic outflows into the IGM. First, we introduce the galaxy formation model, allowing a precise description of both the star formation history and the SN feedback along the galaxy merger tree (MT). These results are then used as inputs to model the evolution of pressuresupported bubbles through the solution of dynamical equations and the turbulence energy spectrum evolution within them.
2.1 Cosmic star formation history
The semianalytic model we introduce in this section to describe the evolution of a galaxy along its MT is similar to the onezone model by Salvadori, Schneider & Ferrara (2007) with the improvements introduced by Salvadori, Ferrara & Schneider (2008). Since in this work we are not interested in following the chemical enrichment in detail, we use a simplified approach described in the following.
In hierarchical models of structure formation, such as Λ cold dark matter (ΛCDM), the formation of a DM halo through accretion and repeated mergers can be described by a MT (Lacey & Cole 1993). MTs, which list the progenitors of a given halo at different redshifts and describe how and when these merge together, contain essentially all the necessary information about the DM content of haloes to build a realistic model of galaxy formation as demonstrated by a longstanding practice in the field (e.g. Kauffmann & White 1993; Somerville & Kolatt 1999).
The MTs can be extracted from Nbody simulations (Springel et al. 2005) or generated by Monte Carlo algorithms (see e.g. Sheth & Lemson 1999). In this work, we compute the mass growth histories of DM haloes using the second method using a public code1 based on the extended Press–Schecter (PS) formalism (Cole et al. 2000), as improved by Parkinson, Cole & Helly (2008), to which we refer the reader for numerical details. As an input, we use the WMAP52 cosmological parameters and power spectrum. As an output, we obtain many realizations of the MT for a given halo mass M_{0} at a certain final redshift z_{f}; each realization lists all progenitors of M_{0} at different redshifts, following the related merging histories down to a mass of M_{res}.
In our model, all the protogalaxies virialize with a gastoDM ratio of 1:5. The neutral gas in these minihaloes cannot cool via atomic hydrogen and relies on the presence of molecular hydrogen, H_{2}, to cool and collapse, ultimately forming stars. Since we assume that feedback effects rapidly suppress star formation in the first minihaloes and that only Lyα cooling haloes (T_{vir} > 10^{4} K) contribute to the star formation history of the galaxy, we choose M_{res}= M_{4} (z = 20), where M_{4}(z) is the mass corresponding to a halo with virial temperature of 10^{4} K at redshift z, given
(for an exact expression, see Barkana & Loeb 2001) and z = 20 is the starting redshift of our simulation. Note that the value of M_{res} sets the limit between progenitors and mass accretion, since the galaxies above M_{res} experience star formation and stellar feedback, ultimately changing their gas content, whereas objects below this threshold retain their original gas content fraction during their evolution, which is finally inherited by another halo in the next hierarchy level.The star formation in gas clouds occurs on a freefall timescale t_{ff} = (3 π/32 G ρ)^{1/2}, where G is the gravitational constant and ρ is the average mass (dark+baryonic) density inside the halo assumed to have a virial radius (r_{vir}) as given in Bryan & Norman (1998). The star formation efficiency of a galaxy is then modelled as a fraction, ε_{*}, of the freefall time:
where ε_{*} represents a free parameter of the model and M_{g} represents the gas mass in the halo which has not yet been converted into stars.According to the standard scenario for galaxy formation, the gas inside the galaxy is depleted by various feedback processes, the most important being gasloss driven by SN explosions. In fact, SN explosions may power a wind which, if sufficiently energetic, may overcome the gravitational pull of the host halo, leading to the expulsion of gas and metals into the surrounding IGM. To model this process, we compare the kinetic energy injected by SNdriven winds with the minimum kinetic energy of a mass M_{w} to escape the galactic potential well (Larson 1974; Dekel & Silk 1986; Ferrara & Tolstoy 2000). The mass of the gas ejected from the galaxy is computed from the equation
where is the kinetic energy injected by SNe and is the escape velocity of the gas from a halo with mass M_{h} and binding energy E_{b} given by (Barkana & Loeb 2001):In equation (4), ε_{w} is a free parameter which controls the conversion efficiency of SN explosion energy in kinetic form, N_{SN} is the number of SNe and 〈E_{SN}〉 = 1.2 × 10^{51} erg is the average SN explosion energy (not in neutrinos). Differentiating equation (3) with respect to time, we find that the gas ejection rate is proportional to the SN explosion rate:
In our model, we only consider Type II SNe (pairinstability SNe are neglected as Population III stars do not contribute significantly to the total SFR, Tornatore, Ferrara & Schneider 2007); hence, we calculate N_{SN} by integrating the adopted initial mass function (IMF) (see below) in the canonical SN range 8–100 M_{⊙}. For any starforming halo of the galaxy hierarchy, we therefore solve the following system of differential equations: The first equation simply defines the SFR. In our model, we have assumed the Instantaneous Recycling Approximation (IRA, Tinsley 1980), according to which stars are divided into two classes: those which live forever, if their lifetime is longer than the time since their formation and those which die instantaneously, eventually leaving a remnant. The transition mass between the two possible evolutions, or turnoff mass (m_{to}), has been computed at any considered redshift. Under the IRA, the returned fraction, that is, the stellar mass fraction returned to the gas through winds and SN explosions, is where ϕ(m) is the IMF of the newborn Population II/I stars and it is assumed to have the form with x =−1.35, m_{cut} = 0.35 M_{⊙} and m in the range [0.1, 100] M_{⊙} (Larson 1998). The quantity w_{m}(m) represents the mass of the stellar remnant left by a star of mass m which explodes as an SN. We have used the grid of models by van den Hoek & Groenewegen (1997) for intermediatemass stars (0.9 < m < 8 M_{⊙}) and by Woosley & Weaver (1995) for massive stars (8 < m < 40 M_{⊙}). With the adopted IMF, the corresponding SN energy per unit stellar mass formed is 1.36 × 10^{49} erg M.The second equation (equation 8) describes the mass variation of cold gas: the latter increases due to accretion [] and it decreases because of both astration and massloss due to SN winds. To model gas accretion, we assume that if a new halo virializes out of the IGM gas, then it forms with a cosmological gas fraction f_{b}=Ω_{b}/Ω_{m} = 1/6, whereas if it results from the merging of two already existing haloes, then its gas content is simply the sum of the gas mass of the progenitors. Similarly, during the merger, the stellar mass of the new galaxy is assumed to be the sum of the stellar masses of the progenitors. The model free parameters (ε_{w} and ε_{*}) are fixed to match the global properties of the MW as we will discuss in more details in the next section.
2.2 Bubble evolution
In our simulation, galactic outflows are treated as pressuredriven bubbles of hot gas emerging from starforming galaxies. They expand by working against the IGM pressure and are driven by the energy injected by multiSN explosions. Most of the sweptup mass, both in the early adiabatic and in the following radiative phases, is concentrated in a dense shell bounding the hot overpressurized interior.
For this reason, galactic bubbles are canonically studied by using the thinshell approximation in which the shell expansion is driven by the internal energy (E_{b}) of the hot bubble gas (Ostriker & McKee 1988; Madau, Ferrara & Rees 2001; Bertone, Stoehr & White 2005; Samui, Subramanian & Srianand 2008).
In this work, we follow the modellization already presented in Madau et al. (2001) (see their section 2). However, in our model, we derive the SFR evolution from a constrained galactic model and from the SFR we derive the mechanical luminosity defined as L(t) = dE/dt, where E is the energy produced by the total contribution of N_{SN} SNe with energy 〈E_{SN}〉 and an efficiency ε_{w}.
In the following, we remind the reader of the relevant equations used to follow the evolution of the interesting bubble properties, namely the bubble radius R_{s}, the shock velocity V_{s} and the internal bubble temperature T_{b}:
where G is the gravitational constant, M(R) is the galactic DM mass assumed to follow a Navarro–Frenk–White profile within a radius R, ρ (P) is the density (pressure) of the ambient medium, is the volume enclosed by the shell, 3 is the pressure of the bubble gas assumed to be isothermal and with adiabatic index 5/3, is the cooling rate per unit volume of the hot bubble gas (whose average hydrogen density and temperature are and , respectively), and C_{1} = 16πμ m_{p}η/25 k and C_{2} = (125/39) πμ m_{p} [where η = 6 × 10^{−7} (cgs units) is the classical Spitzer thermal conduction coefficient, assuming a Coulomb logarithm equal to 30, and μ = 0.59 is the mean molecular weight of a primordial ionized H/He mixture].Finally, the cooling function [represented by Λ(T) in equation 12] depends on the hot bubble density, temperature and metallicity (Z). In principle, one could compute the average metallicity of the bubble, but for simplicity we assume for this quantity a constant value Z = 0.1 Z_{⊙}, noting that the final results are very slightly dependent on this choice within a reasonable range. The cooling rates due to gas radiative processes are taken from Sutherland & Dopita (1993); the other relevant cooling agent is inverseCompton cooling off cosmic microwave background photons (Ikeuchi & Ostriker 1986), which is dominant at higher redshifts.
The first two equations (equations 11 and 12) derive from momentum and energy conservation. In particular, the righthand side of equation (11) represents the momentum gained by the shell from the SNshocked wind which expands against the external pressure and gravitational attraction; the righthand side of equation (12) describes the mechanical energy input, the work done against the shell and the energy losses due to radiation. The third equation (equation 13) is obtained by equating the rate at which the gas is injected from the shell into the cavity with the conductive evaporation rate. Note that, as the bubbles sizes are always much smaller that the horizon scale, the cosmological terms in the above equations can be safely neglected.
Another important quantity to determine the fate of the wind is the total mass in the shell, M_{s}, since the energy required to accelerate it increases with its mass. The shell grow rate is the wind mass deposition rate (equation 3) and the rate at which gas is swept up from the galactic environment (for R_{s}≤ r_{vir}) or from the IGM (R_{s} > R_{vir}):
The ambient gas density is taken to be equal to the halo gas density distribution (assumed to have an isothermal profile) within R_{vir} and to the IGM background density outside the virial radius. To determine the pressure of the ambient medium, we must further assume that the ambient gas behaves as an ideal isothermal gas which, inside haloes, is at the virial temperature (Navarro et al. 1997) and outside R_{vir}, the IGM is taken to be photoheated by the ultraviolet (UV) background radiation to a temperature calculated as in, for example, Choudhury & Ferrara (2006). Note that even prior to reionization, the IGM into which the bubble expands will be locally ionized and heated to roughly the same temperature by the SN progenitor stars.
However, we note here that the bubble reaches very rapidly the virial radius, so that these bubbles are practically expanding in the IGM for most of their time. To see this, we consider as in Madau et al. (2001) a typical object virializing at high redshift (z = 9) with an halo mass M = 10^{8}h^{−1} M_{⊙}. At these epochs, the DM halo of a subgalactic system will be characterized by a virial radius
The Sedov solution (which is a good approximation of the equation systems when the ambient gas pressure, gravity and cooling can be neglected) predicts that the shell radius evolves according to R_{s} = (125/154π)^{1/5} (Lt^{3}/ρ)^{1/5}; assuming a constant luminosity of 10^{38} erg^{−1}, we find that the time taken to reach R_{vir} is ∼10^{7} yr. This time happens to be smaller than the Hubble time at that redshift.The velocity of the shock front decreases during the time because of the adiabatic expansion and cooling of the bubble; it might happen that the shock velocity equals the intergalactic gas sound speed; in this case, similarly to Bertone et al. (2005), we assume that the dynamics of the wind join the Hubble flow and no more mass is accreted.
Since equations (11)–(13) require an initial condition which cannot be 0 (no bubble) to be consistently solved, we used the Sedov solution for a very short time with respect to t_{H} (z = 20) to obtain the primordial bubble to follow the evolution. The evolutionary equations (equations 11–13) can be integrated numerically to follow the evolution of the shell along with the thermodynamic properties of the bubble for a given halo.
A final point concerns the treatment of bubbles when two haloes merge along the hierarchical MT. In that case, we impose mass conservation, that is, the mass of the new shell is the sum of the masses and of the two progenitor bubbles, that is, . Similarly, the final volume is taken to be equal to the sum of the volumes of the two single cavities and by the adopted spherical symmetry, we update the shock radius as . Finally, the shock velocity is given by the momentum conservation, which states that , where , and the new internal energy is the sum of the progenitor bubble internal energies, that is, . The gas temperature in the merged bubble is computed by assuming that the final thermal pressure is the sum of the thermal pressures of the merging bubbles: , where .
2.3 Turbulence evolution
In order to model the evolution of the turbulence developed in the expanding shells, we adopt in the following an approach based on the spectral transfer equation derived by Norman & Ferrara (1996) based on the hydrodynamic Kovasznay approximation. In fact, our aim is to derive averaged properties of the turbulent spectrum.
It is useful to introduce a spectral representation of the velocity field in the fluid. Then the energy density in eddies with wavenumbers between k and k + dk is given by , where ρ is the density of the fluid and is the turbulence spectrum.
The timeevolution of this quantity, under the assumption of local mode interaction in every region of the spectrum, is given by the following equation:
where ν is the kinematic viscosity, is the source function and α is a dimensionless constant that should be fixed experimentally.In order to describe the evolution of the turbulent spectrum, we finally provide the appropriate boundary and initial conditions:
A caveat must be made here. The assumption underlying the above description is one of incompressible, homogeneous, isotropic turbulence. This assumption clearly fails in the highly compressive IGM. However, the bulk of our knowledge on turbulence comes from terrestrial laboratory experiments and most of them deal with liquids; in addition, simulations of compressible turbulence (Lazarian & Cho 2004; Kowal & Lazarian 2007; Schmidt et al. 2009) have shown that nonlinear interactions rapidly transfer most of the energy to noncompressible modes. One can use those results at least as a reasonable guide when discussing compressible turbulence.
To model , we assume that turbulent motions in the outflow are induced by interacting blast waves, that is, we maintain that observed turbulent motions in the shell are ultimately derived from the kinetic energy of SNinduced shock waves which can act as source function for the turbulent cascade with an efficiency of the order of unity. If this is the case, there are some general constraints on the form of the source function. Bykov & Toptygin (1987) have studied in detail the shockinduced turbulence phenomenon and concluded that the turbulent spectrum has, for the McKee & Ostriker (1977) model for the SN shock wave expansions, a k^{−2} dependence in the shortwavelength regime, while for long wavelengths, it is proportional to k^{2}. Thus, the simplest source function retaining this behaviour is
where f(x) = x^{2}/(1 + x^{4}), k_{0} is the wavenumber corresponding to the characteristic length at which the turbulence is injected (and assumed to be k_{0}∼ 1/R_{s}) and S_{0}(t) is a normalization factor obtained by equating the kintegral of the source function to the kinetic energy rate (). Finally, as in Brandenburg, Korpi & Mee (2007), we use a constant kinematic viscosity of ν = 5 × 10^{24} cm^{2} s^{−1}; this is equivalent to scaling the dynamic viscosity coefficient with density.The spectrum evolution is followed for each halo of the MT; again, we need to specify how we treat merger events. To conserve the turbulent energy associated to any MT node, we assume
which entails The turbulent pressure is finally given by3 RESULTS
In this section, we present the results of our model. After a description of the model calibration, useful to fix the model free parameters, we give our prediction for the bubble evolution along a MT extending down to z = 0 along with the IGM turbulence properties.
3.1 Model calibration
Our galaxy evolution model includes several relatively poorly known (albeit important) physical processes that need to be empirically calibrated. To this aim, we have used the observed properties of the MW as a benchmark to fix the best values for the two model free parameters: ε_{*}, the efficiency of star formation (equation 2) and ε_{w}, the SN feedback efficiency (equation 3). This has been accomplished by producing the MT of a DM halo with mass consistent with the MW one (10^{12} M_{⊙}, Xue et al. 2008) at redshift z = 0. We have then compared the resulting properties of the synthetic galaxy with the following ones deduced from MW observations:

Stellar mass. Contributions to the stellar mass come from the disc , the bulge and the halo components (Dehnen & Binney 1998; Brown, Velázquez & Aguilar 2005), yielding a total stellar mass of M_{*}∼ 6 × 10^{10} M_{⊙}.

Gastostellar mass ratio. M_{gas}/M_{*} = 0.13. The mass of the gas has been derived using the observed mass of H i and H ii regions of the Galaxy, (Stahler & Palla 2005).

Current SFR. Murray & Rahman (2010) use the total free–free emission in the WMAP foreground map as a probe of the massive star population and derive a global SFR of 1.3 M_{⊙} yr^{−1}.
For all these observables, we assume a relative uncertainty of 20 per cent. Fig. 1 shows the χ^{2} confidence levels of our model with respect to these observables as a function of ε_{*} and ε_{w}. The distribution presents a clear minimum at (ε_{*}, ε_{w}) = (0.5, 0.012); however, a degeneracy exists with a relatively wide range of choices that can reproduce with sufficient accuracy the global MW properties. Thus, the bestfitting model implies a star formation timescale which is only a factor of ≈2 longer than the freefall time and a relatively inefficient feedback in order not to expel from very first protogalaxies too much gas which will then be crucial to fuel the subsequent star formation in larger objects.
From the same run, we can obtain the redshift evolution of the SFR density along the MW MT for our bestfitting model, after averaging over 100 different realizations of the MT.4 For the sake of comparison only, in Fig. 2 we overplot the cosmic SFR on top of the MW SFR history. Perhaps not coincidentally, the two evolutions match each other quite well, indicating that the Galaxy is a nottoobiased tracer of average cosmic conditions. Finally, in the same figure, we show the SFR density evolution for a model with ε_{w} = 0; you can see that in this case the model is far from the data, which means that, even if the efficiency happens to be small, the feedback by galactic winds is crucial to obtain a consistent model of galaxy evolution.
In the remainder of this paper, we assume that these values for (ε_{*}, ε_{w}) hold for any galaxy present in the MT. This might be a poor approximation, but given the persisting ignorance on star formation and feedback processes, a better performance of more complex choices is not guaranteed to produce a more solid result. In any case, the model success in reproducing key Galactic properties provides at least the first and necessary step towards more refined treatments.
3.2 Galactic outflows
We illustrate the properties of galactic outflows by focusing on the case of the progenitors of a 10^{13}M_{⊙} halo at redshift z = 0 (i.e. a small cluster/group). Since the total number of galaxies at any intermediate redshift is very large, we will deal with mean quantities to describe the global properties of the winds. This precludes the inspection of onetoone relationships between the galaxy and its wind properties, but will make possible to appreciate the correlations among several global quantities, as we will see in the following.
Fig. 3 gives a concise overview of the bubbles we have identified around galaxies at selected redshifts (z = 10, 6, 3, 1) and it can be used to elucidate many physical aspects. Bubbles tend to become older towards lower redshifts with a decreasing age spread; their ages tend to accumulate close to the Hubble time. This behaviour reflects the fact that after an early evolutionary phase in which bubbles grow by number around relatively low mass objects, the average growth of the halo population mass associated with deeper potential wells prevents the formation of new bubbles after z = 3. Hence, most of the bubbles seen around galaxies at intermediate redshifts (z = 3–5), where they are more easily observed, contain material expelled by galaxies at a much earlier time. At the same time, old bubbles merge together to form a fewer large ones. The role of merging events is made clear by the solid curves in the top panels which show the distribution of merging events for haloes that have witnessed at least one merging event. A fraction M/N, indicated by the labels inside the panels in Fig. 3, has evolved in isolation until the considered z, that is, they never suffered a merger. Such ratio decreases steadily with time, going from 81 per cent at z = 10 to 0 at z = 1. The merging activity increases the average ages of bubbles: as discussed above, at high z a large number of young, isolated bubbles exist that are later turned into aged, large ones.
The typical bubble size grows from 5 kpc (physical) at z = 10 to about 100 kpc at z = 1. Merging is definitely an important agent behind this evolution; however, it cannot be the only growth factor for bubbles. This can be realized by combining the information from the age and size distributions. At the lowest redshifts shown in Fig. 3, the bubble age distribution essentially collapses on to the cosmic age value; nevertheless, a considerable size spread, ranging from 1 to a few hundred kpc, exists. Such spread is caused by the corresponding dispersion of galaxy SFRs. Larger galaxies tend to have faster winds, which show up as a tail in the probability distribution function (PDF) of the shock velocity distributions (V_{s} > 30 km s^{−1}). Hence, a minority of the bubbles found at z < 3 are constituted by old bubbles which are still powered by their current galaxy, which is undergoing vigorous star formation. Instead of creating a new bubble, these massive cluster progenitors blow their winds into preexisting, old bubble structures.
From the velocity PDF evolution, it is interesting to note that at early times bubbles (somewhat contrary to naive expectations) expand on average at larger velocities and they tend to decelerate at lower redshifts. This is an obvious implication of the ageing trend described above. As larger galaxies dilute their energy on the larger volume of preexisting cavities, they can sustain the growth of bubbles but cannot reaccelerate them to the high velocities at which they were travelling in the past.
The evolution of the bubble temperature T_{b} (Fig. 3) is governed by competing effects. Adiabatic expansion and radiative cooling (depending on the bubble gas temperature and, to a much limited extent, gas metallicity) lead to a decrease in T_{b}, while the thermalization of the wind energy provides a heat input. Almost all the bubbles have temperatures >10^{5.5} K at z = 10, whereas by z = 3 this fraction decreases to <20 per cent. This is mainly due to the adiabatic expansion of the bubbles.
The expanding outflow is powered by the continuous injection of SN energy both in kinetic (in practice, stored in the shell, where most of the mass is located) and in thermal (stored in the bubble interior) forms. More quantitatively, from Fig. 4, one can appreciate that the kinetic energy is always dominant (by a factor of 2–3) with respect to the internal one. This kinetic energy will be rapidly converted in disordered motions by instabilities, finally resulting in a fully developed turbulent spectrum. The bubble volume filling factor also increases with time, as illustrated in the bottom panel of the same figure, reaching about 10 per cent at z < 2. One has to keep in mind that our study is based on the MT of a cluster/small group and therefore it might reflect a somewhat biased cosmic region, although we do not expect major differences as we extrapolate these results to general IGM. A similar treatment of the outflow evolution has been successfully used by Bertone et al. (2005) and Samui et al. (2008) to investigate the role of galaxies in enriching the IGM. In particular, Samui et al. (2008) used the modified PS formalism to follow the evolution of haloes and porosityweighted averages to investigate the global influence of winds. For a reasonable choice of the model free parameters, they find that outflows can generally escape from the lowmass haloes (M_{h} < 10^{9} M_{⊙}) which then dominate the observed enrichment; nevertheless, MWtype galaxies can create Mpcsize metalenriched bubbles in their surroundings. Bertone et al. (2005) apply the semianalytical prescription for galactic winds to highresolution Nbody simulations of field galaxies. Their main conclusion is that galactic outflows do not perturb the structure of the Lyα forest; mainly, galaxies with 10^{9} < M_{*} < 10^{10} M_{⊙} are responsible for IGM chemical enrichment at z = 3.
The bubble volume filling factor is a useful tool when applying our models to quantify turbulent effects in the quasar absorption line spectra, as we will try to do later on. As a note of caution, the volume filling factor shown in Fig. 4 might be overestimated to some degree as we do not attempt to account for possible, although likely rare, overlaps between bubbles.
3.3 Turbulence evolution
The best way to observationally probe the turbulent content of gas affected by outflows is to measure the Doppler parameter, b, of absorption lines arising from kinetically perturbed gas. For this reason, in the following, we will use b to quantify turbulence evolution.
In general, one can write b as a quadratic sum of the thermal and turbulent components:
where is the thermal contribution and is the contribution due to the turbulent pressure defined in equation (22). Note that the turbulent contribution (contrary to the thermal one) does not depend on the atomic mass m_{a} of the absorbing element. This allows to generalize our results to all species associated with the cold H i shell gas (e.g. C iv). Moreover, since the observed Doppler parameter is the sum of the two contributions and the thermal broadening is smaller for heavier elements, the relative contribution of turbulent broadening should be more important (and hence more easily detectable) for heavier elements.In Fig. 5, we show the distribution of Doppler parameters in simulated galaxies at different redshifts. The b_{t} distribution function at all redshifts shows a similar shape, with a lowenergy cutoff, followed by a maximum and a steady declining tail towards larger values. However, the PDF shifts towards lower b_{t} values with time: while at z ≥ 6, the region affected by the winds can become very turbulent, reaching in a few spots b_{t} > 50 km s^{−1}, as a result of the violent expansion of fresh new bubbles carved by early galactic winds, at later times, winds blowing in preexisting bubbles become more gentle, resulting in decrease in turbulent broadening effects. By z ≈ 1–2 the distribution has stabilized on median values of the order of 1–2 km s^{−1} and even the mostturbulent spots do not exceed 25 km s^{−1}. The trend above is further enhanced by the dissipation of turbulent energy with time, which does not allow to store kinetic energy into eddies indefinitely.
To see this quantitatively, it is useful to estimate the turbulent dissipation timescale. Under the hypothesis of a fully developed, steadystate spectrum, the dissipation timescale is well approximated by the eddy turnover timescale, namely t_{e}∼ R_{s}/V_{s}. If we further assume that to the early stages of bubble evolution the Sedov solution applies, it follows that t_{e}∼ 2.5t_{b}, where t_{b} is the bubble age, is in turn shorter than the Hubble time. This implies that due to their high velocity and relatively small radii, dissipation acts very efficiently at earlier stages, thus decreasing the turbulent energy budget of these objects.
We pause briefly to emphasize an important point. One has to keep in mind though (see Fig. 4) that the filling factor of the turbulent regions is quite small at high redshift and therefore it must be made clear that the above is not the b_{t} parameter applicable to the global IGM, but only to tiny and specific (but welldefined) portions of it. By plotting the volumeweighted evolution of 〈b_{t}〉, shown in Fig. 6, we can find support to the previous statement: at higher redshift, the turbulent pressure contribution is less important because of the small size of outflowaffected regions, whereas in local objects, the bubbles grow enough to occupy a significant fraction of the volume; however, this occurs too late when dissipation has already started to quench eddies on the smallest scales. The tradeoff between injection and dissipation conspires to produce a maximum in 〈b_{t}〉 at z ≈ 1.
Let us turn back to the analysis of the other panels in Fig. 5. It is interesting to determine the level of turbulence in bubbles produced by galaxies of different stellar masses. In general, b_{t} is positively correlated with the stellar mass of a galaxy. This is a consequence of the fact that the SFR (which is ultimately proportional to the energy input rate due to SNe) is more sustained in larger galaxies. However, we see that the onset of a turbulent regime is almost unavoidable in regions that are stirred by galactic winds. Again, the peak at very low M_{*} is interpreted as the imprint of young bubbles in which turbulence has not yet suffered from dissipation: we doublechecked this hypothesis by studying the correlation in different bubble age bins. This analysis is reported in Fig. 7 (for a single MT realization) whose examination clarifies that the youngest bubbles have the largest b_{t} values. Thus, the bursting star formation activity going on in small galaxies momentarily reverses the trend of increasing b_{t} as a function of the stellar mass, clearly indicating that the level of turbulence depends both on the age and on the mass of the object. The spread of the distribution is very limited, suggesting that different mass accretion histories do not influence substantially the final Doppler parameter distribution.
A relatively easily accessible observational quantity when studying the IGM is the neutral hydrogen column density, , of the absorbing shell/system. A reasonable estimate for such quantity can be obtained from the following formula:
In the latter part of the above equation, n_{H} is the mean density in the thin shell, given by the ratio between its mass (M_{s}) and its volume (V_{s}): The neutral hydrogen ionization correction can be obtained by solving the appropriate photoionization equilibrium balance equations given, for example, by Theuns et al. (1998) and assuming the UV background history of the Universe computed by Haardt & Madau (1996). The third row of panels in Fig. 5 show the distribution of b_{t} with . One can note that at redshift z ∼ 3 the turbulent broadening starts making an impact in H i lines with cm^{−2}.Finally, the lowermost panel in Fig. 5 shows the correlation between b_{t} and the distance reached by the outflow. This might have important observational implications. Adelberger et al. (2003) by studying the correlation between galaxies and C iv systems with large velocity spreads showed that most of these systems are associated with starforming galaxies. It is conceivable to interpret such results as a direct evidence of largescale outflows around 2 ≲ z ≲ 3 galaxies. A similar investigation, as discussed in the next section, could directly probe the level of turbulence in these objects and test our predictions.
4 TESTING THE MODEL AGAINST DATA
A direct comparison of our results with the actual level of turbulence in the IGM is not possible since this observable has not been convincingly derived yet. In the meantime, it is possible to perform at least a sanity check for our model as follows.
In order to study the properties of the IGM, a number of highquality QSO (quasar) absorptionline experiments have been carried out during the last decade allowing the buildup of a large sample of observed systems for which both and the Doppler parameter b have been fitted (see an almost complete list of references in Meiksin 2009). The observed H i column densities cover the range 10^{12}–10^{22} cm^{−2} and follow an almost perfect powerlaw distribution with β = 1.5–1.7. The Doppler parameters of these systems fall approximately in the range 10 < b < 100 km s^{−1}, with the vast majority being concentrated between 15–60 km s^{−1}. The observed broadening could be either real, that is, caused by physical mechanisms (e.g. thermal broadening, turbulence and peculiar motions), or artificial, due to the blending of close systems.
Many authors have suggested the use of the lowercutoff in the Doppler parameter distribution versus as a reliable proxy for the gas temperature at a given redshift. Similar approaches have been used by Schaye et al. (2000b) to report evolution in the inferred IGM gas temperature over the range 2.0 < z < 4.5 and by Ricotti, Gnedin & Shull (2000) to extract the equation of state of the IGM from the Doppler parameter distribution of the lowdensity Lyα forest. To make progress, we follow the latter procedure.
In Fig. 8, we show the absorber properties derived in a sample of Lyα forests from 22 highresolution VLT/UVES Large Program QSO spectra (D’Odorico et al. 2008; Saitta et al. 2008), selected at redshift z ∼ 3. To avoid the contamination of absorbers affected by the UV background flux from the QSO, we exclude from our sample absorbers closer than 3000 km s^{−1} to the QSO and we also do not consider in our analysis Doppler parameters >100 km s^{−1} since they are very likely due to errors in the fitting procedure. In the same graph, we report our estimated thermal Doppler parameter (b_{th}) as a function of , obtained from a linear fit of the minimumb values, defined by dividing the range in 50 bins and selecting the smaller b value in each bin.
The sanity check then consists in the comparison of our predicted mean b_{t} distribution with the analogous average residual Doppler parameters from the data after the thermal contribution has been subtracted. Given the reasons for contamination listed above, what we show in Fig. 8 cannot be intended as a direct theory–data comparison as the observed broadening still contains a contribution from IGM peculiar motions, which are not directly related to turbulent motions. Hence, the data points must be considered as strict upper limits which our theoretical Doppler parameters should not exceed. Quite encouragingly, this is indeed the case: our model is consistent with the experimental upper limits in the entire H i column density range. The uncertainties in the exact observational determination of b_{t} do not allow for a more meaningful comparison at this stage. This will become possible in the future if Rauchtype experiments using (lensed) quasar pairs will be carried out. The increased availability of such systems, enabled by the SDSS, has been already shown by Kirkman & Tytler (2008) who have discovered 130 close pairs of QSOs in the SDSS DR5.
As an alternative, the presence of the metal lines in the absorbers could also be used as a possible way to obtain an independent measurement of the IGM turbulent level. The measurements of the Doppler parameters from absorption features corresponding to two or more elements (typically, C iv and S iv) assumed to be cospatial allow the separation between thermal and kinematic broadening contributions. As already mentioned in the previous section, this technique could disentangle the broadening due to the kinematic from the thermal one which is instead massdependent. To achieve this result, highresolution spectroscopy is required in order to properly fit the profile of the absorbers along the QSO line of sight and to determine the redshift of the lines with sufficient precision to make sure that the ions share the same absorber environment. We plan to combine these types of observations with our model in a forthcoming work.
5 SUMMARY AND IMPLICATIONS
Turbulence is an important physical process in essentially all astrophysical environments, from planet formation and atmospheres to (proto)stars and galaxies, but particularly on the scale relevant to cosmological phenomena, it has received a relatively meager attention. This is partly understood from the difficulty to collect reliable data as observations aimed at the study of these aspects are very challenging. Even more surprisingly though, the wellascertained fact that the IGM is polluted on large scales by heavy elements, carried by powerful winds powered by galaxies, should raise the strong suspicion that turbulence must pervade the IGM, albeit at levels yet difficult to understand.
This paper aims at making one of the first steps to clarify the amount of turbulent energy deposited in the IGM and then dissipated through a cascade to smaller scales by galactic winds and to characterize its properties. To this aim, we have developed a simple but physically accurate model based on a standard ΛCDM hierarchical structure formation model to identify haloes hosting galaxy progenitors of a small cluster/group of the total mass M = 10^{13} M_{⊙} at z = 0. In addition, the required star formation history in all their progenitors via a MT technique (carefully calibrated to reproduce the MW) is also followed. As these galaxies progress along their history, their collective SN explosion vents gas, heavy elements and, most importantly for this study, energy into the surrounding IGM. The evolution of these hot intergalactic bubbles is followed by including all the relevant physics, enabling us to determine their thermodynamic properties at best. We then couple this general galaxy–IGM interplay scenario (which can be used and adapted to a number of ancillary applications as e.g. metal enrichment and the Sunyaev'Zel’dovich effect) to a treatment of turbulence evolution using a powerful approach based on the spectral transfer function developed by Norman & Ferrara (1996). In brief, the major advantage of this method is a relatively straightforward derivation of the turbulent IGM turbulent pressure evolution (see equation 22).
The main results can be summarized as follows. At z ≈ 3, the majority of the bubbles around galaxies are old (ages >1 Gyr), that is, they contain metals expelled by their progenitors at earlier times and have a size distribution at that redshift in the range 10–100 (physical) kpc. The velocities reach up to 100 km s^{−1} for larger galaxies but a considerable fraction of the bubbles have already become subsonic and the shock decayed into a sound wave; temperatures in the rarefied bubble cavities are in the range 5.4 < log (T/K) < 6.5. The bubble volume filling factor increases with time, reaching about 10 per cent at z < 2. The energy deposited by these expanding shocks in the IGM is predominantly in kinetic form (mean energy density of 1 μeV cm^{−3}, about two to three times the thermal one), which is rapidly converted in disordered motions by instabilities, finally resulting in a fully developed turbulent spectrum. The corresponding mean turbulent Doppler parameter, b_{t}, peaks at z ≈ 1 at about 1.5 km s^{−1} with maximum b_{t} = 25 km s^{−1}. More informative though is the b_{t} distribution associated with individual bubbles and the relations with the galaxy stellar mass, bubble shell column density and radius (shown in Fig. 5). The shape of the b_{t} distribution does not significantly evolve with redshift, but undergoes a continuous shift towards lower b_{t} values with time, as a result of bubble ageing. We also find a clear trend of decreasing b_{t} with and a more complex dependence on R_{s} resulting from the age spread of the bubbles. We have finally tried to compare our results with available data on the Doppler parameter deduced from QSO absorptionline studies of the Lyα forest, but we are facing the difficulty of accurately removing the contamination due to the thermal component and the more troublesome peculiar contribution. Alternatively, one can extend our study to make predictions for quasar pair systems and/or using observations involving different metal species. We are already working along this direction, which in general requires additional information of the positions of galaxies and therefore the implementation of our turbulent model into a fullscale hydrodynamic simulation.
If turbulence is indeed present in the IGM, as our study seems to suggest, it might have relevant implications for a number of IGM studies. First of all, the longstanding problem of a discrepancy between the observed Doppler parameter and those deduced from numerical simulations of the Lyα forest might be alleviated by a proper inclusion of turbulence broadening. Oppenheimer & Davé (2009) convincingly demonstrated that only a model where subresolution turbulence is explicitly added, and the required turbulence increases with the O vi absorber strength, is able to match all the observations, giving a hint of the correlation between turbulence and enrichment that we are planning to explore in a forthcoming paper. Secondly, the turbulent energy injection might affect the matter power spectrum determination from experiments using the Lyα forest as a probe, particularly on the smallest scales (McDonald et al. 2006). Thirdly, if turbulence is largely present in the IGM, then it might result in a detectable signature in the scintillation of distant quasars (Ferrara & Perna 2001; Lazio et al. 2004). The IGM hosts the necessary conditions for scattering, is largely ionized and is permeated by shocks from largescale formation or galactic feedback. Moreover, the long pathlengths through the IGM can compensate for the lower density with respect to the ISM. Lazio et al. (2008) have considered intergalactic scattering, expanding on previous treatments from intracluster media of Hall & Sciama (1979) to scintillating active galactic nuclei; their findings are broadly consistent with the scenario of a highly turbulent IGM. This possibility was explored so far only as an Ansatz, to which our study now provides a more solid support. Fourthly, turbulence might be carefully modelled in order to interpret the results for experiments aiming at directly measuring the change in the expansion rate of the Universe with time and the variability of fundamental constants (Corasaniti, Huterer & Melchiorri 2007; Liske et al. 2008), as, for example, the proposed CODEXESPRESSO experiment (Cristiani et al. 2007). Last but not the least, the relatively mild level of turbulence we find at z ≈ 3 is consistent so far with the upper limits coming from observations. Moreover, it indicates that the IGM at those epochs had already dissipated a large fraction of the kinetic energy deposited at early times (z > 5) by the galaxies which initiated the metal enrichment and, possibly, the reionization process. All these well agree with the idea that the metals we detect at moderate redshifts were mostly produced during a ‘preenrichment’ phase of the IGM by a population of small galaxy ancestors of the ones we routinely study in detail, the Lyman break galaxies.
We thank V. D’Odorico for useful suggestions and for supplying the Lyα forest data we used in our analysis. Stimulating conversations with P. Petitjean, J. Niemeyer and R. Valdarnini are warmly acknowledged. CE would like to thank Scuola Normale Superiore in Pisa for the warm hospitality during the preparation for this work.