Turbulence and star formation efficiency in molecular clouds: solenoidal versus compressive motions in Orion B^{⋆}
^{1} Univ. Grenoble Alpes, IRAM, 38000 Grenoble, France
email: orkisz@iram.fr
^{2} IRAM, 300 rue de la Piscine, 38406 SaintMartind’ Hères, France
^{3} LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, École normale supérieure, 75005 Paris, France
^{4} LERMA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, 92190 Meudon, France
^{5} Laboratoire d’Astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, Allée Geoffroy SaintHilaire, 33615 Pessac, France
^{6} HarvardSmithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA, 02138, USA
^{7} Joint ALMA Observatory (JAO), Alonso de Cordova 3107 Vitacura, Santiago de Chile, Chile
^{8} ICMM, Consejo Superior de Investigaciones Cientificas (CSIC), 28049 Madrid, Spain
^{9} National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA, 22903, USA
^{10} School of Physics and Astronomy, Cardiff University, Queen’s buildings, Cardiff CF24 3AA, UK
^{11} IRAM, Avenida Divina Pastora, 7, Núcleo Central, 18012 Granada, España
^{12} Maison de la Simulation, CEACNRSINRIAUPSUVSQ, USR 3441, Centre d’Étude de Saclay, 91191 GifSurYvette, France
Received: 30 June 2016
Accepted: 2 January 2017
Context. The nature of turbulence in molecular clouds is one of the key parameters that control star formation efficiency: compressive motions, as opposed to solenoidal motions, can trigger the collapse of cores, or mark the expansion of Hii regions.
Aims. We try to observationally derive the fractions of momentum density (ρv) contained in the solenoidal and compressive modes of turbulence in the Orion B molecular cloud and relate these fractions to the star formation efficiency in the cloud.
Methods. The implementation of a statistical method applied to a ^{13}CO(J = 1−0) datacube obtained with the IRAM30 m telescope, enables us to retrieve 3dimensional quantities from the projected quantities provided by the observations, which yields an estimate of the compressive versus solenoidal ratio in various regions of the cloud.
Results. Despite the Orion B molecular cloud being highly supersonic (mean Mach number ~ 6), the fractions of motion in each mode diverge significantly from equipartition. The cloud’s motions are, on average, mostly solenoidal (excess > 8% with respect to equipartition), which is consistent with its low star formation rate. On the other hand, the motions around the main star forming regions (NGC 2023 and NGC 2024) prove to be strongly compressive.
Conclusions. We have successfully applied to observational data a method that has so far only been tested on simulations, and we have shown that there can be a strong intracloud variability of the compressive and solenoidal fractions, these fractions being in turn related to the star formation efficiency. This opens a new possibility for star formation diagnostics in galactic molecular clouds.
Key words: turbulence / methods: statistical / ISM: clouds / ISM: kinematics and dynamics / radio lines: ISM / ISM: individual objects: Orion B
© ESO, 2017
1. Introduction
The evolution of molecular clouds is controlled by a complex interplay of largescale phenomena and microphysics: chemistry and interaction of the matter with the surrounding farUV and cosmicray radiation control the thermodynamic state of the gas and its coupling to the magnetic field. The medium is highly turbulent, with Reynolds numbers reaching 10^{7} and magnetic Reynolds numbers reaching 10^{4} (Draine 2011). Magnetohydrodynamic (MHD) turbulence is one of the main counteractions to gravity (Hennebelle & Falgarone 2012; Federrath & Klessen 2012; Padoan et al. 2014), as well as the major mechanism that shapes the clouds: their fractal geometry is related to the properties of their turbulent velocity field (Pety & Falgarone 2000; HilyBlant et al. 2008; Federrath et al. 2009). The dissipation timescale for the turbulent energy of a molecular cloud is shorter (~ 1 Myr, Mac Low 1999) than the age of these clouds (~ 20−30 Myr, Larson 1981). Hence, a continuous energy injection must exist (Hennebelle & Falgarone 2012, and references therein). The proposed injection mechanisms may be either external, for instance Galactic shear or nearby supernovae explosions (Kim & Ostriker 2015), or internal, like the expansion of Hii regions and molecular outflows of the recently formed stars (Hennebelle & Falgarone 2012). The nature of turbulence, and notably its solenoidal or compressive forcing, also plays a key role in the star formation efficiency (SFE) of molecular clouds (Federrath & Klessen 2012, 2013). In particular, compressive motions bear the mark of various phenomena related more or less directly to the formation of stars: infall on filaments, collapsing dense cores, expansion around young stars, etc. As a result, we can expect that a more compressive cloud is a more active cloud, probably more likely to be forming stars, as proposed by Federrath & Klessen (2012).
In this work, we propose to measure for the first time (to our knowledge) the relative fractions of momentum in the solenoidal and compressible modes of turbulence in a molecular cloud, following a method devised and tested on numerical simulations by Brunt et al. (2010) and Brunt & Federrath (2014). Our goal is to obtain a quantitative estimation of these fractions from observational data, and to compare their ratio with the star formation efficiency that was derived from independent data. We investigate whether different fractions of compressive and solenoidal motions might provide a diagnostic for the variation of the star formation efficiency among molecular clouds. A different proportion of compressive forcing might be the reason why some molecular clouds form stars at a high rate while others do not (Federrath et al. 2010; Federrath & Klessen 2012; Renaud et al. 2014).
Our object of study is a large region of a nearby giant molecular cloud (GMC), namely the southwestern edge of the Orion B cloud (Barnard 33 or Lynds 1630). Orion B is relatively close to us, at a typical distance of ~ 400 pc (Menten et al. 2007; Schlafly et al. 2014), so that a spatial resolution of 25″ corresponds to 0.05 pc or 10^{4} AU in the cloud. The total mass of Orion B is estimated to be 7 × 10^{4} M_{⊙} (Lombardi et al. 2014), and the average incident FUV radiation field is G_{0} ~ 45 (Pety et al. 2017). Orion B is located in the Orion GMC complex (Kramer et al. 1996; Ripple et al. 2013), east of the famous Orion Belt. Alnitak, the eastmost of the three belt stars shines in the foreground of the cloud. The southwestern edge of the cloud represents an ideal laboratory to study star formation, and features several remarkable regions. First, the cloud is illuminated by the massive star σ Ori that creates an Hii region, the emission nebula IC 434, bounded on its eastern side by an ionization front. Silhouetted against this bright background, a dark cloud can be seen: the famous Horsehead nebula. HD 38087 also creates a small Hii region, IC 435. Still embedded in Orion B, the star forming region NGC 2024, known as the Flame Nebula, hosts several massive Otype young stellar objects, which have created compact Hii regions inside the cloud. NGC 2024 lies just east of the Alnitak star, and is crossed by a filament that is seen in absorption in visible light, and in emission in the radio range. The reflection nebula NGC 2023 is a quieter counterpart of NGC 2024, hosting young Btype stars. It lies northeast of the Horsehead nebula. The rest of the cloud contains extended and quieter areas with strong filamentary structures (Fig. 1). This area has been extensively observed in the 3 mm range with the IRAM30 m telescope (PI: J. Pety). This survey has led to a series of papers (Liszt & Pety 2016; Pety et al. 2017; Gratier et al. 2017) to which this article belongs.
The paper is organised as follows. In Sect. 2, we briefly describe the observations by the ORIONB (Outstanding Radio Imaging of OrioN B) collaboration, and the data we use here. In Sect. 3, we present the concepts and equations of the statistical method and the details of its implementation, from noise filtering to the computation of power spectra and the estimation of velocitydensity correlations. The results are described in Sect. 4, and discussed in Sect. 5 with a special emphasis on the relation of the turbulence properties with the star formation efficiency in Orion B. In the appendix, we present our computation of the Mach number map in the cloud.
2. Observations
2.1. The Orion B project dataset
The Orion B project (PI: J. Pety) has already mapped the southwestern edge of the Orion B molecular cloud with the IRAM30 m telescope over a field of view of 1.5 square degrees in the full frequency range from 84 to 116 GHz at 200 kHz spectral resolution (Pety et al. 2017). The red rectangle in Fig. 1 shows the field of view observed up to now, over an H_{2} column density map produced by the Herschel Gould Belt Survey consortium (André et al. 2010; Schneider et al. 2013). The mapped area covers 56 × 98 arcmin (about 6 × 11 parsecs at the assumed distance of Orion B, 400 pc) in size.
So far, the observations provided about 250 000 spectra over a 32 GHz bandwidth, yielding a positionpositionfrequency cube of 370 × 650 × 160 000 pixels, each pixel covering 9′′ × 9′′ × 0.5 km s^{1} (at Nyquist sampling). The reduced dataset amounts to 100 GB of data. The data reduction is described in details in Pety et al. (2017).
Fig. 1 H_{2} column density map of the southwestern part of the Orion B giant molecular cloud, derived from Herschel Gould Belt Survey observations (André et al. 2010; Schneider et al. 2013). The field observed by the OrionB collaboration and used for this work is overlaid in red. The blue squares mark the nebulae NGC 2024, NGC 2023 and the Horsehead, and the pink star symbols mark the stars Alnitak, HD 38087 and σ Ori (north to south). 

Open with DEXTER 
Fig. 2 Maps of the average brightness temperature of the ^{13}CO(J = 1−0) line in four contiguous velocity ranges. The mainbeam temperature scale is indicated by the color bar on the right. The contour shows the value of 8.9 K km s^{1} in the W_{0} map, corresponding to 0.43 K in the mean temperature map integrated over the −0.5−20.5 km s^{1} velocity range. The set of coordinates used for the observational campaign takes the Horsehead PDR as a reference point, and aligns the IC 434 PDR along the vertical axis (14° counterclockwise rotation with respect to equatorial coordinates). The numbered squares in the first panel show the positions of the spectra presented in Fig. 3, from left to right. 

Open with DEXTER 
This dataset is unique for its extensive coverage, both spatially and spectrally, at a typical resolution of 25′′ and a median noise of 0.1 K[T_{mb}] per channel over a 32 GHz spectral range. From the spatial point of view, the survey gives us access to a large range of scales in the molecular cloud, from 50 mpc to 10 pc. For the analysis of turbulence properties (see Sect. 4.3), it means that we are able to study a large fraction of the inertial range of the turbulence, with a potential view on the injection scale. The dissipation scale, on the other hand, is of the order of a milliparsec (Hennebelle & Falgarone 2012; MivilleDeschênes et al. 2016), and, at a distance of 400 pc, is only accessible using millimetre interferometers, and out of reach for the IRAM30 m telescope.
From the spectral point of view, having such a large bandwidth observed in one go enabled us to image over 20 chemical species (Pety et al. 2017), including those listed in Table 1. As opposed to several small bandwidth mappings, the spectral lines in this survey are observed in the same conditions and are well intercalibrated, which gives an unprecedented spectral accuracy for such a large field of view.
2.2. The ^{13}CO spectral data cube
Most of the work presented here was performed on the ^{13}CO(J = 1−0) datacube^{1}, which covers a velocity range of 40 km s^{1} centred around the source systemic velocity of 10.5 km s^{1} and a rest frequency of 110.201 354 GHz. The datacube presents an root mean square (rms) noise of σ = 0.17 K, and a median signaltonoise ratio of T_{peak}/σ = 7.9.
The ^{13}CO(J = 1−0) line was chosen because it offers one of the highest signaltonoise ratios over the whole map, but it does not feature as much saturation as the ^{12}CO(J = 1−0) line. Signal is present at a S/N greater than 5 in the whole map, except in the Hii regions around σ Ori and HD 38087, where molecular gas is photodissociated. The brightest regions are the NGC 2023 nebula, the center of the NGC 2024 nebula, and its northern edge (Fig. 4, left panel).
^{13}CO is a good tracer of molecular gas, from moderately diffuse and translucent regions (A_{V} = 1−5 mag) up to moderately dense and shielded gas (10^{4} cm^{3}, A_{V} = 10 mag). CO has a small dipole moment (0.11 D), hence the rotational lines have low Einstein coefficients (e.g., Mangum & Shirley 2015). This leads to relatively easy collisional excitation, excitation temperatures approaching the kinetic temperature, and moderate line opacities except for the most abundant species, ^{12}CO. The abundance ratio ^{13}CO/^{12}CO is equal to ^{13}C/^{12}C or about 1/60 when chemical fractionation reactions, which are limited to the most diffuse regions, are inefficient (Wilson & Rood 1994). Therefore, the ^{13}CO abundance relative to H_{2} remains approximately constant at a level of ~ 2 × 10^{6} (Dickman 1978) across most of the cloud volume. ^{13}CO starts to be depleted on dust grains in cold cores, but these represent only a small fraction of the mass and a negligible fraction of the volume of the Orion B molecular cloud (Kirk et al. 2016).
Figure 2 shows the ^{13}CO(J = 1−0) signal integrated over four contiguous velocity ranges. It showcases the complexity of the spectral structure of the cloud, with prominent variations of the line profile with the position, which is a consequence of the strong turbulence at play in the molecular cloud.
Up to four spectral components appear along each line of sight within the field (Fig. 3). A main component is visible around 10 km s^{1}, and a secondary component at lower velocity (about 5 km s^{1}). Sometimes an extra component at higher velocity (about 14 km s^{1}), or secondary peaks around 5 and 10 km s^{1}appear too. The first two components are the most significant ones at the scale of the whole cloud, being the only ones visible in the mean spectrum of ^{13}CO(J = 1−0), and have average velocities of 9.7 km s^{1} and 4.9 km s^{1} respectively. All the components, despite being quite distinct on the spectral axis are, however, thought to be part of the Orion B cloud (see discussion in Pety et al. 2017).
Fig. 3 ^{13}CO(J = 1−0) spectra along selected lines of sights in the Orion B cloud, showing the diversity of velocity components (up to four per spectrum). The coordinates of the various lines of sights are given in arcminutes in our custom set of coordinates (δx, δy). The average spectrum for the whole field, normalized to have the same peak temperature, is superimposed with a dashed line for comparison. The positions of the spectra on the map are shown by white squares in the first panel of Fig. 2. 

Open with DEXTER 
3. Deriving the relative fraction of solenoidal motions from a positionpositionvelocity cube
To measure the fraction of the solenoidal and compressive turbulence modes, we apply the method developed by Brunt et al. (2010) and Brunt & Federrath (2014). In this section, we first recall the method, its assumptions, and the way we implemented it.
3.1. Description of the method
3.1.1. Principles and assumptions
The key point of this method is the fact that the objects we observe are fundamentally 3dimensional (e.g., a molecular cloud), but the observer only has access to a 2dimensional projection along the line of sight of that object. Brunt et al. (2010) developed a method to retrieve properties of the 3dimensional object that we are interested in, and which corresponds to a 3dimensional field F_{3D}, via the properties of the 2dimensional observational F_{2D}, which is a projection of F_{3D} along the z axis. To achieve this, they use the fact that the Fourier transform of the 2dimensional field is proportional to the k_{z} = 0 cut through the Fourier transform of the 3dimensional field. In short, . If these fields are isotropic, i.e., if they are functions of k =  k  alone, with k = (k_{x},k_{y},k_{z}) the wave vector and k the wave number, the 2dimensional field enables us to reconstruct average properties of the 3dimensional field thanks to symmetry arguments.
The Brunt & Federrath (2014) method was developed as an application of the Brunt et al. (2010) method to the case of vector fields. For a vector field, such as a velocity or momentum field, the dimensionality reduction that is due to the projection is made worse by the fact that only one component of the vector (the lineofsight one) can be measured, thanks to the Doppler effect. In this case, the next main tool to retrieve 3dimensional properties is the Helmholtz theorem, which enables to decompose any vector field in its divergencefree (solenoidal) and curlfree (compressive) components, F_{⊥} and F_{}. These components are related via a local orthogonality in Fourier space, . Solenoidal modes can be pictured as the modes of a turbulent incompressible field, made of vertices and eddies. On the other hand, compressive modes, made of compression and expansion motions, are more likely to be generated by phenomena linked to star formation.
The application of these methods implies several requirements on the studied dataset. As mentioned earlier, the statistical isotropy of the cloud is the first necessary point, and enables the use of 2dimensional averages as a means to estimate 3dimensional properties. It means that the method cannot be applied to individual filaments, or to clouds where a strong anisotropy is suspected, e.g., owing to the presence of a strong magnetic field at low Mach numbers.
Second, the field is required to go smoothly to zero on its borders. This property is needed to ensure that the decomposition of the field is unique, since the Helmholtz decomposition is, in theory, defined up to a vector constant. It is also a necessary condition for good behaviours of the Fourier transform, since actual observational data are not periodic fields, unlike hydrodynamical simulations (see discussion in Brunt et al. 2010). This implies that the studied field should be bounded in space like, for example, a gravitationally bound cloud. In the case where the signal extends up to the edge of the observed field, the dataset has to be apodized.
Finally, from a practical viewpoint, Brunt et al. (2010) have shown that their method works best for fields with power spectra that are not too steep. Steep power spectra give measurements that are very sensitive to the low spatial frequencies, which are usually uncertain owing to poor statistics.
The compliance of our dataset with these requirements is discussed in detail in Sect. 5.1.
3.1.2. Equations and notations
The studied quantity is the momentum density field (hereafter momentum), p = ρv, with ρ and v the volume density and the velocity.
The 3dimensional quantity we infer is (1)the ratio of the variance of the transverse (solenoidal) momentum to the variance of the total momentum. For short, R will be referred to as the solenoidal fraction in the rest of this paper.
According to Brunt & Federrath (2014), at hypersonic Mach numbers (M = v/c_{sound} > 5) the solenoidal fraction does not depend any more on the type of forcing, but instead converges towards R ~ 2/3. Brunt & Federrath (2014) note that this behaviour is different from what is observed by Federrath et al. (2011), where the solenoidal fraction converges to different values depending on the type of forcing, but this is due to the fact that Brunt & Federrath (2014) consider the momentum density field, while Federrath et al. (2011) describe the energy density field.
This value of R ~ 2/3 can be simply explained in terms of equipartition of momentum between the compressive and solenoidal modes (see, e.g., Federrath et al. 2008). A value of R lower than 2/3 therefore means that there is more momentum in the compressive modes of the flow, and that the cloud is thus more likely to form stars. The influence of the Mach number is further discussed in Sect. 4.1.
The available observables are a positionpositionvelocity cube, its velocityweighted moments and the power spectra of these moments. We make two major assumptions about this datacube, namely that the ^{13}CO(J = 1−0) line is optically thin and that its emissivity only depends on the ^{13}CO volume density. These assumptions are true within less than 20%, with the exception of a few lines of sight towards the center of NGC 2024, which are more saturated, but represent about 2% of the whole field. Under these assumptions, the positionpositionvelocity cube can be seen as a densityweighted velocity field: the spectrum obtained for each line of sight results from the projection of the emission of the matter present along this line of sight, and moving at various velocities.
The useful moments are the zeroth, first, and second order moments of the momentum field, W_{0}, W_{1}, and W_{2}, which are defined as follows in Brunt & Federrath (2014): (2)The spectral line intensity, I(v), may have contributions from various positions along the line of sight. Given our assumptions on emissivity, and assuming that the natural linewidth is negligible compared to the overall velocity dispersion, we can describe these moments in an alternative way: (3)The solenoidal fraction can be written in terms of these observational quantities as follows (see Brunt & Federrath 2014, for details): (4)The A and B factors are functions of the power spectra of W_{0} and W_{1}, (5)where f(k) and f_{⊥}(k) are the angular averages of the power spectra and , respectively.
Brunt & Federrath (2014) introduce a statistical correction factor, g_{21}, of order unity, which measures the correlations between the variations of the density and the velocity fields, and is defined as (6)Brunt & Federrath (2014) show that this may be written as (7)The 3dimensional variance of the volume density, ⟨ (ρ/ρ_{0})^{2} ⟩, can be derived using the Brunt et al. (2010) method. The exponent ϵ is a small, positive constant, which can be obtained as the exponent of the vs. ρ power law, i.e., the typical velocity dispersion in the cloud as a function of volume density. If the density and velocity fields were uncorrelated, g_{21} would be equal to 1.
3.2. Implementation
Actual data suffer from several limitations that need to be dealt with to apply the method described previously. The field of view, the angular resolution, and the sensitivity are limited. This section describes how these issues were dealt with.
3.2.1. Noise filtering
Computation of line moments is sensitive to noise in the line wings. It is well known that masking the positionpositionvelocity cube where the signal stays undetected improves the determination of the centroid velocity and linewidth. To define the mask containing the pixels detected at high significance, these pixels are first grouped into continuous brightness islands that are made of neighbouring pixels in the positionpositionvelocity space, whose S/N is larger than 2. The noise level, σ, is measured outside the studied velocity range (− 0.5, 20.5 km s^{1}). The list of islands is then sorted by decreasing total flux. The first island contains about 97.8% of the total signal. The following ones are small signal clumps that are spatially or spectrally isolated from this main block and, after about a few hundred islands, we are left with singlecell islands that are just noise peaks.
While it is easy to visually assess that the first few islands correspond to genuine signal, it is more complex to determine the transition to pure noise, as a significant fraction of the total line flux could be hidden in pixels of faint brightness, at low S/N. We thus studied the influence of the number of islands used on R, the solenoidal fraction in the studied cube. This influence is modest, mainly because over 97% of the signal is located in the first island. Using up to about 80 islands yields a very stable value of R with less than 0.1% variation. A steep increase of R is observed when we enter the noisy domain. We thus used the 80 brightest islands for all other calculations.
3.2.2. Moment computation
After selecting the signal islands in the positionpositionvelocity cube, the moments are integrated from −0.5 to 20.5 km s^{1}. The calculations have to be performed in the centerofmass frame of reference of the cloud, which implies determining the centroid velocity of the cloud in the LSR frame. This centerofmass velocity is simply given by (8)where is the first moment field in the observer’s frame of reference. W_{0}, on the other hand, is not velocityweighted and, therefore, does not depend on the frame of reference. For the observed field of view, we obtain V_{c} = 9.16 ± 0.90 km s^{1}.
The velocity scale in the observer’s frame of reference is shifted by V_{c}, before computing W_{1} and W_{2} in the centerofmass frame of reference: (9)The resulting fields are shown in Fig. 4.
Fig. 4 Top: maps of the ^{13}CO(J = 1−0) field moments in the cloud’s frame of reference. Bottom: maps of the physical quantities directly derived from each field, with the centroid velocity being simply W_{1}/W_{0}, and the Full Width at Half Maximum (FWHM) velocity dispersion given by – the normalization of the FWHM corresponds to that of a field with purely Gaussian line profiles. 

Open with DEXTER 
3.2.3. Apodization
Computing the power spectra of the W_{0} and W_{1} fields requires to take the Fourier transform of these fields. Although calculating the FFT of 2dimensional fields is an easy task, several numerical artefacts must be taken care of. In particular, the observed area does not reach the edges of the Orion B molecular cloud in all directions, as illustrated in Fig. 1. This sharp truncation of the ^{13}CO emission will create artefacts in the Fourier transform, owing to the convolution of the true Fourier spectrum by a sinus cardinal function that oscillates at high frequencies. Apodizing the field is required to avoid this behaviour. We have chosen to multiply the intensity by 1−cos(πx/w), where x is the pixel coordinate and w the apodization width. This function goes from 0 to 1 over w pixels. This apodization function is used in Martin et al. (2015). We keep the region affected by apodization as small as possible to minimize signal alteration. An apodization width of about 5% of the smallest dimension of the field (i.e., roughly 25 pixels) was the smallest value that efficiently smoothed out the highfrequency artefacts. This is consistent as well with the width determined in Martin et al. (2015).
Fig. 5 Left: powerlaw fitting of the W_{0} power spectrum. Upper panel: data (blue crosses), fit result (thick solid line) plotted over the fitted domain, power law convolved with the Gaussian beam (dotted line) extrapolated to all spatial frequencies, noise model (dashed line). Lower panel: residuals. The scale in the Fourier space is given in UV distance as for interferometric observations, which enables us to visually relate the resolution and the telescope diameter. Right: same results, except for the W_{1} power spectrum. 

Open with DEXTER 
Once the field is apodized, it must be made square to follow the isotropy requirements of the method. The square was built by padding the right side (west) of the observed area with zeros, since this is the location of the Hii region associated with σ Ori, i.e., no signal is detected past the western edge of the field of view. After this, the Fourier transform is calculated using the FFT implementation of Numpy 1.8 (Cooley & Tukey 1964).
Apodization is a linear filter of the data, and thus has effects on the power spectrum at all frequencies. While apodization allows us to clean the spectrum at high spatial frequencies, it also alters the spectrum at the lowest frequencies. For example, f(0) – the value of the power spectrum at spatial frequency k = 0 – is directly proportional to the spatial integral of , therefore losing some signal because of the apodization will reduce the value of f(0). The method chosen to keep the good parts of the apodized and nonapodized spectra was the following: we first computed the FFT of both the apodized and nonapodized fields, respectively and , then mixed them smoothly around k = 0 using a narrow (about 10 Fourierspace pixels) 2dimensional Gaussian G_{mix}(k). The resulting Fourier field is (10)The power spectrum thus behaves like the nonapodized spectrum at low k, keeping the correct value of the field integral f(0), and like the apodized spectrum at high k, free of the spectral parasites created by the sharp edges of the map.
3.2.4. Power spectra computation and fit
The apodized and corrected FFT of the field needs to be transformed into an angleaveraged power spectrum f(k) or f_{⊥}(k). This is simply done by binning the modulus of the spatial frequencies, and averaging the points found in these radial bins. The resulting discrete function can then be linearly interpolated into a continuous function. A critical element in this exercise resides in the sampling of the spatial frequency axis. On the one hand, the resulting angleaveraged spectrum should be as detailed as possible but, on the other hand, a larger number of bins can lead to empty bins, containing no sampled points at all. As a result, we used a number of bins of S/ 1.45, with S the size in pixels of the square field, so that the size of a bin in the Fourier space corresponds to slightly more than the length of the diagonal of pixels in the Fourier space.
Additional observational constraints (noise, beam shape, etc.) affect the determination of the power spectrum. Following Martin et al. (2015), the power spectra are fitted with a power law, modified to take into account the singledish beam and the noise. In our case, the beam is modelled as the Fourier transform of a Gaussian of FWHM equal to the cube resolution, i.e., 23.5′′. This corresponds to about 2.61 pixels. The convolution in the image space corresponds to a multiplication in the Fourier space that mostly affects the high spatial frequencies.
The noise is not a Gaussian white noise, because of interpixel correlations and systematics. We therefore use the power spectra of 30 signalfree channels and average them to obtain a template of the noise power spectrum. In this case, we use the fully noisy data cube, not the 80 first signal islands, to have the same spatial correlations in noise for each channel (with or without signal), which would not be the case with a masked data cube. The noise template intends to reproduce a systematic behaviour, but we can only use a finite number of channels with random noise. The noise template is therefore smoothed to make this systematic pattern stand out more. Given that W_{1} is, just like W_{0}, a linear combination of the channel maps, the same noise template is used for both power spectra.
The fit is performed in the log(k)log(f) space, so that the straight line of the power law stands out more. The fitting function is therefore the logarithm of (11)where G_{beam} is the Gaussian beam, is its (Gaussian) Fourier transform, and noise is our noise template. The fitted parameters are a, b, and N. During the fit, the data are weighted by the inverse of the variances obtained in each bin when computing the power spectrum.
The choice of fitting the power spectrum with a modified power law implies that the underlying physical processes should produce a power law. This is indeed the case for the inertial range of scales in Kolmogorov turbulence, and can be applied as well to Burgers turbulence (see, e.g., Federrath 2013). However, the power spectrum of turbulence starts to deviate from a power law at scales where the energy is injected (low spatial frequencies) or dissipated (high spatial frequencies). The power spectra computed from our dataset somewhat deviates from a power law at low spatial frequencies (see Fig. 5). The power law can therefore only be fitted and then used above a given spatial frequency. The powerlaw range starts around ~ 5.5′. At high spatial frequencies, no deviations from the power law are detected (the fits render the observations very well). This means that either the dissipation scale happens at lower angular scale than the resolution of the observations or it is hidden in the noise.
Once the fit has been performed, the final version of the power spectra is built using both the fit result and the observational angleaveraged power spectra, and used for further calculations. The final power spectrum is a pure power law (without the beam and noise) above the 1/5.5 arcmin^{1} threshold, and is equal to the linearly interpolated angleaveraged power spectrum below this threshold. In particular, we enforce f_{⊥}(0) = 0, since we are working in the cloud’s rest frame.
3.3. Densityvelocity correlations
We used the information on the mean line profiles to estimate the slope ϵ of the relation between the velocity dispersion and the local density (see Eq. (7)). Among the lines detected in the mean spectrum, we selected lines with different spatial distributions (Pety et al. 2017). To trace the low density gas, we selected ^{12}CO(J = 1−0) and HCO^{+}(J = 1−0) since these lines present very extended emission and have moderate excitation requirements (Pety et al. 2017; Liszt & Pety 2016). We included ^{13}CO(J = 1−0) as our tracer of the bulk of the gas. The somewhat denser and more shielded gas is well traced by C^{18}O(J = 1−0), while we selected N_{2}H^{+}(J = 1−0) for the dense cores. For these five species, we determined the FWHM by fitting a Gaussian line profile to the mean profile of the whole map. Only the 10 km s^{1} component, which is present for all five species, was used for this fit. For N_{2}H^{+} we used the HFS fit method in GILDAS/CLASS^{2}, which makes use of the information on the hyperfine structure.
While these lines are emitted by gas over a wide range of densities, there is a minimum density under which the line is not detected because of lack of excitation or because the molecule is not present in low density gas. It is this minimum density that corresponds to the velocity dispersion of the line. To derive the densities associated with the line emission, we adopted three different methods.
For the low density and extended emission tracers, we derived the volume density by comparing the minimum gas column density where the emission is detected and the resulting size of the emission regions. This leads to gas densities of a few hundred cm^{3} for ^{12}CO(J = 1−0), HCO^{+}(J = 1−0) and ^{13}CO(J = 1−0). The emission of ^{12}CO(J = 1−0) is dominated by the low density regions (Pety et al. 2017). The emission of HCO^{+}(J = 1−0) is also dominated by low to intermediate density regions, and comes from the weak excitation regime (Liszt & Pety 2016). In both cases, opacity broadening is not very significant. We must, however, consider that the line widths for these tracers are upper bounds, owing to the effect of opacity broadening. For these molecules, we are also limited by the sensitivity of the observations, so that the densities should be regarded as upper limits. HilyBlant et al. (2005) analysed the structure and kinematics of the Horsehead nebula and derived the density of the extended region traced by C^{18}O as 3−5 × 10^{3} cm^{3}. We kept this value as the typical density traced by C^{18}O. HilyBlant et al. (2005) show that higher density regions exist with densities significantly larger than 10^{4} cm^{3}. As very few pixels are detected in N_{2}H^{+} in the region of the Horsehead nebula, we used the catalogue of dense cores identified by Kirk et al. (2016) in their SCUBA2 map of the Orion B complex. We found 55 cores associated with N_{2}H^{+}(J = 1−0) emission. The densities were derived using the extracted fluxes and effective radii, a uniform dust temperature of 20 K, and assuming spherical geometry for all cores. The mean density is 10^{4.1} cm^{3} with a scatter of about a factor of two. The temperature of the N_{2}H^{+} cores was not individually checked, but it is likely to be lower than 20 K. In turn, this implies an even higher density of N_{2}H^{+}. Therefore, the derived density should be regarded as a lower limit.
Table 1 presents the resulting data, which are illustrated in Fig. 6. The slope α = −ϵ is derived from a least squares fit of the variation of the FWHM with the density. We derive ϵ = 0.15 ± 0.03. The possible systematics on the densities traced by ^{12}CO(J = 1−0), HCO^{+}(J = 1−0), and N_{2}H^{+}(J = 1−0), as well as on the linewidths of ^{12}CO(J = 1−0) and HCO^{+}(J = 1−0), all tend to make the power law steeper. Therefore, we keep this value of ϵ as an upper limit.
Spectral tracers used in this study, and observed with the IRAM30 m telescope.
Fig. 6 Variation of the FWHM as a function of the gas density. Each point refers to the (J = 1−0) spectral line of a different molecule. The red line shows the least squares fit, with a slope α = −ϵ = −0.15. 

Open with DEXTER 
3.4. Estimation of the uncertainties
The computation of uncertainties implies computing the uncertainty of each element of Eq. (4). We start from the average rms noise level in the data cube, 0.17 K[T_{mb}]. This enables us to compute the noise level for the W_{0}, W_{1}, and W_{2} maps. The computation is straightforward compared to the computation of the uncertainty of the centroid velocity and linewidth because W_{1}, and W_{2} are not normalized by W_{0}, i.e., their noise distributions stay Gaussian, whatever the value of W_{0}. Due to the velocity weighting, the absolute uncertainty increases significantly with the moment order. However, the relative uncertainties on ⟨ W_{0} ⟩, ⟨ W_{1} ⟩ and ⟨ W_{2} ⟩ are similar, with values ranging from about 10% for the whole field to 5% for the deepest zooms on NGC 2023 and NGC 2024 (see Fig. 8). The uncertainties on the sums A and B were explicitly computed according to the error bars described in Sect. 3.2.4. The relative uncertainty on A ranges from 24% for the full field to 11% for the deepest zooms, and stays around 13% for B.
For the overall relative uncertainty, one must not only take into account the errors of the individual terms, but also the correlation between the different variables. In our case, the different variables are strongly correlated, since they are all byproducts of the same data cube. Therefore, we chose to use the average of the various relative errors as a ruleofthumb estimate of the overall error. This approach yielded an overall relative error of 13% for the whole field, and 8% for the deepest zooms. We kept the highest value to allow for a safety margin.
This 13% relative error ΔR/R corresponds to the median noisetosignal ratio in the field, σ/T_{peak}, which is also of the order of 13% – but testing with simulated noise whether this is a coincidence or not is out of the scope of this article.
4. Results
In this section, we first briefly present the derivation of the Mach number in the cloud, then we compare the obtained power spectra with other results in molecular clouds, and finally, we give the results of our computation of the solenoidal fraction R in the Orion B cloud.
4.1. Mach number
According to Brunt & Federrath (2014) and Federrath et al. (2011), at hypersonic Mach numbers (M> 5) the ratio of solenoidal and compressive modes does not depend any more on the type of forcing. It is thus important to also derive the distribution of the Mach number to be able to interpret the results.
Using the maps of the sound speed derived from ^{12}CO(J = 1−0) and dust temperature maps, two different estimations of the Mach number, M_{max} and M_{exc}, were computed (see appendix for details). Figures A.2 and A.3 show their spatial distributions and compare their histograms. The shapes of the histograms of Mach numbers computed with T_{max} and T_{exc} are very similar, and both show a large tail at high Mach numbers. Table A.1 lists several characteristic values of both distributions. The most probable value (~ 3.5) of the Mach number is much smaller than the mean or median values (~ 6), for both distributions.
Schneider et al. (2013) estimate the average Mach number to be of ~ 8, with approximately 30−40% error, deriving this value from Herschel dust temperature and CO linewidth. They also find that Orion B has the highest Mach number of the set of studied clouds. Our results are compatible with this mean value, and they also provide the spatial and statistical distribution at good angular resolution. In particular, the Mach number is much smaller than the average in the star forming regions NGC 2023 and NGC 2024, below the hypersonic regime.
4.2. Power spectra
When the whole field is considered, the fit yields an exponent a_{0} = −2.83 ± 0.02 for the W_{0} field, and a_{1} = −2.50 ± 0.07 for the W_{1} field. When zooming into specific regions of the cloud (NGC 2023 and NGC 2024, see Fig. 8), the values of these exponents range from a_{0} = −2.52 ± 0.08 (widest field) to a_{0} = −3.04 ± 0.05 (smallest field), and from a_{1} = −2.24 ± 0.09 to a_{1} = −2.81 ± 0.06. While the indices of W_{1} power spectra are rarely reported, there are many values of spectral indices for integrated line intensity maps in the literature, and our a_{0} values fall well in the range of spectral indices for observations of CO emission, dust emission, Hi emission and absorption compiled by Hennebelle & Falgarone (2012): the values range from −2.5 to −3.2, with most values around −2.7.
To have a meaningful result for the A and B coefficients, the powerlaw slope of the power spectra must follow the two steepness requirements mentioned by Brunt et al. (2010). On the one hand, the spectra should not be too steep, so that the weight of the low spatial frequencies, for which the available information is scarce, does not become too large in the A and B sums. On the other hand, the slopes should be steep enough to avoid divergence of these sums at large frequencies. A slope of a = −3 is at the limit between these two contradicting constraints, since the divergence of the 3dimensional sum is only logarithmic. In our case, the sums are finite, owing to the finite resolution of the observations, so that with slopes between −2.24 and −3.04, the sums do not grow too quickly and do not give too much weight to the low spatial frequencies.
4.3. Turbulence mode ratio
We determined a relative error of about 13% on the calculation of the ratio R from the positionpositionvelocity cube. The correction factor g_{21} is determined independently, and we assume a range of possible values for g_{21}, the lower limit being given by our calculations of Sect. 3.3, and the upper limit resulting from the minimum estimate of ϵ ≃ 0.05 according to Brunt & Federrath (2014).
For the entire ^{13}CO field, we obtain the following range of values: (12)To gain a deeper understanding of the dynamics at stake in the Orion B cloud, the method was also applied to several subregions of the ^{13}CO map.
First, to check the reliability of the method and the homogeneity of the field, we applied a sliding square window, whose side is equal to the smallest dimension of the mapped area (Fig. 7). This avoids zeropadding the studied field. The results show that, even though the values are in general somewhat lower for these subfields than for the whole ^{13}CO field, they remain marginally compatible with this result, within the estimated uncertainties.
Fig. 7 Solenoidal turbulence fraction for sliding square areas with a width equal to the one of the full map (56 arcmin). The shaded area corresponds to the g_{21} uncertainty, while the error bars show the experimental uncertainties due to observational noise. The vertical line marks the equipartition limit. The map on the left presents the centres of the square areas for which the calculations were performed, superimposed on the ^{13}CO(J = 1−0) T_{peak} map. 

Open with DEXTER 
Second, we searched for systematic variations of the fraction of solenoidal modes when zooming into specific regions of the map. In particular, signs of the solenoidal or compressive forcing are expected to appear mainly in regions of low Mach numbers (see Sect. 4.1). The zooms were thus performed into the NGC 2023 and NGC 2024 star forming regions (Fig. 8), where the Mach number lies between about 3 and 5 (see Fig. A.2), mostly because the speed of sound is higher in regions of higher gas temperature, but also because the velocity dispersion is a bit lower.
Moreover, these regions offer the advantage of presenting a strongly localized emission. One of the requirements of the Brunt & Federrath method is to use isolated fields. This is clearly not the case any more when zooming into these specific regions, and apodization is more necessary than ever to ensure that the signal falls smoothly to zero on the edges of the field. However, by using fields for which most of the signal is concentrated near the center, the effects of apodization are minimized, allowing us to stay as close as possible to the requirement of an isolated field. The smaller the field, the lower the value of R: this indicates an increasing proportion of compressive forcing.
Fig. 8 Solenoidal turbulence fraction for zooms on the NGC 2023 and NGC 2024 star forming regions. The shaded area corresponds to the g_{21} uncertainty, while the error bars show the experimental uncertainties due to observational noise. The upper limit of the plots marks the equipartition limit. The map on the left presents the areas for which the calculations were performed (solid squares: NGC 2024, dashed squares: NGC 2023), superimposed on the ^{13}CO(J = 1−0) T_{peak} map. 

Open with DEXTER 
5. Discussion
5.1. Compliance with the method’s assumptions
To apply the Brunt & Federrath method to real observational data, we were able to overcome several difficulties and sources of uncertainty.
First, the compliance of the dataset with the requirements of the method must be checked, namely the isotropy of the studied cloud, and its isolation. The whole field and the zooms present two opposite situations. The isolation criterion is well met in the case of the whole field: we have almost no signal to the west and to the south of the field, and very diffuse regions to the north and to the east (see Fig. 4, first column), so that there is almost no need for apodization to have the signal going down to zero on the edges of the field. For the zooms, on the other hand, we are well into the cloud, so that there is signal all the way to the edges of the field. However, since we study local maxima of the emission, most of the signal is in the center, which minimizes the effects of the necessary apodization on the final results (only 7.9% of lost signal for the deepest NGC 2023 zoom, 7.6% for NGC 2024).
As far as the isotropy criterion is concerned, the 2D projection is quite obviously isotropic in the case of the zooms, since the considered fields are square, and much less in the case of the whole field, in which the region with signal has an aspect ratio of about 2:1. The third dimension is unknown, and in any case cannot match simultaneously the dimensions of the whole field and those of the deepest zooms. However, for the large diffuse regions like for the bright, compact regions, the dimension along the line of sight is supposed to be of the order of the dimensions in the plane of the sky. In the case of a zoom on a bright region embedded in a diffuse one, the signal and, therefore, the dimension on which it is emitted along the line of sight, is dominated by the bright gas. Therefore, since the zooms are centred on a bright region, the corresponding datacube behaves almost as if the bright region was isolated and isotropic. Besides, the nonangleaveraged power spectra of W_{0} and W_{1} do not show any apparent anisotropy at any scale (except for the windowing effects). If the power spectra are isotropic in two dimensions, then statistically we can expect the third dimension to follow this isotropy as well. Therefore, the isotropy requirement seems to be fulfilled by our dataset and the method can be applied.
Second, as was mentioned by Brunt & Federrath (2014), the method is sensitive to values of the power spectra at large spatial scales (low frequencies of the power spectra) owing to the characteristics of the sums in the parameters A and B. We therefore had to find a way to obtain a smooth and reliable function that would represent an angleaveraged power spectrum at all frequencies. Once the power spectra were binned and fitted, we have two versions of the spectra, each with its flaws. The binned (data only) spectrum suffers from observational effects (noise and beam) but has also larger uncertainties at low spatial frequencies. The fitted spectrum, if extrapolated to all spatial frequencies (fit only), can give unphysical results in the lowest frequencies because they lie outside the power law validity range. For example, if the field had been zeropadded all around to create a very large square, such an extrapolation would give very high values of the spectrum at low frequencies, whereas physically they should be very low, since the field would on average be almost empty. These flaws led us to choose the composite scheme described in Sect. 3.2.4, where the fit result is used only in the powerlaw domain, and the low frequencies keep the angleaveraged power spectrum as it is. Using a different version of the power spectra leads to quite different final results, as illustrated in Fig. 9. However,we note that, even though the absolute value of the solenoidal fraction R varies, the relative variations are consistent across scales, whatever the power spectra computing scheme. Thus, results such as the unusually high solenoidal fraction at the scale of the whole cloud, or the variations of R when zooming out of the star forming regions, stay valid and are further discussed below.
Fig. 9 Comparison of the calculation results in the case of the zoom on NGC 2024 with three different methods of computing the power spectra: either the power spectra resulting from Fouriertransformed data are used directly, with a linear interpolation between the points (data only, blue), or the result of the power law fit at high frequencies is extrapolated to low frequencies to provide a power law throughout the whole range of frequencies (fit only, green), or a composite version of the power spectrum is built using the raw data at low frequencies, and the fit result at high frequencies (data + fit, red). The error bars account for the observational noise as well as for the g_{21} uncertainty. The horizontal line marks the equipartition limit. 

Open with DEXTER 
5.2. Physical interpretation
To our knowledge, this work is the first attempt at applying the Brunt & Federrath method to actual observational data (Brunt & Federrath 2014; Lomax et al. 2015). The results need to be compared with what has been done so far on numerical simulations.
Both Schneider et al. (2013) and our calculations (see Sect. 4.1 and appendix) show that we are in a context of highly supersonic turbulence, with a mean Mach number of about 6. We therefore expect a full mixing of the turbulence modes, so that the momentum equipartition would predict a solenoidal fraction of R = 2/3 and a compressive fraction of 1−R = 1/3. Deviations from this ratio can either be a sign of a specific forcing for the turbulence in the case of transonic Mach numbers, as shown in the case of simulations (Brunt & Federrath 2014), or indicate that an ordered flow is superimposed on top of the turbulent flow.
The global value of R> 0.72 ± 0.09 and the values in Fig. 7 can agree, within the error bars, with the expected value R = 2/3 in the case of equipartition. However, the value R> 0.72 ± 0.09 is still quite high. This can be the sign of a deviation in favour of solenoidal modes. At high Mach numbers, it would imply that an ordered solenoidal flow is superimposed on top of the turbulence. And Fig. 4 (Col. 2) shows that there is a largescale differential motion in the cloud with the southern part receding, while the northern part is approaching. This velocity shift could be the sign of a largescale rotation of the whole cloud, which could dominate the smallerscale motions and be responsible for the large fraction of solenoidal modes.
The fact that the turbulence in the Orion B molecular cloud is, on large scales, mostly solenoidal, is in agreement with the fact that, for its mass, Orion B has a low star formation rate (Lada et al. 2010). Estimations of its SFE vary, with values ranging from 0.4% to 3% (Lada 1992; Carpenter 2000; Federrath & Klessen 2013; Megeath et al. 2016), but all studies show that the SFE in Orion B is about four times lower than in the neighbouring cloud Orion A. In general, Orion B’s SFE is regarded as particularly low, with Megeath et al. (2016) stressing that is has the lowest SFE among all molecular clouds closer than 500 pc. This remarkable feature of the Orion B cloud could be partially explained by the solenoidal flows that drive its velocity field, and hinder collapsing motions that could trigger star formation.
In contrast, Fig. 8 shows a major deviation from the equipartition in favour of compressive motions in two specific regions. When zooming deeply into the star forming regions NGC 2023 and NGC 2024, we obtain solenoidal fractions as low as R = 0.25. The high fraction of compressive flow in these two regions most likely results from the infall of matter onto the star forming region, and/or the expansion of the Hii regions around the young massive stars (Tremblin et al. 2014a; Geen et al. 2015). The Hii regions themselves are not observed in molecular tracers, so that this expansion is detected indirectly through the compression of the molecular gas at the ionization front of the Hii regions.
For both regions, R grows when zooming out, and will eventually reach the average values displayed on Fig. 7. R tends to decrease when reaching a field size of about 20′. This behaviour can be due to the fact that the other star forming region is entering the field of view, since the distance between the two cores is about 22′, but it can also be related to the geometry of each region.
In the case of NGC 2023, this size of 20′ also corresponds to the PDR of the Horsehead Nebula coming into the field and, owing to the pressure at the photodissociation front, it is expected to be a compressive region (WardThompson et al. 2006), which is proven by the detection of at least one young star and one protostar in this region (Megeath et al. 2012).
In the case of NGC 2024, we see that R starts increasing again after 25′, even though NGC 2023 comes more and more into the field. The variations of R around 20′−25′ might therefore be related to the location of the edge of the Hii region around NGC 2024. This Hii bubble exerts a pressure on the surrounding gas (Tremblin et al. 2014a,b), and this edge is therefore a highly compressive region that can be seen in the form of an arc north of NGC 2024 (Megeath et al. 2012). This region has a far lower surface density of stars that the inner part of NGC 2024, but it is likely to be younger that NGC 2024 (since it is a consequence of the expansion of the Hii region), and therefore might have formed only very young protostars, poorly detected by Spitzer (Megeath et al. 2016), or no stars yet – but it might become a very active region in the future.
The sharp contrast between the largescale solenoidal flow and the highly compressive flows in the NGC 2023 and NGC 2024 nebulae is in agreement with the spatial distribution of star formation observed in Orion B: Lada (1992) observes that the largescale SFE in the Orion B cloud is an order of magnitude lower than in the most massive cores, Carpenter (2000) shows that, at his detection level, 100% of young stars in Orion B are located in clusters (NGC 2068 in the north, and NGC 2024 in the south), which is unusual among the studied clouds, and Lada et al. (2013) conclude that Orion B is very ineffective at forming stars at A_{K}< 2.0 mag, compared to other GMCs.
In a broader perspective, not only do these variations of R confirm observationally the intuitive link between compressive motions and star formation, as proposed in simulations (Federrath & Klessen 2012; Padoan et al. 2014), but they also show that there can be an intracloud variability of the solenoidal fraction, in addition to the intercloud one. This shows that the largescale environment of the cloud, although it plays a major role in driving the turbulence of the molecular cloud, cannot alone explain the repartition of solenoidal and compressive motions in the cloud: any denser region created by the density fluctuation in the compressible turbulent gas (e.g., Nolan et al. 2015, and references therein) can lead, under the effect of selfgravity or stellar feedback, to the formation of very localized, strongly compressive regions, even in the context of a mostly solenoidal flow. There is therefore no universal solenoidal fraction that can be applied generally to all clouds, and there are even intrinsic variations of R from region to region within a cloud.
6. Conclusion
From a practical point of view, our work has shown that it is possible to apply the numerical method of Brunt & Federrath (2014) to observational data to determine the fraction of solenoidal and compressive motions in a molecular cloud, using molecular lines as tracers of the density and velocity fields. We were able to pinpoint the observational requirements to apply this method to a dataset.
The spatial dynamic range is an important element, mostly to provide good quality power spectra. The field must be large enough to provide good statistics at low spatial frequencies, but it must also have a good spatial resolution, so that the power spectra have enough points to correct properly for the beam and noise effects. We found our minimum field size to be of at least 50 independent pixels.
In addition to the spatial resolution, having many independent spectral channels is of great help when correcting for the beam and noise effects. The spectral resolution must also be sufficient to resolve the studied spectral line.
The S/N also proved to be a key element during the calculations. An average S/N of at least 5 is desirable: our datacube, which has a mean S/N of 7.8, yielded a relative observational uncertainty of 13% on the fraction of momentum in the solenoidal modes, and Figs. 7 and 8 show that this observational uncertainty contributes significantly to the overall uncertainty, and is even dominant at low solenoidal fractions.
The last point of the computation – the densityvelocity correlation – which also significantly contributes to the overall uncertainty, requires us to use many spectral tracers of various typical densities. In that case, a spectral survey such as the one of the Orion B project is invaluable, in so far as all the needed tracers are available and intercalibrated.
From a physical point of view, the measurements have shown that the motions in the Orion B molecular cloud are highly supersonic, with a mean Mach number of ~ 6. However, the Mach number maps show large variations, with some regions being only moderately supersonic. These variations are due both to the variations of the temperature and to the turbulent velocity distribution in the cloud.
The largest scales of the cloud seem to be dominated by a rotational motion, which can be identified by a high solenoidal fraction in the flow. At smaller scales, we have shown that the motion is largely dominated by compressive (infall and/or outflow) motions in the vicinity of the NGC 2023 and NGC 2024 star forming regions. The northern edge of NGC 2024 and the photodissociation front of the Horsehead nebula are also likely to be highly compressive regions, according to our results. This is in agreement with the observations of the star formation efficiency in Orion B, which is unusually low at the scale of the whole cloud, and exclusively concentrated in clusters (NGC 2023 and NGC 2024).
The example of Orion B also shows that the star formation efficiency in a molecular cloud does not only depend on its overall fraction of momentum in the solenoidal modes of turbulence, but also on the local variations of this fraction, which can be driven by internal phenomena such as selfgravity and stellar feedback.
This method could be applied in the future to study the variations of the solenoidal fraction between different molecular clouds, or between different regions or different chemical tracers within a given cloud. In the case of Orion B, we intend to analyse other data cubes for tracers such as C^{18}O or HCO^{+}, which trace, respectively, more compact and more diffuse regions of the cloud, to probe different environments in terms of density and temperature.
The data products associated to this paper are available at http://www.iram.fr/~pety/ORIONB
See http://www.iram.fr/IRAMFR/GILDAS/ for more details on the GILDAS software.
Acknowledgments
We wish to thank Patrick Hennebelle, Edith Falgarone, and Pierre Lesaffre for fruitful conversations on hydrodynamics and magnetohydrodynamics, as well as our referee for enlightening remarks on compressible turbulence in the interstellar medium. This research also made use of data from the Herschel Gould Belt survey (HGBS) project (http://gouldbeltherschel.cea.fr). The HGBS is a Herschel Key Programme jointly carried out by SPIRE Specialist Astronomy Group 3 (SAG 3), scientists of several institutes in the PACS Consortium (CEA Saclay, INAFIFSI Rome and INAFArcetri, KU Leuven, MPIA Heidelberg), and scientists of the Herschel Science Center (HSC). This work was supported by the CNRS/CNES program Physique et Chimie du Milieu Interstellaire (PCMI). We thank the CIAS for its hospitality during the two workshops devoted to this project. V.V.G. thanks for support from the Chilean Government through the Becas Chile program. P.G.’s postdoctoral position was funded by the INSU/CNRS. P.G. thanks ERC starting grant (3DICE, grant agreement 336474) for funding during this work. NRAO is operated by Associated Universities Inc. under contract with the National Science Foundation.
References
 André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brunt, C. M., & Federrath, C. 2014, MNRAS, 442, 1451 [NASA ADS] [CrossRef] [Google Scholar]
 Brunt, C. M., Federrath, C., & Price, D. J. 2010, MNRAS, 403, 1507 [NASA ADS] [CrossRef] [Google Scholar]
 Carpenter, J. M. 2000, AJ, 120, 3139 [NASA ADS] [CrossRef] [Google Scholar]
 Cooley, J. W., & Tukey, J. W. 1964, Math. Comp., 19, 297 [CrossRef] [MathSciNet] [Google Scholar]
 Dickman, R. L. 1978, ApJS, 37, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton University Press) [Google Scholar]
 Federrath, C. 2013, MNRAS, 436, 1245 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., & Klessen, R. S. 2013, ApJ, 763, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., Klessen, R. S., & Schmidt, W. 2008, ApJ, 688, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., Klessen, R. S., & Schmidt, W. 2009, ApJ, 692, 364 [NASA ADS] [CrossRef] [Google Scholar]
 Federrath, C., RomanDuval, J., Klessen, R. S., Schmidt, W., & Mac Low, M.M. 2010, A&A, 512, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Federrath, C., Chabrier, G., Schober, J., et al. 2011, Phys. Rev. Lett., 107, 114504 [NASA ADS] [CrossRef] [Google Scholar]
 Geen, S., Hennebelle, P., Tremblin, P., & Rosdahl, J. 2015, MNRAS, 454, 4484 [NASA ADS] [CrossRef] [Google Scholar]
 Goldsmith, P. F. 2001, ApJ, 557, 736 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Gratier, P., Bron, E., Gerin, M., et al. 2017, A&A, 599, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennebelle, P., & Falgarone, E. 2012, A&ARv, 20, 55 [NASA ADS] [CrossRef] [Google Scholar]
 HilyBlant, P., Teyssier, D., Philipp, S., & Güsten, R. 2005, A&A, 440, 909 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 HilyBlant, P., Falgarone, E., & Pety, J. 2008, A&A, 481, 367 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kim, C.G., & Ostriker, E. C. 2015, ApJ, 802, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Kirk, H., Di Francesco, J., Johnstone, D., et al. 2016, ApJ, 817, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Kramer, C., Stutzki, J., & Winnewisser, G. 1996, A&A, 307, 915 [Google Scholar]
 Lada, E. A. 1992, ApJ, 393, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Lada, C. J., Lombardi, M., & Alves, J. F. 2010, ApJ, 724, 687 [NASA ADS] [CrossRef] [Google Scholar]
 Lada, C. J., Lombardi, M., RomanZuniga, C., Forbrich, J., & Alves, J. F. 2013, ApJ, 778, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, R. B. 1981, MNRAS, 194, 809 [NASA ADS] [CrossRef] [Google Scholar]
 Liszt, H., & Pety, J. 2016, ApJ, 823, 124 [NASA ADS] [CrossRef] [Google Scholar]
 Lomax, O., Whitworth, A. P., & Hubber, D. A. 2015, MNRAS, 449, 662 [NASA ADS] [CrossRef] [Google Scholar]
 Lombardi, M., Bouy, H., Alves, J., & Lada, C. J. 2014, A&A, 566, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mac Low, M.M. 1999, ApJ, 524, 169 [NASA ADS] [CrossRef] [Google Scholar]
 Mangum, J. G., & Shirley, Y. L. 2015, PASP, 127, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, P. G., Blagrave, K. P. M., Lockman, F. J., et al. 2015, ApJ, 809, 153 [NASA ADS] [CrossRef] [Google Scholar]
 Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2012, AJ, 144, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2016, AJ, 151, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Menten, K. M., Reid, M. J., Forbrich, J., & Brunthaler, A. 2007, A&A, 474, 515 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MivilleDeschênes, M.A., Duc, P.A., Marleau, F., et al. 2016, A&A, 593, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nolan, C. A., Federrath, C., & Sutherland, R. S. 2015, MNRAS, 451, 1380 [NASA ADS] [CrossRef] [Google Scholar]
 Padoan, P., Federrath, C., Chabrier, G., et al. 2014, in Protostars and Planets VI, 77 [Google Scholar]
 Pety, J., & Falgarone, É. 2000, A&A, 356, 279 [NASA ADS] [Google Scholar]
 Pety, J., Guzmán, V. V., Orkisz, J. H., et al. 2017, A&A, 599, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Renaud, F., Bournaud, F., Kraljic, K., & Duc, P.A. 2014, MNRAS, 442, L33 [NASA ADS] [Google Scholar]
 Ripple, F., Heyer, M. H., Gutermuth, R., Snell, R. L., & Brunt, C. M. 2013, MNRAS, 431, 1296 [NASA ADS] [CrossRef] [Google Scholar]
 Rohlfs, K., & Wilson, T. L. 2004, Tools of radio astronomy (Berlin: Springer) [Google Scholar]
 Schlafly, E. F., Green, G., Finkbeiner, D. P., et al. 2014, ApJ, 786, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, N., André, P., Könyves, V., et al. 2013, ApJ, 766, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Tremblin, P., Anderson, L. D., Didelon, P., et al. 2014a, A&A, 568, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tremblin, P., Schneider, N., Minier, V., et al. 2014b, A&A, 564, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 WardThompson, D., Nutter, D., Bontemps, S., Whitworth, A., & Attwood, R. 2006, MNRAS, 369, 1201 [NASA ADS] [CrossRef] [Google Scholar]
 Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Mach number estimation
In this appendix, we describe the computations that yielded the Mach number maps that enable us to estimate if the studied regions of the molecular clouds were below or above the hypersonic limit, which is characterized by the equipartition of momentum in the solenoidal and compressive modes of turbulence. We first compute the speed of sound using two different temperature maps, then derive the speed of the flow from the ^{13}CO(J = 1−0) positionpositionvelocity cube, and finally obtain statistics on the 3dimensional Mach number.
Appendix A.1: Speed of sound
The speed of sound in the gas needs to be derived from a temperature map, by means of , where γ is the adiabatic index, R_{gas} the universal gas constant and m_{mol} the molar mass of the gas.
We have direct access to two temperature maps: the dust temperature, computed by the Herschel Gould Belt Survey consortium (Schneider et al. 2013), and the peak temperature of the ^{12}CO map, from our IRAM30 m observations. The dust temperature is expected to be a lower bound for the kinetic temperature, because dust radiates more efficiently than gas. The gas and dust become coupled only for densities larger than 10^{4} cm^{3} (Goldsmith 2001).
For optically thick gas, the excitation temperature T_{exc} is determined (e.g., in Rohlfs & Wilson 2004) by (A.1)where ν = 115.271202 GHz is the frequency of the ^{12}CO(J = 1−0) transition, and T_{CMB} = 2.728 K is the cosmic microwave background temperature. We assume that the excitation temperature is close to the kinetic temperature of the gas (local thermodynamic equilibrium).
We assume that at very low column density (and therefore low T_{peak}), the kinetic temperature is underestimated (since the computation is only valid for optically thick gas). The dust temperature is supposed to be a lower limit close to the kinetic temperature of the gas, except when it is not deemed reliable any more: above a threshold of 60 K, we deem that Herschel and Planck observations do not yield the best temperature estimate due to their limited wavelength ranges (50−600μm for Herschel). At these temperatures, ^{12}CO(J = 1−0) is usually saturated enough to give a good result for the excitation temperature, as can be seen from the ^{12}CO(J = 1−0)/^{13}CO(J = 1−0) line ratio which diverges significantly from the ^{12}C/^{13}C abundance ratio (Pety et al. 2017; Ripple et al. 2013). We can therefore construct a third temperature map T_{max} using the dust temperature and the
excitation temperature of ^{12}CO(J = 1−0): we use the excitation temperature whenever it is above the 60 K threshold, and in other regions we use the maximum of the gas excitation temperature and the dust temperature (Fig. A.1, panels 1 to 3).
Fig. A.1 First three panels: speed of sound computed for the estimated excitation temperature of the gas, the dust temperature, and the maximum of the two previous ones. Last panel: turbulent flow velocity dispersion. We assume that having one or several spectral components along the line of sight is irrelevant, since one component gives a turbulent velocity dispersion, and two components show a highvelocity shock. 

Open with DEXTER 
When computing the speed of sound from the temperature map, the nature of the gas comes into play. We considered a 75−25% mixture in mass of molecular hydrogen and helium, which yields a molecular mass of 2.513 kg mol^{1}. Since we are not simply dealing with a monoatomic or diatomic ideal gas, the adiabatic index γ of this mixture has to be determined. To compute γ, we used tabulated NIST values of the calorific capacities C_{P} and C_{V} of the two gases at an average temperature of 24 K. The resulting value, γ = 1.66674, is very close to the value 5/3 that would be expected for a monoatomic ideal gas, which could be expected since, at such low temperatures, the rovibrational modes of H_{2} are frozen.
Appendix A.2: Flow velocity
The turbulent velocity dispersion was computed using the velocity dispersion along each line of sight, which can be derived from the W_{2} map. We determined the flow velocity dispersion as (Fig. A.1, last panel).
This computation implies two assumptions. First, we assume that the natural and thermal widths of the lines are negligible, compared to the turbulent broadening. Second, we neglect the opacity broadening as well. While it is correct that thermal broadening has less than 1% effect on the line width, the assumption for the opacity broadening is less obvious. We determined that the expected correction would be of the order of 10% for the brightest lines of sight, reducing the width of the lines and therefore the Mach number. However, given that the correction would have been difficult to implement for nonGaussian lines, as is the case in our ^{13}CO(J = 1−0) field, and that the correction would be significant only for a small fraction of the lines of sight, we left the turbulent velocity field uncorrected.
Appendix A.3: Results
The ratio of the turbulent velocity to the sound speed gives us the Mach number along the z axis. The total Mach number is , when the turbulence is isotropic. The Mach number is determined twice, using two maps: the excitation (gas) temperature, and the maximum of dust and excitation temperatures, as described above. We then compare the results of the two computations (Figs. A.2 and A.3).
To estimate the typical Mach number at the scale of the whole cloud, we plot the histogram of the obtained maps (Fig. A.3). The shapes of the histograms are very similar for T_{max} and T_{exc}. Owing to the shape of both distributions, with a large tail at high Mach number, the most probable value of the Mach number for both distributions (M = 3.5), is significantly smaller than the mean and the median values, as shown in Table A.1.
Fig. A.2 Mach number maps computed using the two different temperature maps, left using T_{exc}, right using T_{max}. 

Open with DEXTER 
Fig. A.3 Histogram of the Mach number map, for T_{max} ant T_{exc}. The vertical bars show the position of the mean value of each distribution. 

Open with DEXTER 
Typical values of the temperature and the Mach number distribution, for the two computed temperature maps.
All Tables
Typical values of the temperature and the Mach number distribution, for the two computed temperature maps.
All Figures
Fig. 1 H_{2} column density map of the southwestern part of the Orion B giant molecular cloud, derived from Herschel Gould Belt Survey observations (André et al. 2010; Schneider et al. 2013). The field observed by the OrionB collaboration and used for this work is overlaid in red. The blue squares mark the nebulae NGC 2024, NGC 2023 and the Horsehead, and the pink star symbols mark the stars Alnitak, HD 38087 and σ Ori (north to south). 

Open with DEXTER  
In the text 
Fig. 2 Maps of the average brightness temperature of the ^{13}CO(J = 1−0) line in four contiguous velocity ranges. The mainbeam temperature scale is indicated by the color bar on the right. The contour shows the value of 8.9 K km s^{1} in the W_{0} map, corresponding to 0.43 K in the mean temperature map integrated over the −0.5−20.5 km s^{1} velocity range. The set of coordinates used for the observational campaign takes the Horsehead PDR as a reference point, and aligns the IC 434 PDR along the vertical axis (14° counterclockwise rotation with respect to equatorial coordinates). The numbered squares in the first panel show the positions of the spectra presented in Fig. 3, from left to right. 

Open with DEXTER  
In the text 
Fig. 3 ^{13}CO(J = 1−0) spectra along selected lines of sights in the Orion B cloud, showing the diversity of velocity components (up to four per spectrum). The coordinates of the various lines of sights are given in arcminutes in our custom set of coordinates (δx, δy). The average spectrum for the whole field, normalized to have the same peak temperature, is superimposed with a dashed line for comparison. The positions of the spectra on the map are shown by white squares in the first panel of Fig. 2. 

Open with DEXTER  
In the text 
Fig. 4 Top: maps of the ^{13}CO(J = 1−0) field moments in the cloud’s frame of reference. Bottom: maps of the physical quantities directly derived from each field, with the centroid velocity being simply W_{1}/W_{0}, and the Full Width at Half Maximum (FWHM) velocity dispersion given by – the normalization of the FWHM corresponds to that of a field with purely Gaussian line profiles. 

Open with DEXTER  
In the text 
Fig. 5 Left: powerlaw fitting of the W_{0} power spectrum. Upper panel: data (blue crosses), fit result (thick solid line) plotted over the fitted domain, power law convolved with the Gaussian beam (dotted line) extrapolated to all spatial frequencies, noise model (dashed line). Lower panel: residuals. The scale in the Fourier space is given in UV distance as for interferometric observations, which enables us to visually relate the resolution and the telescope diameter. Right: same results, except for the W_{1} power spectrum. 

Open with DEXTER  
In the text 
Fig. 6 Variation of the FWHM as a function of the gas density. Each point refers to the (J = 1−0) spectral line of a different molecule. The red line shows the least squares fit, with a slope α = −ϵ = −0.15. 

Open with DEXTER  
In the text 
Fig. 7 Solenoidal turbulence fraction for sliding square areas with a width equal to the one of the full map (56 arcmin). The shaded area corresponds to the g_{21} uncertainty, while the error bars show the experimental uncertainties due to observational noise. The vertical line marks the equipartition limit. The map on the left presents the centres of the square areas for which the calculations were performed, superimposed on the ^{13}CO(J = 1−0) T_{peak} map. 

Open with DEXTER  
In the text 
Fig. 8 Solenoidal turbulence fraction for zooms on the NGC 2023 and NGC 2024 star forming regions. The shaded area corresponds to the g_{21} uncertainty, while the error bars show the experimental uncertainties due to observational noise. The upper limit of the plots marks the equipartition limit. The map on the left presents the areas for which the calculations were performed (solid squares: NGC 2024, dashed squares: NGC 2023), superimposed on the ^{13}CO(J = 1−0) T_{peak} map. 

Open with DEXTER  
In the text 
Fig. 9 Comparison of the calculation results in the case of the zoom on NGC 2024 with three different methods of computing the power spectra: either the power spectra resulting from Fouriertransformed data are used directly, with a linear interpolation between the points (data only, blue), or the result of the power law fit at high frequencies is extrapolated to low frequencies to provide a power law throughout the whole range of frequencies (fit only, green), or a composite version of the power spectrum is built using the raw data at low frequencies, and the fit result at high frequencies (data + fit, red). The error bars account for the observational noise as well as for the g_{21} uncertainty. The horizontal line marks the equipartition limit. 

Open with DEXTER  
In the text 
Fig. A.1 First three panels: speed of sound computed for the estimated excitation temperature of the gas, the dust temperature, and the maximum of the two previous ones. Last panel: turbulent flow velocity dispersion. We assume that having one or several spectral components along the line of sight is irrelevant, since one component gives a turbulent velocity dispersion, and two components show a highvelocity shock. 

Open with DEXTER  
In the text 
Fig. A.2 Mach number maps computed using the two different temperature maps, left using T_{exc}, right using T_{max}. 

Open with DEXTER  
In the text 
Fig. A.3 Histogram of the Mach number map, for T_{max} ant T_{exc}. The vertical bars show the position of the mean value of each distribution. 

Open with DEXTER  
In the text 