Characterizing the structure of diffuse emission in HiGAL maps
Abstract
We present a study of the structure of the Galactic interstellar medium through the variance technique, related to the power spectrum and the fractal properties of infrared/submm maps. Through this method, it is possible to provide quantitative parameters which are useful to characterize different morphological and physical conditions, and to better constrain the theoretical models. In this respect, the Herschel^{1}^{1}1Herschel is an ESA space observatory with science instruments provided by Europeanled Principal Investigator consortia and with important participation from NASA. Infrared Galactic Plane Survey carried out at five photometric bands from 70 to 500 m constitutes an unique database for applying statistical tools to a variety of regions across the Milky Way. In this paper, we derive a robust estimate of the powerlaw portion of the power spectrum of four contiguous HiGAL tiles located in the third Galactic quadrant (, ). The low level of confusion along the line of sight testified by CO observations makes this region an ideal case. We find very different values of the power spectrum slope from tile to tile but also from wavelength to wavelength (), with similarities between fields attributable to components located at the same distance. Thanks to the comparison with models of turbulence, an explanation of the determined slopes in terms of the fractal geometry is also provided, and possible relations with the underlying physics are investigated. In particular, an anticorrelation between ISM fractal dimension and star formation efficiency is found for the two main distance components observed in these fields. A possible link between the fractal properties of the diffuse emission and the resulting clump mass function is discussed.
Subject headings:
Methods: data analysis — Methods: statistical — ISM: clouds — ISM: structure — Infrared: ISM — Stars: formation1. Introduction
One of the most intriguing tasks in the observational study of the interstellar medium (ISM) is to extract information about the 3dimensional structure of the clouds, starting from the 2dimensional maps of these objects, generally taken at different wavelengths and with different techniques and resolutions. Although a certain degree of selfsimilarity of the ISM maps over a given range of spatial scales can be in many cases perceived by eye, there are numerous more solid arguments suggesting this can be the case, starting from the work of sca90.
In this respect the phenomenon mainly responsible of selfsimilar morphologies is turbulence. This is a largely recognized fact in molecular clouds, being a typical scalefree phenomenon inducing fractality (see, e.g., sre89). It is indeed characterized by the lack of a specific length scale, then it can produce a fractal distribution of matter in a molecular cloud over a wide range of scales. Therefore, to determine the starting and the ending point of these ranges is generally considered a tentative way to get an estimate of the turbulence injection and dissipation scales. An extensive and detailed review of the observational evidences of the presence of turbulence in molecular clouds and its role in shaping their structure in fractal sense can be found in vaz99 and sch11.
It is noteworthy that the ISM clouds belong to the category of the stochastic fractals, whose structure does not appear perfectly self similar, but rather selfaffine: although a given set and a part of it have not exactly the same appearance, they have the same statistical properties and it is still possible to use a fractal description for them.
There are many observational grounds supporting the fractal scenario. The observations of the low CO and CO emission lines in several starforming molecular cloud complexes (e.g., fal96; sch98; wil99b) show that the measured line intensities, shapes and ratios cannot be produced in clouds of uniform gas temperature and density, suggesting the idea that these interstellar objects are far from being homogeneous, being instead organized in small clumps with a filling factor lower than unity (elm97a). Interestingly, such a structure is also able to justify further observed characteristics of the investigated region, as for example the clump mass function (sha11) and the stellar initial mass function (elm02). These remain meaningful observables although in the last years the picture of the ISM has changed with the recognition of filaments as intermediate structures (e.g., ros96; wil99b), which have definitely been found ubiquitous in the recent Herschel observations (e.g., mol10b; sch13). In any case, the cloud description based on a hierarchical decomposition in recognizable substructures (hou92) is not incompatible with the fractal approach. Indeed stu98 have shown that these two points of view are consistent: an ensemble of clumps with a given mass and size spectrum can give rise to a fractal structure of the cloud.
Statistical descriptors which, in general, can be related to the fractal properties of a cloud are powerful methods to characterize its structure. The techniques initially used to estimate the fractal dimension of the interstellar clouds were based on the isocontours of the images, as for example the perimeterruler and the areaperimeter relation (see e.g. san05, and references therein). Subsequently, statistical tools have been applied, namely descriptors based on the value and the spatial distribution of the single pixels, providing quantitative information on one or more aspects of the investigated morphology (a relevant part of them is summarized in elm04). The direct estimate of the power spectrum (e.g. ing04; miv07; mar10; gaz10) can be used to infer the fractal structure of the ISM, although to deal with real observational sets other algorithms have been demonstrated to be more adequate (stu98). Other statistical estimators are the structure function (padc02; pad03; kri04; cam05; gus06; kow07; row11), the variance (see below), the autocorrelation function (cam05), and the adapted correlation length (car06), whereas a further development of these monofractal descriptors is represented by the multifractal spectrum (cha01; vav01). In particular, the variance method was introduced by stu98 and subsequently improved by ben01 and oss08a to analyze the drift behavior of observed scalar functions such as the intensity distribution in molecular clouds, real or synthesized. It has been applied not only to maps of line emission (see also ben01; oss01; oss08a; sch11; row11) and dust extinction (cam04; sch11) or emission (rus13), but also to the recovered velocity field (oss06; fed10), or to 3dimensional density fields of turbulence simulations (fed09).
The aim of this paper is to contribute both to the enlargement of the sample of the regions whose structural properties have been studied by means of fractal techniques, and to the improvement in characterizing the response of statistical tools to different observing conditions. The HiGAL survey (Herschel Infrared GALactic plane survey, mol10a) represents an extraordinary resource for carrying out statistical studies of the ISM. Indeed a large coverage is obtained in five different bands, so that a large variety of morphologies and physical conditions can be investigated at unprecedented spatial resolution. Moreover, these large Herschel maps offer the chance to probe a wide range of spatial scales, since the number of available pixels is very important for the reliability of the statistical descriptors.
Galactic plane observations suffer from confusion due to the superposition of different components along the line of sight, especially in the first and fourth Galactic quadrants. To minimize the problem of confusion, the first available observations of the third Galactic quadrant (in the range ) are studied as a first test case, in which we are more confident that the observed ISM emission corresponds to a morphology which is quite coherent from the spatial point of view. These observations have been presented by eli13 (hereafter eli13, ), and are briefly summarized in Section 2.
As a paradigm of synthetic cloud images used for testing the statistical tools used in this work we consider the class of socalled fractional Brownian motion images (hereinafter fBm). They have been already used, for example, by stu98; ben01; kha06; miv07; sha11 for testing their algorithms. We briefly discuss the properties of this class of images in Section 3.
In this paper we adopt the variance algorithm to derive a robust estimate of the power spectrum slope of the maps. In Section 4 this method is briefly described, and its application to synthetic maps is discussed to characterize the response of the algorithm in case of departure of the analyzed image from the ideal fBmlike behavior.
In Section 5 we preent the results of our variance analysis, and discuss the obtained power spectrum slopes and selfsimilarity ranges, searching for crosscorrelations among different maps and observational wavebands. Moreover, links with turbulence and observables related to star formation (as star formation efficiency and mass functions) are investigated. Finally, the results are summarized in Section 6.
2. Observational datasets
The Herschel (pil10) open time key project HiGAL (mol10a) is a fiveband photometric survey initially aimed at studying the stellar life cycle in the inner Galaxy () and subsequently extended to the whole Galactic plane.
The first HiGAL observations available for the outer Galaxy, presented in eli13, consist of four adjacent tiles of the Galactic third quadrant. We will denote these with , , , , respectively, according to the Herschel Data Archive nomenclature. These farinfrared maps of the outer Galaxy represent an ideal case for studying the structure of the ISM, for two main reasons: the lower occurrence, in general, of compact bright sources and of star forming regions, and the lower degree of confusion along the line of sight.
The reduction procedure and the main characteristics of these observations are described in eli13, so here we resume only those features which are useful for the discussion in this paper. The observed wavebands are centered around 70 and 160 m (PACS, pog10) and 250, 350 and 500 m (SPIRE, gri10), with nominal resolutions of about , , , , , respectively. The pixel sizes of these maps are , , , , , at 70, 160, 250, 350 and 500 m, respectively. The fields were observed by simultaneously acquiring PACS and SPIRE images in the five aforementioned photometric bands. This observing mode generally implies that the areas imaged by the two instruments are not exactly the same; in this paper, in particular, we consider only the common area of each tile, because we are interested in comparing the statistics of the same regions of the sky seen at different wavelengths.
The final images extracted at each wavelength for each tile are shown in Figure 1, where the color coding used to identify the different Herschel bands throughout the article is introduced. Since for a given tile the same area of the sky is considered in each band, the total number of pixels depends on the band: the size of each image in pixels is reported in Table 1.
Tile  70 m  160 m  250 m  350 m  500 m ^{a}^{a}The size of the column density maps is the same as that of the corresponding 500 m maps. 

217  1971  1403  1051  789  549 
220  1951  1387  1041  781  545 
222  1961  1395  1047  785  547 
224  1977  1405  1055  791  551 
The column density maps of each tile are analyzed here as well. They have been derived from a pixeltopixel modified black body fit as explained in more detail in eli13. To this end, first the maps were absolutely calibrated correcting their zero level by means of offset values derived from the comparison with the Planck/IRAS data ber10. Then, the 160, 250, and 350 m maps were reprojected onto the grid of the 500 m ones. The 70 m maps were not involved in this calculation because they can contain emission coming from the socalled Very Small Grains (com10), not at thermal equilibrium, and/or from warmer spectral components, such as the protostellar content of the clumps, reflection nebulae, etc. For these reasons the column density maps we obtained are more suitable to describe the cold component of dust in these regions. At a first glance, they look very similar to the SPIRE maps, and in particular to the 500 m ones, having also the same resolution: both these aspects turn out to be important for the spatial analysis reported in the following sections.
The four tiles show the presence of star forming activity, but 217 and 224 contain the brightest and most extended regions, mainly arranged in large filamentary shapes.
The knowledge of the kinematic distances of the clouds is somehow fundamental for identifying really coherent nebular structures, instead of superpositions of different distance components, to which the statistical descriptors can be applied. The distances for the HiGAL tiles analyzed here are obtained from CO(10) line observations carried out at the NANTEN submm telescope (miz04), and presented in detail in eli13. In detail, gas emission in the 217 field is dominated by a component located at a kinematic distance of kpc (see Table 1 of eli13, their component II). CO emission is fainter in the 220 tile, yet the component corresponding to the 2.2 kpc distance is still predominant.
Looking at the second raw of Figure 1, the 222 and 224 tiles are found to be dominated by a bright component corresponding to an average distance of kpc, i.e. the socalled component I of eli13. Two further distance components, III and IV, can be considered negligible in the present analysis. The former, located at an average distance of 3.3 kpc (then likely in the Perseus arm) coincides with a bright region in the southeastern corner of the original 222, which however lies outside the region considered in this paper as a consequence of the tile cropping described above. The latter is constituted only by a bright but small portion of the 217 field.
3. Power spectrum and fBm images
The statistical descriptors we use in this paper revolve around the central concept of power spectrum of the observable , defined as the square modulus of its Fourier transform . The variable is defined in a dimensional space, with in the case of image analysis. A convenient choice is to study the shellaveraged power spectrum , where . Because a power law is a distinctive experimental signature seen in a large variety of complex systems, frequently a search for a powerlaw behavior of the spectrum is carried out. An emblematic case is that of the Kolmogorov’s power spectrum of turbulence (kol41), calculated for the velocity field of incompressible fluids; a dependence as is expected, with in three dimensions, in two and in one, respectively. Different models, as for example the classic Burgers’ turbulence (e.g. bec07), still show a similar power spectrum, although with different slope (, bis03). The presence of such a powerlaw behavior (in the full range of values of the wave number, or over a limited part of it) can be interpreted as an indication of turbulence and can suggest the characteristic scales at which both energy injection and dissipation take place.
The power spectra (or portions of them) of the ISM maps often exhibit a powerlaw behavior. The link with turbulence is quite natural: it is commonly accepted that a turbulent velocity field in the ISM can also shape the density field (cf. bol02; pad03; row11). To derive the slope of the power spectrum, however, requires some care. Indeed, the Fourier transform of a noninfinite mapped signal inevitably introduces unwanted frequencies due to the spatial sampling and to the limited size of the image, since the Fourier transform implicitly assumes wraparound periodicity) leading to aliasing (stu98; ben01). This has lead, in the past, to the use of statistical tools more robust than the direct determination of power spectrum. The improved performance of the new generation of observing instruments have allowed the production of maps with large numbers of pixels (as in the case of the HiGAL tiles), which might mitigate the aforementioned issues. Nevertheless in this paper we prefer to keep exploiting one of these indirect methods, namely the variance technique, briefly discussed in Section 4. In this way we can both make possible a direct comparison with the literature, and exploit further information that this technique can provide about the structure of the maps (see Section 4).
3.1. Fractional Brownian motion images
There is a category of stochastic fractals that can be helpful to test the statistical algorithms for structure analysis, since their Fourier transform has specific analytic properties which make them easy to be generated: the fractional Brownian motion (pei88). In the twodimensional case, these show a good similarity with molecular cloud maps (stu98; ben01; miv03, see also Figure 2). A detailed description of the fBm image properties is given in stu98; here we summarize only the most relevant ones, which will prove to be useful in the analogy with ISM maps we will make in next sections. First, they are characterized by having a powerlaw power spectrum with exponent , where is the Euclidean dimension of the considered space (for cloud maps, ), and is the socalled Hurst exponent, ranging from 0 to 1. Therefore can take values from 2 to 4. Second, the phases of their Fourier transform are random. Based on these two constraints, it is quite easy to obtain fBm images once is assigned and a random phase distribution is generated^{1}^{1}1The fBm images generated in such a way are periodic, i.e. it is possible to place side by side and to connect with continuity the image to itself. Clearly, this can not be the case of real ISM maps.. In Figure 2, six pixel fBm images are shown. They have the same phase distribution and differ only by , which ranges from 2.0 to 4.0 in steps of 0.4. It can be seen that the phase distribution determines the overall appearance of the “cloud”, but as increases the structure becomes smoother and smoother due to the transfer of power from high to low spatial frequencies.We recall that the case would correspond to white noise.
In terms of the fractal description it has been shown that the fractal dimension of an dimensional fBm image is given by:
(1) 
so that the direct relation between and is:
(2) 
An important property of the fBm images is that the power spectrum of the ()dimensional projection on an dimensional fBm set is again a powerlaw with the same spectral index (stu98). This result turns out to be important for establishing a link between the observed 2dimensional column density and the real 3dimensional cloud density field for clouds that can be considered isotropic, as suggested through theoretical arguments (bru10) and empirical evidences (fed09; gaz10).
Thus, invoking Equation 2, it is found that
(3) 
i.e. the fractal dimension of a fBm object changes, under projection, by 1.5, and not by 1 as one could erroneously expect.
4. variance analysis
The variance method is a generalization of the Allan variance (all66), elaborated and characterized in detail by stu98. For a twodimensional observable field , the variance at the scale is defined as the variance of the convolution of with a filter function :
(4) 
where
(5) 
is the downupdown cylinder (or French hat) function and . The two nonzero terms of the above definition represent the core and the annulus component of the filter, respectively.
oss08a recommended, as a possible alternative to the French hat function to obtain a more reliable estimate of the spectral index , to use the smoother Mexican hat, defined as
(6) 
where the two main terms in the right side of the equation represent the core and the annulus components, respectively, and is the diameter ratio between them. To speed up calculations, the same authors suggest to perform the operation in Equation 4 as a multiplication in the Fourier domain,
(7) 
where is the power spectrum of , and is the Fourier transform of the filter function.
The fundamental relation that relates the slopes of variance and of the power spectrum () was shown by stu98:
(8) 
Given the expression above, one can derive the power spectrum slope by performing a linear fit over the range of spatial scales for which the logarithm of variance manifests a linear behavior. In this work, we adopted this procedure following the prescriptions of oss08a^{2}^{2}2The IDL package for calculating the variance can be found at http://hera.ph1.unikoeln.de/ ossk/Myself/deltavariance.html]http://hera.ph1.unikoeln.de/ossk/Myself/deltavariance.html., i.e. using a Mexican Hat filter with . Furthermore, we do not adopt any strategy based on assigning different weights to the pixels involved in the variance calculation, which is recommended by oss08a in case of maps characterized by a variable data reliability. The portions of the HiGAL maps we chose, indeed, being far from the tile boundaries, are characterized by a quite uniform coverage (see also tra11) and, consequently, by a stable rms noise.
4.1. The contribution of compact sources
From Equation 8 the variance behavior of a fBm image is expected to be a perfect power law. Although the ISM maps generally exhibit a fBmlike behavior (see Sect. 3.1), it is important to identify and characterize all the signatures in the power spectrum ascribable to possible departures from the ideal fBm case, first of all the presence of bright compact sources. For this purpose, we performed a test by simulating a PACS m map of a portion of the Galactic plane, using some typical parameters of this band, such as the pixel size of . The steps of the recipe can be also followed in Figure 3, together with the effects they achieve on the corresponding variance curve:

A pixel fBm background has been generated with a “typical” power spectrum slope (see, e.g., sch11). Therefore, to avoid dealing with a periodic image (see Section 3.1), this has been truncated, extracting a subimage of pixels.

To simulate the presence of very bright small regions, another pixel fBm set has been generated with a significantly higher power spectrum slope () and a different phase distribution. The resulting image has then been exponentiated in order to enhance the highsignal regions. Again, a subimage of pixels has been extracted, making sure that the extracted image still contains the maximum of the original image. Finally, this has been added to the image obtained in i), resulting in regions of enhanced brightness.

To reproduce the luminosity decrease off the Galactic plane, a modulation through a Gaussian profile has been applied. To make this profile more realistic, both the FWHM of the Gaussian and position of the peak slowly float as a function of the longitude, following a Gaussian distribution and a longperiod sinusoid, respectively. A relevant decrease of emission moving away from the plane is more pronounced in the HiGAL observations of the inner Galaxy with respect to those considered in this work. However, the goal of this test is to identify qualitatively the effect on the variance curve of peculiar structures , thus it is instrumental to exacerbate these contributions.

A population of 500 compact sources has been spread across the map, generating random 2dimensional Gaussians whose size and peak flux distributions follow those found for the HiGAL field (eli10). The probability of displacing a source in a given position of the map has been weighted with the intensity of the image in that position, to obtain a more realistic concentration of compact sources in regions with bright diffuse emission. Before coadding the sources, the background image has been scaled by a given amount, such that its dynamical range gets similar to that of the background at m. Finally the image has been convolved with the PACS m beam, and a lowlevel white noise is summed over, using an additional fBm image with .
In the right panels of Figure 3, the variance of the simulated images displayed on the left is plotted. All the images have been normalized between 0 and 1, to prevent overflows in calculation. As a consequences, the units of variance are arbitrary. We note that the variance slope is not affected by this rescaling. In step i), the extracted subimage curve (black line) shows a slightly steeper slope than that of the original set, plotted as reference (green line). While the latter looks as a line as expected, the former exhibits a linear behavior only over a limited range of scales (the grey area). The flattening of the curve at is just due to the truncation of the original fBm set. Moving to step ii), a slope similar to the original one is still found on the limited range of scales (darker grey area ), while the steeper slope at is caused by the coaddition of the bright spots in the image. A further steepening is produced in step iii) by modulating emission with a low spacial frequency profile. This is particularly evident at the largest scales, where it compensates the flattening of the variance seen in the previous steps. Finally, the insertion of compact sources in step iv) is responsible for the appearance of a bump for , which is the typical size range of the injected sources. This clearly corresponds, in light of the correspondence between variance and power spectrum, to the component of the power spectrum discussed by miv07 and mar10. Here the diffuse emission behavior can be recovered only for a limited range of scales, since at the largest scales the effect due to the Galactic plane shape is obviously still present. At the smallest scales, instead, a flattening of the variance curve is seen, due to white noise (cf. ben01). In any case, no physical information can be extracted at scales smaller than the instrumental beam.
Notice, however, that the significance of the effects described in this section depends on the analyzed map. For example, for the third Galactic quadrant maps we don’t expect a strong influence of the Galactic plane shape. Furthermore, the bump in variance due to compact sources is likely not as sharp as seen in Figure 3, because in the real maps the transition between a compact source and the surrounding cirrus emission is smoother than in our simulations. Larger clumps, Hii regions and filaments present in the real maps are additional intermediate structures between the two scale regimes of compact sources and the diffuse emission, contributing to enlarge the bump toward larger spatial scales, and to connect it more smoothly with the linear portion of the variance. The importance of the filaments in the scenario of star formation and their ubiquity in the ISM have been highlighted by Herschel observations (e.g., mol10b; arz11; rus13; sch13), and they are certainly the main responsible for the departure of the images from a fbmlike behavior at intermediate scales between compact sources and cirrus.
5. variance of the third Galactic quadrant HiGAL fields
After normalizing, as described in Section 4.1, the HiGAL maps shown in Figure 1, the variance vs spatial lag curves have been computed. They are plotted in Figure 4 using the colorband coding introduced in Figure 1.
5.1. General results
Some general considerations can be drawn from the global trends exhibited by these curves. The lack of a sufficient level of diffuse emission in the m maps influences the corresponding variance spectrum, which shows peculiar trends compared with other wavelengths. Therefore, the curves at this wavelength are plotted only for completeness and the corresponding slopes are not shown. Only as a general remark, we notice that in the tiles 217 and 224, namely those with a significant emission of compact emission at m, a bump peaking at about 40″and 25″, respectively, is present. At larger spatial scales, the m curves show in all cases slopes smaller those of the larger wavelengths, even negative in three cases out of four. This behavior is expected in presence of low signaltonoise ratio in the maps, keeping in mind Equation 8 and the white noise (=0) borderline case for the power spectrum slope.
For the remaining wavebands and for the column density maps, only one common linear range has been identified for each tile. Noteworthy, scales shorter than 100 have been not considered because of the possible contamination by compact sources. The curvature of each line, defined as has been estimated, by allowing only ranges where reasonably low values of are found^{3}^{3}3In reality, only few variance curves exhibiting no linear behavior in the common interial range, e.g. the one of m of 224, have not been considered in the evaluation of the linear range. However, despite these exceptions, the constraint we apply on the linear range appears more robust than the qualitative criteria adopted in the literature.. This procedure keeps the extremes of the range far from possible peaks of the curve (e.g. in 222), where the most relevant departure from linearity is expected.
The estimated fitting ranges and the corresponding slopes of the power spectrum (obtained from the variance ones through Equation 8) are reported in Figure 4, and summarized in Table 2, together with the corresponding fractal dimensions of the maps (derived through Equation 2).
Field  Distance  Fit range  

[kpc]  [pc]  160 m  250 m  350 m  500 m  Col. dens.  160 m  250 m  350 m  500 m  Col. dens.  
217  2.2  1.3–4.2  1.02  1.45  1.65  1.90  2.15  3.49  3.27  3.17  3.05  2.93 
220  2.2  1.3–5.3  2.12  2.21  2.25  2.35  2.35  2.94  2.89  2.87  2.83  2.83 
222  1.1  0.5–3.8  2.65  2.72  2.74  2.79  2.65  2.67  2.64  2.63  2.61  2.68 
224  1.1  0.5–1.7  2.17  2.27  2.41  2.61  2.77  2.91  2.86  2.79  2.70  2.61 
As a general remark on some evident trends found in all the four tiles, we notice that, from the qualitative point of view, for each tile the three SPIRE bands have quite similar spectral behaviors, whereas the general shape of the PACS 160 m curves appear relatively different. More importantly, a systematic increase of the slope with the wavelength is found from 160 to 500 m.This suggests that not only the emission morphology changes when observed at different wavelengths, but also its statistical properties are found to be different. In particular, at 160 m the contribution of warmer Very Small Grains can be still relevant (com10) and seems to be responsible for a more uniform distribution of the power of the image through the different spatial scales, resulting in a shallower . In other words, the warmer dust turns out to be more diffuse and spread around than the cold dust, which is expected to be preferentially concentrated in denser environments like filaments (e.g., pad06; cam07), producing a global smoothing of the observed emission features.
This effect can represent a possible explanation for the systematic discrepancy, found by sch11, between the variance slope of and CO maps (shallower and steeper, respectively) of the same areas of the sky. Most likely, these two tracers do not describe exactly the same components of the ISM in the same manner (e.g. goo09). Moreover, the ISM is optically thinner at the SPIRE wavelengths than at 160 m, so at long wavelengths one expects a more enhanced contrast between emission from highdensity and lowdensity regions, resulting in a possible steepening of the power spectrum (and, equivalently, of the variance). However this reasoning could be too simplistic, because moving towards SPIRE wavelengths a number of cold small scale filaments can manifest themselves thus contributing to the power spectrum and counterbalancing the effect described above.
Another consideration concerns the variance of the column density maps: although the maps look quite similar to the SPIRE ones, the power spectrum behavior is generally found to be slightly different from those. Anyway, all the slopes we found lie in the typical range of values found for clouds studied in previous works (ben01; sch11; row11; rus13).
5.2. Results of individual maps
Going into detail of single tiles, we start from the westernmost field, 217, associated with distance component I (see Figure 1). As in the case of 224 discussed in the following, the abundance of compact sources and filamentd in this field produces a visible bump at (see Section 4.1) at all wavelengths, whose peak and upper endpoint shift towards larger scales with increasing wavelength. It is generally followed by a descending trend which stops around 1000″, corresponding to a physical scale of pc at a distance of 2.2 kpc. The presence of the large filament associable with Sh 2287 (sha59) in the northern part of the map is probably the responsible of this slope change, and of the departure of the maps from selfsimilarity at the smaller scales.
The 220 tile exhibits a more extended range of linearity. In fact this tile is, in the data set we consider, the poorest of bright features, so that the cirrus component can be really probed. The inertial range (1.35.3 pc) partially overlap those found by sch11 for some lowmass star forming clouds. A good correspondence is also found with the Herschelbased analysis of the NGC 6334 star forming region () of rus13, where three out of four separate subregions exhibit inertial ranges similar that of 220. At longer scales, around 15 pc, the variance curves flatten, which does not seem to correspond to any visible structure in the maps, such as filaments, ridges or bubbles (and corresponding cavities). In this case it is likely that the upper limit of the selfsimilarity range is really correlated with the injection scale of turbulence.
In tile 222 the steepest slopes are found, over an inertial range of 0.52.7 pc. As in the case of 220, the small number of bright features in the maps translates into a wide range of linearity of the variance and into a weak compact source bump. A peak is present at ( pc at kpc), corresponding approximately to one sixth of the map size. Although we do not have sufficient information to claim that this feature is associated with the bubble located on the East side of the tile (see Figure 1), we believe that this structure is responsible for the high values of we find: indeed it generates a certain degree of segregation between relatively empty and bright regions (cf. for example panels and of Figure 2), hence a transfer of power towards large scales.
Finally for 224, similarly to 217, we find that the contribution of compact sources and filaments introduces features in the variance curves which make it difficult to identify a possible inertial range. The one we find between 0.5 and 1.7 pc (neglecting the 160 m curve) is compatible with that of 222. The slopes are generally steeper than those of the component II tiles, but shallower than those of 222.
This comparison suggests that the region of the plane covered by the eastern tiles (222 and 224), which is quite coherent from the kinematical point of view being associated with distance component I, show some common global statistical properties, different from those of the western tiles (217 and 220, component II), generally characterized by shallower slopes. Furthermore, within the two distance components, the tiles containing bright features (217 and 224, respectively) show slopes shallower than those of the corresponding lowemission tiles (220 and 222, respectively).
Anyway, the variety of powerspectrum slopes we find in different tiles reinforces the scenario of the nonuniversality of the ISM fractal properties. From the morphological point of view, this corresponds to a different distribution of the power of the image on a spatial frequency range. From the point of view of the underlying physics, different power spectrum slopes are related to different conditions of the compressible turbulence which is widely considered the most realistic situation in the ISM, especially in presence of star formation (hen84; padn02). Unlike the “rigid” results of the kol41 incompressible turbulence (see Section 3), the compressible one is able to produce, in the ISM, a variety of morphologies and hence of power spectrum profiles (fed09).
5.3. The fractal dimension of the images
The fractal dimension is another important observable, useful to characterize the ISM morphology. As mentioned in Section 1, several computational approaches have been adopted to derive it. Here we use Equation 2, assuming that the analyzed maps have a fBmlike behavior in the recognized inertial ranges. This makes the descriptions based on and completely equivalent, however to speak in terms of fractal dimension here allows us to make further comparisons with observational and theoretical results present in literature. A check on possible alternative methods for calculating the fractal dimension is provided in Appendix A.
The linear relation between the power spectrum slope and the fractal dimension contained in Equation 2 clearly expresses the intuitive concept that a “smoother” texture (high , see Figure 2) must correspond to a lower degree of fractality (i.e. low ). The factor appearing in the equation can half the perception of the variation of , which in fact is expected to vary only between 2 and 3. Although in the literature a variation of between two values of is presented as negligible, it actually corresponds to significant structural differences of the maps. In this respect, the fractal dimensions reported in Table 2 reassert again i) the dicrease of at increasing observed wavelength (already discussed as increase of in Section 5.1), and ii) the significant differences between the structure observed in the western and in the eastern tiles.
In the three tiles showing fractal behavior (220, 222 and 224), all the values we derived range from (500 m of 222 column density of 224) to (160 m of 220), and most of them are compatible with the typical range of variability found through statistical techniques (row11; sch11, and references therein), but also in most part with the found by san07; san09 through the perimeterarea relation.
Instead, the values we obtained are significantly larger than the average value found on IRAS 100 m maps by miv07 (, which is in the fBm approximation), in particular those found at 160 m, which is the closest band to the IRAS 100 m one.
In any case, the values we find are far away from the first estimates of the fractal dimension of interstellar clouds (, fal91; elm96; elm97b), which initially described a quite constant and universal behavior of the ISM structure, and were found also to be in good agreement with the predictions of the Kolmogorov incompressible turbulence ( for the density field, then , fal94). Instead, the values found here, as well as the others obtained through the variance (sch11, and references therein), are more compatible with the power spectra slopes of compressible turbulence, (e.g. fed09), and in particular for the case of solenoidal driving of turbulence (, ).
We can also compare with the theoretical predictions of gaz10 on the power spectrum of both 3dimensional and column density in turbulent thermally bistable flows^{4}^{4}4A gas in which temperatures above and below a given instability range can coexist in thermal pressure equilibrium is called bistable (see gaz05, and references therein). The atomic ISM is generally believed to be bistable.. Their slope (calculated between 7 and 25 pc, i.e. a range of scales larger than ours) gets shallower as the Mach number increases: from at to at (a similar trend can be clearly found also in kri06; kow07). Although these simulations essentially reproduce the Hi emission and despite the different investigated spatial ranges, the behavior we find suggests that the two eastern fields are characterized by smaller Mach numbers. In turn, the star formation efficiency of a cloud is recognized to inversely depend on the Mach number, whose contribution, however, can be found to be concomitant and then surpassed by that of other physical factors (e.g., ros10).
To investigate this aspect in our data, we derived the 3D turbulent Mach numbers from the NANTEN data cube for the velocity components I and II, as the ratio between the local values of the deconvolved CO(10) line width and the sound speed :
(9) 
Here km s, where is taken from the dust temperature maps obtained in eli13 simultaneously with the column density maps. The amount of CO spectra suitable for this analysis is limited by few constraints. First, only lines with a reasonable SNR are considered (see details in eli13, ); second, the line full width at half maximum must be larger than the spectral resolution element ( km s) to allow deconvolution; third, the location CO spectra must lie within one of the four fields we investigate in this paper. These are shown in Figure 5 (left panel), overlaid to the distribution of the Mach number across the entire NANTEN data set.
We note that we can derive meaningful statistical information only for the brightest regions, while the observations are not sensitive enough to trace the cirrus component, which is the one showing the selfsimilar behavior that justifies such a test. For this reason, we stress that observations with better sensitivity and spectral resolution are required for confirming this result.
We recall that in eli13 the distance components I and II have been found to differ in star formation efficiency by a factor 2 (0.008 against 0.004, respectively). Interestingly, a direct correspondence between a larger fractal dimension and a higher star formation efficiency is found in this case. Unfortunately, the inconclusive argument of the Mach number can not support the aforementioned theoretical picture, so that further checks are required and will be feasible when star formation efficiencies will be available for other portions of HiGAL.
From this discussion it clearly emerges that the estimation of the fractal dimension of farinfrared dust emission from interstellar clouds (or portions of them) constitutes a further constraint to be put on theoretical models of cloud structure. In this sense, the vast amount of data provided by the multiwavelength HiGAL survey represents a real minefield which can be searched for such observables.
5.4. Cloud structure and clump mass function
The quantitative indications obtained in the previous sections allow us to probe a relation derived by stu98, linking the power spectrum of the ISM to the mass function of the clumps (CMF hereinafter) which originate from such a structure. Indeed, as mentioned in Section 1, these authors showed that the fractal description of the ISM and the approach based on a hierarchical decomposition (e.g., hou92) are compatible. stu98 demonstrated that, for a fBmlike cloud characterized by a powerlaw power spectrum with slope , and assuming a clump mass spectrum in the cloud (in its highmass end) and a relation between the clump mass and size as , then
(10) 
This theoretical relation contains the intuitive concept that, once is fixed (where the cases and indicate clumps with constant volume and column density, respectively), a steeper power spectrum corresponds to a power transfer towards larger scales, hence to the presence of bigger clumps, which makes the CMF shallower.
So far Equation 10 has been tested on very few data sets, due to the lack of observational constraints, namely the simultaneous knowledge of the power spectrum and of the CMF of a cloud. stu98 combined the “average” slopes and found for a sample of molecular clouds ( and , respectively), deriving . More recently, sha11 explored this relation analyzing 3dimensional fBm synthetic clouds, starting from the assumption that molecular clouds have on average (different, however, from the shallower slopes we find in and ).
Here we have the possibility to check this relation using the estimates of for the investigated fields, and exploiting the clump masses derived in eli13 to build a statistically significant CMF by selecting the sources as follows. First, it makes sense to consider only prestellar cores, namely starless sources gravitationally bound (i.e. those exceeding their corresponding BonnorEbert mass, see e.g. gia12). Second, the considered source samples should be coherent from the spatial point of view, namely all the sources should belong to a physically connected region, following the steps described in eli13 for 398 compact sources associated with the closest velocity component (I), which is the component dominating the and fields. Here we consider also component II for comparison, identifying 131 sources suitable for this analysis. The amounts of sources associated with components III and IV, instead, are not statistically relevant to extend the analysis. The slope of the CMF for component I has already been estimated in the eli13 as , using an algorithm independent of the histogram binning (olm13, and references therein) which also allows us to determine the lower limit of the validity range of the power law. Similarly, for component II here we obtain (Figure 6, panel ).
For the power spectrum slopes, we decided to use the values obtained for the column density maps of the two tiles characterized by low contamination of compact sources, namely 220 for component I and 222 for component II, respectively. Thus we assume and for components I and II, with error bars of 0.1, and we get and , respectively.
A comparison between these results and the observational mass vs radius relation can be performed using the same sources identified for building the CMFs of components I and II (Figure 6, ). However, the distributions are affected by a significant dispersion and lack of a clear powerlaw trend, as noted in previous analysis of this kind (e.g., and10; gia12), so no clearcut conclusion can be reached. Indeed, the lack of a definite scaling relation is indicative of a departure from pure selfsimilarity. Moreover, despite such clumps have been selected in order to build samples as homogeneous as possible, residual differences of physical and evolutionary conditions can be encountered. This might increase the scattering in the mass vs size plot, making it rather difficult to confirm the theoretically expected value for .
It is interesting, instead, to discuss the two very different values of we find, especially in the context of the applicability of Equation 10 to various cases. The large discrepancy between these values depends on the different average found for the two components, but mostly on the different values of appearing at the denominator of Equation 10. The discrepancy between the CMF slopes of the two components can be justified, in turn, with the relatively different intrinsic nature of the observed sources. eli13 discussed the fact that the detected compact sources, characterized by the same range of angular sizes (from one to few instrumental point spread functions) correspond however to a variety of physical sizes, depending on the source distance. In this case the typical factor between distances of component I and II sources is responsible of the clear segregation in size between the two populations, visible in Figure 6, panel , such that the former are clumps, while the latter are a mixture of clumps and cores, according to the typical definition of these categories (e.g., ber07). It is widely acknowledged that the mass function slope of the cores is similar to that of the stellar initial mass function (e.g., and10) and steeper than that of larger clumps (e.g., , kra98), as a consequence of the underlying fragmentation process. Samples characterized by contamination generally present intermediate slopes (e.g., gia12). In this scenario, the CMF slope of the component II looks more similar to the case of the Polaris flare CO clumps (, hei98) used by stu98 to derive a value of . The component I source sample, on the contrary, consists of sources intrinsically smaller (as better resolved), closer to spatial scales where gravity dominates the morphology of the ISM, definitely stopping its scalefree behavior. As a consequence, in this case the applicability of Equation 10 turns out to be dubious.
These indications must be confirmed in future by extending this kind of analysis to a larger variety of cases, offered by further HiGAL fields characterized by favorable observing conditions, typically found in the outer Galaxy.
6. Summary
We have analyzed the first four fields of the outer Galaxy (, ) observed by Herschel as part of the HiGAL survey, to give a quantitative description of the structure of the far infrared diffuse emission. We exploited the power spectrum slope as descriptor, and the variance algorithm as a tool to derive it in a robust way. The low degree of confusion along the line of sight, revealed by CO line observations, and both the resolution and the large size of the Herschel maps represent an ideal case to study the structure of the diffuse ISM on the Galactic plane. In this respect, the results of our analysis have a double value, because on one hand we characterize the morphology of the ISM in this portion of the plane, on the other hand we provide a set of general prescriptions, considerations and caveats for a possible extension of this kind of analysis to other Herschel maps. It is difficult in some cases to separate global from local results. Accordingly, in the following we try to summarize the conclusions going from the more general ones to those which have a more specific value:

The presence of compact sources in the maps affects the variance curves in terms of a bump up to scales of . Other relevant effects are produced by the interruption of the selfsimilarity due to the Galactic plane latitude emission profile, and, on smaller spatial scales, to the presence of bright filaments and star forming regions. In our sample, however, these effects (especially the former) are less relevant than, for example, in the HiGAL fields of the inner Galaxy.

The maps obtained in four PACS/SPIRE wavebands (160500 m) for the investigated fields, together with the column density maps derived from them, show common features of the variance (or, equivalently, of the power spectrum), i.e. spatial scales corresponding to peaks and turn over points. However, the slope of the power spectrum can remarkably change from one band to another, generally increasing from 160 m to 500 m, probably due to a different spatial displacement of small against large grain dust components.

The obtained power spectrum slopes strongly vary from one tile to another, but remain in the typical range of slopes calculated for different phases of the ISM by sch11. The power spectrum, and the corresponding fractal dimension, are then far from being constant and universal as initially suggested by the analysis of the boundaries of the interstellar clouds (fal91). In two of the four investigated fields, the presence of regions of intense emission is responsible of a peak in the power spectrum, followed by a linear trend with slope shallower than the case of the other two tiles, dominated by the cirrus.

None of these slopes is consistent with the Kolmogorov’s incompressible turbulence case. We find that the model of supersonic isothermal turbulence with solenoidal forcing of fed09 seems more realistic, as already found by sch11 using different tracers.

The power spectrum slopes of the two eastern fields, and are steeper than those of the western fields, and . Each of these two pairs is dominated by the emission of a different distance component ( kpc and 2.2 kpc, respectively). Morphological differences of the density field suggest a higher degree of turbulence in the second component, also suggested by a lower star formation efficiency (yet not confirmed by the Mach number analysis). This represents an interesting evidence of the connection between fractality and star formation efficiency predicted by the theory.

We tested the linear relation, predicted for fBmlike clouds, between the mass function of the over dense structures (clumps) and the cloud power spectrum slope, using as a probe the exponent of the clump mass vs radius power law relation. However the mass vs radius distribution we obtain is affected by significant dispersion and can not be used to draw strong conclusions. Also, this relation seems to work better on the clumps at 2.2 kpc, while its predictions at 1.1 kpc are unrealistic, probably because the sample of clumps associated to this component is contaminated by cores.
Appendix A Box counting fractal dimension
In this paper we followed the fBm approach to derive the fractal dimension of the analyzed images from the analysis of their power spectrum, through Equation 2. Several other approaches might be used, however, to derive for grayscale images. These can be seen as surfaces defined on a 2dimensional plane and embedded in a 3dimensional space (”mountain surfaces”). A fractal image, in particular, is a complex surface whose dimension is larger than 2 (as the Euclidean geometry would expect).
Methods aimed at estimating the fractal dimensions of the image intensity contour levels, as the perimeterruler and the areaperimeter relation, are based on the distorted assumption that the image fractal dimension can be derived from the contour dimension simply by adding 1. This appears immediately wrong in the case of highly anisotropic fractals, but also in more isotropic cases, as the fBm sets (see Equation 3) and simulated IS clouds (san05).
Instead, considering the fractal set itself, it is possible to use approaches based on directly sampling the set with grids having different spacings. For example, the box counting approach is one of the most intuitive ways to derive the fractal dimension, and to understand its sense as well. Following man83, a 3dimensional cube containing the fractal surface is partitioned into grids of cubes of variable size . The number of volume elements occupied by the fractal set is a function of the element size: , where designates the boxcounting fractal dimension, which can be obtained from the least square linear fit of the logarithms. However, whereas the estimate of is rather trivial for binary fractal sets, it requires additional care in the case of grayscale images as ours. One possible approach consists in partitioning the plane in a grid of squares (boxes) of size , and in estimating the difference between the maximum and the minimum values (normalized to integer multiples of ) achieved by the image in each box (sar94). This is needed to estimate the volume spanned by the image in such box. At each explored value of , the average of the volume over all boxes is computed. Again, the socalled mass fractal dimension can be derived from the powerlaw scaling of :
(A1) 
(man83).
For each tile and at each wavelength between 160 and 500 m, we evaluated the linear fit in the same spatial scale range reported in Table 2. In Figure 7, left panel, an example of the calculation is shown for the 222 tile at four different wavelengths.
The comparison of the boxcounting vs. the variancebased fractal dimension of all the HiGAL maps analyzed in this article is shown in Figure 7, right panel. In most cases, a good agreement is found between the two methods at the SPIRE wavelengths, whereas strong discrepancies are found for PACS 160 m, in at least three tiles out of four. It is possible that, being the 160 m band the one exhibiting the most evident departures from a linear (i.e. fractal) behavior (see Figure 4), to estimate an univocal fractal dimension for maps at this wavelength is particularly difficult and plagued by strong uncertainties. In fact, it is well known that the variance is more sensitive to variations of the fractal behavior over different spatial scale ranges (stu98; oss08a; sch11), highlighted by distortions or local maxima/minima of the curves, so that a fractal dimension can be derived in a limited inertial range which is integral part of the information. Instead, the curves obtained through box counting algorithms show linearity over larger scale ranges, as in Figure 7, left panel, but the corresponding fractal dimension represents a coarser description of the global statistical properties of the images.