text:
# Turbulent dissipation in the ISM: the coexistence of forced and decaying regimes and implications for galaxy formation and evolution

## 1 Introduction

The large-scale star formation (SF) cycle in disk galaxies plays a key role in modeling galaxy formation and evolution. This cycle crucially depends on the dissipative properties of the turbulent interstellar medium (ISM). In normal, noninteracting disk galaxies, SF is believed to be statistically stationary, and self-regulated111Throughout the paper, by self-regulation we only mean that the SF rate feeds back on itself, generally in a negative fashion, maintaining a statistically stationary value over long time scales, but of course allowing for significant spatial and temporal fluctuations. by a balance between the energy injection and dissipation in the turbulent ISM in the disk (e.g., Scalo & Struck-Marcell 1984; Vázquez & Scalo 1989; Dopita 1990; Firmani & Tutukov 1992, 1994; Dopita & Ryder 1994; Wang & Silk 1994; Firmani, Hernández, & Gallagher 1996; Avila-Reese 1998; Struck & Smith 1999; Avila-Reese & Firmani 2000). A self-regulating SF mechanism has also been used in semi-analytical models of galaxy formation (e.g., White & Frenk 1991; Kauffmann, White, & Guiderdoni 1993; Cole et al. 1994; Somerville & Primack 1999; van den Bosch 1999), but in this case it has been applied to a hypotethic intrahalo medium in hydrostatic equilibrium with its own gravitational potential plus that of the cosmological dark matter halo, which extends up to approximately 15-20 times the optical radius of the galaxy. In these models the disk ISM is virtually ignored and the feedback from stars is assumed to efficiently reheat and drive back the cooled gas into the intrahalo medium, the reheated gas fraction being a strong function of the halo mass. Unfortunately, in most galaxy formation and evolution models, the detailed thermo-hydrodynamics of the ISM has not been treated explicitly. Therefore, several assumptions and approximations have had to be made about the SF feedback and the dissipative properties of the ISM.

Stars —mainly massive OB stars and the SN explosions they produce— are sources of thermal and turbulent kinetic energy in the ISM. Results from thermo-hydrodynamical simulations in different contexts have shown that the kinetic energy $`E_\mathrm{k}`$ injected by SNe and massive stars to the ISM deeply affects the dynamics of the gas in disk galaxies causing it to be highly turbulent (e.g., Bania & Lyon 1980; Chiang & Prendergast 1985; Chiang & Bregman 1988; Navarro & White 1993; Mihos & Hernquist 1994,1996; Rosen & Bregman 1995; Vázquez-Semadeni, Passot & Pouquet 1995; Friedly & Benz 1995; Passot, Vázquez-Semadeni & Pouquet 1995; Gerritsen 1997; Avillez 1999; Korpi et al. 1999; see also the reviews by Scalo 1987 and by Vázquez-Semadeni et al. 2000). Cloud formation and disruption, the structure and dynamics of the vertical gaseous disk (in particular its scale height), fountains and chimneys, and even a huge hot gas corona (the “intra-halo medium”) in virial equilibrium (including the potential of a huge cosmological dark matter halo), might be some of the phenomena and processes related to the turbulence produced by stellar kinetic energy sources. A key ingredient in all these phenomena and processes is the kinetic energy dissipation rate in the turbulent ISM, since it determines the effectiveness of large-scale SF feedback on galaxy formation and evolution.

Dissipation in turbulent incompressible fluids is mainly controlled by a turbulent cascade from large to small scales, since it is in the latter where kinetic energy is dissipated. Several analytical estimates and numerical simulations have shown that incompressible MHD turbulence decays as $`t^\alpha `$ with $`0.7\genfrac{}{}{0pt}{}{_<}{^{}}\alpha \genfrac{}{}{0pt}{}{_<}{^{}}1.0`$ (Biskamp 1994; Hossain et al. 1995; Galtier, Politano, & Pouquet 1997). In the case of the large-scale ISM, turbulence involves supersonic MHD compressible non-isothermal flows, which most likely are dominated by shocks. High-resolution 3D simulations have been used recently to investigate dissipation in compressible MHD isothermal flows with parameters corresponding to Galactic molecular clouds (Mac Low et al. 1998; Stone, Ostriker & Gammie 1998; Padoan & Nordlund 1999; Mac Low 1999a). For the decaying regime, Mac Low et al. (1998) and Stone et al. (1998) have found $`E_\mathrm{k}`$ to decay as $`t^\alpha `$, with $`\alpha 0.81.0`$. For the forced regime, Stone et al. (1998) and Mac Low (1999a) concluded that the characteristic turbulent dissipation time $`t_d`$ is of the order of or smaller than the crossing time for the driving (or forcing) length at the rms turbulent velocity $`\upsilon _{\mathrm{rms}}`$, even in the presence of strong magnetic fields. In those simulations, turbulence has been driven in Fourier space, with a fixed kinetic energy injection rate, $`\dot{E}_k^{\mathrm{in}}`$, generating random large-scale velocity fluctuations with a constant $`u_{\mathrm{rms}}`$. As a result, kinetic energy is injected ubiquitously (i.e., everywhere in space), and typically at large spatial scales, although Mac Low (1999a) also considered intermediate injection scales.

The way in which $`E_\mathrm{k}`$ is injected into the ISM certainly differs from the prescription used in the studies mentioned above. The scales at which kinetic energy injection occurs (the regions directly heated or accelerated by the stellar activity, of sizes of several pc to a few hundreds of pc), are small compared to the scales of interest in the large-scale ISM (up to a few kiloparsecs).222Strictly speaking, since such forcing is confined to a limited range of (small) spatial scales, in Fourier space it contains components of all spatial frequencies. In reference to the stellar driving we consider, in this paper the “scale of the forcing” will always refer to its characteristic size in physical space. Besides, the sources are of short duration and located at discrete, generally isolated sites. Moreover, the global energy input rate is not expected to be constant in general. This situation implies that in the ISM both forced- and decaying-turbulence regimes should coexist, the latter being the regions located far from the energy injection sources.

In this paper we use two-dimensional (2D) numerical MHD simulations of turbulent compressible fluids resembling the ISM of normal disk galaxies (Vázquez-Semadeni et al. 1995, 1996; Passot et al. 1995) to explore the ability of the ISM to dissipate the turbulent kinetic energy injected by stellar sources. In order to do this, we use an ISM-like turbulent kinetic energy injection mechanism, based on “stellar winds” applied randomly in the medium; thus, we study how the dissipative properties of the flow depend on the parameters of the injection sources.

It should also be noted that the ISM at large scales is far from isothermal (see, e.g., Myers 1978), to the extent that this non-isothermality has constituted the starting point for multi-phase models of the ISM (e.g., Field, Goldsmith & Habing 1969; Cox & Smith 1974; McKee & Ostriker 1977). Therefore, in this paper we consider turbulent flows in the presence of parameterized radiative cooling and background heating (mimicking the heating from a background UV field and cosmic rays).

In §2 we describe the numerical method, the implementation of the kinetic energy injection mechanism, and the parameters and initial conditions of the simulations. In §3 we discuss the dissipation time scales for driven (forced) turbulence simulations using various energy injection regimes, finding that most of the injected energy is dissipated locally; we also present a simple model to calculate $`t_\mathrm{d}`$ as a function of observable quantities of the ISM. In §4 we explore the decay and propagation of the “residual” turbulent energy that “survives” the strong local dissipation, by considering the simulations after turning off the stellar energy input sources. In §5 we discuss possible caveats of our simulations, such as the dimensionality and the low resolution, and present higher resolution runs which support our conclusions. Then §6.1 compares our results with previous work, and §§6.2 and 6.3 discuss the implications of our results for models of galaxy evolution and formation, respectively. Finally, in §7 we present a summary and our conclusions.

## 2 Numerical method and simulations

### 2.1 The model

We use a slight variation of the numerical model presented in a series of previous papers (Vázquez-Semadeni et al. 1995, 1996; Passot et al. 1995) and the reader is referred to those papers for details; here we just sketch the method. The MHD equations, including the internal energy conservation equation, are solved in two-dimensions (2D) in the presence of model terms for radiative cooling, background heating, and stellar winds. While in previous papers self-gravity, the Coriolis force and large-scale shear were included, here we do not consider them because of various reasons. First, self-gravity causes runaway gravitational collapse when the stellar energy injection is not applied at the densest regions, as we do here (see §2.2 below). Second, the Coriolis force modifies the dissipation rate at long times in decaying regimes, but we consider this to be an unrealistic situation (see §4). Third, the large-scale shear introduces a non-zero lower bound to the kinetic energy, which we feel is irrelevant for the purposes of the present paper.

The numerical technique is pseudo-spectral with periodic boundary conditions, and we use a combination of regular second-order viscosity of the form $`\mu [^2𝐮+(1/3)(𝐮)]`$, and hyperviscosity of the form $`\nu ^8𝐮`$ in the momentum equation. The hyperviscosity operator is a very steep function of wavenumber $`k`$, and, upon adequate choice of the coefficient $`\nu `$, confines viscous effects to a much smaller range of scales than the regular second-order (or lower-than-eigth-order in general) viscosity, allowing a more extended inertial (self-similar) range in the kinetic energy spectrum. Our simulations show clear power-law inertial ranges over more than one decade in wavenumbers (§6.1). It should be pointed out that finite-difference schemes necessarily introduce numerical dissipation which is typically second-or third-order. On the other hand, hyperviscosity introduces spurious oscillations in the vicinity of strong shocks (Passot & Pouquet 1988). In order to “filter” these out we add a small amount of second-order viscosity, which, however, can be kept much smaller than what would be necessary without the hyperviscosity. Note that the use of hyperviscosity should not affect the global dissipation rate of the flow, as this rate adjusts itself to counterbalance the energy injection by the stars, independently of the specific form of the dissipation operator. What changes from one operator to another at a given dissipation rate is the dissipation scale (e.g., Frisch 1995) and the saturated energy content (Stone et al. 1998).

Additionally, the code uses second-order mass diffusion (with coefficient $`\mu _\rho `$) with the purpose of smoothing out shocks even further, as the spectral scheme we use cannot handle exceedingly steep gradients. We do not expect the mass diffusion to affect the dissipation rates in the code since its only purpose is to prevent the formation of excessively small density structures. The dissipative coefficients are chosen in such a way that the dissipation operators are already strongly active a safe margin before the shocks become too steep to be resolved.

In the $`128^2`$ simulations we use $`\nu =10^{11}`$, $`\mu =0.003`$ and $`\mu _\rho =0.06`$. We take zero electrical resistivity, and a thermal diffusivity $`\eta =\gamma [\mu +\nu (n/4)^6]`$, where $`\gamma `$ is the ideal-gas ratio of specific heats, and $`n`$ is the resolution. All coefficients are expressed in the code’s non-dimensional units, and are, in the momentum and magnetic flux equations, essentially inverse Reynolds numbers.

We have chosen physical units in such a way that the simulations resemble the ISM in the plane of the Galaxy near the solar neighborhood at the 1 kpc scale. We adopt as normalization values a number density $`n_0=1`$ cm-3, a temperature $`T_0=10^4`$ K, a velocity $`u_0=11.7`$ km s-1, and a length $`L_0=1`$ kpc. The velocity and length units imply a time which (divided by $`2\pi `$ since in non-dimensional units, the length of the box is $`2\pi `$) gives the code time unit, $`t_0=(L_0/u_0)/2\pi =1.36\times 10^7`$ yr. The magnetic field is written as $`𝐁=B_0e_x+𝐛`$, where $`B_0=1.5\mu `$G represents the uniform azimuthal component of the field and $`𝐛`$ is a superposed fluctuating component of rms amplitude $`=5\mu `$G.

### 2.2 “Star” formation and kinetic energy injection

In previous papers, an “OB star” (or, in general, a source) was turned on at grid point $`𝐱`$ whenever the density exceeded a certain (arbitrary) threshold value and the velocity field had finite negative divergence. A “star” consisted of a local source of heat, which caused the nearby gas to expand, forming warm bubbles resembling HII regions. After some experimentation, in this paper we have chosen to use “winds” (local outward radial accelerations, see Vázquez-Semadeni et al. 1996) rather than heating, and to place the stars randomly with a given probability, rather than at the density peaks.

The reasons for our choice are that the heating scheme deposited unspecified amounts of kinetic energy into the flow, the exact amount depending on the local cooling rate, which in turn depends on the local density and temperature, rendering the measurement of the kinetic energy injection rate impossible. Moreover, the density-directed placement of the sources caused three phenomena which obstruct the study of how the results depend on the the initial parameters of the energy injection mechanism. First, the OB star formation rate (SFR) fluctuated strongly, with alternating periods of dormancy and of strong activity (see Vázquez-Semadeni et al. 1995). This complicates the determination of the characteristic energy injection time, $`E_\mathrm{k}/\dot{E}_\mathrm{k}^{\mathrm{in}}`$, where $`\dot{E}_\mathrm{k}^{\mathrm{in}}`$ is the global energy injection rate, when the latter goes to zero. Second, the density-directed SFR produces strong clustering of the sources, because generally more than one grid point satisfy the density criterion within a density peak (a “cloud”). This causes their outputs to interfere destructively (since oppositely directed accelerations acting on the same grid point cancel out) in such a way that the total energy input of a cluster is smaller than the sum of the inputs from its constituent stars taken separately. Besides, due to clustering, large “supershells” form; these supershells are energy input sources which cannot be easily related to the initial injection parameters. In fact, destructive interference between the sources and “supershell” formation still occur in the random-star placement scheme at the largest SFRs, but to a much lesser extent. Finally, the density-directed scheme involves the presence of an initial negative velocity divergence at the forcing sites, making it impossible to determine the final velocity attained by the flow around the sources, since the acceleration acts on a pre-existing velocity field of fluctuating magnitude. In summary, the random placement of stars we use here, although certainly less realistic than the density-directed placement, makes for a much more controlled energy input, and we adopt it in our forced simulations.

We add a fixed number of sources at each code timestep by specifying the probability with which a source is placed in each grid point. However, since the timestep is variable in the code, the instantaneous number of sources $`n_\mathrm{s}`$ in the simulation is not perfectly constant in time. This is why we prefer to use the SFR per unit of area, $`\dot{\mathrm{\Sigma }}_{\mathrm{OB}}`$, measured directly from the simulations, in such a way that at any one time $`n_\mathrm{s}=\dot{\mathrm{\Sigma }}_{\mathrm{OB}}\mathrm{\Delta }t_{\mathrm{OB}}`$, where $`\mathrm{\Delta }t_{\mathrm{OB}}=6.8`$ Myr is a typical lifetime of OB stars. Once a source has been turned on at a given grid point x, it stays on for a time $`\mathrm{\Delta }t_{\mathrm{OB}}`$. During this time the gas around the source receives an acceleration $`𝐚`$ directed radially away from x, with a magnitude

$$a(r)=\frac{a_{\mathrm{max}}r}{\sigma }\mathrm{exp}\left[(r^2\sigma ^2)/2\sigma ^2\right],$$
(1)
where $`a_{\mathrm{max}}`$ is the amplitude of the acceleration, $`r`$ is the distance from the position of the star, and $`\sigma `$ is a free parameter of the order of a few pixels. At a resolution of 128 grid points per dimension, 1 pixel $`7.8`$ pc. The acceleration $`𝐚`$ produces an evolving velocity profile $`v(r,t)`$ around the source. We show an example of the acceleration and velocity profiles at the end of a stellar lifetime in fig. 1 for the fiducial simulation called run 30 (see Table 1). It is easily seen that $`a(r)`$ has its maximum at $`r=\sigma `$, and that by $`r=2.2\sigma `$, it has decreased to $`0.2a_{\mathrm{max}}`$. We choose this as a characteristic radius of the forcing region, which we denote $`l_\mathrm{f}`$. It is important to note that $`l_\mathrm{f}`$ is not the final radius of the expanding shells that result, which are actually a response of the flow to the energy injection mechanism, and expand to sizes quite larger than $`l_\mathrm{f}`$. This is particularly noticeable when stars cluster, in which case the shells may reach sizes of a few hundred parsecs, becoming “supershells”.

The characteristic power $`\dot{e}_\mathrm{k}`$ and energy $`e_\mathrm{k}`$ injected to the flow by one source can be roughly estimated as

$`\dot{e}_\mathrm{k}={\displaystyle _{S_\mathrm{f}}}\mathrm{\Sigma }_\mathrm{g}𝐮𝐚d^2x\mathrm{\Sigma }_\mathrm{g}_\mathrm{f}u_\mathrm{f}a_\mathrm{f}\pi l_\mathrm{f}^2=M_\mathrm{f}u_\mathrm{f}a_\mathrm{f}`$
$`e_\mathrm{k}=\dot{e}_\mathrm{k}\mathrm{\Delta }t_{\mathrm{OB}},`$ (2)
where $`𝐮`$ is the flow velocity vector, $`S_\mathrm{f}=\pi l_\mathrm{f}^2`$ is the area of the forcing region, $`\mathrm{\Sigma }_\mathrm{g}_\mathrm{f}`$ is its initial average surface density, and $`u_\mathrm{f}`$ and $`a_\mathrm{f}`$ are characteristic values of the velocity and acceleration, respectively. Note that, since the simulations are two-dimensional, it is necessary to specify some effective disk height $`H_\mathrm{f}`$ for the forcing regions in order to define their surface density and mass. In the Galaxy, the regions of highest young cluster density oscillate above and below the plane by $`50`$ pc (Alfaro, Cabrera-Cano, & Delgado 1991). Therefore, we take $`H_\mathrm{f}=100`$ pc. Then, $`\mathrm{\Sigma }_\mathrm{g}_\mathrm{f}=\rho _0H_f`$ with $`\rho _0=m_\mathrm{H}n_0`$, and $`M_\mathrm{f}=S_\mathrm{f}\mathrm{\Sigma }_\mathrm{g}_\mathrm{f}`$; we use the simulation mean density $`n_0`$ as representative of the mean density in the forcing regions because the sources are seeded randomly in the simulation.

The product $`a_\mathrm{f}u_\mathrm{f}`$ may be defined as the characteristic power per unit of mass injected per forcing region, $`\dot{ϵ}_\mathrm{k}`$. We define $`u_\mathrm{f}u(\sigma ,t=\mathrm{\Delta }t_{\mathrm{OB}})<u_{\mathrm{max}}(t=\mathrm{\Delta }t_{\mathrm{OB}})`$ and $`a_\mathrm{f}`$ as the value of $`a(r)`$ where the area contained under the acceleration profile is divided horizontally in two equal regions. This roughly happens at $`r=2\sigma `$, where $`a(r)=0.45a_{\mathrm{max}}`$; therefore, we take $`a_\mathrm{f}=0.45a_{\mathrm{max}}`$. These values were decided upon empirically, but we have checked that they apply reasonably well to all of our simulations. Note that $`u_\mathrm{f}`$ is not trivially related to $`a_\mathrm{f}`$, $`\mathrm{\Delta }t_{\mathrm{OB}}`$ and $`l_\mathrm{f}`$, due to the shape of the acceleration profile and to the fixed time over which the acceleration is active. In particular, we have checked that the simple formula $`v^22a_\mathrm{f}l_\mathrm{f}`$ does not apply very accurately.

Together, the three parameters $`a_\mathrm{f}`$, $`l_\mathrm{f}`$ and $`u_\mathrm{f}`$ completely define the local energy injection process.333We wish to stress here that while the former two parameters are controlled by us, the latter is already a response of the flow to the first two, and has to be measured from the simulations, rather than being user-defined. In this sense, our model here is still semi-empirical. However, $`a_\mathrm{f}`$ and $`u_\mathrm{f}`$ have meaning mostly within the context of our simulations only, while $`\dot{e}_\mathrm{k}u_\mathrm{f}a_\mathrm{f}l_\mathrm{f}^2H_\mathrm{f}`$ (or, equivalently, $`e_\mathrm{k}`$, since we take $`\mathrm{\Delta }t_{\mathrm{OB}}=6.8`$ Myr=const in all cases) is a more physically and observationally relevant quantity. Thus, in what follows, we will consider $`e_\mathrm{k}`$ and $`l_\mathrm{f}`$ as the determining parameters of the local properties of the injection process. Additionally, there is a third, global, parameter, namely $`\dot{\mathrm{\Sigma }}_{\mathrm{OB}}`$. In the following sections we explore how the dissipative properties of the turbulent flow depend on the forcing parameters: $`l_\mathrm{f}`$, $`e_\mathrm{k}`$, and $`\dot{\mathrm{\Sigma }}_{\mathrm{OB}}`$.

### 2.3 The simulations

We have performed a sizeable number of simulations, summarized in Table 1. Most of the simulations have been done at a low resolution of 128 grid points per dimension, in order to achieve reasonable coverage of parameter space. The dependence of the dissipation rate on resolution is further discussed in §5, where test runs at higher resolution ($`512^2`$) are presented as well, for which the results obtained from the low-resolution simulations continue to hold.

We adopt the simulation labeled run 30 as the “fiducial” case. The set of simulations was initially chosen in such a way as to have at least two different values of each of the three parameters $`l_\mathrm{f},u_\mathrm{f},`$ and $`n_\mathrm{s}`$, for (roughly) constant values of the others (see Table 1). However, as our study progressed, we realized that the most relevant parameters are not exactly those we considered initially and, as a result, this philosophy is not clearly reflected in the set of derived parameters $`l_\mathrm{f}`$, $`e_\mathrm{k}`$, and $`\dot{\mathrm{\Sigma }}_{\mathrm{OB}}`$. In any case, the three derived input parameters span a reasonable range of values and, as will be described in §3, we have constructed a simple model which combines them, to which our simulations will be seen to conform to within $`\genfrac{}{}{0pt}{}{_<}{^{}}20`$%.

In the fiducial run, the parameters of the kinetic energy injection were chosen to represent a compromise between values typical of expanding HII regions and SN remnants at late evolutionary stages (unfortunately, our code cannot handle very strong shocks). For run 30, the power and the total kinetic energy injected by each source is $`\dot{e}_\mathrm{k}=1.35\times 10^{35}`$ erg s-1 and $`e_\mathrm{k}=2.84\times 10^{49}`$ erg, respectively (see Table 1). The estimated radius of a typical HII region before the SN explosion is $`5060`$ pc; at this time the kinetic energy deposited by the region –expanding with a velocity of $`10`$ km s-1 on average– into the ISM is $`12\times 10^{49}`$ erg (e.g., Elmegreen 1991). The cooled shell swept up by a SN (type II) ejecta has roughly the same amount of kinetic energy, although this energy is injected in a time period much smaller than in the case of the HII region stage.

We start the simulations with uniform density and temperature and nearly zero velocity everywhere ($`u_{\mathrm{rms}}0.1`$km s-1), and begin placing sources randomly as described in §2.2. After one source lifetime, the number of sources in the simulation reaches an approximately stationary number. It is not precisely constant because of the probability-based algorithm we use and the variable timestep. The sources generate expanding shells which interact in a complicated fashion, producing a network of filaments and distorted shells, in a statistically stationary regime. Figure 2 shows the density and $`x`$-component of the acceleration at a typical instant in the evolution of the fiducial simulation, run 30.

We have also performed a simulation with the density-directed scheme of source placement (run 6) in order to compare our results to this more realistic simulation. The density threshold used for turning on a source was $`4\rho _0`$. The input parameters for this run were $`l_\mathrm{f}=30`$ pc and $`u_\mathrm{f}=7`$ km s-1 (estimated from the final expansion velocity of the shells). Finally, in order to study a strictly decaying regime, we performed a simulation without SF activity (run 12 in Table 1). This run used a late evolutionary stage of run 6 as its initial condition.

## 3 Dissipation in driven turbulence

### 3.1 The global dissipation rate

For all the forced runs described in §2.3 we measure several global quantities of the flow: the kinetic energy content, $`E_\mathrm{k}`$, the kinetic energy injection rate, $`\dot{E}_\mathrm{k}^{\mathrm{in}}`$, the total rate of change of $`E_\mathrm{k}`$, $`dE_\mathrm{k}/dt`$, and the rms velocity dispersion, $`u_{\mathrm{rms}}`$. Figure 3 shows the evolution of these quantities for the fiducial run except $`u_{\mathrm{rms}}`$, which we do not plot because, as shown in Fig. 4, it scales very closely as $`E_\mathrm{k}^{1/2}`$ (see also Mac Low 1999a). These global quantities are seen to be statistically stationary, without exceedingly strong fluctuations in time. The injection rate $`\dot{E}_\mathrm{k}^{\mathrm{in}}`$ is seen to fluctuate moderately, typically some 50% around its mean value. This is not only because the number of sources varies slightly in time, but also because of the complex interaction between the acceleration from the sources and the pre-existing velocity field (see eq. ). Similar behaviour is observed in all forced runs with random source placement. However, a much more irregular behavior is present (not shown) in the run with density-directed SF, in which the SFR alternates between zero and strong bursts. Table 2 gives the average of these quantities between $`t=40.8`$ and 272 Myr for all runs.

From the measured values of $`E_\mathrm{k}`$ and $`\dot{E}_\mathrm{k}^{\mathrm{in}}`$, we can define the energy injection timescale as

$$t_{\mathrm{in}}E_\mathrm{k}/\dot{E}_\mathrm{k}^{\mathrm{in}}.$$
(3)
However, the actual quantity of interest is $`t_\mathrm{d}`$, the kinetic energy dissipation timescale, defined as

$$t_\mathrm{d}=E_\mathrm{k}/\dot{E}_\mathrm{k}^{\mathrm{diss}},$$
(4)
where $`\dot{E}_\mathrm{k}^{\mathrm{diss}}`$ is the energy dissipation rate. In a perfectly stationary regime, $`\dot{E}_\mathrm{k}^{\mathrm{in}}`$ and $`\dot{E}_\mathrm{k}^{\mathrm{diss}}`$ coincide, as in the well-known Kolmogorov (1941) theory in the incompressible case, and in the simulations of Stone et al. (1998) and Mac Low (1999a) in the compressible case. However, in the more realistic case of a regime stationary only in an averaged sense, the two rates need not coincide. We thus estimate $`\dot{E}_\mathrm{k}^{\mathrm{diss}}`$ from the equation:

$$\frac{dE_\mathrm{k}}{dt}=\dot{E}_\mathrm{k}^{\mathrm{in}}\dot{E}_\mathrm{k}^{\mathrm{diss}},$$
(5)
where $`dE_\mathrm{k}/dt`$ can be measured directly in the simulation from the evolution of $`E_\mathrm{k}`$. Afterwards, $`t_\mathrm{d}`$ can be computed from eq. (4). The evolution of $`\dot{E}_\mathrm{k}^{\mathrm{diss}}`$ and $`t_\mathrm{d}`$ is also shown in fig. 3. The latter quantity is shown averaged over periods of 20 successive outputs of the code, because of its high sensitivity to fluctuations in $`\dot{E}_\mathrm{k}^{\mathrm{in}}`$. The values of $`\dot{E}_\mathrm{k}^{\mathrm{diss}}`$ and $`t_\mathrm{d}`$ averaged between $`t=40.8`$ and 272 Myr are given also in Table 2.

It should also be noted that eq. (5) is actually incomplete, and should include transfer terms from (to) kinetic to (from) other energy forms, such as magnetic or thermal. Thus, such transfer terms are effectively included in $`\dot{E}_\mathrm{k}^{\mathrm{diss}}`$. In particular, this includes kinetic energy losses due to first converting it to those other forms and then dissipating these by resistivity or radiative cooling. On the other hand, this implies that $`\dot{E}_\mathrm{k}^{\mathrm{diss}}`$ may also act as a source of kinetic energy, in which case it may become negative. However, it is seen that $`\dot{E}_\mathrm{k}^{\mathrm{diss}}<0`$ only during a brief period at the beginning of the simulation, which, together with the similarity between $`\dot{E}_\mathrm{k}^{\mathrm{in}}`$ and $`\dot{E}_\mathrm{k}^{\mathrm{diss}}`$ discussed below, reassure us that the error incurred by omitting the transfer terms in eq. (5) is not too severe.

Remarkably, fig. 3 shows that $`\dot{E}_\mathrm{k}^{\mathrm{in}}`$ and $`\dot{E}_\mathrm{k}^{\mathrm{diss}}`$ are always very close to each other. Since the energy injection occurs at very localized and scattered sites, while the dissipation rate is computed over the whole flow, their highly similar evolutions suggest that most of the dissipation must occur at or near the injection regions, and dominates the global dissipation rate444A posteriori, this is actually not surprising; in stationary forced incompressible regimes, it is well known that the dissipation rate is proportional to $`E_\mathrm{k}^{3/2}`$ (see, e.g., Frisch 1995; see also Stone et al. 1998; Mac Low 1999a for the compressible case). Since the kinetic energy density is greatest at the injection sites, most of the dissipation must occur there as well. Thus, the dissipation rate is actually a function of position (through the local kinetic energy density) and may fluctuate strongly from one position to another, as is also the case in incompressible turbulence (see, e.g., Frisch 1995).. Only in this case can the global dissipation rate track the locally-originated injection so closely. In turn, this implies that only a small fraction of the energy injected at the sources “escapes” to more distant regions. We refer to this as the “residual” turbulence, which coexists in the flow with the forced regions. We discuss the propagation of turbulent kinetic energy from the injection to the quiescent sites in §4. As a consequence, in the remainder of this paper we refer to the injection and the dissipation times indistinctly.

Analyzing all the forced runs, we find that $`t_\mathrm{d}`$ ($`t_{\mathrm{in}}`$) does not significantly depend on two of the energy injection parameters, namely the energy input rate of the source, $`\dot{e}_\mathrm{k}`$ (or just its total energy input $`e_\mathrm{k}`$) and the source formation rate, $`\dot{\mathrm{\Sigma }}_{\mathrm{OB}}`$ (see §2.2 for definitions). The relevant forcing parameter seems to be the forcing scale, $`l_\mathrm{f}`$. We obtain the following approximate empirical correlation of $`t_\mathrm{d}`$ on the input parameters of the forcing:

$$t_\mathrm{d}\frac{l_\mathrm{f}^{0.60}}{(\dot{ϵ}_\mathrm{k}\dot{\mathrm{\Sigma }}_{\mathrm{OB}})^{0.12}}\frac{l_\mathrm{f}^{0.84}}{(e_\mathrm{k}\dot{\mathrm{\Sigma }}_{\mathrm{OB}})^{0.12}}.$$
(6)
### 3.2 Relating $`t_{\mathrm{in}}`$ to measurable properties of the ISM

In order to estimate theoretically $`t_\mathrm{d}(t_{\mathrm{in}})`$ as a function of frequently used astronomical quantities of the ISM, we can use the definition of $`t_{\mathrm{in}}`$:

$$t_{\mathrm{in}}E_\mathrm{k}/\dot{E}_\mathrm{k}^{\mathrm{in}}=\frac{_A\mathrm{\Sigma }_\mathrm{g}u^2𝑑A^{}}{_A\mathrm{\Sigma }_\mathrm{g}𝐚𝐮𝑑A^{}}.$$
(7)
The integral in the numerator can be estimated as

$$_A\mathrm{\Sigma }_\mathrm{g}u^2𝑑A^{}\mathrm{\Sigma }_\mathrm{g}u_{\mathrm{rms}}^2A,$$
(8)
where $`\mathrm{\Sigma }_\mathrm{g}`$ is the average gas surface density and $`A`$ is the total area where the flow is contained. Similarly to what was done in eq. (2), the integral in the denominator can be estimated as follows:

$`{\displaystyle _A}\mathrm{\Sigma }_\mathrm{g}𝐚𝐮𝑑A^{}={\displaystyle _{A_\mathrm{f}}}\mathrm{\Sigma }_\mathrm{g}𝐚𝐮𝑑A^{}\mathrm{\Sigma }_\mathrm{g}_\mathrm{f}a_\mathrm{f}u_\mathrm{f}A_\mathrm{f}`$
$`=\mathrm{\Sigma }_\mathrm{g}_\mathrm{f}\dot{ϵ}_\mathrm{k}A_\mathrm{f},`$ (9)
where the first equality follows because $`𝐚=0`$ outside of the forcing regions, whose total area is $`A_\mathrm{f}`$. Therefore, from eqs. (7), (8), and (9) we obtain

$$t_{\mathrm{in}}\frac{u_{\mathrm{rms}}^2}{\dot{ϵ}_\mathrm{k}f},$$
(10)
where $`f=A_\mathrm{f}/A`$ is the filling factor of forcing regions. Naively, one could take $`f=\pi l_\mathrm{f}^2n_\mathrm{s}/A`$. However, one has to take into account the fact that some of the forcing regions overlap; the probability of this occurring is larger at larger $`l_\mathrm{f}`$. Thus, $`f\genfrac{}{}{0pt}{}{_<}{^{}}\pi l_\mathrm{f}^2n_\mathrm{s}/A`$ in general, and in practice we prefer to measure $`f`$ directly from the simulation as the fraction of the simulation area in which the stellar acceleration is larger than 20% of $`a_{\mathrm{max}}`$ (see eq. ), corresponding to the area within a distance $`l_\mathrm{f}`$ from the centers of the forcing.

Table 2 gives the values of $`u_{\mathrm{rms}}`$ and $`f`$ measured in all the forced runs ($`\dot{ϵ}_\mathrm{k}`$ is an input parameter given in Table 1), as well as the values of $`t_{\mathrm{in}}`$ resulting from them. These “theoretical” injection times are seen to be reasonably close to those measured directly in the simulations according to eq. (3). For most of the runs, the agreement is better than 20%, with the largest deviations being $`45`$% and $`30`$%.

A posteriori, it is not surprising that eq. (10) gives estimates for $`t_{\mathrm{in}}`$ consistent with the values obtained using eq. (3), since, after all, we have done nothing more than estimating the integrals intervening in eq. (7) in terms of averaged quantities. However, the relevance of this result is twofold. First, it shows that reasonable estimates can be obtained using global averaged quantities even when the dissipation is extremely localized, and second, it provides a means of relating the dissipation time to frequently used quantities of the ISM. Indeed, using the expression $`f\genfrac{}{}{0pt}{}{_<}{^{}}\pi l_\mathrm{f}^2n_\mathrm{s}/A`$ for the filling factor in eq. (9) we obtain

$`{\displaystyle _A}\mathrm{\Sigma }_\mathrm{g}𝐚𝐮𝑑A^{}{\displaystyle \genfrac{}{}{0pt}{}{_<}{^{}}}\mathrm{\Sigma }_\mathrm{g}_\mathrm{f}\dot{ϵ}_\mathrm{k}\pi l_\mathrm{f}^2n_\mathrm{s}=`$
$`M_\mathrm{f}\dot{ϵ}_\mathrm{k}\mathrm{\Delta }t_{OB}\dot{\mathrm{\Sigma }}_{\mathrm{OB}}A=e_\mathrm{k}\dot{\mathrm{\Sigma }}_{\mathrm{OB}}A,`$ (11)
where $`M_\mathrm{f}=\mathrm{\Sigma }_\mathrm{g}_\mathrm{f}\pi l_\mathrm{f}^2`$, $`e_\mathrm{k}=M_\mathrm{f}\dot{ϵ}_\mathrm{k}\mathrm{\Delta }t_{OB}`$, and $`\dot{\mathrm{\Sigma }}_{\mathrm{OB}}`$ has been defined in §2.2. Using eqs. (8) and (11) in eq. (7), one obtains

$$t_{\mathrm{in}}\frac{\mathrm{\Sigma }_\mathrm{g}u_{\mathrm{rms}}^2}{e_\mathrm{k}\dot{\mathrm{\Sigma }}_{\mathrm{OB}}},$$
(12)
which is our desired connection to measurable quantities of the ISM. For example, taking $`u_{\mathrm{rms}}=8`$ km s-1, $`\dot{\mathrm{\Sigma }}_{\mathrm{OB}}=3\times 10^5`$ yr-1 kpc-2 (Tammann, Loffler & Schroeder 1994), $`e_\mathrm{k}=3.5\times 10^{49}`$ erg, and the gas surface density typical for the solar vicinity, $`\mathrm{\Sigma }_\mathrm{g}12M_{}`$pc-2, we find that $`t_{\mathrm{in}}1.5\times 10^7`$ yr. Note that the value of $`e_\mathrm{k}`$ we used reflects a contribution of roughly $`1.5\times 10^{49}`$ ergs per HII region plus $`2\times 10^{49}`$ ergs of mechanical energy per SN event. The resulting value of $`t_{\mathrm{in}}`$ is in rough agreement with our fiducial simulation.

It is worthwhile to note that $`t_{\mathrm{in}}`$ depends on $`\mathrm{\Sigma }_\mathrm{g}`$ and $`\dot{\mathrm{\Sigma }}_{\mathrm{OB}}`$, and since both of these vary with Galactocentric distance, one expects $`t_{\mathrm{in}}`$ (and consequently the dissipation time) to exhibit such a dependence. This question will be discussed in further detail in section §6.2.

One more point should be remarked regarding eq. (10) (or, equivalently, relation ), namely that the numerical simulations clearly show that the dissipation time increases with the size of the forcing regions (see eq. ), even though relation (12) has no explicit dependence with $`l_\mathrm{f}`$. This is equivalent to asking the question: how does the empiric expression for $`t_{\mathrm{in}}`$ (eq. ) connect with relation (12), which was derived directly from the definition of $`t_{\mathrm{in}}`$? These two relations imply that $`u_{\mathrm{rms}}`$ must scale approximately as $`l_\mathrm{f}^{0.42}(e_\mathrm{k}\dot{\mathrm{\Sigma }}_{\mathrm{OB}})^{0.44}`$. Recalling that in relation (12) the inequality accounts for overlapping of the regions, one expects that an extra factor depending on $`l_\mathrm{f}`$ must be included in the RHS of relation (12) in order for it to become an equality. In the simulations we empirically find that roughly $`u_{\mathrm{rms}}l_\mathrm{f}^{0.40}(e_\mathrm{k}\dot{\mathrm{\Sigma }}_{\mathrm{OB}})^{0.42}`$. A full MHD theory of the problem should account for this dependence.

### 3.3 A more realistic simulation

In the runs presented above, we have not included self-gravity nor the Coriolis force, and the energy input sources were placed randomly with a fixed probability, in order to have good control of the forcing parameters, and to uncover the dependences presented in $`\mathrm{\S }\mathrm{\S }`$ 3.1 and 3.2, at the expense of some realism. Somewhat more realistic simulations should include those additional agents, and involve a density-directed placement of the sources. In particular, the latter scheme causes the energy injection rate to fluctuate strongly and the input sources to be highly clustered; large supershells form and significant cancellation of the energy input from neighboring sources occurs (see §2.2). It is thus interesting to see if our results still hold in this case.

To this effect, we have also performed a simulation with these ingredients (run 6, see §2.3). In fig. 5 we present images of density field of this simulation at times $`t=14.3`$ Myr (left) and $`t=144.3`$ Myr (right). At the earlier time, small shells are seen, whose sizes are comparable to the sphere of influence of the “winds” themselves (of radius $`l_\mathrm{f}=30`$ pc). At the later time, the shells are seen to be larger, but their larger sizes are not only a consequence of the inertial motions induced by the forcing, but also of induced-SF events in the shells themselves. Even though for this run we cannot estimate $`u_\mathrm{f}`$ precisely, and consequently $`\dot{ϵ}_\mathrm{k}`$, we can measure the total injection rate and, from simple visual inspection, we see that the average radius of the shells in the simulation is at least twice as large as $`2l_\mathrm{f}`$. Moreover, on average we measure $`E_\mathrm{k}0.18`$ (in code unities), $`u_{\mathrm{rms}}7`$ km s-1, and $`t_\mathrm{d}18`$ Myr. We thus see that this run is in a very similar regime as runs 30 and 35, and in particular the dissipation time is within 15% of that measured in these runs, and is 30% larger than the lower-bound estimate for the solar neighborhood given in §3.2. Thus, this run reassures us that the results obtained with the random placement of the stellar sources are applicable in more realistic situations.

## 4 Dissipation in decaying turbulence

In order to study the decay of turbulence far from the injection sites, we now consider a simulation in a decaying (i.e., without energy injection) regime. To this end, we consider a continuation of run 6, named run 12, but with the SF turned off. In fig. 6 we present images of the density field from this run at two times, one immediately after turning off the SF ($`t=0`$) and the other at a much more advanced stage ($`t=180.7`$ Myr). The former time shows a large number of expanding shells, while the latter contains no shells anymore. The evolution of the kinetic energy for this simulation is shown in fig. 7. We find that $`E_\mathrm{k}`$ decays roughly as

$$E_\mathrm{k}(t)=E_0\left(1+\frac{t}{t_0}\right)^\alpha $$
(13)
with $`\alpha 0.8`$, in reasonably good agreement with previous studies of decaying MHD turbulence for isothermal flows (Mac Low et al. 1998, Stone et al. 1998), and $`t_0`$ is the code time unit. A characteristic decay time can be defined in this case as the time $`E_\mathrm{k}`$ to decay by a factor of two at the beginning of the simulation. For $`\alpha 0.8`$, this occurs at $`t=1.38\times t_0=18`$ Myr.

We can explore the decay of turbulent motions somewhat further. This analysis is greatly simplified in the case of long times ($`t_0`$), which allows us to drop the term 1 inside the parenthesis in eq. (13), writing

$$E_\mathrm{k}(t)E_0\left(\frac{t}{t_0}\right)^\alpha $$
(14)
This limit is reasonable, since we are interested in the behavior of the turbulence far from the injection sites. Then, assuming that the rms velocity $`u_{\mathrm{rms}}`$ scales roughly as the square root of the total kinetic energy (see Fig. 4), from eq. (13) we can write:

$$u_{\mathrm{rms}}(t)u_1\left(\frac{t}{t_0}\right)^{\alpha /2},$$
(15)
where $`u_1`$ is the rms velocity dispersion at $`t=t_0`$ ($`7`$ km s-1). Similarly to the forced case, here we may formally define the dissipation time scale as the ratio between $`E_\mathrm{k}`$ and $`|\dot{E}_k|`$. Thus, we obtain:

$$t_\mathrm{d}(t)\frac{E_\mathrm{k}}{|\dot{E}_k|}\frac{t_0}{\alpha }\left(\frac{t}{t_0}\right).$$
(16)
Thus, the characteristic dissipation time increases linearly with time in the decaying case, implying that the instantaneous dissipation time of the flow depends on the net energy content of the flow, being smaller for larger energy contents.

Identifying now the decaying case with the regions far from the injection sites in a forced case with small-scale, spatially intermittent forcing, the temporal decay of the energy can be translated into a decay with distance. We can estimate the distance “reached” by the turbulent motions in a time $`t`$ as

$$\mathrm{}(t)u_{\mathrm{rms}}(t)t\mathrm{}_0\left(\frac{t}{t_0}\right)^{1\alpha /2},$$
(17)
where $`\mathrm{}_0=u_1t_010^2`$ pc is the typical length traveled by a fluid parcel at the characteristic time and velocity at the beginning of the decaying regime. This length does not exceed the size scale of the expanding shells produced by the forcing. Solving for $`t`$ in eq. (17) and substituting in eqs. (14) and (15), we obtain the dependence of $`E_\mathrm{k}`$ and $`u_{\mathrm{rms}}`$ with $`\mathrm{}`$:

$$E_\mathrm{k}(\mathrm{})E_0\left(\frac{\mathrm{}}{\mathrm{}_0}\right)^{\frac{2\alpha }{2\alpha }}$$
(18)
$$u_{\mathrm{rms}}(\mathrm{})u_1\left(\frac{\mathrm{}}{\mathrm{}_0}\right)^{\frac{\alpha }{2\alpha }}$$
(19)
For $`\alpha 1`$ this implies that $`E_\mathrm{k}`$ and $`u_{\mathrm{rms}}`$ decay with distance as $`(\mathrm{}/\mathrm{}_0)^2`$ and $`(\mathrm{}/\mathrm{}_0)^1`$, respectively. Thus, “residual” turbulent motions with an rms velocity of roughly 6-10 km s-1, as observed in our and other disk galaxies, would decay to roughly 1/10 of that value on distances $`0.61`$ kpc, in the absence of energy sources. Note however that, since the decay law (14) is a power law (rather than, say, an exponential), the time required for $`E_\mathrm{k}`$ to decrease by consecutive constant factors increases as the the kinetic energy content of the flow decreases (cf. eq. ). This applies also to the decay with distance, and thus the distance required for the rms velocity to decrease by a factor of 10 depends on the initial energy content of the medium. This is also why energy is dissipated so efficiently near the injection sites.

## 5 Effects of resolution and dimensionality

For economy, and to maximize our ability to explore parameter space, in this paper we have limited ourselves to low-resolution, 2D numerical simulations. However, even though in this paper we are not attempting to give detailed quantitative results, but only first-order estimates, it is important to test that even these can indeed be attained with the simulations we have employed.

The two-dimensionality poses an apparently serious problem, since it is well known that, in the incompressible case, the energy dissipation rate has a fundamentally different behavior in 2D with respect to 3D. While in 3D this rate is believed to remain finite even in the limit of vanishing viscosity, in 2D it does tend to zero in this limit (see, e.g., Lesieur 1990, sec. IX.3.5). Moreover, in 2D the kinetic energy cascades in an “inverse” way, from small to large scales (always in the incompressible case). All of this questions the validity of using 2D simulations for estimates of the dissipation rate and, in particular, implies that convergence tests may be meaningless in 2D since the kinetic energy dissipation rate does not converge as the resolution is increased and the dissipation coefficients are decreased.

Fortunately, this problem may not be present in the highly compressible case, in which shocks dominate the dissipation rate, since it is well known that shocks cause direct transfer from all scales to the dissipation scales (Kadomtsev & Petviashvili 1973). This process is independent of the dimensionality of the flow since the shocks are essentially one-dimensional structures. Moreover, the dissipative nature of the shocks suggests that the problem of vanishing dissipation in the limit of vanishing viscosity in 2D is eliminated in their presence. This suggestion is supported by the fact that in our forced simulations the dissipation rate follows the injection rate closely. In summary, we are confident that in the highly compressible 2D case, the dissipation rate can still be meaningfully measured, and that convergence tests can be performed.

To this effect, we have produced a higher-resolution simulation in the decaying regime, at $`512^2`$ grid points, still in 2D, labeled run 17. In this run, we have reduced the hyperviscosity coefficient by a factor of $`4^7=16384`$, and the second-order viscosities and diffusivities by a factor of $`4^1=4`$, in an attempt to cause the dissipative and diffusive scales to be a factor of 4 smaller than in the $`128^2`$ runs, i.e., roughly the same size in grid points. The factors chosen take into account the order of the hyperviscous and second-order viscous operators (8 and 2, respectively), and the fact that the mean Fourier amplitude of the velocity scales as $`k^1`$ if the kinetic energy spectrum $`E(k)`$ scales as $`k^2`$, as expected in a flow dominated by shocks (see, e.g., the review by Vázquez-Semadeni et al. 2000 and references therein). As is the case with the decaying run at $`128^2`$ resolution (run 12), run 17 is again a restart of a forced simulation which used the old scheme of placing the stars at the highest density sites.

Figure 8 shows the evolution of the kinetic energy $`E_\mathrm{k}`$ for run 17. It is seen that over the early epochs of the simulation ($`\xi \mathrm{log}(1+t/t_0)\genfrac{}{}{0pt}{}{_<}{^{}}0.4`$), the energy varies nearly as a power law with time, with exponent $`0.9`$, in rough agreement with the $`0.8`$ exponent observed for run 12, and also with the exponents reported by Mac Low et al. (1998) and Stone et al. (1999). However, at later times, a departure from this behavior is observed, with the decay rate of $`E_\mathrm{k}`$ monotonically decreasing for $`0.4\genfrac{}{}{0pt}{}{_<}{^{}}\xi \genfrac{}{}{0pt}{}{_<}{^{}}0.9`$, indicating that at the higher resolution of run 17 the decay rate indeed becomes much smaller after the shocks have subsided. To support this interpretation, in the same figure we also show the evolution of the most negative value of the divergence of the velocity field (shown in log of its absolute value). It is seen that the departure from a power-law decay rate coincides with a drop of the divergence by nearly a factor of 3 at $`\xi 0.5`$. Moreover, a sharp drop in $`E_\mathrm{k}`$ at $`\xi 0.95`$ coincides with a large peak in the divergence, after which the latter decreases even further and $`E_\mathrm{k}`$ returns to a slow decay rate ($`\xi 1.1`$). In summary, Fig. 8 supports the view that in 2D the dissipation rate follows roughly the same decay rate as 3D simulations (Mac Low et al. 1998; Stone et al. 1999) as long as the shocks dominate the dissipation, giving us confidence that our 2D simulations provide realistic measurements of the dissipation rate in the forced regime, which contains large numbers of shocks, and in the early epochs of the decaying regimes.

Concerning the forced case, a potential problem is that for some of the low-resolution runs, the injection scale overlaps with the dissipation scale (J. Brasseur 2000, private communication), and thus there is the risk that this energy never reaches the larger scales of the flow. We believe this is not really a problem, because, as mentioned above, shocks “short-circuit” the energy directly from the energy-containing scales to the dissipations scales, so the cascade through intermediate scales is eliminated anyway. In physical space, what happens is that, even if the forcing has characteristic scales larger than the dissipation scale, its effect is to induce shocks into the medium, which are small-scale structures. Indeed, runs 40 and 41, which have a value of $`l_\mathrm{f}`$ twice as large as that of run 30, still agree with our model prediction within $`45`$% and $`30`$%, respectively.

Nevertheless, it is important to verify that the agreement between our semi-empirical model of §3.2 and the low resolution simulations is maintained in higher-resolution runs. To this effect, we have performed a $`512^2`$ forced run, labeled run 44 (Table 1) in which the physical size of the box (1 kpc) is maintained, but the smallest resolved scale is now 1.95 pc. This run uses the same viscosities and difussivities as run 17, and its parameters are chosen so that it is a rescaling of run 30 (the fiducial run): $`l_\mathrm{f}`$ is 4 times larger in pixels, so it keeps the same physical size in parsecs. However, the dissipation scale is four times smaller in physical size, becoming safely smaller than the injection scale. All other parameters are the same as in run 30. From Table 2, it can be seen that run 44 conforms nicely to our semi-empirical model, its predicted and measured injection times agreeing within $`15`$%.

## 6 Discussion

### 6.1 Comparison with previous work

The main difference between our MHD simulations of supersonic compressible flows and previous ones aimed at exploring turbulence dissipation resides in the method we use to drive the turbulence. As pointed out in §2.2, we have tried to represent the way in which kinetic energy is injected into the ISM of disk galaxies by stellar sources. Moreover, since our interest in this paper focuses in the large-scale ISM, we have taken parameters of the flow corresponding to the solar neighborhood, as done in previous papers modeling the ISM at large (Chiang & Bregman 1988; Rosen & Bregman 1995; Vázquez-Semadeni et al. 1995, 1996; Passot et al. 1995), and have not assumed the flow to be isothermal. All of this contrasts with recent studies of turbulence dissipation in molecular clouds (Mac Low et al. 1998; Stone et al. 1998; Padoan & Nordlund 1999; Mac Low 1999a), which have considered isothermal flows with large-scale, ubiquitous forcing. On the other hand, those works have been based on 3D simulations, while here we have used 2D simulations. However, large differences on the dissipation rates and timescales between 2D and 3D simulations are not expected (see e.g., Stone et al. 1998 and §5 above). Specifically, the kinetic energy dissipation rates in 3D are expected to be slightly higher than in 2D, particularly in the regime of decaying turbulence, because at higher dimensionality, more shocks are produced due to the larger number of degrees of freedom (e.g., Mac Low et al. 1998). For example, Stone et al. (1998) reported that the energy decay times found for $`2+1/2`$D MHD simulations are a factor of 1.50—1.75 times larger than those for 3D.

At the results level, in spite of the very different physical regimes we consider, we agree with previous works in that compressible MHD flows dissipate their turbulent kinetic energy very efficiently. More specifically, for their molecular cloud simulations, Stone et al. (1998) and Mac Low (1999a) find that the dissipation time is of the order to the flow crossing time at the energy-containing scale $`\lambda _\mathrm{E}`$, i.e. $`t_\mathrm{d}\lambda _\mathrm{E}/u_{\mathrm{rms}}`$. However, in their work, this scale is very close to the injection scale because of the large-scale and ubiquitous nature of the injection. In our case, on the other hand, $`\lambda _\mathrm{E}`$ is in general quite larger than the driving scale $`2l_\mathrm{f}`$, because of the expansion of the shells, which is a way of transferring the energy to larger scales, and possibly also because energy that is not dissipated in shocks may escape towards the largest scales through an inverse cascade due to the two-dimensionality (to be investigated elsewhere). Indeed, fig. 9a shows the energy spectrum for run 30 at times $`t=1.36`$, 6.8 and 136 Myr. For reference, fig. 9b shows the $`x`$-component of the velocity field at $`t=6.8`$ Myr. It can be seen that at $`t=1.36`$ Myr, the spectrum contains only one peak at $`\mathrm{log}k1.4`$ (the bulk of the spectrum at lower $`k`$ at this time corresponds to the very mild initial conditions with $`u_{\mathrm{rms}}0.1`$ km s-1). A scale (in pixels) corresponding to wavenumber $`k`$ can be defined as $`\lambda =128/k`$. For $`\mathrm{log}k=1.4`$, $`\lambda 5`$ pixels, and indeed at this time there is only one expanding shell in the simulation with “radius” $`4`$ pixels, measured as the distance from the center to the velocity maximum of the shell (this distance has been measured from cuts through the $`x`$-component of the velocity passing through the shell center, similar to those in fig. 1, not shown).

At $`t=6.8`$ Myr, this peak has moved to $`\mathrm{log}k1.22`$ ($`\lambda 7.7`$ pixels) due to the expansion of the shell, whose radius has grown to $`6`$ pixels. New peaks at larger $`k`$ have appeared corresponding to the formation of new sources. The maximum of the spectrum (roughly, the “energy-containing scale”) at $`t=6.8`$ Myr, at $`\mathrm{log}k0.8`$, corresponds to the diameter of the largest shell at that time. Energy at even larger scales must correspond to the size of the region spanned by the various sources at that time, since the shells themselves have not grown further than $`\mathrm{log}k=0.8`$ yet. By $`t=136`$ Myr, the spectrum is very close to a power law in the range $`0.6\genfrac{}{}{0pt}{}{_<}{^{}}\mathrm{log}k\genfrac{}{}{0pt}{}{_<}{^{}}1.5`$, suggesting the existence of fully developed turbulence at that time.

Since the results above suggest a clear separation between the injection and energy-containing scales in our simulations, in order to compare our results to the works of Stone et al. and of Mac Low, we can quantitatively estimate the energy-containing scale in our simulations as the scale corresponding to the centroid of the energy distribution, given by (see, e.g., Novikov 1978; Santangelo et al. 1989; McWilliams 1990)

$$\lambda _\mathrm{E}2\pi \left(\frac{_𝐤k|𝐮_𝐤|^2}{_𝐤|𝐮_𝐤|^2}\right)^1,$$
(20)
where $`𝐮_𝐤`$ is the Fourier amplitude of the mode with (vector) wavenumber k. Table 2 gives the values of $`\lambda _\mathrm{E}`$ for all runs at $`t=136`$ Myr, a time at which the spectrum has become essentially stationary, and the value of the large-scale crossing time $`t_{\mathrm{cr}}\lambda _\mathrm{E}/u_{\mathrm{rms}}`$, calculated with the rms velocity at that same time. It is seen that the crossing times do not correlate well with the measured dissipation times, and are in general quite larger. We speculate that this is due to the fact that $`t_{\mathrm{cr}}`$ measures essentially the dissipation time at the largest scales, and fails to capture the local dissipation near the sources. In particular, the local and large-scale dissipations contribute in different proportions to the total rate depending on the size, number and specific energy injection rate of the sources, and thus the crossing time is not a good estimator of the true injection and dissipation times in the case of non-ubiquitous, small-scale forcing.

In conclusion, we do not expect that a simple law of the type $`t_\mathrm{d}=\lambda _\mathrm{E}/U`$, with $`U`$ either the rms of the characteristic forcing velocity, can hold in the case with non-ubiquitous forcing: another parameter, the filling factor of the forcing regions, is crucial, as it determines the total amount of energy injected by sources of given local parameters, the proportion in which $`u_\mathrm{f}`$ contributes to $`u_{\mathrm{rms}}`$ and the amount of overlapping between the regions.

### 6.2 Implications for models of galaxy evolution

From the point of view of galactic evolution, the global ISM processes, which take place at large space and time scales, are the most relevant. As mentioned in the Introduction, several galactic evolution models are based on the idea that the SF rate is controlled by a disk vertical balance in the ISM between the rates of energy injection due to SF and energy dissipation. Since turbulent kinetic energy is believed to be one of the dominant sources of pressure in the ISM, the relevant dissipation rate for these models is that of turbulence, motivating our study. We have shown that $`E_\mathrm{k}`$ is dissipated very efficiently near the input sources (locally). This implies that very large “active” turbulent zones can exist in galaxies only if the forcing regions are themselves large, as a consequence of either very powerful sources, or strong clustering and self-propagating SF. This might be the rule for starburst galaxies, but generally the exception in normal disk galaxies. Besides, as results from our simulations have shown (§3), the dissipation time $`t_\mathrm{d}`$ indeed increases as the source size and filling factor increase.

In the energy balance in the vertical direction assumed in the galaxy evolution models mentioned above (taking a given energy input rate, given by the SF rate, and estimating the dissipation rate as $`u_{\mathrm{rms}}^2/t_\mathrm{d}`$), the equilibrium gas velocity dispersion is determined by $`t_\mathrm{d}`$. Hence, assuming vertical (one-zone) hydrostatic equilibrium, the disk HI thickness will depend basically on $`t_\mathrm{d}`$. According to the results of these galaxy evolution models (e.g., Firmani et al. 1996; Avila- Reese & Firmani 2000; Firmani & Avila-Reese 2000), a realistic HI gas thickness at the solar neighborhood is obtained with values of $`t_\mathrm{d}`$ roughly twice as large as the values found here (all the other relevant quantities as gas surface density, circular velocity and star formation rate, are self-consistently calculated and agree well with observations). From the hydrostatic equilibrium condition, indeed it is easy to see that a turbulent gas layer supported by a velocity dispersion $`9`$ km s-1 would have a thickness roughly 3 times smaller than observed, suggesting that probably other kinds of pressure, such as the magnetic field and cosmic-rays, are also important for supporting the galactic gas layer.

On the other hand, it is important to have in mind that the vertical galactic disk is actually stratified; as the gas density decreases, the survival and propagation of turbulent motions could be easier, and therefore $`t_\mathrm{d}`$ would be longer. In fact, shells are well known to accelerate when they transit from a high- to a low-density medium. In this sense, the dissipation times calculated here for a fixed gas average density can be interpreted as a lower limit; hence, the disk thicknesses estimated with our turbulent dissipation times will be also a lower limit. In any case, according to the results of galaxy evolution models, with self-regulated SF, the dissipation time obtained here is not too far from the values necessary to reproduce several observational data. Thus, our results do not disagree with the idea that turbulent motions induced by SF feedback contribute significantly to support the vertical HI disk (e.g., Lockman & Gehman 1991; Firmani & Tutukov 1992, 1994; Ferrara 1993; Dopita & Ryder 1994), although some extra support may be necessary to keep the HI and gas ionized layers up at the observed heights.

Finally, it is important to remark that the turbulent dissipation time may change along the galactic disk since the gas densities and dynamical conditions change with radius. According to Firmani et al. (1996), $`t_\mathrm{d}R`$ when the rotation curve is flat. In §3.2 we obtained a theoretical expression for $`t_\mathrm{d}`$ (actually $`t_{\mathrm{in}}`$; eq. ) which in particular implies a dependence on galactic radius through $`\mathrm{\Sigma }_\mathrm{g}`$ and $`\dot{\mathrm{\Sigma }}_{\mathrm{OB}}`$. The OB-star formation rate is directly proportional to the total formation rate $`\dot{\mathrm{\Sigma }}_{}`$ and, in most disk galaxies, $`\dot{\mathrm{\Sigma }}_{}\mathrm{\Sigma }_\mathrm{g}^n`$, with $`n1.42.0`$ (e.g., Kennicutt 1998 and references therein). Therefore, from eq. (12) one obtains that $`t_{\mathrm{in}}(r)\mathrm{\Sigma }_\mathrm{g}^{n1}(r)`$. The gas surface density typically decreases with radius. For example, in the Galaxy, from 4 to 16 kpc, $`\mathrm{\Sigma }_\mathrm{g}`$ decreases roughly as $`r^1`$ (Dame 1993). Hence, the dissipation time should decay as the Galactocentric radius $`R`$ to some power, probably less than unity.

### 6.3 Implications for models of galaxy formation

Several approaches of galaxy formation within the context of the hierarchical CDM-based scenario —in particular the so-called semi-analytical models (e.g., Kauffmann et al. 1993; Cole et al. 1994; Somerville & Primack 1999; van den Bosch 1999)— require the feedback from SF to be able to reheat the cooled gas up to the virial temperature of the whole dark matter halo. This requirement implies that a significant fraction of the released SN energy (probably more than 10%) remains as kinetic energy and that this energy is not dissipated within the disk.

The physical picture invoked in the semi-analytical models implies a self-regulated SF not at the level of the disk ISM, but at the level of the hypothetical intra-halo medium. The SNe formed in the disk are assumed to reheat the cold disk gas up to the virial temperature of the cosmological dark matter halo, driving it back into the intra-halo medium, and occasionally expelling it completely from the system (White & Frenk 1991). Thus, a crucial question for models of galaxy formation is whether the energy released by SNe and stars in the disk is able not only to maintain the warm and hot phases and the stirring of the disk ISM, but also to maintain a huge hot gas corona in quasi-hydrostatic equilibrium with the cosmological dark halo (the intra-halo medium); the estimated sizes of these dark halos are several tens and hundreds of kiloparsecs for dwarf and giant galaxies, respectively. Note that this intra-halo medium should not be confused with the diffuse ionized gas and high-velocity-dispersion HI gas at scale heights of one or a few kiloparsecs above the disk plane, also frequently referred to as a “halo” in the ISM community (see e.g., Reynolds 1997; Kalberla & Kerp 1998; Mac Low 1999b).

The question is then whether the intra-halo gas can be dynamically heated by the turbulent energy input due to stellar winds produced by SNe and OB stars. Indeed, numerical simulations show that up to $`510`$% of the energy produced by SNe can be transformed into kinetic energy in the extreme case of multiple SN explosions forming supershells (e.g., Silich et al. 1996; Korpi et al. 1999); for isolated SNe, this fraction is lower than $`1`$%.

However, according to the results obtained here, most of kinetic energy in these “active” turbulent regions is dissipated locally and the $`u_{\mathrm{rms}}`$ of the “residual” turbulent motions quickly decays with distance (roughly $`u_{\mathrm{rms}}\mathrm{}^1`$). In order for the turbulent gas to be driven back into the intra-halo medium, the typical sizes of most of the forcing shells should exceed the gaseous disk height (supershells). This could be a common situation in dwarf starburst galaxies, but it is not the rule in normal disk galaxies. In fact, even in dwarf starburst galaxies, gas ejection from the disk seems to be not very efficient, according to recent hydrodynamical simulations by Mac Low & Ferrara (1999). Besides, numerical simulations have shown that the disk magnetic field is an efficient shield able to confine most superbubbles preventing immediate blowout (e.g., Slavin & Cox 1992; Franco et al. 1995; Tomisaka 1992, 1998). Therefore, we conclude that it is very unlikely that the turbulent kinetic energy injected by SNe and massive OB stars into the disk ISM can survive, propagate and drive back the gas into the intra-halo medium as the semi-analytic models of galaxy formation require.

We should note that our results are for densities typical of the ISM in the disk plane at the solar neighborhood. One may think that the disk stratification might reduce the dissipation rate so much that the turbulent kinetic energy could actually reach the intra-halo medium. However, since the HI disk scale height is larger than the typical dissipation lengths we have found, we believe any escaping energy will be too weak to provide support to the intra-halo medium. Simulations in a stratified medium will be presented elsewhere, with a code that does not require periodic boundary conditions.

An important question not treated in this paper is where the turbulent dissipated energy goes. It is well known that shocks efficiently transform kinetic energy into thermal energy. Thus, the heated gas at relatively large heights could be a source of high-energy photons able to ionize the low density extraplanar medium or even the hypothetical intra-halo medium. For example, Slavin (2000) has showed that hot gas in cooling SN remnants can be an important source of photoionization for the diffuse ionized gas. Dissipation in turbulent mixing layers has been also proposed as a mechanism to produce energetic radiation (Slavin, Shull, & Begelman 1993). Minter & Spangler (1997) have found that turbulent dissipation heating like that due to the ion-neutral collisional damping in the fluid-like turbulence model is able to provide a substantial contribution to the energy budget of the diffuse ionized gas. More studies on dissipation of turbulence and the thermo-hydrodynamics of the extraplanar and intra-halo medium are certainly necessary.

## 7 Concluding remarks

### 7.1 Summary

We have used 2D simulations of MHD, compressible, self-gravitating flows in the presence of parameterized heating and cooling and wind-like stellar kinetic energy input in order to study turbulence dissipation in the ISM. The “winds” are described by the typical radius $`l_\mathrm{f}`$ of their region of influence and the total kinetic energy $`e_\mathrm{k}`$ they inject into the medium. These parameters and the OB star formation rate $`\dot{\mathrm{\Sigma }}_{\mathrm{OB}}`$ completely characterize the forcing mechanism. It is important to emphasize that $`l_\mathrm{f}`$ is not the size ultimately reached by the expanding shells, but the scale over which the force is directly applied. The remainder of the shell expansion is inertial and strictly speaking does not constitute a forcing. Our main conclusions and implications are:

1. The non-ubiquitous, small-scale nature of the forcing implies that the forced regions have small filling factors, giving rise to the coexistence of both forced and decaying turbulent regimes within the same flow.

2. The kinetic energy injection and dissipation rates are always very close to each other, strongly suggesting that most of the dissipation occurs at or near the sources, where shocks are common, and that, for practical purposes, the injection and dissipation times are equal. In general, the flow dissipates its turbulent energy rapidly, in roughly 15–20 Myr.

3. In the forced regime, the global dissipation timescale $`t_\mathrm{d}`$ is empirically found to depend mainly on the forcing scale $`l_\mathrm{f}`$, and only very weakly on $`e_\mathrm{k}`$ (or $`\dot{ϵ}_\mathrm{k}`$, the specific power injected per source) and $`\dot{\mathrm{\Sigma }}_{\mathrm{OB}}`$.

4. Expressing $`t_{\mathrm{in}}`$ in terms of average quantities in the flow and parameters of the forcing, we showed that $`t_\mathrm{d}u_{\mathrm{rms}}^2/(\dot{ϵ}_\mathrm{k}f)`$, where $`f`$ is the filling factor of the forcing regions. In terms of measurable properties of the ISM, $`t_\mathrm{d}\mathrm{\Sigma }_\mathrm{g}u_{\mathrm{rms}}^2/(e_\mathrm{k}\dot{\mathrm{\Sigma }}_{\mathrm{OB}})`$, where $`\mathrm{\Sigma }_\mathrm{g}`$ is the average gas surface density. For the solar neighborhood, $`t_\mathrm{d}1.5\times 10^7`$ yr. Since $`\dot{\mathrm{\Sigma }}_{\mathrm{OB}}\dot{\mathrm{\Sigma }}_{}\mathrm{\Sigma }_\mathrm{g}^n`$, $`t_\mathrm{d}\mathrm{\Sigma }_\mathrm{g}^{n1}(R)`$, $`t_\mathrm{d}`$ is expected to vary with Galactocentric radius.

5. In the decaying regime, the turbulent kinetic energy decays as $`E_\mathrm{k}(t)(1+t/t_0)^\alpha `$, with $`\alpha 0.8`$. For this value of $`\alpha `$, the characteristic time for decay to half the initial energy is $`18`$ Myr, but this time depends on the total kinetic energy content of the flow, because of the power-law dependence. We associate the decaying case to the regime existing in regions distant from the energy sources, which contains the “residual” turbulence left over from the strong local dissipation operating in the neighborhoods of the sources.

6. From dimensional arguments, we suggest that the kinetic energy $`E_\mathrm{k}`$ and the rms velocity dispersion $`u_{\mathrm{rms}}`$ of the “residual” turbulence (far from the forcing regions) should decay with distance $`\mathrm{}`$ as $`\mathrm{}^{2\alpha /(2\alpha )}`$ and $`\mathrm{}^{\alpha /(2\alpha )}`$, respectively. For $`\alpha 1`$, $`E_\mathrm{k}\mathrm{}^2`$ and $`u_{\mathrm{rms}}\mathrm{}^1`$.

7. Our results imply that energy is dissipated much more efficiently at larger kinetic energy contents in the flow, in such a way that the dissipation time is rather insensitive to the total energy content.

8. If our results are applicable to the stratified vertical direction in the Galactic disk, then turbulent motions produced near the disk plane will propagate up to distances not too much smaller than the observed HI disk semi-thickness. This is consistent with models of galaxy evolution where the HI disk thickness is determined mainly by the turbulent kinetic energy content of the medium, which in turn results from the balance between kinetic energy injection by stars and its dissipation rate. However, some extra support is probably necessary to keep the HI and ionized layers up at the observed heights.

9. Our result of mostly local dissipation, again if applicable in the vertical direction, is in conflict with models of galaxy formation where the turbulent kinetic energy injected by SNe and stellar winds is assumed to reheat and drive back the gas from the disk into the cosmological dark matter halo (of sizes 15-20 times the disk size) in such a way that SF is self-regulated at the level of the whole intra-halo medium. Only in cases of non-stationary runaway SF (starbursts), most of the superbubbles might be able to blowout of the disk expelling large amounts of gas and energy into the intra-halo medium.

### 7.2 Importance of small-scale, intermittent forcing

The results of this paper imply that ISM turbulence is essentially different from ideal homogeneous turbulence due to the nature of the intervening energy injection (forcing) mechanisms. Although several authors (e.g., Scalo 1987; Fleck 1983; Norman & Ferrara 1996; see also the reviews in Vázquez-Semadeni 1999; Vázquez-Semadeni et al. 2000) have already discussed the fact that in “traditional” turbulence energy is injected at large scales while in the ISM it is injected over a range of scales (and preferentially at small ones), the consequences of the spatially- and temporally-discrete nature of interstellar energy injection have barely been discussed to our knowledge.

The effects of intermittent forcing have recently begun to be considered in other frameworks as well. Bec, Frisch & Khanin (2000) have studied the statistics and solvability properties of Burgers turbulence in the presence of “kick” forcing (i.e., the forcing is a series of delta functions in time), pointing out that such a flow combines features of both purely decaying and continuously forced cases. Léorat, Passot & Pouquet (1990) and Vázquez-Semadeni, Gazol & Scalo (2000) have shown that small-scale forcing is able to prevent the development of the gravitational and thermal instabilities, respectively. Finally, Kornreich & Scalo (2000) have considered the exposure of starless interstellar clouds to the intermittent passage of shock waves originated at distant injection sites, to show that the typical time between passages should be comparable to the decay time within the clouds, explaining their apparently turbulent state even in the absence of local energy injection sources.

In conclusion, the discrete, small-scale nature of interstellar kinetic energy injection seems to be responsible for a new, rich kind of turbulent flow, whose dissipation properties seem to pose clear constraints on models of galaxy formation and evolution.

This work has benefitted extensively from comments from J. Brasseur, A. Brandenburg, E. Ostriker and J. Stone while one of us (E.V.-S.) attended the “Astrophysical Turbulence Program” of the Institute of Theoretical Physics of the University of California at Santa Barbara, and from a thorough referee’s report by M.-M. Mac Low. This work was supported by CONACYT grant 27752-E to E.V.-S.  REFERENCES

Alfaro, E.J., Cabrera-Cano, J., & Delgado, A.J. 1991, ApJ, 378, 106

Avila-Reese, V. 1998, PhD. Thesis, U.N.A.M.

Avila-Reese, V., & Firmani, C. 2000, RevMexAA, 36, 23

Avillez, M. A. 1999, MNRAS, in press

Bania, T. M. & Lyon, J. G. 1980, ApJ, 239, 173

Bec, J., Frisch, U. & Khanin, K. 1999, J. Fluid Mech., submitted (chao-dyn/9910001)

Biskamp , D. 1994, Nonlinear Magnetohydrodynamics (Cambridge Univ. Press, Cambridge, England)

Chiang, W., & Bregman, J.N. 1988, ApJ, 328, 427

Chiang, W., & Prendergast, K. H. 1985, ApJ, 297, 507

Cole, S., Aragon-Salamanca, A., Frenk, C.S., Navarro, J., & Zepf, S. 1994, MNRAS, 271, 781

Cox, D. & Smith, 1974, ApJ, 189, L105

Dame, T. 1993, in Back to the Galaxy, eds. S. Holt, F. Verter, (AIP), p. 267

Dopita, M.A. 1990, in The Interstellar Medium in Galaxies, eds. H.A. Thronson, Jr. & J.M. Shull (Dordrecht: Kluwer), p. 437

Dopita, M.A., & Ryder, S.D. 1994, ApJ, 430, 163.

Elmegreen, B.G. 1991, in The Galactic ISM, eds. W.B. Burton, B.G. Elmegreen, R. Genzel (Springer-Verlag), p. 157

Ferrara, A. 1993, ApJ, 407, 157

Field, G.B., Goldsmith, D.W., & Habing, H.J. 1969, ApJ, 155, L149

Firmani, C., & Avila-Reese, V. 2000, MNRAS, 315, 457

Firmani, C., & Tutukov, A.V. 1992, A&A, 264, 37

Firmani, C., & Tutukov, A.V. 1994, A&A, 288, 713

Firmani, C., Hernández, X., & Gallagher, J. 1996, A&A, 308, 403

Fleck, R. C. Jr. 1983, ApJ, 272, L45

Franco, J., Santillán, A., & Martos, M. 1995, in The formation of the Milky Way, G.Tenorio-Tagle, M.Prieto, & Sánchez, F. eds. (Cambridge Univ. Press), p. 515

Friedli, D., & Benz, W. 1995, A&A, 301, 649

Frisch, U. 1995, Turbulence (Cambridge: Cambridge Univ. Press, Cambridge, England)

Galtier, S., Politano, H., & Pouquet, A. 1997, Phys.Rev.Lett, 79, 2807

Gerritsen, J.P.E. 1997, PhD. Thesis, Groningen University

Gerritsen, J.P.E., & Icke, V. 1997, A&A, 325, 972

Hossain, M., Gray, P., Pontius, D., Matthaeus, W., & Oughton, S. 1995, Phys.Fluids 7, 2886

Kadomtsev, B.B., & Petviashvili, V.I. 1973, Sov. Phys. Doklady, 18, 115

Kalberla, P.M.W., & Kerp, J. 1998, A&A, 339, 745

Kauffmann, G., White, S.D.M., & Guiderdoni, B. 1993, MNRAS, 264, 201

Kennicutt, R.C. 1998, ApJ, 498, 541

Kolmogorov, 1941, Dokl. Akad. Nauk 30,301

Kornreich, P. & Scalo, J. 2000, ApJ, 531, 366

Lèorat, J., Passot, T. & Pouquet, A. 1990, MNRAS, 243, 293

Lesieur, M. 1990, Turbulence in Fluids, 2nd ed. (Dordrecht:Kluwer)

Lockman, F.J, & Gehman, C.S. 1991, ApJ, 382, 182

Mac Low, M.-M. 1999a, ApJ, 524, 169

Mac Low, M.-M. 1999b, in New Perspectives on the Interstellar Medium, eds. A. R. Taylor and T. L. Landecker (ASP: San Francisco), in press

Mac Low, M.-M., & Ferrara, A. 1999, ApJ, 513, 142

Mac Low, M.-M., Klessen, R. S., & Burkert, A., & Smith, M.D. 1998, Phys. Rev. Lett., 80, 2754

McWilliams, J. C. 1990, Phys. Fluids A, 2, 547

McKee, C.F., & Ostriker, J.P. 1977, ApJ, 218, 148

Mihos, J.C., & Hernquist, L. 1994, ApJ, 437, 611

Mihos, J.C., & Hernquist, L. 1996, ApJ, 464, 641

Minter, A.H., & Spangler, S.R. 1997, ApJ, 485, 182

Myers, P.C. 1978, ApJ, 225, 380

Navarro, J.F., & White, S.D.M. 1993, MNRAS, 265, 271

Norman, C.A., & Ferrara, A. 1996, ApJ, 467, 280

Novikov, Y. A. 1978, Izv. Atmos. Ocean. Phys., 14, 474

Padoan, P., & Nordlund, A. 1999, ApJ, 526, 279

Passot, T., Vázquez-Semadeni, E., & Pouquet A. 1995, ApJ, 455, 533

Passot, T. & Pouquet, A. 1988, J. Comput. Phys. 75, 300

Reynolds, R.J. 1997, in Proc. of 156 WE-Heraeus-Seminar on The Physics of Galactic Halos, eds. H. Lesch et al. (Berlin: Akademie Verlag), p. 57

Rosen, A., & Bregman, J.N. 1995, ApJ, 440, 634

Santangelo, P., Benzi, R. & Legras, B. 1989, Phys. Fluids A, 1, 1027

Scalo, J.M. 1987, in Interstellar processes, ed. D.J.Hollenbach & H.A.Thronson (Dordrecht: Reidel), 349

Scalo, J.M., & Struck-Marcell, C. 1984, ApJ, 276, 60

Silich, S.A. Franco, J., Palous, J., & Tenorio-Tagle 1996, ApJ, 468, 722

Slavin, J.D., & Cox, D.P. 1992, ApJ, 392, 131

Slavin, J.D. 2000, in Astrophysical plasmas: codes, models, and observations, eds. J. Arthur, N. Brickhouse, & J. Franco, RevMexAA (Serie de Conferencias), v. 9, 246

Slavin, J.D., Shull, J.M., & Begelman, M.C. 1993, ApJ, 407, 83

Somerville, R.S., & Primack, J.R. 1999, MNRAS, 310, 1087

Stone, J.M., Ostriker, E.O., & Gammie, C.F. 1998, ApJ, 508, L99

Struck, C., & Smith, D.C. 1999, ApJ, 527, 673

Tammann, G.A., Loffler, W., & Schroeder, A. 1994, ApJS, 92, 487

Tomisaka, K. 1992, PASJ, 44, 177

Tomisaka, K. 1998, MNRAS, 298, 797

van den Bosch, F.C. 2000, ApJ, 530, 177

Vázquez, E.C. & Scalo, J.M. 1989, ApJ, 343, 644

Vázquez-Semadeni, E. 1999, in Millimeter & Submillimeter Astronomy: Chemistry and Physics in Molecular Clouds, Proc. of the 1996 INAOE Summer School of Millimeter-Wave Astronomy, eds. W. F. Wall, A. Carramiñana, L. Carrasco, and P. F. Goldsmith (Dordrecht: Kluwer), 161

Vázquez-Semadeni, E., Gazol, A. & Scalo, J. 2000, ApJ, 540, 271

Vázquez-Semadeni, E., Passot, T., & Pouquet A. 1995, ApJ, 441, 702

Vázquez-Semadeni, E., Passot, T., & Pouquet A. 1996, ApJ, 473, 881

Vázquez-Semadeni, Ostriker, E.C., Passot, T., Gammie, C.F., & Stone, J.M. 2000, in Protostars and Planets IV, eds. V.Mannings, A.Boss, S.Russell (Tucson: Univ. of Arizona Press), p. 3

Wang, B., & Silk, J. 1994, ApJ, 427, 759

White, S.D.M, & Frenk, C.S. 1991, ApJ, 379, 52


format: text

text:
# 1 Introduction

## 1 Introduction

The effects of hadron density and/or temperature on the mass and lifetimes of vector mesons have been the subject of numerous investigations, experimental as well as theoretical . Under conditions of chiral symmetry restoration, the quark condensate vanishes, $`<\overline{q}q>0`$. The condensate is not, however, an observable but it is related to the meson mass, which is. Thus, the behaviour of vector meson masses under high hadron densities and/or temparatures has become a testbed of our understanding of hadron or quark dynamics in nuclear matter. A renormalization of vector meson masses in hot and dense matter is of significant importance to the pursuit of quark gluon plasma research, has relevance to the equation of state for nuclear matter in supernovae explosions, and could affect lower energy nuclear physics, as well. It has been argued that the mass modification of vector mesons should be measurable even at normal nuclear matter densities .

Even though a reduction of vector meson mass in hadronic matter is supported by most of the more recent theoretical models, the spectral shape (manifested in the width) is more of an open question. If there is substantial broadening as reported in , then extracting meaningful mass change information from $`\rho `$ production and decay experiments becomes problematic, especially so for massive nuclei. In this respect, production of either the $`\omega `$ in massive nuclei or the $`\rho `$ from light nuclei may be a better solution. In any case, production of vector mesons, by whatever method, is not enough to probe medium modifications. The production of mesons must be accompanied by a substantial decay fraction within the nuclear volume, a fact which favours either low beam energies or restricted phase space in which the vector mesons are essentially at rest with respect to the recoiling nuclear system.

## 2 The TAGX Results from $`𝐀(\gamma _𝐭,\pi ^+\pi ^{})𝐗`$ Reactions

### 2.1 The Subtreshold Energy Regime

Photo- or lepto-production of vector mesons from nuclei has been and is being explored at many facilities. Given the quantum numbers of the vector mesons and those of the photon (real or virtual), the primary interaction of the latter proceeds via the interaction of the former with the nucleons in nuclei, the well established VMD model of $`\gamma `$-N interactions. At incident photon energies of approximately 1.2 GeV and higher, the diffractive mechanism dominates, in which the photon flactuates into a $`\rho ^0`$ or $`\omega `$, which then scatters off the nucleon field with sufficient four-momentum transfer to bring the meson on the mass shell. In the case of nuclei, the single $`\gamma `$-N amplitudes add coherently and large production cross sections result. Since the momentum transfer $`t`$ is less than 0.1 $`\mathrm{GeV}^2/\mathrm{c}^2`$, coherent or diffractive cross sections are very forward peaked. As a result, they also conserve helicity and the reaction is characterized by a large impact parameter of approximately 10 fm or more.

While coherent photoproduction produces vector mesons prolifically, the probability of the vector mesons penetrating the nuclear interior and decaying within it is rather small. The observation, then, of medium modifications in such reactions is diminished by “geometry” (including trajectory of the produced meson with respect to the nucleus) and the large Lorentz boost of the meson with respect to the nucleus. This is particularly severe for the $`\omega `$ due to its long $`c\tau `$ compared to that of the $`\rho ^0`$. One is thus forced to depend on the large number of mesons produced to observe the small fraction that penetrates and decays within the nuclear volume, and which is perharps overwhelmed in the invariant mass signature from the low tail of the unmodified mesons decaying freely. Where the fraction is larger, like in the $`\rho `$-meson case, the free $`\rho ^0`$ width is also large, thus effectively masking the modified mass. In the case of the $`\omega `$, the narrow width makes it easier to separate the modified from the unmodified masses, but the fraction of the former to the latter is very small.

An alternative is to produce the $`\rho ^0`$ below the free $`\gamma +\mathrm{N}\rho ^0`$ nominal threshold of 1083 MeV. In this case, the reaction is characterized by large $`t`$ and the Fermi four-momentum of the struck nucleon is utilized to bring the $`\rho ^0`$ on the mass shell. A penalty one pays for such subthreshold $`\rho ^0`$ production is the small cross section compared to coherent production. There are numerous advantages, however. First, the $`\rho ^0`$ production vertex is forced to be within the nuclear volume (small impact parameter) by the interaction of the photon with a bound nucleon. The reaction then is quasi-free $`\gamma `$-N interaction in nature. Second, the $`\rho ^0`$ is produced with small relative momentum and Lorentz boost with respect to the struck nucleon or the nucleus, thus maximizing the probability of decay within the hadronic volume. Third, the large opening angle (of $`180^\mathrm{o}`$ in the $`\rho ^0`$ rest frame) is essentially preserved in the lab frame due to the small Lorentz boost and this provides a characteristic signature of the $`\rho ^0`$ decay, which can be used effectively in the suppression of competing processes.

### 2.2 The TAGX Results on $`{}_{}{}^{\mathrm{𝟑}}\mathrm{𝐇𝐞}`$

The choice of $`{}_{}{}^{3}He`$ was made based on several reasons. (a): Its core nuclear density is near saturation density for heavier nuclei, while the small number of nucleons reduces final state interactions (FSI) between the pions and the recoiling nucleus or nucleons. (b): The kinematical analysis is much simpler due to the small number of spectator nucleons and the Monte Carlo (MC) simulations carry more confidence. And, (c): extrapolations of cross sections for background processes on proton and $`{}_{}{}^{2}H`$ are more accurate for $`{}_{}{}^{3}He`$ than more complex nuclei. Thus, it was felt that any mass modification, even if small, could carry higher confidence due to the simplicity of the probe-target system.

The analysis was based on a detailed and kinematically constrained fit to several distributions simultaneously for all processes that can lead to a $`\pi ^+\pi ^{}`$ final state. Details are published in , , and . The results were consistent with a large mass modification for the $`\rho ^0`$, which was also found to depend on incident photon energy. At photon energies as low as $`660\pm 40`$ MeV, the best fit was obtained for a $`m_{\rho ^0}^{}=490\pm 40`$ $`\mathrm{MeV}/\mathrm{c}^2`$, while for $`840\pm 40`$ MeV the mass of the $`\rho ^0`$ was $`m_{\rho ^0}^{}=640\pm 40`$ $`\mathrm{MeV}/\mathrm{c}^2`$. For the highest energy bin available to the experiment (which was centered around 1080 MeV), the data were shown to be insensitive to any $`\rho ^0`$ mass modification. The conclusions survived a number of different analyses, and even the data sets for two of the analyses were different and independent of each other.

The analysis was based on the simultaneous fitting of a number of “key” variables, such as missing momentum ($`p_{miss}`$), missing mass ($`m_{miss}`$), invariant mass ($`m_{\pi \pi }`$), pion opening angle ($`\theta _{\pi \pi }`$), and emission angle of the di-pion system ($`\theta _{IM}`$). All these are kinematical observables and follow phase space distributions, applicable to each individual reaction assumed. For the 80 MeV wide photon energy bins in the analysis, the assumption was made that the matrix element is constant. The yield then is the product of several factors:

$$Y=|M^{}(s,t)|^2\mathrm{\Phi }_{{}_{}{}^{3}He}(p)\mathrm{\Psi }(s,I,J)[MBPSF(m,\mathrm{\Gamma })]$$
(1)
where Y = Yield, $`|M^{}|`$ = invariant matrix element for the reaction, $`\mathrm{\Phi }`$ = the single nucleon momentum distribution, $`\mathrm{\Psi }`$ = function of spin, isospin, and total angular momentum, and $`[MBPSF]`$ = the multi-body phase space for this reaction with the mass and width of any resonances folded in. Therefore, even though a variation of $`|M^{}|`$ within the 80 MeV bin is not taken explicitly into account, the width $`\mathrm{\Gamma }`$ of the states involved has a large part of this energy dependence folded in. Each background and foreground reaction was simulated using Monte Carlo (MC) generators, which incorporated equation 1.

Cross sections for known processes, such as $`\mathrm{\Delta }\pi ,N^{}\pi ,\mathrm{\Delta }\mathrm{\Delta }`$ and $`\pi ^+\pi ^{}\pi ^0`$, were extracted and compared to established measurements in the literature, if such existed. The cross sections for $`\rho ^0`$ production with modified masses were quite small (a few $`\mu `$b), while that of competing processes are of the order of 200 $`\mu `$b . The small $`\rho ^0`$ signal could be extracted due to the features of the TAGX acceptance. The trigger requirement of left-right (with respect to the photon beam) $`\pi ^+\pi ^{}`$ coincidences, and the small out of plane acceptance of the spectrometer, result in substantially larger acceptance for the $`\rho ^0\pi ^+\pi ^{}`$ decay than the other competing two-step processes leading to two-pion emission. Furthermore, the requirement of simultaneous fitting of so many different kinematical variables provided severe contraints to the strength of these background reactions in different regions of the various distributions.

One of the difficulties of the above analysis is the model dependency and the smallness of the $`\rho ^0`$ production cross sections compared to all the other competing processes. The analysis, however, has been shown to be insensitive to reasonable variations in assumed cross sections for the dominant background processes because the $`\rho ^0`$ events occupy quite distinct regions of phase space. Thus, variations in $`\mathrm{\Delta }\pi `$ production cross sections, for example, affected $`\mathrm{\Delta }\mathrm{\Delta }`$ fits and resulting cross sections, but they did not affect the value of $`m_{\rho ^0}^{}`$. Another conclusion was that $`m_{\pi \pi }`$ is not a very sensitive parameter and good fits can be obtained for large differences in the input parameters. On the other hand, $`\theta _{\pi \pi }`$ and $`p_{miss}`$ are difficult (and sensitive) distributions to fit.

In the case of coherent $`\rho ^0`$ production and at very small angles of emission in the helicity frame of reference, the angular distribution for $`\theta _\pi ^{}`$ of one pion with respect to the di-pion center of mass momentum in the helicity frame, has a $`\mathrm{sin}\theta _\pi ^{}`$ -like distribution for unpolarized photons. This is a result of $`s`$-channel helicity conservation. In the case of helicity non-conservation, for example at large angles of $`\rho ^0`$ emission due to large $`t`$ exchange, a random distribution of polarization would result in an isotropic distribution for $`\theta _\pi ^{}`$. When, however, such a $`\theta _{\pi ^+}^{}`$ distribution was formed, an enhancement in the regions of $`0^o`$ and $`180^o`$ was observed in the “raw” event sample. This has been interpreted as a substantial polarization component in the $`\rho ^0`$ event sample which corresponds to the $`l=1,m=0`$ substate. Such a $`Y_1^0`$ -like distribution was most prominent in the same $`m_{\pi \pi }`$ bins where the independent kinematical analysis placed the mean of the modified $`\rho ^0`$ mass.

### 2.3 The TAGX Results on $`{}_{}{}^{\mathrm{𝟐}}𝐇,^{\mathrm{𝟏𝟐}}𝐂`$

The analysis of the $`{}_{}{}^{2}\mathrm{H}`$ and $`{}_{}{}^{12}\mathrm{C}`$ results is still ongoing and a final value for the modified masses has not been extracted yet. The philosophy and methodology is quite different than the analysis of the $`{}_{}{}^{3}\mathrm{He}`$ data. The simpler $`{}_{}{}^{2}\mathrm{H}`$ nucleus, coupled with photoabsorption cross sections in the literature for most of the channels leading to di-pion emission, allow the use of experimental amplitudes and cross sections and a likelyhood method of fitting the kinematical observables. Unlike the analysis of the $`{}_{}{}^{3}\mathrm{He}`$ data, in this analysis the $`\theta _{\pi ^+}^{}`$ distributions are explicitly involved among the kinematical variables fitted. The cross section results are not yet final, however, the analysis based on the $`\theta _{\pi ^+}^{}`$ distribution is quite interesting. In Figure 1, the $`\mathrm{cos}\theta _{\pi ^+}^{}`$ distributions are shown for $`{}_{}{}^{12}\mathrm{C}+^2\mathrm{H}`$ data, (a) without any restrictions on the data, and (b) with a $`\theta _{\pi \pi }120^o`$ cut. Several observations can be drawn:

* In the $`M_{\pi \pi }`$ bin of 300-400 $`\mathrm{MeV}/\mathrm{c}^2`$, the shape of the distribution and the small number of events near the -1 and +1 limits is characteristic of acceptance effects due to trigger and geometry of TAGX. As such, the population in these two regions for the invariant mass bins in the 500-600 $`\mathrm{MeV}/\mathrm{c}^2`$ regions, signify substantial $`\rho ^0\pi ^+\pi ^{}`$ -like production.
* The application of the opening angle ($`\theta _{\pi \pi }`$) cut, removes events form the central regions of the distribution, but leaves the population around the two limits unchanged.

The $`\mathrm{cos}\theta _{\pi ^+}^{}`$ distributions for all three nuclei exhibite the same overall behaviour. The preliminary results from the cross section fitting of the $`{}_{}{}^{2}\mathrm{H}`$ data are consistent with a few $`\mu `$b cross section for polarized $`\rho ^0`$ production in the $`l=1,m=0`$ substate and no strength in the $`l=1,m=1`$ substate. This is qualitatively in agreement with the $`{}_{}{}^{3}\mathrm{He}`$ data which also showed no statistically significant population in the small angle regions in plots such as shown in Figure 1. This was verified by successively applying increasingly tighter opening angle cuts and comparing the behaviour as a function of $`\mathrm{cos}\theta _{\pi ^+}^{}`$ values.

### 2.4 The Effects of $`\theta _{\pi \pi }`$ Cuts

The $`l=1,m=0`$ signature is model independent and difficult to argue away. Nevertheless, there have been a number of possible explanations, other than a modified and polarized $`\rho ^0`$. One is that there may be other reactions which lead to two-pion emission in a relative $`p`$-state and with a $`m=0`$ substate. However, no such specific reaction has ever been identified or proposed as a candidate; in any case, if there is indeed one, this by itself is a new and interesting result. Another concern put forward is that the application of ever tighter $`\theta _{\pi \pi }`$ requirements may artificially produce such $`l=1,m=0`$ signatures. In order to address this valid concern, a detailed investigation of the behaviour of two of the most prolific background processes has been made as a function of $`\theta _{\pi \pi }`$ cuts. MC simulations of three production reactions have been made and analyzed in the same fashion as the real data are analyzed in the TAGX detector. These are shown in Figure 2 for $`{}_{}{}^{2}\mathrm{H}`$, plotted as a function of $`\mathrm{cos}\theta _H`$ in the helicity frame.

The $`\mathrm{\Delta }\pi `$ and $`\mathrm{\Delta }\mathrm{\Delta }`$ background channels are the two strongest channels leading to $`\pi ^+\pi ^{}`$ production. As can be seen in Figure 2, whether the $`\mathrm{\Delta }`$ is produced polarized (indicated by the arrow) or not, the application of successively more demanding $`\theta _{\pi \pi }`$ cuts does not alter the shape of the distributions. In the case of unpolarized $`\rho ^0`$ production, the distribution is essentially flat with a slight enhancement at the limits due to detector acceptance. Only in the case of polarized $`\rho ^0`$ production do the MC simulations in Figure 2 agree with the observed behaviour of the data in Figure 1.

## 3 Conclusions

The analyses of $`{}_{}{}^{3}\mathrm{He}`$ data using MC and model dependent, as well as model independent variables, agree on a substantial mass modification of the $`\rho ^0`$. The preliminary results on $`{}_{}{}^{2}\mathrm{H}`$ and $`{}_{}{}^{12}\mathrm{C}`$, using primarily the model independent $`\theta _{\pi ^+}^{}`$ analysis, agree with the earlier results and strongly support the production of longitudinally polarized (zero helicity) and mass modified $`\rho ^0`$ -mesons. The similarities among the results for these three different nuclei is suggestive of a $`\rho ^0`$-N based interaction effect rather than a nuclear effect. Since the relative $`\rho ^0`$-N momenta at such low incident photon energies are also very low, the mean decay length of the $`\rho ^0`$ from its production vertex is less than the nucleon radius. Hadronic densities at the $`interior`$ of a nucleon are very high, thus qualitatively accounting for the large mass modifications reported so far from TAGX results. The excellent agreement on the $`\theta _{\pi \pi }`$ dependence for data and MC simulations conclusively shows that only a $`l=1,m=0`$ substate of the $`\rho ^0`$ is present in the data. Such a production can be accounted for by a spin-flip transition for the target nucleon, consistent with a $`N^{}`$ excitation in the $`\rho ^0`$ production process.