Turbulence in the Interstellar Medium

Turbulence is ubiquitous in the insterstellar medium and plays a major role in several processes such as the formation of dense structures and stars, the stability of molecular clouds, the amplification of magnetic fields, and the re-acceleration and diffusion of cosmic rays. Despite its importance, interstellar turbulence, alike turbulence in general, is far from being fully understood. In this review we present the basics of turbulence physics, focusing on the statistics of its structure and energy cascade. We explore the physics of compressible and incompressible turbulent flows, as well as magnetized cases. The most relevant observational techniques that provide quantitative insights of interstellar turbulence are also presented. We also discuss the main difficulties in developing a three-dimensional view of interstellar turbulence from these observations. Finally, we briefly present what could be the the main sources of turbulence in the interstellar medium.


Introduction: the basics of turbulence
Turbulence is characterized by chaotic motions in a fluid (Rempel et al., 2004;He & Chian, 2005;Chian et al., 2006Chian et al., , 2007Chian et al., , 2010) that lead to diffusion of matter and dissipation of kinetic energy. It is to be stressed that not all chaotic motions in a fluid may be called "turbulent". Because of its chaotic nature turbulence can only be studied and modelled in terms of statistical quantities. Long-term deterministic local properties of a turbulent fluid are unpredictable.
For nearly incompressible and unmagnetized fluids, the temporal evolution of the fluid velocity field is given by the Navier-Stokes equation: ∂u (x, t) ∂t +u(x, t)·∇u(x, t) = − ∇p(x, t) ρ(x, t) +ν∇ 2 u(x, t)+F (x, t), (1) where u(x, t) represents the velocity field, p the pressure, ν the kinematic viscosity, and F an external force normalized by the local density. ρ is the gas mass density and is set constant in the incompressible case (with ∇ · u = 0). Even in this simplified mathematical description the fluid dynamics is not a trivial solution. Equation 1 is non-linear, as seen from the advective term in the left hand side, and non-local -in the sense that the local properties of the fluid are related to all the other regions -, through the pressure term. The incompressibility condition results in an infinite sound speed, and in an instantaneous propagation of any perturbation in the fluid. Burgers (1939) modeled the time evolution of the simplified version of the Navier-Stokes equation by considering ∇p = 0. This equation has exact solutions, which may sound interesting, but it results in non-universal "turbulence". Even-though Burgers turbulence models have gained increasing interest due to their ability to describe the statistics of shock induced structures, and many other applications (see review by Bec & Khanin, 2007).
In the full Navier-Stokes equation, perturbations in u(x, t) are expected to have their distribution changed due to nonlinear terms. These instabilities may drive local vorticity and result in the fragmentation of large amplitude eddies into smaller ones, creating a turbulent pattern. As imagined by Richardson (1922), big whirls have little whirls that feed on their velocity, and little whirls have lesser whirls, and so on to viscosity. This statement represents one of the first conceptual descriptions of the energy cascade in turbulent flows. The shear drives unstable motions at large scales, which are broken and fragmented into smaller vortices, down to the smallest scales where they are damped, e.g. due to viscosity. In an incompressible viscous fluid this damping scale is that at which the timescale for viscous damping is of the order of the turnover dynamical time. At that scale, the eddy kinetic energy is transferred to internal energy due to viscosity. Turbulence is naturally developed over larger range of scales if viscosity is small, i.e. with large Reynolds number (Re = UL/ν ≫ 1), being the characteristic velocity U injected at a lengthscale L. Kolmogorov (1941, hereafter K41) realized that it would be possible to solve the Navier-Stokes equation for a turbulent flow if u(x, t) is considered a stochastic distribution. One of the key assumptions in the K41 theory is that the energy transfer rate ǫ should be constant at all scale. It is defined as ǫ ≃ δu 2 l /τ l , where δu l is the velocity fluctuation amplitude at lengthscale l, and τ l = τ eddy = l/δu l its dynamical timescale 1 . Therefore, one obtains: δu l ≃ (ǫl) 1 3 .
(2) Equation 2 means that turbulence can be modeled by scaling laws. This would be true within the so called inertial range of scales, i.e. the scales where the energy transfer rate is constant, generally between the energy injection and the viscous damping scales. The velocity power spectrum P u (k) is defined 2 here by ∞ k=1/l P u (k ′ )dk ′ = δu 2 l , from which we obtain the standard Kolmogorov power spectrum for the velocity field: In other words, it is possible to reinterpret Kolmogorov's idea in Fourier space in terms of non-linear interaction between similar wavenumbers. This theory is a result of the so-called locality, i.e. similar wavenumbers, k = 2π/λ, of the non-linear wave-wave interaction that result in the energy cascade through smaller scales (Kraichnan, 1965a). From the spectral form of the Navier-Stokes equation, the three-wave interactions follow the selection rule k 3 = k 1 +k 2 . The extrema are found at k 3 → 0 and k 1 = k 2 , which is the locality assumed in Kolmogorov's theory, resulting in k 3 = 2k 1 .
The theory also predicts the scaling laws for the moments of velocity spatial increments, known as velocity structure functions, defined as: where p is a positive integer representing the moment order and l is the spatial increment vector. In incompressible fluids, if the turbulence is considered homogeneous, isotropic and self-similar, i.e. scale invariant, then: where C(p) was initially assumed by Kolmogorov to be constant with p.
One of the main successes of the Kolmogorov-Obukhov turbulence theory is the explanation of the empirical determination of the diffusion coefficient by Richardson (1926), done more than a decade before K41. The diffusion coefficient is related to the time evolution of the separation between Lagrangian points (e.g. particles dragged by the flow) in a turbulent medium. The probability distribution function Φ of pairs of points separated by a distance r may be described as: where K(r) represents the diffusion coefficient. It is easy to determine, from dimensional analysis, that ifṙ = u(r) ∝ r 1/3 as in the Kolmogorov scaling, the diffusion coefficient for the inertial range will be K(r) = k 0 ǫ 1/3 r 4/3 , the scaling proposed by Richardson (1926). This diffusion coefficient for the inertial range substituted in Equation 6 then results in: where A is a normalization coefficient. The Richardson distribution is therefore non-Gaussian. Several experiments and numerical models have shown the validity of the turbulent diffusion scaling (Elliott & Majda, 1996;Fung & Vassilicos, 1998;Zouari & Babiano, 1994;Boffetta & Sokolov, 2002), as has also been recently used in the predictions of stochastic magnetic reconnection 3 (Lazarian, Eynk & Vishniac, 2012). This theory of turbulence has been quite successful in reproducing most of experimental data, and there is a flourishing literature with hundreds of works available e.g. Armstrong, Rickett & Spangler (1995); Leamon et al. (1998); Bale et al. (2005); Koga et al. (2007); Bourras et al. (2009);Chian & Miranda (2009) ;Chepurnov & Lazarian (2010); Sahraoui et al. (2010); Chian & Muñoz (2011);Chang, Bewley & Bodenschatz (2012); Hurricane et al. (2012);Miranda, Chian & Rempel (2013), just to mention a few. Naturally, many authors criticized the fact of C(p) is a constant in Kolmogorov's initial theory, given the breakdown of self-similarity at small scales and the possible non-universality of turbulence (given its "memory" related to the energy injection). These criticisms have been later addressed in the Kolmogorov-Obukhov turbulence theory (Kolmogorov, 1962;Obukhov, 1962), including the effects of intermittency. Intermittency results from rare and large local fluctuations in the velocity field which break the similarity condition (Frisch, 1995). One of the effects of intermittency is observed in the probability distribution function (PDF) of velocity longitudinal increments δu l = [u(r + l) − u(r)] ·l, which shows large deviations from the Gaussian distribution at small scales, with large amplitude tails and peaked distributions at δu l ∼ 0 (see Figure 1). Kraichnan (1991) pointed that sharp shocks could, for intance, result in more regions with smooth fluid flows and also more regions with sharp transitions in velocities, compared to the standard picture of the self-similar K41 turbulence. We would then expect non-Gaussian PDFs at both small and large scales.
Many authors attempted to theoretically determine the scalings of turbulence with intermittency. One of the most successful approaches is the multifractal description for the energy dissipation field proposed by She & Leveque (1994). This theory results in S p (l) ∝ l ζ(p) , with: where D ′ represents the dimensionality of the dissipation structures. In the Kolmogorov-Obukhov theory, structures of highest dissipation are filamentary, better described then by D ′ ∼ 1, while recent numerical simulations reveal a dominance of two dimensional intermittent structures at small scales (e.g. Moisy & Jimenez, 2004;Kowal, Lazarian & Beresnyak, 2007;Boldyrev & Perez, 2012), what is also supported 3 this term accounts for the magnetic reconnection that is induced by turbulent motions near the current sheet -separation layer between fields with components of opposed directions -, which would then result in reconnection rates as a function of the stochastic motions of the fluid. by experimental data (e.g. Fredriksen et al., 2003;Thess & Zikanov, 2007). Multifractal analysis of Voyager 1 and 2 in situ data have also showed intermittent features on the magnetic turbulence at the solar wind and the termination shock (Macek & Szczepaniak, 2008;Macek, Wawrzaszek & Carbone, 2011. On the theoretical side, Birnir (2013) derived a statistical solution of the stochastic Navier-Stokes equation from the linear Kolmogorov-Hopf differential equation, accounting for the She-Levêque intermittency corrections. His results satisfactorily reproduce the PDFs built on observations and numerical simulations of turbulent flows. Compressibility and coupling between magnetic fields and the plasma flow -both present in the dynamics of the interstellar medium (ISM) -make the description of the interstellar turbulence even more complex.

Supersonic turbulence
Compressible plasmas are of great interest in astrophysics, and particularly in the case of interstellar turbulence. Compressibility in turbulent flows results in the formation of a hierarchy of density structures, viewed as dense cores nested in less dense regions, which are in turn embedded in low density regions and so on. Such a hierarchical structure was described by von Weizsäcker (1951) as: where ρ ν represents the average density of a structure at hierarchical level ν, at a lengthscale l, and α the compressibility degree, assumed to be the same at each level. The dimensionality of the system is obtained by D ′ = 3 −3α. Therefore, www.jn.net J. Name the average mass within each substructure must follow the relation M l ∝ l 3−3α . The density hierarchy as described above must then be coupled to the local turbulent motions. The energy density transfer rate must now be rewritten as ǫ l = ρ l δu 3 l /l to account for the density changes at different scales (Lighthill, 1955). If, once again, one assumes the constancy of the energy transfer rate across scales within the inertial range (Fleck, 1996), one obtains the scaling of the amplitude of the velocity fluctuations: and the velocity power spectrum is then given by: Note that for stationary energy distribution solutions in compressible turbulence, α > 0 which results in steeper velocity power spectra, compared to the standard K41 scaling. The density power spectrum, on the other hand, instead of following the velocity field as a passive scalar would do, presents a distinct power spectrum given by: i.e. for α ∼ 1/6, the power spectrum of the density field becomes flat in the inertial range. One of the most striking results of the hierarchical model for the density field in compressible turbulence is its ability to recover the standard Kolmogorov scalings for the density weighted velocity field v ≡ ρ 1/3 u (Fleck, 1996). Numerical simulations of compressible turbulence have confirmed the scalings described above for α ≃ 0.15 (Kritsuk et al., 2007;, close to α = 1/6 for which the density power spectrum becomes flat. The velocity power spectrum on the other hand becomes P u (k) ∝ k −2 . Remarkably, this is the exact slope obtained for Burger's turbulence, despite the different framework of that theory.

Magnetized turbulence
Magnetic fields introduce further complexity in the plasma dynamics that can be described by the magnetohydrodynamic (MHD) equations in the fluid approximation and assuming perfect coupling between the field and the plasma: where B is the magnetic field and η the plasma resistivity (η = 0 for ideal plasmas). Let us first consider an external uniform magnetic field B 0 . Any perturbation in the fluid velocity field will be coupled to the magnetic field. The magnetic tension/pressure results in a decrease of the non-linear growth of perturbations, but only of those perpendicular to the magnetic field lines. This complex coupling between the flow and magnetic field makes the modelling of turbulence in magnetized plasmas an interesting task 4 .

The Iroshnikov-Kraichnan model
An useful simplification to the equations above is made by considering B = B 0 +δB, and using the Elsässer variable z ± = u ± δB, whereB = B/(4πρ) 1/2 . This has been independently derived by Iroshnikov (1963) and Kraichnan (1965a,b) (IK hereafter). From this change of variables, Eqs 13 and 14 result in (see Schekochihin & Cowley, 2007): where v A = B 0 / 4πρ is the Alfvén velocity and ∇ || is the spatial derivative parallel to the direction of the mean magnetic field.
In their model, Iroshnikov and Kraichnan proposed that incompressible magnetized turbulence results from the nonlinear interactions of counter propagating waves packets. The timescale for the two wave packets to cross each other is of order of the Alfvén time τ A ∼ l || /v A , where l || is the lengthscale of the wave packet parallel to the mean magnetic field. In their phenomenological description of the MHD turbulence, the interactions between the wave packets are weak, i.e. |z ± | ≪B 0 or the field perturbations are much smaller than B 0 . Notice that, superimposed to the magnetic fluctuations, the fluid is also perturbed and the dynamical timescale of a fluid "eddy" is τ eddy ≡ l/δu l . The different wave modes (mechanical and magnetic perturbations) thus interact with each other. For the interaction between modes to be weak the Alfvén time must be much smaller than the dynamical timescale, i.e. τ A ≪ τ eddy . The non-linear decay of the wave packets in such weak interactions, and subsequently the turbulent cascade, can only occur after several interactions. Since interactions are random, the wave packet amplitude changes in a random walk fashion, i.e. N = (τ l /τ A ) 1/2 interactions are needed for the wave packet to significantly change. At the same time, N is also defined by the number of crossings in a decay timescale N = τ l /τ eddy , which results in: 4 More details on MHD turbulence may be found in Biskamp (2003) Therefore, the turnover time at scale l is longer by a large factor and, as expected, the non-linear cascade proceeds much more slowly.
The second major assumption in the IK theory of weak turbulence is its isotropy, i.e. l || ∼ l. Substituting this scaling into the relation ǫ = δu 2 l /τ l , one obtains: There is evidence for an IK cascade in the solar wind and interplanetary medium (e.g. Bamert et al., 2008;Ng et al., 2010). However, many observations of the solar wind turbulence also suggest a more Kolmogorov-like turbulence, i.e. ∝ k −5/3 (e.g. the early studies of Coleman, 1968;Matthaeus & Goldstein, 1982; and the more recent papers by Alexandrova, Lacombe & Mangeney, 2008;Chian & Miranda, 2009;Sahraoui et al., 2010;Li et al., 2011;Chian & Muñoz, 2011;Kozak et al., 2012;Hellinger et al., 2013). It is possible though that a mix of both cascades may occur, as pointed by e.g. Salem et al. (2009) and Alexakis (2013), which showed a mix of K41 and IK cascades for the magnetic and velocity field fluctuations, respectively. Moreover, most of these data also reveal the solar wind turbulence to be highly anisotropic (i.e. δu || l δu ⊥ l ) with respect to the local magnetic field (Horbury et al., 2008(Horbury et al., , 2012. As pointed by Goldstein (2001), one of the main issues raised by the solar wind is why is the power spectrum of this anisotropic, compressible, magnetofluid often Kolmogorovlike? 1.2.2 The Goldreich-Sridhar model Marsch (1990) remarked that if, instead of an Alfvén time, the timescale for the waves to non-linearly interact with each other was the regular eddy turnover time, i.e. τ l ≃ τ eddy ∼ l || /δu l , one would get a K41 cascade for the magnetized turbulence. This would be true also for the case of strong turbulence, |z ± | >B 0 . The isotropy condition was retained, which was raising a problem, most of the observational data mentioned above revealing strongly anisotropic turbulence.
Goldreich & Sridhar (1995, GS95 hereafter) proposed a turbulent model based on anisotropic fluctuations, with strong coupling between the wave modes. Strictly speaking the GS95 model assumes a critical balance between mechanical and Alfvénic modes in such a way that l ⊥ /δu l ≃ l || /v A . Therefore: From Eq. 18, not only the magnetized turbulence is anisotropic but it is local in the sense that the anisotropy is measured in the reference frame of the local magnetic field. Such an anisotropy is expected to occur in both the dispersion of velocity (δv) and wave vectors k, though it is easier to observe velocity dispersion anisotropies from the interstellar medium, as discussed below. Therefore, statistically, a large number of eddies with local fields randomly distributed in space result in an average zero anisotropy (even at small scales). In the strong magnetized cases though, the anisotropy would be more clearly detected in experiments and observations. Several direct numerical simulations of magnetized turbulence in a quasi-incompressible regime have been performed in the past decade. Many numerical experiments reveal that MHD turbulence indeed has a large part of its energy cascade close to a K41 distribution. However, as shown by , 2003; Cho, Lazarian & Vishniac (2002) and , the decomposition of the different modes in MHD turbulence actually reveals that, although Alfvén and slow modes behave as K41 type of turbulence and are anisotropic, the fast modes are isotropic and follow IK statistics (see Figure 2). Effects of imbalanced (or cross-helicity) turbulence in the cascade and statistics of the local fields have also been addressed in the past few years (Lithwick, Goldreich & Sridhar, 2007;Beresnyak & Lazarian, 2008Wicks et al., 2011;Markovskii & Vasquez, 2013, and references therein). Imbalanced turbulence occurs when waves traveling in opposite directions along the mean magnetic field are of unequal amplitudes, i.e. carry different energy fluxes to small length scales, so that z + l /z − l 1 and ǫ + l /ǫ − l 1. The imbalance may arise in MHD turbulence since the interaction timescales between the waves z + l and z − l are different, and the cascade generally occurs faster for z − l . This is understood as the number of interactions (N) is much larger for counter-propagating wave packets, resulting in ǫ + l /ǫ − l > 1. In such a scenario, numerical simulations show that the anisotropy is not equal for the different wave modes.
Locality of scales for wave-wave interactions has also been the subject of recent studies in turbulence (Carati et al., 2006;Alexakis, 2007;Mininni, Alexakis & Pouquet, 2008;Aluie & Eyink, 2009;Beresnyak & Lazarian, 2010). Magnetic fields are responsible for long range interactions, from the Lorentz force acting over the whole fluid frozen to it. Therefore, different wavelengths may interact with each other non-linearly. Bi-spectra of fluctuations of density are discussed in Burkhart et al. (2009), and the nonlocal interactions appear to be important in MHD and supersonic turbulence models. A similar approach is used for studying the non-local interactions of Elsässer modes (Cho, 2010), resulting in a substantial fraction of nonlocal interactions in MHD turbulence. The role of the non-local interactions in the turbulent cascade is still not clear though. Turbulence in magnetized collisionless plasmas has been also studied in the past few years (e.g. Hellinger et al., 2006;Schekochihin et al., 2008;Bale et al., 2009, and others) in order to determine the role of collision-  2013), reveal that the statistics are still dominantly Kolmogorov-like, though strong asymmetries may also arise due to instabilities (firehose, mirror and cyclotron instabilities) are small scales.

Signatures of a turbulent ISM
In the previous section some theoretical aspects of turbulence have been presented. Its direct comparison to the dynamics of the interstellar gas is not trivial, as we discuss in the following. However, we will present here some observational evidences for a turbulent ISM, and discuss the possible turbulent regimes that may be inferred from these. The recognition of a turbulent interstellar medium dates back to 1950's with the work of von Weizsäcker (1951) on the spatial distribution of dense structures in the plane of the sky. He recognized the hierarchy of structures and suggested its turbulent origin. The identification of turbulent motions was provided shortly after it was measured from velocity dispersions (von Hoerner, 1951). Later on, the observational and theoretical supports for a turbulence dominated ISM have grown considerably (see reviews by Elmegreen & Scalo, 2004;Mac Low & Klessen, 2004;Hennebelle & Falgarone, 2012, and references therein), causing a major shift in the uderstanding of the ISM nature, from a thermal pressure dominated system, as thought before, to a very dynamic multi-phase system.

Density distributions
As mentioned above, one of the main signature of the turbulent character of the ISM is related to the density distribution of its contents. Up to now, tracers of the gas density distributions of the ISM at large scales have been dominantly indirect 5 . They rely on spectral lines and continuum emissions from the different phases of the ISM: the hot and fully ionized (HIM), the warm and fully/partially ionized medium (WIM/WNM), and the cold weakly ionized (CNM). These emissions being integrated along lines of sight and projected in the plane of the sky, sophisticated inversion methods have to be implemented. It is the statistical approach of the temporal and spatial variability of these emission fluxes that are the readily accessible observational techniques for studying interstellar turbulence.
With hydrogen being the most abundant element in the Universe, the λ 21cm line of neutral hydrogen is a key diagnostic. Its line integrated emission is proportional to the bulk of the hydrogen column density, since its opacity remains low over most of the ISM. Statistics of the HI intensity spatial distributions have therefore been used to probe interstellar turbulence but the results are far from homogeneous. Green (1993) studied the power-law of the spatial power spectrum of the HI emission from different fields in our Galaxy. He obtained power spectra with slopes between -2.2 and -2.8 at a scale range between 35 and 200 pc. From the HI 21 cm absorption towards Cas A. Roy et al. (2009) derived a power law with index -2.7, consistent with Kolmogorov turbulence in the diffuse interstellar medium. However, (Miville-Deschênes et al., 2003) find an impressive power-law in the nearby ISM at high galactic latitude with the same slope of -3.6 over two orders of magnitude in scales (between 0.1 and 25 pc). Similar studies have been performed since then, including other density tracers such as the CO and 13 CO line emission of molecular clouds and power-laws have also been inferred (e.g. Bensch, Stutzki & Ossenkopf, 2001;Hill et al., 2008). A review of the scatter of the power-law slopes measured is given in (Hennebelle & Falgarone, 2012). The scatter of the slope values is certainly affected by projection effects: one would expect a 2D power spectrum k −8/3 for an intrinsic Kolmogorov scaling. However, the integration along lines of sight crossing often large amounts of turbulent ISM with different properties tends to blur such a simple law. Moreover, the different tracers originate in truly different phases of the ISM with varying amounts of small scale structure that may affect the power spectrum of the density distributions (i.e. in many cases, like supersonic turbulence, density fluctuations are not simply advected by turbulence as passive scalars, see Audit & Hennebelle (e.g. 2005)). Indeed, as seen in Fig. 10 of (Hennebelle & Falgarone, 2012), many studies (including the power spectrum of the dust thermal emission) give power-laws indices close to -2.7. It is not possible though to presume that a Kolmogorov-like cascade operates in the ISM, with scalings given by Equations 2 and 3. Even though compressibility seems to play little effect on the statistics of the ISM, except for small scales (∼ pc scales) and cold and dense regions, magnetization effects may be important, as we discuss further below. Armstrong, Rickett & Spangler (1995) used another tracer of density fluctuations, the scintillation of the background radiation (i.e. changes in the refraction index due to the turbulent motions in the ionized components of the ISM) in order to obtain the density spectrum along the line-of-sight. As a complementary method, fluctuations of the Faraday rotation measurements (RM) in the plane of sky are also used to estimate density fluctuations (once the magnetic field is known) on the line-of-sight (e.g. Minter & Spangler, 1996). The combined data provide the density fluctuations along the line-of-sight, but for different lengthscales, as seen in Figure 3. The turbulence probed by both methods (scintillations and RM) present a most impressive spectrum, with a unique Kolmogorov-like slope across more than ten orders of magnitude in wavenumber.
Similar works have been done for external galaxies. Turbulence has been characterized based on similar techniques for the Small Magellanic Cloud (see Stanimirovic et al., 1999;Stanimirovic & Lazarian, 2001;Chepurnov et al., 2008;Burkhart et al., 2010) and revealed spatial variations of HI morphology. Dutta et al. (2013) calculated HI intensity fluctuation power spectrum for a sample of 18 spiral galaxies and found slopes in the range of -1.9 to -1.5. Shallower spectra, compared to K41, could be evidence for two-dimensional eddy dominated turbulence at scales larger than the disk thicknesses.

Direct statistical analysis
Spectral lines of several species observed with high spectral resolution may be used to infer the turbulence velocity distributions in the different phases of the ISM, such as hydrogen lines (mostly) and some ions for the diffuse ISM (e.g. Bowen et al., 2008), and molecular spectral lines ( 12 CO and 13 CO in most surveys) for the molecular clouds. The early surveys of Larson (1981) and Solomon et al. (1987) revealed the universal line-width and mass distribution scalings among molecular clouds. Notably, both works pointed to a velocity dispersion relation σ v ∝ l α ν , with α ∼ 0.5 (see Figure 4, left). Many similar studies were carried out to study the velocity distribution in molecular clouds, such as the work by Goldsmith et al. (2008); Yoshida et al. (2010); Qian, Li & Goldsmith (2012); Heyer & Brunt (2012) in the Taurus Molecular Cloud; Gustafsson et al. (2006) and Liu, Wu & Zhang (2012) for the Orion Complex, and many others.
More recent studies confirmed the same scaling relation although with slopes varying significantly (Heyer & Brunt, 2004;Qian, Li & Goldsmith, 2012). Qian, Li & Goldsmith (2012) for instance used the variance of the velocity difference of cores in molecular clouds, instead of the line width, and obtain α ν ∼ 0.7. On the other hand, massive cores are known to exhibit shallower slopes compared to what is frequently assumed (i.e. α ν < 0.5).
Recently, Ballesteros-Paredes et al. (2011) compiled different observational surveys and concluded that while in general terms, the typical CO clouds observed by Heyer et al. (2009) lie close to Larson's relation, this is clearly not the case for the dense and massive cores, which exhibit large velocity dispersions for their relatively small sizes (Figure 4). Those authors propose that the large dispersion observed at small scales are related to increased velocities as the clouds become gravitationally bound. However the increased dispersion at small scales has already been reported though, based on numerical simulations, in Falceta-Gonçalves et al. (2010a) without self-gravitating objects. For these authors the large dispersion observed at small scales is an intrinsic feature at the turbulent gas. The broad dispersion of the scaling relation indicates a turbulent regime dominated by compressible motion at small scales, as discussed in Section 1.1, though regular incompressible turbulence dominates at larger scales. Compressibility, as described in Fleck's model (see Equation 10), naturally give larger slopes for the dispersion relation, with a value of α ∼ 0.16 favoured by observations. It is not clear though what is the actual role of gravity in the statistics of the molecular cloud emissions.
At the large scale end of the cascade, the apparent uniqueness of the scaling of the velocity dispersion with size scale suggests an universal source (or mixture of sources) of energy for the molecular gas turbulence in our Galaxy. (Chepurnov & Lazarian, 2010) presented statistical analysis of high-latitude HI turbulence in the Milky Way based on velocity coordinate spectrum (VCS) technique. They found a velocity power spectrum P u (k) ∝ k −3.8 and an injection scale of ∼ 140±80pc. The alightly steeper slope, compared to K41, can be the result of shock-dominated (compressible) turbulence, with averaged sonic Mach numbers ∼ 7 − 8 (see Section 1.1 above).
Two-point statistics are also used but, since in situ measurements are not yet available, one easily accessible observable turns out to be the variations in the plane-of-the-sky of the line-of-sight centroid velocity of spectral lines. Lis et al. (1996) showed that they trace the planeof-the-sky projection of the vorticity. Using a sample of about one million independent CO spectra in a diffuse field, Hily-Blant & Falgarone (2009) identified, on statistical grounds, the ensemble of positions at which vorticity departs from a Gaussian distribution. These form coherent elongated structures at the parsec-scale that are found to harbor substructures of most intense velocity shears down to the milliparsec scale (Falgarone, Pety & Hily-Blant, 2009). These coherent structures are proposed to be the manifestations of the intermittency of turbulent dissipation in diffuse molecular clouds (see the review of (Hennebelle & Falgarone, 2012)), which may be compared to Equation 8 above. Li & Houde (2008) studied the scaling relations of the velocity dispersions from different neutral and ionized molecular species, namely HCN and HCO + , in the region of M17. As it occurs in many other star forming regions, the ionized molecules systematically present smaller dispersion of velocity compared to the neutral. Such a difference arises as turbulent energies dissipate differently for the species due to ambipolar diffusion. Falceta-Gonçalves, Lazarian & Houde (2010b) showed that the dispersion for ions is typically smaller than that for the neutral species basically due to the damping of the ion turbulence at the ambipolar diffusion scales (≃ 0.01pc).
The direct comparison between statistics of observational data and the theory must be done with caution. Column density projections, or in other sense emission maps, are influenced by projection effects. Different structures projected on the same line-of-sight, but decorrelated at a given lengthscale, may be observed as a single structure in the projected emission map. Some deconvolution is possible though once the velocity profile is known, with high spectral resolution.

Indirect access to the velocity field via maser emission
The low surface brightness of the above tracers and projection effects make the direct analysis of turbulent flows in the ISM difficult. Maser spots that are bright point sources and are transported by turbulence as passive scalars (because they are tiny and low mass structures), turn out to be powerful tracers of the turbulent velocity field. Maser radiation in molecular lines appear in dense regions where population levels inversion can be generated by radiative pumping, for instance, e.g. in the dense molecular gas of starforming regions (SFRs) associated with ultra-compact HII regions, embedded IR sources, hot molecular cores, Herbig-Haro objects, and outflows (Litvak, 1974;Reid & Moran, 1981;Elitzur, 1992;Lo, 2005). Maser emissions are often characterized by high brightness temperatures and highdegrees of polarization. Intense maser emission is detected in the molecular lines of hydroxyl (OH), water (H 2 O), sili- con monoxide (SiO), ammonia (NH 3 ), methanol (CH 3 OH), among others. Walker (1984) used the Very Long Baseline Interferometry (VLBI) maps of the H 2 O maser source in W49N to demonstrate that both two-point velocity increments and two-point spatial correlation functions exhibit power-law dependencies on the maser spot separation, which is indicative of a turbulent flow. Gwinn (1994) performed statistical analysis of VLBI data for W49N to confirm the power-law dependence of the velocity dispersion and spatial density of masing spots on spatial scale, and interpreted this observation as an evidence of turbulence. Imai, Deguchi & Sasao (2002) reported sub-milliarcsec structures of H 2 O masers in W3 IRS 5. A cluster of maser spots (emission spots in individual velocity channels) displays velocity gradients and complicated spatial structure. Two-point spatial correlation functions for the spots can be fitted by the same power laws in two very different spatial ranges. The Doppler-velocity difference of the spots as a function of spot separation increases as expected in Kolmogorov-type turbulence. Strelnitski et al. (2002) used VLBI data to investigate the geometry and statistical properties of the velocity field traced by H 2 O masers in five starforming regions. In all sources the angular distribution of the H 2 O maser spots shows approximate self-similarity over almost 4 orders of magnitude in scale. The lower order structure functions for the line-of-sight component of the velocity field can be fitted by power laws, with the exponents close to the Kolmogorov value. Similar results were also obtained for other regions (e.g. Richards et al., 2005;Strelnitski, 2007;Uscanga et al., 2010).

Turbulent magnetic fields
The magnetic fields in our Galaxy is modelled as a superposition of different components: i) a large scale field, following a spiral structure in the plane of the galactic disk, and extending high above the plane into the Galactic halo, and ii) a complex component of locally disturbed magnetic fields, which are related to molecular clouds and star formation regions. The spiral pattern in the disk aligns with the spiral arms (e.g. Han, 2006). This is expected since the shear of gaseous motion around the center of the galaxy stretches the field lines in this direction (see review of mean field dynamo by Beck et al., 1996).
There are four main methods to study the fluctuations in the ISM magnetic field, namely the polarization of dust thermal emission (both in emission in the far-infrared (FIR) and absorption in the visible and near-IR), Zeeman effect of spectral lines, Faraday rotation and polarization of the synchrotron emission. Polarized synchrotron emission can also be mapped in order to provide the geometry of the field lines in the plane of the sky. Faraday rotation and synchrotron polarization measurements excel in probing the magnetic field of the diffuse ionized medium of the ISM, i.e. they are excellent tools to study the large scale fields of galaxies in general. More extensive reviews both on magnetic fields in star formation regions and galactic scale magnetic fields are given in Crutcher (2012) and Han (2006), respectively.
As mentioned earlier, synchrotron emission polarization can be used for mapping the large scale structure of the magnetic fields in galaxies (see review by Beck, 2009). The fields traced by the polarized synchrotron emission present intensities of the order of ∼ 10 − 15µG. However, the synchrotron emission probes the ionized medium only, which is less use-ful in determining the turbulence properties of the star formation regions, dominated by the dense and neutral components of the ISM. Therefore, a magnetic field with intensity ∼ 10 − 15µG is supposed to thread most of galactic disk, except the dense regions of the arms where the local properties of the plasma and stellar feedback may dramatically change the field properties. Oppermann et al. (2012) compiled an extensive catalog of Faraday rotation measure (RM) data of compact extragalactic polarized radio sources in order to study the angular distribution of the all-sky RMs. The authors found an angular power spectrum P(k) ∝ k −2.17 for the Faraday depth, which is given by the product of the line-of-sight magnetic field component B LOS and the electron number density n e . The combination of the RM and polarization vectors of the synchrotron emission allows to reconstruct the three-dimensional structure of galactic magnetic fields. Such angular fluctuations of the Faraday depth is thought to be related to the turbulent ISM. However, the relationship between the fluctuations of the RM and the local fluctuations of electron density and magnetic fields is not clear yet. This, for instance, is an interesting subject for further comparisons with simulations (as in Gaensler et al. (2011)).
Possibly the most direct method for estimating the magnetic field intensity in the dense and cold ISM relies on the detection of Zeeman effect (see Robishaw, Quataert & Heiles, 2008, for details). For instance, Sarma et al. (2002) detected and studied the Zeeman effect in H 2 O masers in several SFRs and determined line-of-sight magnetic field strengths ranging from 13 to 49 mG. They found a close equilibrium between the magnetic field energy and turbulent kinetic energy in masing regions. Alves et al. (2012) showed that shock-induced H 2 O masers are important magnetic-field tracers of very highdensity gas in low-mass protostellar core IRAS 16293-2422. They investigated whether the collapsing regime of this source is controlled by magnetic fields or other factors such as turbulence, and concluded that the magnetic field pressure derived from data is comparable to the ram pressure of the outflow dynamics. This indicates that the magnetic field is energetically important for the dynamical evolution of the protostellar core.
Due to its brightness, maser emission is better for probing magnetic fields, but they are rare and limited in extent. The Zeeman effect in non-masing regions has been detected for HI, OH, and CN lines for which the turbulent broadening is typically larger than the Zeeman splitting in frequency. The compilation by Crutcher (1999) and recent CN Zeeman observations in SFRs by Falgarone et al. (2008) show that the turbulent motions within the SFRs and molecular clouds are supersonic but sub-Alfvénic. The upper limit magnetic field intensity scales with density, estimated from a Bayesian analysis, as B ∝ n κ , with κ ∼ 0.47 (Crutcher et al., 2010). Collapsed structures along the mean field would produce κ → 0, while shock compressions per-pendicular to the field lines result in κ → 1. The observed relation with κ ∼ 0.47 is expected, for instance, in Alfvénic perturbations and is in agreement with MHD simulations (e.g. Burkhart et al., 2009). It was also claimed in that work that, despite its relative importance in the overall dynamics of clouds, the uniform magnetic fields in these clouds are in general not strong enough to prevent gravitational collapse based on the mass-to-flux (M/Φ) ratios observed. Other major compilations of Zeeman measurements in molecular clouds are given, for example, in Bourke et al. (2001) and  with similar results.

Polarization maps of molecular clouds
Radiation may be polarized due to a prefered direction for emission/absorption from aspherical dust grains, as well as by some molecules and atoms. The ISM in known to be populated by a complex distribution of grain sizes and shapes. Depending on its composition an aspherical rotating dust particle may align with the magnetic field line. The orientation of the polarization of radiation is then linked to the orientation of the magnetic field itself (see review by Lazarian, 2007). Many observational data has been made available in the past decade both on absorption and emission dust polarization (e.g. Heiles, 2000;Chapman et al., 2011).
The strengths of magnetic fields can be estimated from polarization maps by the Chandrasekhar & Fermi (1953) (CF) technique. The CF method is based on the assumption that the magnetic and turbulent kinetic pressures are the dominant ones within the cloud, and that the fluid motions are coupled to the magnetic field lines. In this sense, any perturbation from the fluid turbulence will result in a change in the orientation of the field lines. Major improvements on the CF technique are given e.g. by Falceta-Gonçalves, , Hildebrand et al. (2009) andHoude et al. (2009). If the velocity dispersion δv los is known, e.g. from spectral lines, the mean magnetic field in the plane of sky can be estimated as (Falceta-Gonçalves, : where δφ represents the dispersion in the polarization angle. From the equation above, the ratio δB/B sky -assumed to be ∼ tan −1 δφ -is directly related to the Alfvénic Mach number of the turbulence. Notice that the dependence of the projected δB/B sky with the actual 3D MHD turbulence may be removed from higher order statistical analysis, as proposed in Falceta-Gonçalves, .
The left image of Figure 5 presents the polarization map of the Muska Dark Cloud (Pereyra & Magalhães, 2004) in the optical wavelengths, as a result of dust absorption. Vectors represent the magnetic field orientation. The filamentary morphology of the dark cloud is perpendicular to the external field, which is very uniform indicating a sub-Alfvénic turbulent regime. At the right hand side of Figure 5, the polarization map overplotted on the column density projection of a 3D MHD numerical simulation of sub-Alfvénic turbulence (Falceta-Gonçalves, . Such comparisons between MHD numerical simulations and measurements of magnetic fields in the ISM are important in unveiling the physics of MHD turbulence, and its role on other phenomena such as star formation. Spatial dispersion of magnetic fields in molecular clouds from polarization maps may be used to characterize the power spectrum of magnetized turbulence in the inertial and dissipation ranges. Houde et al. (2011) found a power law inertial range for the magnetic field spatial distribution that is ∝ k −2.9±0.9 , and a cutoff at scales ∼ 0.009pc, which is claimed by the authors to be related to the ambipolar diffusion scales.
Again, as another issue in a proper modelling of the statistics of velocity, gravity is claimed to interfere in the statistics of the observed polarization maps (Koch, Tang & Ho, 2012a,b). Gradients in emission towards the cores of molecular clouds have been shown to be associated to gradients in the polarization angles. A transition from a magnetically subcritical to a supercritical state 6 could then explain the trend, and this technique could provide an independent way to estimate the local magnetic force compared to gravity. Heyer & Brunt (2012) showed that the turbulence in the densest regions of Taurus molecular cloud is super-Alfvénic, while the reverse is true in the surrounding lower density medium, threaded by a strong magnetic field. This observational result is in agreement with the transition expected between scales as dense structures are formed, e.g. by shocks, in a supersonic but sub-Alfvénic large scale turbulence (see discussion in Falceta-Gonçalves, Heyer et al., 2008;Burkhart et al., 2009).
Similar to the synchrotron radiation case, by combining dust polarization maps with Zeeman measurements in molecular clouds one can determine the three-dimensional structure of the magnetic field. Poidevin et al. (2013) recently succeeded in testing this approach for a number of objects of the SCUBA Polarimeter Legacy (SCUPOL) data catalog. The authors were able to determine the orientation of the mean field with respect to the line-of-sight, as well as to estimate the turbulence regime within several molecular clouds. The authors also claimed that all observed clouds seem to present a universal large scale turbulence that is supersonic (M s ∼ 6 − 8) and sub-Alfvénic (M a ∼ 0.5 − 0.9), at scales as large as 50pc.
In terms of comparing these data with basic theories of magnetized turbulence, most observations point towards a magnetically dominated turbulence at scales larger than few tenths of parsecs. Heyer et al. (2008) also showed one of the first evidences for anisotropic turbulence in molecular clouds, with respect to the large scale magnetic field orien-tation. The observations of the Taurus Molecular Cloud revealed a significant anisotropy in the dispersion of velocity (δv), being larger for lags perpendicular to the mean large scale field lines. Even though a Goldreich-Sridhar similarity relation is not obtained, the anisotropy observed is a strong indication for strong coupling between MHD wave modes in the insterstellar turbulence, as predicted by GS95 model. We could extrapolate a bit and say that a GS95 model combined with fractal density distributions, as given in Fleck (1996), is favoured for the ISM turbulence based on current observations.

Origins of interstellar turbulence
Surveys of different atomic and molecular line emissions have shown us that the diffuse ISM is turbulent at scales > 150pc, with δv ≥ 50km s −1 . This results in a specific energy transfer rate 7 of ǫ ≃ m H n H δv 3 L /L ∼ 10 −25 −10 −24 erg cm −3 s −1 . Brunt, Heyer & MacLow (2009) estimated the driving scales of turbulence for molecular clouds by comparing observed and synthetic CO velocity dispersions from numerical simulations. They found that only models with large scale sources of turbulence, such as supernovae-driven outflows (SNe) and galactic dynamics, fit well to the observed data.
One issue is that numerical simulations of SNe driven turbulence create superbubbles that are far too hot and diffuse (see Avillez & Breitschwerdt, 2005;Joung & Mac Low, 2006;Hill et al., 2012). Other critical arguments disfavour SNe as a main driver mechanism as well, at least for our Galaxy. Zhang et al. (2001) analysed the CO emission lines from the Carina Complex and obtained a turbulent energy flux per mass density unit cascading over scales ∼ 10 −7 (km s −1 ) 2 yr −1 , that could not be explained from stellar feedback, but is in rough agreement with the injection rate of energy from the gravitational interaction of the ISM gas and the galactic spiral arms. Sánchez-Salcedo, Santillán & Franco (2007) also showed that HI mapping of our Galaxy is consistent with a turbulence injection rate that is not directly related to the star formation rate, but is about constant with respect to the galactocentric radius. Also, the correlation lengths related to SNe turbulence is strongly dependent on local properties (such as local density and temperature) (see  Pereyra & Magalhães, 2004). Right: simulated polarization map from a three dimensional simulation of MHD sub-Alfvénic turbulence (extracted from Falceta-Gonçalves, . Leão et al., 2009). Such local dependence also occurs with respect to the height related to the galactic plane, since SNe energy is easily released outwards (e.g. Melioli et al., 2009;Marasco, Marinacci & Fraternali, 2013). The universality of the observed properties of turbulence in our Galaxy, together with the extremely large injection scales (> 100pc) suggest a Galactic scale driving source, which is later amplified, as second order effects, by local stellar feedback. (Qian, Li & Goldsmith, 2012), for instance, obtained similar core and ambient turbulent statistics, which suggested that molecular cores condense from more diffuse gas, and that there is little (if not none) additional energy from star formation into the more diffused gas.
Turbulence driven by galactic dynamics models, such as driven by velocity-shears in the galactic disk, have been posed long ago (e.g. Fleck, 1981). Instabilities such as the magneto-rotational instability have also been proposed (e.g. Sellwood & Balbus, 1999;Kim & Ostriker, 2002). Interactions between the arms of the Galaxy and the disk gas also generate perturbations, as large as 20km s −1 (Gómez & Cox, 2002), that could explain most of the injection of energy into turbulent motions. It is not clear yet which of these mechanisms (SNe or galactic dynamics) is more important for the observed turbulence in the ISM. Certainly, it is a promising subject for studies in the upcoming years, both from theoretical and observational sides.

Conclusions
In this work, we briefly reviewed part of the current understanding on incompressible, compressible and/or magnetized turbulence, which can be applied to characterize the interstel-lar medium. There is a vast literature available for each of these and a complete review on turbulence is out of the scope of this work. We discussed the recent theoretical improvements made on the modelling and characterization of the different turbulent regimes. Multifractal description, statistics of probability functions, and spectral analysis are just a few that have been currently employed to characterize spatial and temporal variations of plasma properties associated to turbulent motions.
Phenomenological descriptions of turbulence in Fourier space, such as that of Kolmogorov-Obukhov, are particularly simple and still very useful on the diagnostics of interstellar turbulence. Since scaling relations for compressible, incompressible and magnetized turbulence of these theories may differ among each other, observations can be used to determine the turbulent regime of the ISM.
Spectroscopy has been long used to probe the velocity distributions along the line of sight. The observed amplitudes of the turbulent motions indicate that the ISM transits from a supersonic turbulent regime at scales of tens to hundreds of parsecs, at which the turbulence is driven, to subsonic at smaller scales. The scales where turbulence is subsonic depend on the "phase"' of the ISM plasma. Dense molecular clouds present lower temperatures, which result in subsonic turbulence only at very small scales (≪ 1pc). The warmer and more diffuse media, such as warm neutral medium and the warm diffuse medium, present subsonic flows at scales of few parsecs due to the larger local sound speeds. It is interesting to mention that this transition is deeply related to the origin of the dense molecular clouds. These objects are either originated due to the large scale compressible motions of the gas (e.g. Williams et al., 2000), or at small scales due to other mechanisms, such as thermal instabilities. Current observations favour the first, given their lengthscales and internal dynamics (see Poidevin et al., 2013). Spatial gas distributions over the plane of sky are also provided observationally. The filamentary structure observed reveals a compressible dominated turbulent regime, at least at most of the observed scales. Observations also reveal a magnetized ISM. All these ingredients combined result in compressible and magnetized ISM turbulence, challenging theorists to provide a phenomenological description of the combined effects of supersonic flows and strong magnetic fields. Despite the good agreement between observations and the Goldreich-Sridhar model for magnetized turbulence, and Fleck's model for compressible one, such as spectral slopes, scaling relations and multifractal analysis applied to emission maps, a complete unified theory is yet to be developed.
One of the major problems in comparing statistics of observed quantities to theories of turbulence relies on projection effects. Observations are spatially limited in the sense that all statistics are done either along the light-of-sight (e.g. scintillation, velocity dispersion from spectral lines, Faraday rotation) or in the plane-of-the-sky. In addition, even the plane-of-sky maps are related to integrated quantities (e.g. emission lines, column density, Stokes parameters for the polarization maps). One must therefore be careful when comparing these with theories of three dimensional turbulence.
Other effects may also make the direct comparison between theory and observations challenging. Self-gravity of dense gas and stellar feedback, for instance, have been neglected in this paper. These processes are responsible for extra sources of energy and momemtum, but are not easily linked to the turbulent cascade. Despite their obvious importance on the process of e.g. star formation, their role on the statistics of the turbulence is not completely clear. Naturally, fragmentation and clumping would be enhanced if self-gravity is considered (Vázquez-Semadeni et al., 1996;Ballesteros-Paredes et al., 2011;Cho & Kim, 2011), however its role on the cascade itself and on intermittency is unknown.
Future studies from the theoretical side are possibly to be focused on the understanding of combined processes, such as magnetic fields, gravity, compressibility and radiation, on the energy transfer among scales. Formation of coherent structures and how their statistics relate to the bulk of the fluid are vital for theories of star formation. New data is also expected for the upcoming years. Although the Herschel mission ended in early 2013 its data are not yet fully explored. Other major observational facilities, such as the Planck 8 satellite and the Atacama Large Millimeter Array (ALMA), will provide complementary data at radio and microwave frequencies with very large sensitity, therefore going "deeper" than reached by other instruments. Also, (Gurnett et al., 2013) recently presented the first in situ measurement of the interstellar plasma as Voyager 1 has crossed the heliopause and started to probe the nearby interstellar plasma. This opens new possiblities in studying interstellar turbulence locally. It is clear that the future of interstellar turbulence science is going to be very exciting.