© E. Jaupart and G. Chabrier 2022

1 Introduction

Observations of molecular clouds (MCs) show that their main properties (velocity and column-density) exhibit large fluctuations. These fluctuations are at the heart of the star formation process (Padoan & Nordlund 2002; Mac Low & Klessen 2004; Hennebelle & Chabrier 2008; Hopkins 2012), implying that knowledge of their statistical characteristics is of prime importance. The accurate determination of the statistics of any quantity must rely on either a large enough number of samples or a large enough sample, so that a natural question arises: can we derive accurate statistical properties of MCs from observations, and if so how can we evaluate the level of accuracy? The relevance of a general statistical analysis of the global properties of MCs (e.g., mass, density PDF, temperature and velocity dispersion) deduced from observations and numerical simulations for studies of star formation processes must be assessed properly. For example, all of the theories that aim at determining the mass spectrum, that is the initial mass function (IMF) or the star formation rate (SFR) in a molecular cloud, rely on the assumption that a restricted number of observations or numerical simulations are representative of any MCs with similar properties (see e.g., Hennebelle & Chabrier 2008; Hopkins 2012; Vázquez-Semadeni et al. 2019 and reference therein). This key assumption must be tested.

Indeed, in studies of star formation based on observations or numerical simulations, one only has access to a small number of samples (and in reality only one most of the time). Therefore, in order to evaluate the statistics of the various stochastic fields of interest, one makes the basic assumption, sometimes called the “fair-sample hypothesis”, that the available sample is large enough for volumetric (or time) averages to be meaningful (see e.g., Peebles 1973 for a discussion in the context of cosmology). This assumption is only valid for stochastic fields that are statistically homogeneous and ergodic (Papoulis & Pillai 1965). Here, one should note that statistical homogeneity must not be confused with spatial homogeneity (we come back to this point below). The assumption of statistical homogeneity has been adopted by many authors, for example in studies of turbulent flows with or without self gravity (Chandrasekhar 1951a,b; Batchelor 1953; Pope 1985; Frisch 1995; Pan et al. 2018, 2019a,b; Jaupart & Chabrier 2020, 2021) and in cosmology for studies of the dynamics of structures in the Universe (Peebles 1973; Heinesen 2020). This assumption, however, provides no information on the magnitude of fluctuations around the average.

Quantifying (i) whether the “fair-sample hypothesis” is correct and, if so, (ii) what are the statistical error bars that derive from it, is an important issue when addressing PDF determinations in star-forming clouds. This is related to the completeness of the observations which we reformulate in this study in terms of statistical accuracy. In the context of the PDF of column-densities, (Alves et al. 2017), for instance, made an attempt to illustrate this problem by introducing the concept of open and closed density contours. These authors suggest that complete observations, that is considered as statistically significant and worth studying, correspond to closed contours. Although interesting, such an approach, however, can not be considered as a robust statistical determination of the bias and the statistical errors corresponding to incomplete observations of a cloud PDF. It is one of the very aims of this article to provide such a robust statistical analysis, using standard tools of random field theory and signal processing. Notably, one of the goals of the paper is to identify the statistical properties of the density field in a cloud inferred from column density data, and to derive procedures, based on the aforementioned tools, to accurately estimate these properties, whatever the PDF (lognormal or not). As such our approach does not make any assumption about the initial driving mechanism of the random motions responsible for the PDF of a cloud: turbulence or gravity.

The proper and standard way of addressing the previous issue relies on Ergodic theory. Ergodic Theory allows one to circumvent the problem of dealing with a single sample and to derive a robust measure of the accuracy of field statistics derived from the available data. In the present context, it also enables us to assess and quantify the relevance of a statistical approach on the evolution of star forming MCs. The key quantity is the correlation length, which is defined in terms of the integral of the auto-covariance function (see e.g., Papoulis & Pillai 1965). The fundamental result is that ergodic estimates are accurate if the dimensions of the sample, i.e. a whole cloud or part of it, are large enough compared to the correlation length. A proper determination of the correlation length in MCs is therefore of prime importance.

In pioneer works, (Scalo 1984); (Kleiner & Dickman 1985) studied the correlations of centroid velocities of the ρ-Oph and Taurus complex respectively and only found evidence of weak correlations at short scale on the order of their resolution. (Kleiner & Dickman 1984) studied the correlation of the column density field of the Taurus complex in search of a statistically significant length scale characterizing the separation of condensations within the complex but did not perform an evaluation of the correlation length, as defined above.

The objectives of this article are twofold. First, our main objective aims at examining the relevance and validity of a statistical approach based on ergodic theory to study the stochastic fields of star-forming MCs. Second, we seek to identify which statistical properties of the density field can be inferred from column density data and we derive procedures to obtain accurate estimates of these properties. The article is organized as follows. In Sect. 2, we outline the mathematical framework that yields the definition of the auto-covariance function (ACF) and correlation length of any statistical sample. In Sect. 3 we derive ways to determine the correlation length of any stochastic field without having to compute the ACF. In Sect. 4, we examine the case of astrophysical stochastic fields induced by compressible turbulent motions. In Sect. 5, we focus on the case of star-forming clouds and on the ways to infer the statistics of these fields from observations of column-densities. In Sect. 6, we apply our calculations to the typical star-forming cloud Polaris. We identify artifacts that are generated when one uses the statistical properties of the column-density field to infer those of the real density field; we show how to reduce these biases. In particular we derive a procedure to obtain proper error bars for the column density PDF. In Sect. 7, we examine the case of Orion B. Section 8 is devoted to the conclusion.

2 Methods: Mathematical Framework for a Statistical Approach

As mentioned in the introduction, a statistical approach of the properties of a cloud (or part of) is valid if this latter is large enough, compared to the correlation length of the quantity of interest, for the measured statistical quantities to be representative with high confidence of the genuine quantities. How to measure this confidence level and thus the relevance of a statistical approach is given by the ergodic theory, as commonly used in statistical physics or in the study of dynamical systems. Indeed, ergodicity implies by definition that different observations and realizations of a given statistical quantity yield results comparable enough for each of them to be representative of the average real quantity. It is described in the next section.

2.1 Ergodic Theory

We rederive here some ergodic theorems that lead to the definition of the correlation length. Let us consider a (scalar) stochastic field X(y), which depends on a D-dimensional position vector y (D = 1, 2 or 3). For a specific and fixed y, X(y) is a random variable of which we want to accurately determine the statistics.

2.1.1 Frequency Interpretation and Repeated Trials

The usual way to estimate the statistical average or expectation $\mathbb{E}\left(X\left(y\right)\right)$ of the random variable X(y) is to observe N samples X(y, ω_i), 1 ≤ i ≤ N, of X(y) and to build the unbiased estimator $${\widehat{X}}_y{,}_N=\frac{1}{N}{\displaystyle \sum _{i=1}^{N}X\left(y,{\omega }_i\right),}$$(1)

with variance $$\mathrm{var}\left({\widehat{X}}_{y,N}\right)=\sigma {\left({\widehat{X}}_{x,N}\right)}^2=\frac{\mathrm{var}\left(X\left(y\right)\right)}{N}=\frac{\sigma {\left(X\left(y\right)\right)}^2}{N},$$(2)

where σ is the standard deviation (std). From Bienaymé-Tchebychev inequality (Papoulis & Pillai 1965), we know that, for any real number m > 0, $$\mathbb{P}\left(\left|{\widehat{X}}_{y,N}-\mathbb{E}\left(X\left(y\right)\right)\right|\le m\sigma \left({\widehat{X}}_{y,N}\right)\right)\ge 1-\frac{1}{m^2},$$(3)

where ℙ denotes the probability of a given event. We note that this inequality is valid for any random field, whether it is Gaussian or not. Although Tchebychev inequality gives a lower limit for the probability, it allows to give a confidence interval to measure the accuracy of the estimator given by Eq. (1). The larger the number of samples N, the smaller the std $\sigma \left({\widehat{X}}_{y,N}\right)$ and the more accurate the estimate in Eq. (1).

In case of statistical homogeneity for field X, the expectation $\mathbb{E}\left(X\left(y\right)\right)$ and std σ (X(y)) are no longer functions of the positions and one can drop the reference to y in Eqs. (2) and (3).

2.1.2 Ergodic Theorems, Autocovariance Function, Correlation Length

In the context of the study of MCs, one usually has only a single sample of X(y). As mentioned in the introduction, if one wants to be able to describe the stochastic fields at play, one assumes statistical homogeneity and build the unbiased estimator $${\widehat{\text{x}}}_L=\frac{1}{L^D}{\displaystyle {\int }_{\Omega }X\left(y\right)\text{d}y\text{,}}$$(4)

where $\Omega ={\left[-\frac{L}{2},\frac{L}{2}\right]}^D$ is a control volume of linear size L and volume L^D, which we want to be as large as possible ¹. The ergodic estimator ${\widehat{X}}_L$ has a variance $$\mathrm{var}\left({\widehat{X}}_L\right)=\frac{1}{{\left(L\right)}^D}{\displaystyle {\int }_{{\left[-L,L\right]}^D}{C}_X\left(y\right){\displaystyle \coprod _{k=1}^D\left(1-\frac{\left|yk\right|}{L}\right)\text{d}y,}}$$(5)

where $C_X\left(y\right)=\mathbb{E}\left(X\left(y^{\prime }+y\right)X\left(y^{\prime }\right)\right)-\mathbb{E}{\left(X\right)}^2$ is the autocovariance function (ACF) of X at a lag y. The stochastic field X is said to be mean ergodic if the estimator ${\widehat{X}}_L$ converges toward $\mathbb{E}\left(X\right)$ as L → ∞ either in the mean square (MS) sense, meaning: $$\mathbb{E}\left({\left|{\widehat{X}}_L-\mathbb{E}\left(X\right)\right|}^2\right)=\mathrm{var}\left({\widehat{X}}_L\right)\underset{L\to \infty }{\to }0,$$(6)

or in probability meaning that, for every ϵ > 0, $$\mathbb{P}\left(\left|{\widehat{X}}_L-\mathbb{E}\left(X\right))\right|>\in \right)\underset{L\to \infty }{\to }0.$$(7)

Bienaymé-Tchebychev inequality (Eq. (3)) not only shows that if X is MS mean ergodic ${\widehat{X}}_L$ also converges in probability, but it also provides a confidence interval for the estimate ${\widehat{X}}_L$. Slutsky’s theorem allows to write an equivalence for the ergodicity of X in a more convenient form: indeed, X is MS mean ergodic if and only if $$\frac{1}{{\left(L\right)}^D}{\displaystyle {\int }_{{\left[-L,L\right]}^D}{C}_X\left(y\right)\text{dy}\underset{L\to \infty }{\to }0.}$$(8)

From there we obtain two sufficient (physical) conditions for X to be mean ergodic. Either: $${\displaystyle {\int }_{{\text{R}}^D}{C}_X\left(y\right)\text{d}y<\infty ,}$$(9)

or $$c_x\left(y\right)\underset{\left|y\right|\to \infty }{\to }0$$(10)

We assume both, and use the common definition of the correlation length l_c(X) of the field X as a function of the ACF (Papoulis & Pillai 1965): $${\left(l_c\left(X\right)\right)}^D=\frac{1}{2^D{C}_X\left(0\right)}{\displaystyle {\int }_{{ℝ}^D}{C}_X\left(y\right)dy}.$$(11)

This definition generalizes the usual definitions for 1D fields: $$l_c\left(X\right)=\frac{1}{Cx\left(0\right)}{\displaystyle {\int }_{\left[0,+\infty \right]}{C}_x\left(y\right)}\text{d}y=\frac{1}{2C_X\left(0\right)}{\displaystyle {\int }_{ℝ}{c}_X\left(y\right)}\text{dy.}$$(12)

For l_c(X) «L we then have from Eq. (5) $$Var\left({\widehat{X}}_L\right)\simeq \text{Var}\left(X\right)\left(\frac{2l_c{\left(X\right)}^D}{L}\right)=\mathrm{var}\left(X\right){\left(\frac{l_c\left(X\right)}{R}\right)}^D,$$(13)

where R = L/2. If we compare Eq. (13) with Eq. (2), we see that instead of having the number of samples, N, we now have the ratio (R/l_c)^D, where R (or L) is usually an observationally accessible quantity. We can thus interpret the ratio (R/l_c)^D as an effective number of “independent” samples. This result is of prime importance in the analysis of fluctuations within any stochastic field.

Furthermore, the correlation length is linked to the value of the power spectrum Ƥ_X(k) of X at k = 0. Indeed $$\begin{array}{c}{\left(l_c\left(X\right)\right)}^D=\frac{1}{2^D{C}_X\left(0\right)}{\displaystyle {\int }_{{ℝ}^D}{C}_X\left(y\right)}\text{d}y,\\ =\frac{1}{2_D{c}_X\left(0\right)}{\mathcal{P}}_X\left(0\right).\end{array}$$(14)

2.1.3 Ergodic Hypothesis and Ergodic Theory

The results of ergodic theory derived above enable us to define under which conditions volumetric averages correspond to statistical averages and to provide a confidence interval for the estimate ${\widehat{X}}_L$ of the expectation of X. However, these results rely on the knowledge of the statistical properties of X and more precisely of its ACF, which is in general not known. To apply this theory to the study of a real field, such as the density field for example, one must use ergodicity as an assumption. The above results can then be used to test the validity of this assumption and the accuracy of the estimates that are derived from it in a self-consistent way.

2.2 Estimates of the Autocovariance Function and Correlation Length

As shown in the previous section, the knowledge of the ACF of X (or of the value of the power spectrum of X at k = 0) is of crucial importance to measure the relevance of a statistical approach in studies of the properties of large (astrophysical) systems. In practice, however, the ACF of X must be evaluated from data.

2.2.1 Reliability of the Estimators of the Auto-Covariance and the Power Spectrum

In most cases, data are drawn from a finite size sample so that the ACF is not reliable at large lag (large scales). To simplify the notation, we now introduce the variable $X_{\mu }=X-\mathbb{E}\left(X\right)=X-\mu$ and define the estimate, for a sample of size L, $${\widehat{C}}_X^L=\frac{1}{{\displaystyle \prod {}_i\left(L-\left|y_i\right|\right)}}{\displaystyle \iint {\displaystyle {\int }_{-R+\frac{\left|y_i\right|}{2}}^{R-\frac{\left|y_i\right|}{2}}{X}_{\text{μ}}\left(u+\frac{y}{2}\right)\text{d}}u}$$(15) $$=\frac{1}{{\displaystyle \prod {}_i\left(L-\left|y_i\right|\right)}}{\displaystyle \iint {\displaystyle {\int }_{-L+\left|y_i\right|}^{L-\left|y_i\right|}{X}_{\text{μ}}\left(\frac{u-y}{2}\right)X_{\text{μ}}\left(\frac{u+y}{2}\right)\frac{\text{d}u}{2^D}.}}$$(16)

This is an unbiased estimate of C_X(y) but its variance is increasing as $\left|y_i\right|\to L$ and eventually becomes very large due to poor sampling. We thus introduce the biased estimate $${\widehat{C}}_{X,L}\left(y\right)=\frac{{\displaystyle \prod {}_i\left(L-\left|y_i\right|\right)}}{L^D}{\widehat{C}}_X^L\left(y\right),$$(17)

which is still a good estimate at small scales compared to L and has a reduced variance. We note however that it is an unbiased estimator of the quantity entering the integral in Eq. (5). Finally, it is also the Fourier Transform of the periodogram S_L which is defined as: $$S_L\left(K\right)=\frac{1}{L^D}{\left|{\displaystyle {\int }_{{\left[-\frac{L}{2},\frac{L}{2}\right]}^D}X\left(y\right)e^{ik.y}dy}\right|}^2.$$(18)

It is the usual estimate of the power spectrum of X, Ƥ_X. It is, however, a biased estimator of the power spectrum Ƥ_X and is only unbiased asymptotically, in the limit L → ∞. Moreover, the variance of the estimator S_L does not vanish as L → ∞ (Papoulis & Pillai 1965), which makes it quite unreliable.

We thus see that, because of the finite size of the sample, one cannot obtain a reliable estimate of the ACF (or of the power spectrum) for all lag values. Furthermore, in many cases, the mean value of X is not known and is replaced in Eq. (16) by its estimate ${\widehat{X}}_L=\widehat{\mu }$, which introduces further, but reasonable, bias (see Papoulis & Pillai 1965 for a more complete discussion).

2.2.2 Periodic Estimators

To get rid of the effect of finite sampling, one may perform simulations in periodic calculation boxes or may artificially add some periodicity to the available data to obtain the following estimate: $${\widehat{C}}_{X\text{,per}}\left(y\right)=\frac{1}{L^D}{\displaystyle {\int }_{{\left[-\frac{L}{2},\frac{L}{2}\right]}^D}\left(X_{\widehat{\mu }}\left(y+u\right)X_{\widehat{\mu }}\left(u\right)\right)}\text{d}u,$$(19)

where $X_{\widehat{\mu }}=X-\widehat{\mu }=X-{\widehat{X}}_L$ and where one makes the identification $X_{\widehat{\mu }}\left(y+nL\right)=X_{\widehat{\mu }}\left(y\right)$. However, in such cases, the spatial average of the estimated ACF is necessarily 0. Indeed $$\begin{array}{lll}{\displaystyle {\int }_{{\left[-\frac{L}{2},\frac{L}{2}\right]}^D}\frac{{\widehat{C}}_{X,\text{per}}\left(y\right)}{L^D}dy}\hfill & =\hfill & {\displaystyle {\int }_{{\left({\left[-\frac{L}{2},\frac{L}{2}\right]}^D\right)}^2}\frac{\left(X_{\widehat{\mu }}\left(y+u\right)X_{\widehat{\mu }}\left(u\right)\right)}{L^{2D}}du\text{\hspace{0.17em}d}y}\hfill \\ \hfill & =\hfill & {\left({\displaystyle {\int }_{{\left[-\frac{L}{2},\frac{L}{2}\right]}^D}\frac{X_{\widehat{\mu }}\left(u\right)}{L^D}\text{d}u}\right)}^2=0,\hfill \end{array}$$(20)

due to the assumption that $X_{\widehat{\mu }}$ is periodic.

Therein lies a significant problem: as the correlation length is defined as an integral over all possible lags, it is not easy to evaluate the reliability of estimates that are obtained in this manner.

Therefore, one traditionally produces an estimate for l_c (or the integral scale l_i) in either of the next two ways. Either one searches for the e^–1 value of the reduced ACF ${\widehat{C}}_X/\text{Var}\left(X\right)$ to obtain an estimate of the correlation length (see e.g. Kleiner & Dickman 1984, 1985), assuming some exponential envelop for the ACF. Or, if the ACF decays fast enough at scales larger than l_c(X), as is the case in turbulence (see previous section), the ACF is then generally extrapolated with a decaying exponential in regions where it becomes non monotonic (see e.g. Batchelor 1953; Reinke et al. 2016, 2018) so one can perform the integral and give a reliable estimate of l_c(X) if l_c(X) «L.

3 Fluctuations and Estimation of the Correlation Length

3.1 Expected Fluctuations in Repeated Trials

Be it for (numerical) experiments that can be repeated several times or for a statistically homogeneous and stationary field, volume averaged quantities fluctuate around their true expectations. In the former case, the volume averaged quantities fluctuate between the different samples while in the latter case they fluctuate in time. As we show in the following, these fluctuations depend on the ratio (l_c/R). By studying these, we thus aim to obtain an accurate estimate of (l_c/R), without having to calculate the ACF.

We consider here the case where one can reproduce several times the same experiment, as can be done for instance with numerical experiments or as can be approximated for clouds with similar conditions. We wish to determine the expected amplitudes of fluctuations of volume averaged quantities between samples.

To each experiment i of the N trials corresponds a value of the estimate ${\widehat{X}}_{L,i}$ defined by Eq. (4). From Bienaymé-Tchebychev inequality, we know that ${\widehat{X}}_{L,i}$ lies around the true expectation $\mathbb{E}\left(X\right)$ within a distance such that, in probability, $$\mathbb{P}\left(\left|{\widehat{X}}_{L,i}-\mathbb{E}\left(X\right)\right|\le m\sigma \left(X\right){\left(\frac{l_C\left(X\right)}{R}\right)}^{D/2}\right)\ge 1-\frac{1}{m^2}.$$(21)

The average over the N trials (or sample average) $${\widehat{X}}_L^{\left(N\right)}=\frac{1}{N}{\displaystyle \sum _{I=1}^N{\widehat{X}}_{L,i}}$$(22)

is obviously a better estimate of $\mathbb{E}\left(X\right)$ as $$\mathbb{P}\left(\left|{\widehat{X}}_L^{\left(N\right)}-\mathbb{E}\left(X\right)\right|\le m\frac{\sigma \left(X\right)}{\sqrt{N}}{\left(\frac{l_c\left(X\right)}{R}\right)}^{D/2}\right)\ge 1-\frac{1}{m^2}.$$(23)

We then expect the ${\widehat{X}}_{L,i}$ to fluctuate around the sample average ${\widehat{X}}_L^{\left(N\right)}$ with variance $$\text{Var}\left({\widehat{X}}_{L,i}-{\widehat{X}}_L^{\left(N\right)}\right)=\sigma {\left(X\right)}^2{\left(\frac{l_c\left(X\right)}{R}\right)}^D\left(1-\frac{1}{N}\right),$$(24)

such that $$\mathbb{P}\left(\left|{\widehat{X}}_{L,i}-{\widehat{X}}_L^{\left(N\right)}\right|\le m\sigma \left(X\right){\frac{L_C\left(X\right)}{R}}^{D/2}{\left(1-\frac{1}{N}\right)}^{1/2}\right)\ge 1-\frac{1}{M^2}.$$(25)

If l_c(X) is known, Eqs. (21), (24) and (25) allow one to give statistical error bars. Conversely, if l_c(X) is not known, these equations give information on the product σ(X)(l_c(X)/R)^D/2 by performing the same statistical experiment several times and studying the dispersion of the ${\widehat{X}}_{L,i}$ around ${\widehat{X}}_L^{\left(N\right)}$.

Indeed, the half length l_50%, of the segment centered on ${\widehat{X}}_L^{\left(N\right)}$ within which 50% of the estimate ${\widehat{X}}_{L,i}$ falls, verifies² $$\sqrt{2}\sigma \left(X\right){\left(\frac{l_c\left(X\right)}{R}\right)}^{D/2}{\left(1-\frac{1}{N}\right)}^{1/2}\underset{˜}{<}{l}_{50\%.}$$(26)

This gives a quick and easy qualitative determination of σ(X)(l_c(X)/R)^D/2. More quantitatively, the empirical variance of the sample of the N trials $${\text{Var}}^{\left(N\right)}=\frac{1}{N-1}{{\displaystyle \sum _{i=1}^N\left({\widehat{X}}_{L,i}-{\widehat{X}}_L^{\left(N\right)}\right)}}^2,$$(27)

is an unbiased estimator of $$\text{Var}\left({\widehat{X}}_{L,i}\right)=\sigma {\left(X\right)}^2{\left(l_c\left(X\right)/R\right)}^D.$$(28)

Computing the variance Var^(N) thus yields an easy and rigorous method to determine (l_c(X)/R)^D and the correct error bars for statistical experiments.

3.2 Expected Temporal Fluctuations for a Statistically Homogeneous and Stationary Field

We consider here the case of a statistically homogeneous and stationary field (i.e., whose statistical properties are invariant under space and time translations). This can for example describe a steady turbulent flow (e.g. as simulated in a periodic box) without gravity.

As before, to each time t corresponds a value of the estimate ${\widehat{X}}_L\left(t\right)$ defined by Eq. (4) which lies around the true expectation $\mathbb{E}\left(X\right)$ within a distance such that, in probability, $$\mathbb{P}\left(\left|{\widehat{X}}_L\left(t\right)-\mathbb{E}\left(X\right)\right|\le m\sigma \left(X\right){\left(\frac{l_c\left(X\right)}{R}\right)}^{D/2}\right)\ge 1-\frac{1}{m^2}.$$(29)

The time average over a large timescale T $${\widehat{X}}_{L,T}=\frac{1}{T}{\displaystyle {\int }_{to}^{to+T}{\widehat{X}}_L\left(t\right)\text{d}t,}$$(30)

where t₀ is a time at which the steady state is reached, is a better estimate of $\mathbb{E}\left(X\right)$ as it has a variance $$\mathrm{var}\left({\widehat{X}}_{L,T}\right)\simeq \sigma {\left(\widehat{X}\right)}^2{\left(\frac{l_c\left(X\right)}{R}\right)}^D\frac{{\tau }_{c\left(X_L\right)}}{2T},$$(31)

where τ_c(X_L) « T is the correlation time. Then, the signal ${\widehat{X}}_L\left(t\right)$ will fluctuate around ${\widehat{X}}_{L,T}$ with an empirical (temporal) variance $${\mathrm{var}}_T=\frac{1}{T}{\displaystyle {\int }_{to}^{to+T}{\left({\widehat{X}}_L\left(t\right)-{\widehat{X}}_{L,t}\right)}^2\text{d}t}$$(32)

which, providing that τ_c(X_L) « T, is an accurate estimate of $$\mathrm{var}\left({\widehat{X}}_L\left(t\right)\right)=\sigma {\left(X\right)}^2{\left(l_c\left(X\right)/R\right)}^D.$$(33)

Hence, computing the variance Var_T yields an easy and robust estimate of (l_c(X)/R)^D and the correct error bars for a statistically stationary field.

3.3 Fluctuations of Integrated Fields Overa Column

The previous cases work for experiments that can be repeated or for statistically stationary fields. However, in some situations, the two conditions cannot be fulfilled either because it is impossible to reproduce the experiment a large number of times or because the fields are not stationary (for instance in the presence of gravity).

In that case an estimate of ratio (l_c/R) can be obtained if one has access to the integral of field X over a column of fixed length L = 2R: $${\displaystyle {\sum }_x\left(r\right)}={\displaystyle {\int }_{\left[0,L\right]}X\left(r,z,\right)}\text{dz},$$(34)

where r is a vector of D − 1 dimension (typically 2). The column must have a constant length to avoid creating spurious biases (see Sect. 5).

As ∑_X/L corresponds to averaging X along one direction, we thus expect that its fluctuations will be reduced in comparison of those of X. The longer the length L of the column, the smaller the fluctuations of ∑_X/L are expected to be. More quantitatively, the variance of ∑_X/L is $$\mathrm{var}\left(\frac{{\displaystyle \sum {}_X}}{L}\right)=\frac{\mathrm{var}\left({\displaystyle \sum {}_X}\right)}{L^2}=\frac{\mathbb{E}\left({\displaystyle \sum {}_X\left(r\right)}{\displaystyle \sum {}_X\left(r\right)}\right)}{L^2}$$(35) $$=\frac{1}{L^2}{\displaystyle {\int }_{{\left[0,L\right]}^2}{C}_X\left(0,z-z^{\prime }\right)}\text{d}z\text{d}{z}^{\prime }$$(36) $$=\frac{1}{L}{\displaystyle {\int }_{\left[-L,L\right]}{C}_{\rho }\left(0,u\right)}\left(1-\frac{\left|u\right|}{L}\right)\text{d}u.$$(37)

A similar equation for the centroid velocities was given in (Scalo 1984). Now, if the ACF of X is isotropic at short lags and l_c(X) « L, one can make the approximation $$\mathrm{var}\left(\frac{{\displaystyle \sum {}_X}}{L}\right)\simeq \frac{1}{L}{\displaystyle {\int }_{\left[-L.L\right]}C\rho \left(\left|u\right|\right)\text{d}u,}$$(38) $$\simeq \mathrm{var}\left(X\right)\frac{2l_i\left(X\right)}{L},$$(39)

where l_i(X) is the integral scale (Batchelor 1953) which, in most cases, verifies l_i(X) ≃ l_c(X) (see Jaupart & Chabrier 2021). Thus a quick and easy estimate of ratio l_c(X)/R is given by $$\frac{\text{Var}\left({\displaystyle \sum {}_X}\right)}{L^2\mathrm{var}\left(X\right)}\simeq \frac{lc\left(X\right)}{R}.$$(40)

This method was applied to the density field (ρ = X) in (Jaupart & Chabrier 2021). The above estimate (Eq. (40)) was shown to produce the trends predicted analytically and thus to be a good approximation of the actual ratio l_c(ρ)/R.

4 Application to Astrophysical Fields

The general results derived in Sect. 2 can be applied to many physical and astrophysical systems. They have been used extensively in cosmology but have somehow been overlooked in studies of star formation. Today, it is generally accepted that star formation is triggered by density fluctuations generated by compressible turbulence injected at a large scale in MCs (see e.g. McKee & Ostriker 2007 and reference therein). In this context, the density field ρ (or the logarithmic density field $s=\ln \left(\rho /\mathbb{E}\left(\rho \right)\right)$) is of prime interest and its cumulative distribution function (CMF) and probability density function (PDF) must be determined accurately.

Each of these statistical quantities is associated with a stochastic field X to which the results of Sect. 2 can be applied. For instance, the CMF F_ρ(ρ₀) at ρ₀ is linked to the stochastic field h_ρ0(y) = Θ (ρ₀ − ρ(y)) (where Θ is the Heavyside function), as $$F_{\rho }\left({\rho }_0\right)\Delta \rho =\mathbb{E}\left(h_{{\rho }_0}\left(y\right)\right)-\mathbb{E})\left(h_{{\rho }_0}\right),$$(41)

while the PDF f_ρ(ρ0) is given by $f_{\rho }\left({\rho }_0\right)=\mathbb{E}\left({\delta }_{{\rho }_0}\left(y\right)\right)=\mathbb{E}\left({\delta }_{{\rho }_0}\right)$ where δ_ρ0(y) = δ (ρ0 − ρ(y)) is Dirac’s distribution. Usually the PDF is rather deduced from histograms with some bin size ∆ρ such that $$f_{\rho }\left({\rho }_0\right)\Delta \rho \simeq F_{\rho }\left({\rho }_0+\Delta \rho \right)-F_{\rho }\left({\rho }_0\right)=\mathbb{E}\left(h_{{\rho }_0+\Delta \rho }-h_{\rho 0}\right).$$(42)

In principle, knowledge of the ACF of all these fields is required to establish the accuracy of the estimations. Fortunately, it can be shown that sometimes, with a few simplifying assumptions, one can proceed with the ACF of ρ, only, in some situations. This is explained in detail in Appendix B.

4.1 Exact Results Regarding the Properties of the Auto-Covariance Function (ACF) ofρ

For a statistically homogeneous field, the ACF of ρ, the density field, can be expressed in term of the second order structure function: $$S_{\rho }^{\left(2\right)}\left(y\right)=\mathbb{E}\left({\left\{\rho \left(u+y\right)-\rho \left(u\right)\right\}}^2\right),$$(43)

as $S_{\rho }^{\left(2\right)}\left(y\right)=2\left(C_{\rho }\left(0\right)-C_{\rho }(y\right)$. A similar statement can be made for the logarithmic density field $s=\ln \left(\rho /\mathbb{E}\left(\rho \right)\right)$. This helps us to grasp some key features of the ACF. At very short scale (below the viscous scale), the density field is supposed to be differen-tiable and hence C_ρ must possess second-derivatives at y = 0. Then, due to the parity of the ACF and because it is maximal at y = 0, its gradient must exist and be equal to 0 at y = 0.

Furthermore, (Jaupart & Chabrier 2021), generalizing the work of (Chandrasekhar 1951a), show that for a statistically homogeneous density field, the quantity $$\mathbb{E}\left(\rho \right)\mathrm{var}\left(e^s\right)l_c{\left(\rho \right)}^3=\frac{\mathrm{var}\left(\rho \right)}{\mathbb{E}\left(\rho \right)}{l}_c{\left(\rho \right)}^3$$(44)

is an invariant of the dynamics.

4.2 Phenomenology of (Compressible) Turbulence

The phenomenology of compressible turbulence (Kritsuk et al. 2007) can be derived, with some adjustments, from that of incompressible turbulence (Frisch 1995). Thus, we use the latter to derive some expected features of the density ACF in star-forming MCs that can be described by such phenomenology.

In isotropic turbulence, the second order structure function is observed to be a monotonic increasing function of separation distance, at least in the inertial range, and to converge rapidly toward 2Var (ρ) at scales that are larger than the integral scale l_i. This integral scale (not to be confused with the injection scale of turbulence, see Sect. 4.3) is defined in the same manner as the correlation length (Batchelor 1953) $$l_i=\frac{1}{c\left(0\right)}{\displaystyle {\int }_0^{\infty }C\left(r\right)}\text{d}r.$$(45)

In many situations, l_c ~ l_i, as shown in (Jaupart & Chabrier 2021). Thus, at small scales (short lags) and in the inertial range, the ACF must be a monotonically decreasing function. Above the inertial range, it is often assumed that the structure function and the ACF are still monotonic and the ACF is usually approximated by a decaying exponential, even though density fluctuations are likely to generate oscillations of the observed and estimated ACF as it tends to zero (Batchelor 1953; Reinke et al. 2016, 2018).

In compressible isothermal and stationary turbulence, the density field ρ is found to be approximately lognormal (Kritsuk et al. 2007; Federrath et al. 2010), implying that the logarithmic density field $s=\ln \left(\rho /\mathbb{E}\left(\rho \right)\right)$ is Gaussian with variance $\sigma {\left(s\right)}^2=\ln \left(1+{\left(bℳ\right)}^2\right)$. In such Gaussian conditions, the ACFs of ρ and s are linked by the following equation: $$C_{\rho }\left(y\right)=\mathbb{E}{\left(\rho \right)}^2\left(e^{c_s\left(y\right)}-1\right).$$(46)

As a consequence, if C_ρ (or C_s) is monotonically decaying toward 0, we deduce that: $${\left(\frac{\sigma {\left(s\right)}^2}{e^{\sigma {\left(s\right)}^2}-1}\right)}^{1/3}{l}_c\left(s\right)\le l_c\left(\rho \right)\le l_c\left(s\right),$$(47)

where we have used the following inequalities: e^ax − 1 ≤ x(e^a − 1) for 0 ≤ x ≤ 1 and ax ≤ e^ax − 1 ∀x. For typical star forming conditions, σ(s)² ≲ 4, implying that: $$0.4l_c\left(s\right)\underset{˜}{<}{l}_c\left(\rho \right)\le l_c\left(s\right),$$(48)

or, $$l_c\left(s\right)\sim l_c\left(\rho \right).$$(49)

This shows that under Gaussian conditions, for the two lengths l_c(s) and l_c(ρ), knowledge of one of them is sufficient to characterize the other one within an order of magnitude.

4.3 Large Injection Scale but small Correlation Length

Density and velocity fluctuations in MCs are thought to originate from turbulent motions driven at large scale (McKee & Ostriker 2007; Brunt et al. 2009), i.e. at scale comparable to the cloud scale L. This means that the energy of these turbulent motions is injected at an injection scale l_inj ~ L, below which the turbulent cascade eventually occurs.

The injection scale l_inj, however, is not the correlation length l_c of either the velocity, kinetic energy or density fields. Otherwise, if l_inj = l_c, every estimate produced from volumetric averages of the former fields would be inaccurate and far from the actual statistical values. This could result in large fluctuations of these averaged quantities either between different simulations (samples) or at different times for steady turbulent flows (see Sect. 3). What is instead observed in numerical simulations of compressible and steady turbulent flow is that volume averaged quantities of these fields display rather small fluctuations around their mean values. This is the case, for example, for the rms Mach number $ℳ={\upsilon }_{\text{rms}}/c_s$, where υ_rms is the root of the volume average square velocity υ² and c_s the sound speed. To be more explicit, $${ℳ}^2=\frac{1}{L^3}{\displaystyle {\int }_{{\left[-\frac{L}{2}.\frac{L}{2}\right]}^3}\frac{v^2\left(y\right)}{c_s^2}}\text{d}y,$$(50)

so the result of Sect. 3 can be applied with $X={\upsilon }^2/c_s^2$ and ${ℳ}^2={\widehat{X}}_L$.

In (Federrath 2013), a series of numerical simulations of isothermal compressible turbulence driven to $ℳ\simeq 17$ with an injection scale l_inj = L/2 = R is presented. Once a statistical steady state is reached, the volume averaged Mach number $ℳ$ (or ${ℳ}^2$ which is a measure of the volume averaged specific kinetic energy) displays fluctuations that are rather small compared to their average values (their Fig. 1). Would we have l_c(υ²) = l_inj = R, fluctuations of the signal ${ℳ}^2\left(t\right)$ would have yielded a temporal variance (Eq. (32)): $$Va{r}_T\simeq {\sigma }^2\left(v^2\right)\left(\frac{l_c{\left(v^2\right)}^3}{R}\right)={\sigma }^2\left(v^2\right),$$(51)

from Eq. (33). Since the statistics of υ are close to being Gaussian (their Fig. A1), this would imply $$Va{r}_T\simeq 2{ℳ}_T^4,$$(52)

where ${ℳ}_T$ is the average of the signal $ℳ\left(t\right)$ over a time T. This would thus yield large fluctuations incompatible with their Fig. 1. The actual temporal variance Var_T of signal $ℳ{\left(t\right)}^2$ in (Federrath 2013) rather yields a ratio $$\frac{l_c\left(v^2\right)}{R}=\frac{l_c\left(v^2\right)}{l_{inj}}\simeq 0.1,$$(53)

which shows that l_c(υ²) « l_inj = R.

Furthermore, (Jaupart & Chabrier 2021) used the estimate produced by Eq. (40) to compute the correlation length of the density field ρ. They found that within a factor of order unity, $$l_c\left(\rho \right)\simeq {\lambda }_s=L/{ℳ}^2\ll L=2_{\text{linj}},$$(54)

where λ_s is the sonic length which is found to be close to the average width of filamentary structures in isothermal turbulence (Federrath 2016).

The above results show that a large injection scale does not imply a large correlation length and that, on the contrary, correlation lengths in star-forming MCs are small compared with the injection scale (see above and Jaupart & Chabrier 2021).

4.4 Practical Assumptions Regarding the ACF

From the above results, we thus assume that the ACF decays rapidly at scales larger than the correlation length l_c (l_c ~ l_i, the integral scale) and then that the defining integral Eq. (11) can be calculated only up to a few l_c. Moreover, we assume that the ACF can be bounded by a decaying exponential exp(−|y|/λ), where λ ~ l_c above and in the inertial range to allow the computation of the correlation length (we note that such an exponential behavior is prohibited at very small scales due to the differentiability of ρ).

5 Star-Forming Clouds. Column Densities as Tracers of the Underlying Density Field

We now turn to observations of star-forming molecular clouds. Measurements provide values of the column-density ∑(x, y), which is the integral of density along the line of sight (l.o.s.(x, y)): $$\begin{array}{lll}{\displaystyle \sum {}_{\left(x,y\right)}}\hfill & =\hfill & {\displaystyle {\int }_{\text{l,o,s}\left(x,y\right)}\rho \left(x,y,z\right)\text{d}z}\hfill \\ \hfill & =\hfill & \mathbb{E}\left(\rho \right)l\left(x,y\right)+{\displaystyle {\int }_{\text{l.o,s}\left(x,y\right)}\delta \rho \left(x,y,z\right)\text{d}z\text{,}}\hfill \end{array}$$(55)

where l(x, y) is the thickness of the cloud along the l.o.s. at (x, y) and $\delta \rho =\rho -\mathbb{E}\left(\rho \right)$ is the density fluctuation. Column densities are the only data that depend directly on the density field and one must determine how to retrieve reliable information from them.

5.1 Inhomogeneity and Anisotropy due to Integration over the Line of Sight

Star forming clouds are shaped by turbulent motions conferring statistical properties to their geometrical characteristics, and hence to the area projected in a plane perpendicular to the line of sight and to the thickness projected along the line of sight. This is responsible for difficulties in evaluating exactly the statistical average of ∑(x, y). However, provided that the cloud thickness is much larger than the correlation length, that is if l(x, y) » l_c(ρ), we can reasonably assume that (see Eq. (55)): $$\mathbb{E}\left({\displaystyle \sum \left(x,y\right)}\right)\simeq \mathbb{E}\left(\rho \right)\mathbb{E}\left(x,y\right)).$$(56)

One must note here that we are dealing with the statistical average and not with the spatial average. This equation shows that ∑(x, y) may not be statistically homogeneous even if the density field ρ is, just because of integration effects. To illustrate this important point, let us imagine two idealized situations. In one of them, the cloud is a sphere of radius R. In the other one, the cloud is a “cube” of side L misaligned with the line of sight and seen from one of its edges such that the projected surface is of size $\sqrt{2}L\times L$ (see Fig. 1). For the sphere, the thickness along the line of sight is: $$\mathbb{E}\left(l_{s\left(\text{R}\right)}\right)\left(x,y\right)=2R\left(1-\frac{x^2+y^2}{R^2}\right),x^2+y^2<R^2,$$(57)

whereas for the cubic cloud it is: $$\mathbb{E}\left(l_{c\left(\text{L}\right)}\right)\left(x,y\right)=\sqrt{2}L\left(1-\frac{\sqrt{2}\left|x\right|}{L}\right),\left|x\right|\le \frac{L}{\sqrt{2}},\left|y\right|\le \frac{L}{2}.$$(58)

Even though they are very simple, these two examples demonstrate that the column-density field may exhibit large-scale gradients and hence may not be statistically homogeneous, even if the density field is. Furthermore, as seen with the example of the cube, integration effects can also generate some anisotropy in the column-density field.

To reduce these effects, one can use a low pass filter to filter out large-scale gradients (Kleiner & Dickman 1984) and then treat the column density field as if it were homogeneous. Furthermore, most of the integration effects are expected to be produced by the first term of the r.h.s of Eq. (55). Thus they are expected to affect column densities that are around the (surface) average of the column density map 〈∑〉. In contrast, high column density (∑(x, y) > 〈∑〉, the regions of interest for star formation) are most likely to originate from the second term of the r.h.s of Eq. (55) and be produced by dense pockets along the line of sight. They are thus expected to be less affected by the integration effects and by the low pass filter. Studying the statistics of these high column density is thus expected to bear insights on the bias introduced by integration effects.

In Sect. 6, we will apply the above considerations to the observations of the Polaris cloud.

5.2 Column-Density Field in a Simulation Box

For a cubic simulation domain of size L, projecting the density field along one of the 3 principal directions of the cube leads to a statistically homogeneous column density field such that: $$\mathbb{E}\left({\displaystyle \sum \left(x,y\right)}\right)=\mathbb{E}\left(\rho \right)\times L.$$(59)

The results of Sects. 2 and 3 can thus be applied with X = ∑ and the ACF of ∑ in this case is $$\begin{array}{lll}{C}_{\sum }\left(r\right)\hfill & =\hfill & \mathbb{E}\left(\left({\displaystyle \sum \left(u+r\right)}-\mathbb{E}\left(\rho \right)L\right)\left({\displaystyle \sum \left(u\right)}-\mathbb{E}\left(\rho \right)L\right)\right)\hfill \\ \hfill & =\hfill & {\displaystyle {\int }_{{\left[-L/2,L/2\right]}^2}C\rho \left(r,z-z^{\prime }\right)\text{dzd}{\text{z}}^{\prime }}\hfill \\ \hfill & =\hfill & {\displaystyle {\int }_{\left[-L,L\right]}C\rho \left(r,u\right)}\text{d}u{\displaystyle {\int }_{-L+\left|u\right|}^{L-\left|u\right|}\frac{\text{d}v}{2}}\hfill \\ \hfill & =\hfill & L{\displaystyle {\int }_{\left[-L,L\right]}{C}_{\rho }\left(r,u\right)\left(1-\frac{\left|u\right|}{L}\right)\text{d}u.}\hfill \end{array}$$(60)

Its variance is $$\text{Var}(\sum )=C_{\sum }\left(0\right)=L{\displaystyle {\int }_{\left[-L,L\right]}{C}_{\rho }\left(0,u\right)\left(1-\frac{\left|u\right|}{L}\right)\text{d}u}.$$(61)

Then, if the density field is statistically isotropic at small scales (i.e. the ACF is isotropic at short lag) and l_c(ρ) « L, $$\text{Var}(\sum )\simeq 2L{l}_c\left(\rho \right)\text{Var}\left(\rho \right).$$(62)

Then Eq. (62) yields: $$\text{Var}(\sum )\simeq \text{Var}\left(\frac{\rho }{\text{E}(\sum )}\right)\simeq \text{Var}\left(\frac{\rho }{\text{E}\left(\rho \right)}\right)\frac{2l_c\left(\rho \right)}{R},$$(63)

which is a reformulation of Eq. (40). This is an important result because it gives a measure of l_c(ρ)/R without having to compute the ACF. (Brunt et al. 2010); (Federrath et al. 2010) for example found a ratio Var $\left(\text{Σ}/\mathbb{E}\left(\text{Σ}\right)\right)/\text{Var}\left(\rho /\mathbb{E}\left(\rho \right)\right)$ between 0.03 and 0.15.

(Vázquez-Semadeni & García 2001) were the first to study the impact of the l_c(ρ)/R ratio on the statistics of column-density fields. Based on a crude interpretation of the central limit theorem (CLT), they proposed that, for l_c(ρ)/R → 0, the column-density PDF should appear to be Gaussian instead of lognormal. This is not consistent with the apparent lognormality of the observed column-density PDFs, which led these authors to conclude that l_c(ρ)/R cannot be vanishingly small and that it must be on the order of 10^–1. However, the CLT only applies to independent variables and can hardly be valid for the sum ofcor-related variables, even if correlations decay. This casts doubt on the conclusions of (Vázquez-Semadeni & García 2001). More recently, (Szyszkowicz & Yanikomeroglu 2009) and (Beaulieu 2011) have shown that, for some special types of correlations, the sum of a large number N of lognormal variables tends to a lognormal distribution as N → ∞. We conclude that knowledge of the l_c(ρ)/R value does not allow robust conclusions on the shape of the column-density PDF. However, as shown by Eq. (63), the variance Var $\left(\frac{\text{Σ}}{\mathbb{E}\left(\text{Σ}\right)}\right)$ does become vanishingly small as l_c(ρ)/R tends to zero. In that case, one can show with high probability that: $$\ln \left(\frac{\sum }{\text{E}(\sum )}\right)\simeq \frac{\sum -\text{E}(\sum )}{\text{E}(\sum )}.$$(64)

Thus, in the limit of vanishing values of l_c(ρ)/R, the distributions of $\text{Σ}/\mathbb{E}\left(\text{Σ}\right)$ and its logarithm are both Gaussian if one of them is.

Fig. 1 Projection of the two idealized situation. Left panel: case of a sphere. Right panel: case of a cuboid mis-aligned with the line of sight.

5.3 Decay Length of Correlations

We now examine how the decay of correlations of ρ impacts the decay of correlations of ∑. For sake of simplicity, we again consider the case of a cubic box in order to avoid unncessary complications. For the 2D field ∑, the correlation length is given by: $$\begin{array}{c}{l}_c{(\sum )}^2=\frac{1}{4}\frac{1}{\text{Var}(\sum )}{\displaystyle \iint C_{\sum \left(r\right)}\text{d}r,}\\ =\frac{1}{4}\frac{1}{\text{Var}(\sum )}{\displaystyle \iint L}{\displaystyle {\int }_{\left[-L,L\right]}{C}_{\rho }\left(r,u\right)\left(1-\frac{\left|u\right|}{L}\right)\text{d}u\text{\hspace{0.17em}d}r},\\ \simeq 2\frac{L\text{Var}\left(\rho \right)}{\text{Var}(\sum )}{l}_c{\left(\rho \right)}^3.\end{array}$$(65)

Using Eq. (62), this implies that: $$l_c{(\sum )}^2\simeq l_c{\left(\rho \right)}^2.$$(66)

This shows that correlations of the column-density fields are decaying over a characteristic length close to l_c(ρ), the correlation length of the underlying density field. In general, we can thus assume that l_c(∑) ~ l_c(ρ), so that information gathered from the column-density yields an estimate of the characteristic decay length of correlations of the underlying density field ρ.

6 Application to the Observations of the Polaris Cloud

As mentioned in Sect. 5, observations trace back the column-density (Kleiner & Dickman 1984; Schneider et al. 2015; Ossenkopf-Okada et al. 2016). These observations of the column-density PDFs in MCs show that regions where star formation has not occurred yet exhibit lognormal PDFs while regions with numerous prestellar cores develop power-law tails (PLTs) at high column densities (Kainulainen et al. 2009; Schneider et al. 2013). In addition to the integration effects yielding potentially the observed column-density to be anisotropic and inhomogeneous, observational data suffer further biases due to line of sight (l.o.s) contamination and noise (Schneider et al. 2015; Ossenkopf-Okada et al. 2016). L.o.s contamination causes two important biases. The observed power-law tail appears to be steeper than its corrected and uncontaminated counterpart while the observed variance in the lognormal part appears to be smaller than its corrected counterpart (Schneider et al. 2013). The overall effect of l.o.s contamination is to produce an underestimation of the total variance of the column-density.

6.1 Polaris

As a typical example of initial conditions of star formation in MCs, we focus on the Polaris flare, where line of sight contamination appears to be negligible (André et al. 2010; Miville-Deschênes et al. 2010; Schneider et al. 2013). Furthermore, most of the stellar cores in this cloud are still unbound (André et al. 2010), showing that star formation activity is very recent. Polaris is therefore a good candidate to probe the statistics of initial phases of star formation in MCs.

Data from Herschel Gould Belt survey extend across part of this cloud over approximately a 10 square degrees region with a linear size L ~ 10 parsecs (pc) (André et al. 2010). The cloud total mass and area above an extinction A_v ≥ 1 are $M_{c,A_{\text{v}}\ge 1}=1.21\times {10}^3{M}_{\odot }$ and $A_{c,A_{\text{v}}\ge 1}=3.9{\text{pc}}^2$, respectively. Dust temperatures are in a narrow T_dust = 13 + 1K interval, indicating fairly isothermal conditions with an average Mach-number $ℳ\simeq 3$ (Schneider et al. 2013).

The Polaris logarithmic column-density field ŋ, where ŋ = ln(∑/ 〈∑〉), has a PDF with an extended Gaussian part and two emerging power-law tails, a first one with exponent α_ŋ,1 ≃ −4 followed by a shallower one with exponent α_ŋ,2 ≃ −2 (Fig. 2). (Jaupart & Chabrier 2020) have shown that the first steep PLT is due to gravity beginning to affect turbulence in parts of the cloud and hence records an early stage of (local) collapse. Furthermore the authors outlined a procedure to reconstruct the underlying logarithmic volume density PDF, noted s-PDF, where $s=\ln \left(\rho /\mathbb{E}\left(\rho \right)\right)$, from data on the η-PDF. The underlying s-PDF displays a Gaussian part and two PLTs with exponents α_S,1 = −2 and α_S,2 = 3/2, respectively (see Fig. 2).

Fig. 2 Column and volume density PDF of the Polaris cloud. Left: Observed logarithmic column-density (ŋ = ln(∑/ 〈∑〉)-PDF (Schneider et al. 2013; Jaupart & Chabrier 2020). Right: Estimated and reconstructed underlying logarithmic density $\left(s=\ln \left(\rho /\mathbb{E}\left(\rho \right)\right)\right)$-PDF with the procedure from (Jaupart & Chabrier 2020).

6.2 Filtering Large-Scale Gradients

As seen in Sect. 5, integration effects can produce large-scale gradients and break statistical homogeneity as well as isotropy in the column density field.

Filtered and unfiltered column-density maps of the Polaris flare are displayed in Figs. 3 and 4. The low pass filter does not alter qualitatively the intricate structures that exist, while the high-pass filter reveals a large-scale gradient likely due to integration effects. In order to partially reduce measurement artifacts, we use a low pass filter that screens out structures larger than L/2 in the column-density contrast (∑ − 〈∑〉), where L is the linear size of the observed region and we recall that 〈∑〉 is the (surface) average of the column density map. We can then treat the column-density field as if it was homogeneous.

The low pass filter slightly diminishes the variance Var(∑/〈∑〉) which is ≃0.20 and ≃0.17 for the unfiltered and low pass filtered data, respectively. It barely affects structures with a positive column-density contrasts but increases the occurrence of highly negative column-density contrasts. This is seen in Fig. 5, which portrays the η-PDFs of the unfiltered and low pass filtered column-density maps. As mentioned in Sect. 5, high column density are not expected to be sensitive to integration effects.

6.3 Estimated ACF and Correlation Length

6.3.1 Correlation Length of Ŋ from the ACF

We now estimate the ACFs of the logarithmic column-density field ŋ = ln(∑/ 〈∑〉) for the three data sets (unfiltered, low and high pass filtered), using Eq. (17). The 2D heat-maps of the reduced ACFs ${\widehat{C}}_{\eta }/\text{Var}\left(\eta \right)$ are given in Fig. 6. We only display the top-right quadrant of possible lags (x > 0, y > 0) which amounts to half of the space useful to study the ACF due to its symmetry. The high pass filtered ACF illustrates the bias that can be introduced by integration effects.

The three ACFs all seem to be fairly isotropic at very short lags (scales) but are anisotropic at large ones. The low pass filtered ACF seems to decay more rapidly at short lags with a reduce anisotropy than the unfiltered one. Figure 7 displays the reduced ACF of the low pass filtered map in three different directions, x (θ = 0), x = y (θ = π/4) and y (θ = π/2). As can be seen from the heat maps but also from Fig. 7, a strong anisotropy is detected at large scales in the x direction (x/L ≥ 2 × 10^–2), while the ACF if fairly isotropic at shorter lags. From the y-direction to the π/4-direcťıon, the data seem to be fairly isotropic and bounded by an exponential with λ/L ≃ 5× 10^–2. Anisotropy is most pronounced along the x-direction and the resulting estimated correlation length ${\widehat{l}}_c\left(\eta \right)$ is: $${\widehat{l}}_c\left(\eta \right)\simeq 6\times {10}^{-2}L\simeq \frac{1}{2}{\left(2\pi \right)}^{1/2}\lambda ,$$(67)

or ${\widehat{l}}_c\left(\eta \right)/R\simeq 1.2\times {10}^{-1}$, thus ${\widehat{l}}_c\left(\eta \right)/R\sim {10}^{-1}$. We then use ${\widehat{l}}_c\left(\eta \right)$ as an estimate of l_c(ρ) to within an order of magnitude, such as $$l_c\left(\rho \right)~{10}^{-1}R.$$(68)

In fact, we expect Eq. (67) to provide upper bounds for ratios l_c(ŋ)/R and l_c(ρ)/R, because integration artifacts are only partially cancelled by the low pass filter.

6.3.2 Correlation Length of ŋ from Eq. (40)

To obtain an estimate of l_c(ŋ) we could also apply the results of Sect. 3 to the low pass filtered map. By integrating the column density map along the x (θ = 0) or y (θ = π/2) direction and computing the variance of the resulting integrated field we can obtain with Eq. (40) two estimates ${\widehat{l}}_{c,x}\left(\eta \right)$ and ${\widehat{l}}_{c,y}\left(\eta \right)$ of l_c(ŋ) within a factor of order unity. However, the estimated ACF displays a strong anisotropy in the x direction at large scales (x/L ≥ 2 10^–2), so we expect the estimates to give rather different results. Computing the estimates yields $${\widehat{l}}_{c,x}\left(\eta \right)\simeq 2.6\times {10}^{-1}R$$(69) $${\widehat{l}}_{c,y}\left(\eta \right)\simeq 3.5\times {10}^{-2}R,$$(70)

with ${\widehat{l}}_{c,x}\left(\eta \right)>{\widehat{l}}_{c,y}\left(\eta \right)$, as can be expected from the anisotropy of the observed ACF. Since the anisotropy of the ACF starts before it is significantly smaller than the variance, it is not clear which of the two estimates gives the better approximation of the actual l_c(ŋ). However, they are both within a factor 3 of the estimate produced by Eq. (67) from the ACF which yielded ${\widehat{l}}_c\left(\eta \right)\sim {10}^{-1}R$.

Using l_c(ŋ) as an estimate of l_c(ρ) yields again, to within an order of magnitude, l_c(ρ) ≲ 10^–1 R.

6.3.3 Correlation Length of ρ from Eqs. (40) or (63) and the Variance of ∑

As discussed in Sect. (5.2) and Eq. (63), one can also estimate the ratio l_c(ρ)/R by (1) computing the variance Var $\left(\text{Σ}/\mathbb{E}\left(\text{Σ}\right)\right),\left(2\right)$ giving an estimate of Var $\left(\rho /\mathbb{E}\left(\rho \right)\right)$ and (3) giving an estimate of the average thickness of the cloud (the length of the line of sight), for example by assuming that the cloud has roughly the same dimension in the three directions.

In pure isothermal turbulence, Var $\left(\rho /\mathbb{E}\left(\rho \right)\right)\approx {\left(bℳ\right)}^2$, which is ≃ 1 for the Polaris case $\left(b\simeq 0.3-0.4,ℳ\simeq 3\right)$. However, when gravity starts generating power-law tails in the density PDF, the variance becomes larger than ${\left(bℳ\right)}^2$ (Jaupart & Chabrier 2020). For Polaris, the column-density PDF displays a power-law tail with exponent α_ŋ≃−4, which is linked to an underlying density PDF with a power law tail exponent α_s≃−2 (Federrath & Klessen 2013; Jaupart & Chabrier 2020). Using the reconstructed s-PDF of Fig. (2) from the procedure described in (Jaupart & Chabrier 2020), we can derive an estimate of Var$\left(\rho /\mathbb{E}\left(\rho \right)\right)$. In principle, for such a model PDF, the variance is infinitely large due to the power-law tails exponents α_s,1=−2 and α_s,2=−3/2. However, we expect a cut-off at high (column)-density, which is indeed visible in the data. This cutoff may be due to a change of thermodynamic conditions of the cloud, e.g. from isothermal to adiabatic conditions. For a typical cut-off number-density n_ad = 10¹⁰ cm^–3 (Masunaga & Inutsuka 2000; Machida et al. 2006; Vaytet et al. 2013, 2018) and for a cloud of average density $\overline{n}={10}^3{\text{cm}}^{-3}$, the cutoff occurs at s_ad ≃ 16. However, there may be other causes for a high density cut-off. In order to assess this possibility, we thus determine three different estimates of the variance Var $\left(\rho /\mathbb{E}\left(\rho \right)\right)$ from the reconstructed s-PDF of Fig. 2: one densities up to 6.3 (s ≤ 6.3), which corresponds to the onset of the 2nd PLT, a second one for s ≤ 8 in order to include contributions from the 2nd PLT, and a third one for s ≤ 16 ≃ s_ad in order to include all the data up to the adiabatic limit. We obtain Var $\left(\rho /\mathbb{E}\left(\rho \right)\right)\simeq 5,7,227$, respectively, such that: $${\widehat{f}}_{X,L}\left({\xi }_0\right)\Delta \xi \simeq {\widehat{F}}_{X;L}\left({\xi }_0+\Delta \xi \right)-{\widehat{F}}_{X;L}\left({\xi }_0\right),$$(71)

This provides us with the conservative estimate l_c(ρ)/R ~ 10^–2, which is an order of magnitude smaller than the value estimated from the ACF (Eq. (67)) but closer to the estimate Eq. (70).

It is thus important to understand whether most of the anisotropy in the ACF originates from some integration artifacts and whether it causes or not an overestimation of the correlation lengths of ŋ or ρ.

Fig. 3 Column-density maps of the Polaris cloud. Left panel: without filter. Middle panel: with a high-pass filter filtering scales smaller than L/2. Right panel: with a low pass filter filtering scales larger than L/2. The low pass filter does not alter qualitatively the richness of structures found in the Polaris flare, while the high-pass filter shows a large-scale gradient that can be produced by an integration effect.

Fig. 4 Same as Fig. (3), but for the binary map Θ(log(∑/ 〈∑〉))where Θ is Heaviside’s step function. Regions where ∑ > 〈∑〉 appear darker than regions where ∑ < 〈∑〉.

Fig. 5 η-PDFs. Blue round and purple triangular symbols represent the PDFs of the unflltered and low pass filtered maps, respectively. The filter does not alter regions with η > 0 but increases the occurrence of regions with η < −1. Horizontal errorbars represent bin spacing.

6.4 Ergodic Estimate of the Observed PDF, Real Error Bars, and Reduced Integration Artifacts

As mentioned earlier, column-density PDFs serve as tracers of the statistics of the underlying density field. The various forms of these PDFs can be attributed to the various processes that are operating in MCs, from a fully lognormal distribution when purely turbulent motions dominate to a lognormal distribution with high density PLTs when gravitational effects become significant (Vazquez-Semadeni 1994; Passot & Vázquez-Semadeni 1998; Kainulainen et al. 2009; Schneider et al. 2013). This calls for a precise determination of the statistical uncertainty on the observed PDF, especially at high-density values.

The empirical PDF ${\widehat{f}}_X\left({\xi }_0\right)$ of stochastic field X (here X will be the column density ŋ) is deduced from histograms with some bin size ∆ξ. Error bars are usually estimated from Poisson statistics (using the number of points per bin) and can therefore be very small (Schneider et al. 2013). It is worth delving deeper into this issue. A histogram yields the following estimate: $${\widehat{f}}_{X,L}\left({\xi }_0\right)\Delta \xi \simeq {\widehat{F}}_{X;L}\left({\xi }_0+\Delta \xi \right)-{\widehat{F}}_{X;L}\left({\xi }_0\right),$$(72)

where ${\widehat{F}}_{X;L}$ is the empirical cumulative distribution function. Formally, this amounts to the ergodic estimate of the average of the following field, noted g_ξ0(y): $$g_{{\xi }_0}\left(y\right)=h_{{\xi }_0+\Delta \xi }\left(X\left(y\right)\right)-h_{{\xi }_0}\left(X\left(y\right)\right),$$(73)

where $$\text{where\hspace{0.17em}}{h}_{{\xi }_0}\left(X\left(y\right)\right)=\Theta \left({\xi }_0-X\left(y\right)\right),$$(74)

(see Sect. 2.1.2 and 4). Thus, proper statistical error levels must be calculated using the results of Sect. 2 and in general are not given by Poisson statistics.

In Appendix B.2, we study in detail ergodic estimates of average quantities. In general, the correlation length of g_ξ0 is a function of ξ₀itself. For Gaussian (or lognormal) distributions, an important result is that the confidence interval becomes quite large for values $\left|{\xi }_0-\mathbb{E}\left(X\right)\right|\ge \sigma \left(X\right)$, resulting in large errors if the sample size is too small. Thus, a reliable evaluation of the statistics of rare events (away from the average) requires very large sample sizes.

6.4.1 Reduced Integration Effects at High Density Contrasts

In this study, we focus on the column-density field X = ŋ and its PDF, noted p(ŋ). Using g_ŋ0(y) and its ACF for various values η₀, we are able to determine the appropriate statistical error bars and to get rid of some of the artifacts that are due to integration along the line of sight. In practice, we expect that such artifacts are not significant in high column-density regions (see Sect. 6.2). For example, anisotropy of the Polaris column density ACF in the x-direction is likely due to integration effects (see Sect. 6.3).

We expect, however, that the ACF of field g_ŋ0 for ŋ₀ > 0 is expected to show a reduced anisotropy at short scales. We thus obtain an empirical ACF of g_ŋ0 using Eq. (17). Figure 8 displays the estimated PDF of g_ŋ0 for the low pass filtered column-density map. At low column-density (η₀ = −1.06), a strong anisotropy is observed in the x-direction starting at x/L ≥2× 10^–2, as for the ACF of ŋ (see Fig. 7). For positive column density contrasts (η₀ > 0), this anisotropy is reduced and the ACFs are fairly isotropic at small scales in both the x and θ = π/4 directions, up to x/L ~ 10^–1 and r/L ~ 10^–1, respectively, where r denotes separation distance in the θ = π/4 direction. At larger separation distances, the data become quite noisy. This is consistent with the fact that the low path filtering procedure does not modify the PDF significantly in regions where ŋ > 0 (see Fig. 5).

This suggests that most of the ŋ ACF anisotropy in the x-direction at scales in the 10^–2−10^–1 range is due to integration effects. The peak of the correlation in the π/4-direction at high column-densities (ŋ₀ = 1.58) is probably due to the presence of the “Saxophone”-shaped filamentary structure that may be seen at the top of Fig. 3, which hosts most of the Polaris high density regions (Schneider et al. 2013).

6.4.2 Statistical Errorbars

Using the statistics of ɡ_ŋ0(y) has several advantages. One is that it reduces the impact of l.o.s. integration artifacts. In addition, it leads to proper error estimates for the empirical PDF.

Introducing some function of ŋ₀ noted ϕ(ŋ₀) which is expected to increase for increasing values of |ŋ₀|, the confidence interval above (1 − 1/m²) can be written as follows (see Bienayme-Tchebychev inequality, Eq. (21)): $$p\left({\eta }_0\right)\equiv f_n\left({\eta }_0\right)={\widehat{f}}_L\left({\eta }_0\right)\left(1\pm m{\left(\phi \left({\eta }_0\right)\right)}^{1/2}{\left(\frac{l_c\left(\eta \right)}{R}\right)}^{D/2}\right),$$(75)

with D = 2 and where $$\phi \left({\eta }_0\right)=\frac{\text{Var}\left(g_{{\eta }_0}\right)}{{\widehat{f}}_L{\left({\eta }_0\right)}^2{\left(\Delta \eta \right)}^2}\times {\left(\frac{l_c\left(g_{{\eta }_0}\right)}{l_c\left(\eta \right)}\right)}^2,$$(76) $$=\frac{1}{{\widehat{f}}_L\left({\eta }_0\right)\left(\Delta \eta \right)}\times {\left(\frac{l_c\left(g_{{\eta }_0}\right)}{l_c\left(\eta \right)}\right)}^2,$$(77)

because $\text{Var}\left(g_{{\eta }_0}\right)\simeq {\widehat{f}}_L\left({\eta }_0\right)\left(\text{Δ}\eta \right)$ and where l_c(ɡη₀)² ∝ Δη; so that the error bars on the PDF do not depend on the choice of bin size (for small Δη, see Appendix B.2).

From the empirical ACF $C_{g_{{\eta }_0}}$, one can then estimate the correlation length of ɡ_ŋ0 and thus ϕ(ŋ₀) for every ŋ₀. Unfortunately, this procedure is hampered by the fact that the ACF becomes increasingly noisy at high contrasts |η₀| > 1, due to sample sizes that are too small.

In principle, to determine ϕ(ŋ₀) and its variation one must calculate the complete integral that defines l_c(ɡ_ŋ0). This may be avoided as follows. The growth of ϕ(ŋ₀) may be obtained by looking at the short scale behavior of the ACFs of ɡ_ŋ0. In Fig. 8, it appears that the values of ${\widehat{C}}_{g_{{\eta }_0}}$ for positive column density contrasts (ŋ > 0) are isotropic and close to being ∝ |y|^–1/2 at short scales. We thus write that: $${\widehat{C}}_{g_{{\eta }_0}}\left(y\right){\left(f_{\eta }\left({\eta }_0\right)\Delta \eta \right)}^2\times c_{\sqrt{}}{\left|y/L\right|}^{-1/2}$$(78)

where $c_{\surd }$ is a constant of proportionality that depends on ŋ₀. Values of $c_{\surd }$ as a function of ŋ₀ are given in Fig. 9. We have only studied ɡ_ŋ0 for −0.7 ≤ ŋ₀ ≤ 1.58, because the ACFs are extremely noisy at high positive density contrasts (ŋ ≥ 1.58) due to poor sampling. At negative density contrasts (ŋ ≤ −0.7), where integration artifacts are the largest (see Sect. 6.4.1), the ACFs are no longer sufficiently isotropic and do not conform to a scaling in |y|^–1/2. Figure 9 shows that $c_{\surd }\left({\eta }_0\right)$ is an increasing function of |η₀| for large |η₀|, illustrating the fact that ϕ(ŋ₀) is expected to be large compared to f_ŋ(ŋ₀) for large contrasts |η₀| > 1.

Fig. 6 Reduced ACF function of $\eta \left({\widehat{C}}_{\eta }/\text{Var}\left(\eta \right)\right)$ for the Polaris flare. Left panel: without filter. Middle panel: with a high pass Alter filtering scales smaller than L/2. Right panel: with a low pass Alter Altering scales larger than L/2. Contours from black to purple to blue to light blue give the value of the reduced ACF at 0.5, e–1 ≃ 0.37, 0.1, −0.1.

Fig. 7 Reduced ACF of the low pass Altered map in three different directions. Blue line: x-direction (y = 0). Purple line: y-direction (x = 0). Red line: π/4 or x = y-direction. Green line: exponential pro Ale giving an estimate of the rate of decay (here λ/L 5 × 10–2). A strong anisotropy is present in the x direction at large scales (x/L ≥ 2 × 10–2).

6.4.3 Correlation Length of ɡ_η0 from Eq. (40)

To obtain an estimate of l_c(ɡ_ŋ0) and thus of ϕ(ŋ₀), we can use the results of Sect. 3. For each ŋ₀ we compute the field ɡ_η0(x, y) from the column density map ŋ(x, y). We then produce the two fields ${\displaystyle {\sum }_{g_{{\eta }_0},x}\left(y\right)}$ and ${\displaystyle {\sum }_{g_{{\eta }_0},y}\left(x\right)}$ obtained from the integration of ɡ_η0(x, y) along the x and y direction, respectively: $${\sum }_{g_{{\eta }_0,}x}\left(y\right)={\displaystyle \int g_{{\eta }_0}\left(x,y\right)\text{d}x,}$$(79) $${\sum }_{g_{{\eta }_0,}y}\left(y\right)={\displaystyle \int g_{{\eta }_0}\left(x,y\right)\text{d}y.}$$(80)