© ESO, 2013

1. Introduction

In this paper, we address a general statistical problem: the fit of density models to discrete data. Consider a non-negative function ρ(x | θ), where x ∈ R^d is a point of a Euclidean space of dimension d (in typical situations d = 1, 2, or 3), and θ is a vector of model parameters. Suppose then that we observe a realization of a non-homogeneous Poisson point process with density¹ρ(x | θ): that is, suppose that we observe a random set of points x_n ∈ R^d distributed such that the number of points N_A,B in two disjoint sets A and B of R^d are independent, and both N_A and N_B follow a Poisson distribution P(N) = e^− μμ^N/N ! with means $\begin{array}{ccc}{\mathit{\mu }}_{\mathit{A}}\mathrm{=}{\mathrm{\int }}_{\mathit{A}}\mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\theta }\mathrm{)}\hspace{0.17em}{\mathrm{d}}^{\mathit{d}}\mathit{x,} {\mathit{\mu }}_{\mathit{B}}\mathrm{=}{\mathrm{\int }}_{\mathit{B}}\mathit{\rho }\mathrm{(}\mathit{x}\mathrm{\left|}\mathit{\theta }\mathrm{\right)}\hspace{0.17em}{\mathrm{d}}^{\mathit{d}}\mathit{x}\mathit{.}& & \end{array}$(1)We seek a way to infer the parameters θ from the points { x_n }.

The framework just introduced can be used to describe many different situations. For example, we could use it to model the stellar distribution in globular clusters, using the King model for the stellar positions; or we could use it to model the number counts of galaxies. Sarazin (1980) considered the same framework for fitting the distribution of galaxies in galaxy clusters. In this paper, we use it to investigate the local Schmidt law of star formation in molecular clouds. In all cases we would be interested in inferring the relevant parameters of the density models used: center, richness, and concentration parameter for the King model; normalization, slope, and completeness limit for the number counts. In general, the relevant parameters θ can be constrained using a statistical frequentist approach, or a Bayesian approach. In both cases, we need to evaluate the likelihood function, i.e. the conditional probability to observe a given set of data points { x_n } given the value of the parameters θ.

2. The log-likelihood

Suppose that we observe a set of N points { x_n } that we know are a realization of a random Poisson point process with a given density ρ₁(x | θ), which we assume to be normalized to unity: $\mathrm{\int }{\mathit{\rho }}_{\mathrm{1}}\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\theta }\mathrm{)}\hspace{0.17em}{\mathrm{d}}^{\mathit{d}}\mathit{x}\mathrm{=}\mathrm{1.}$(2)Since this function, from a statistical point of view, is a probability density, we can immediately write the likelihood of the configuration { x_n } given the parameters θ as (Sarazin 1980) $\mathrm{ℒ}(\mathrm{\left\{}{\mathit{x}}_{\mathit{n}}\mathrm{\right\}}\mathrm{|}\mathit{\theta })\mathrm{=}\underset{\mathit{n}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}{\mathit{\rho }}_{\mathrm{1}}\mathrm{(}{\mathit{x}}_{\mathit{n}}\mathrm{)}\mathit{.}$(3)More generally, if we deal with a spatial density ρ(x | θ) that is not normalized to unity, we can always consider the function $\begin{array}{ccc}{\mathit{\rho }}_{\mathrm{1}}\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\theta }\mathrm{)}\mathrm{=}\mathit{\rho }\mathrm{(}\mathit{x}\mathrm{\left|}\mathit{\theta }\mathrm{\right)}\mathit{/}\mathit{\mu ,} \mathit{\mu }\mathrm{\equiv }\mathrm{\int }\mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\theta }\mathrm{)}\hspace{0.17em}{\mathrm{d}}^{\mathit{d}}\mathit{x,}& & \end{array}$(4)which by construction is normalized, and apply the likelihood (3)to ρ₁. However, this procedure prevents us from inferring any information on the normalizing constant μ; additionally, it has the drawback of forcing us to deal with normalized densities, which in many contexts is not a natural choice.

To clarify these points, consider a simple model of star formation in nearby molecular clouds, where we postulate that ${\mathrm{\Sigma }}_{\mathrm{YSO}}\mathrm{=}\mathit{\kappa }{\mathrm{\Sigma }}_{\mathrm{gas}}^{\mathit{\beta }}$(5)where Σ_YSO is the surface density of protostars (stars pc^-2) in a molecular cloud and Σ_gas is the mass surface density (M_⊙.pc^-2) of the molecular gas. This relation is similar to the version of the well-known Schmidt Law (Schmidt 1959) introduced by Kennicutt (1998) which relates the globally averaged surface densities of the star formation rate (Σ_SFR) and the gas (Σ_gas) in galaxies in a power-law fashion. Suppose that we know the locations { x_n } of N protostars and the gas column density Σ_gas(x) in a region of the sky, and that we intend to infer the parameters κ and β. In this case the likelihood discussed above does not provide any information on κ, and this is unfortunate since this parameter is directly linked to a critical piece of information, the star formation rate of the molecular cloud.

In order to solve this problem, we note that the probability density to observe N points { x_n } from ρ(x | θ) can be factorized as the probability to observe N points (a Poisson distribution with mean μ), and the probability that these N points be located at the positions { x_n } (as provided by Eq. (3)): $\mathrm{ℒ}(\mathrm{\left[}{\mathit{x}}_{\mathit{n}}\mathrm{\right]}\mathrm{|}\mathit{\theta })\mathrm{=}{\mathrm{e}}^{\mathrm{-}\mathit{\mu }}\frac{{\mathit{\mu }}^{\mathit{N}}}{\mathit{N}\mathrm{!}}\underset{\mathit{n}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}{\mathit{\rho }}_{\mathrm{1}}({\mathit{x}}_{\mathit{n}}\mathrm{|}\mathit{\theta })\mathrm{=}\frac{\mathrm{1}}{\mathit{N}\mathrm{!}}{\mathrm{e}}^{\mathrm{-}\mathrm{\int }\mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right)}\hspace{0.17em}{\mathrm{d}}^{\mathit{d}}\mathit{x}}\underset{\mathit{n}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}\mathit{\rho }({\mathit{x}}_{\mathit{n}}\mathrm{|}\mathit{\theta })\mathit{.}$(6)In this expression the data points [x_n] are taken to be ordered; if instead we consider the unordered set { x_n }, we just have to multiply this result by the number of permutations of N elements, i.e. N !. With this choice, the log-likelihood takes the simple form $\ln \mathrm{ℒ}(\mathrm{\left\{}{\mathit{x}}_{\mathit{n}}\mathrm{\right\}}\mathrm{|}\mathit{\theta })\mathrm{=}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}\ln \mathit{\rho }\mathrm{(}{\mathit{x}}_{\mathit{n}}\mathrm{\left|}\mathit{\theta }\mathrm{\right)}\mathrm{-}\mathrm{\int }\mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\theta }\mathrm{)}\hspace{0.17em}{\mathrm{d}}^{\mathit{d}}\mathit{x}\mathit{.}$(7)

3. Parameters inference

Using Bayes’ theorem $\mathit{P}(\mathit{\theta }\mathrm{\left|}\mathrm{\right\{}{\mathit{x}}_{\mathit{n}}\mathrm{\}})\mathrm{=}\frac{\mathrm{ℒ}(\mathrm{\left\{}{\mathit{x}}_{\mathit{n}}\mathrm{\right\}}\mathrm{\left|}\mathit{\theta })\mathit{p}\mathrm{\right(}\mathit{\theta }\mathrm{)}}{\int \mathrm{ℒ}(\mathrm{\left\{}{\mathit{x}}_{\mathit{n}}\mathrm{\right\}}\mathrm{\left|}{\mathit{\theta }}^{\mathrm{\prime }})\mathit{p}\mathrm{\right(}{\mathit{\theta }}^{\mathrm{\prime }}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{\theta }\prime }\mathit{,}$(8)we can infer the posterior probability distribution of the parameters, P(θ | { x_n }) given the prior distribution p(θ). In our specific case, if we write ρ(x | θ) = μρ₁(x | θ), then the likelihood factorizes into the product of a term that depends only on μ and a term that depends on the other density parameters θ. As a result, if the prior p(μ,θ) = p(μ)p(θ) also factorizes, then the posterior distribution factorizes and μ is independent of θ (which implies that the two are also uncorrelated). In particular, if we use the uninformative improper prior p(μ) = 1/μ for μ > 0, we find that the posterior for μ is a simple gamma distribution $\mathit{P}\mathrm{\left(}\mathit{\mu }\mathrm{\right|}\mathrm{\left\{}{\mathit{x}}_{\mathit{n}}\mathrm{\right\}}\mathrm{)}\mathrm{=}\frac{{\mathrm{e}}^{\mathrm{-}\mathit{\mu }}{\mathit{\mu }}^{\mathit{N}\mathrm{-}\mathrm{1}}}{\mathrm{\Gamma }\mathrm{\left(}\mathit{N}\mathrm{\right)}}\mathrm{·}$(9)Below, as often the case in astronomy, we deal with an unnormalized density ρ, and therefore we do not expect the various parameters θ to be independent.

4. Frequentist description

Alternatively, we can use to (log-)likelihood to obtain a maximum-likelihood point estimate of the parameters θ (Fisher 1922), i.e. find .

Bounds on the errors on the parameters can be evaluated from the Fisher information matrix (see, e.g., Feigelson & Babu 2012), that we recall is defined as ${\mathit{I}}_{\mathit{ij}}\mathrm{=}\mathrm{E}\left[\frac{\mathit{\partial }\ln \mathrm{ℒ}}{\mathit{\partial }{\mathit{\theta }}_{\mathrm{i}}}\frac{\mathit{\partial }\ln \mathrm{ℒ}}{\mathit{\partial }{\mathit{\theta }}_{\mathrm{j}}}\right]\mathrm{=}\mathrm{-}\mathrm{E}\left[\frac{{\mathit{\partial }}^{\mathrm{2}}\ln \mathrm{ℒ}}{\mathit{\partial }{\mathit{\theta }}_{\mathrm{i}}\hspace{0.17em}\mathit{\partial }{\mathit{\theta }}_{\mathrm{j}}}\right]\mathit{,}$(10)where the symbol E denotes the operation of ensemble average. The Fisher information matrix is related to the minimum covariance matrix that can be attained by an unbiased estimator, as provided by the Cramér-Rao bound: $\mathrm{Cov}\mathrm{\left(}\mathit{\theta ̂}\mathrm{\right)}\mathrm{\ge }{\mathrm{I}}^{-1}\mathit{.}$(11)Since the maximum-likelihood estimator is asymptotically efficient (i.e. it attains the Cramér-Rao bound when the sample size tends to infinity) and the resulting errors on tend to a multivariate Gaussian distribution, it is interesting to obtain an analytic result for the information matrix. In Appendix A.2, we show that $\begin{array}{ccc}{\mathit{I}}_{\mathit{ij}}& \mathrm{=}& \int \frac{\mathrm{1}}{\mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\theta }\mathrm{)}}\frac{\mathit{\partial \rho }\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\theta }\mathrm{)}}{\mathit{\partial }{\mathit{\theta }}_{\mathit{i}}}\frac{\mathit{\partial \rho }\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\theta }\mathrm{)}}{\mathit{\partial }{\mathit{\theta }}_{\mathit{j}}}\hspace{0.17em}{\mathrm{d}}^{\mathit{d}}\mathit{x}\\ & \mathrm{=}& \mathrm{\int }\mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\theta }\mathrm{)}\frac{\mathit{\partial }\ln \mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\theta }\mathrm{)}}{\mathit{\partial }{\mathit{\theta }}_{\mathit{i}}}\frac{\mathit{\partial }\ln \mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\theta }\mathrm{)}}{\mathit{\partial }{\mathit{\theta }}_{\mathit{j}}}\hspace{0.17em}{\mathrm{d}}^{\mathit{d}}\mathit{x}\mathit{.}\end{array}$(12)Finally, we can evaluate the goodness of the fit obtained by comparing the maximum likelihood value with what we are expected to obtain. The expected (average) likelihood can be estimated analytically by performing an ensemble average over the points { x_n }. As show in Appendix A.2, the average value of the log-likelihood calculated at the true value of the parameters is $\mathrm{E}\mathrm{\left[}\ln \mathrm{ℒ}\mathrm{\right]}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathrm{\int }\mathit{\rho }\mathrm{\left(}\mathrm{x}\mathrm{\right|}\mathit{\theta }\mathrm{\left)}[\ln \mathit{\rho }\mathrm{\right(}\mathrm{x}\mathrm{\left|}\mathit{\theta }\mathrm{\right)}\mathrm{-}\mathrm{1}]\hspace{0.17em}{\mathrm{d}}^{\mathrm{d}}\mathrm{x}\mathit{.}$(13)For our purposes it is more relevant to obtain the expectation value of the likelihood at its maximum value : this value can then be compared with the observed maximum of the log-likelihood to provide a goodness-of-fit proxy. The result obtained is $\mathrm{E}\mathrm{\left[}\ln \mathrm{ℒ}\mathrm{\right]}\mathrm{\left(}\mathit{\theta ̂}\mathrm{\right)}\mathrm{=}\mathrm{E}\mathrm{\left[}\ln \mathrm{ℒ}\mathrm{\right]}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{+}\frac{\mathrm{J}}{\mathrm{2}}\mathit{,}$(14)where J is the number of the free parameters in the density, i.e. the dimensionality of the space of θ.

Finally, the variance of the log-likelihood is provided by $\mathrm{Var}\mathrm{\left[}\ln \mathrm{ℒ}\mathrm{\right]}\mathrm{\left(}\mathit{\theta ̂}\mathrm{\right)}\mathrm{=}\mathrm{\int }\mathit{\rho }\mathrm{\left(}\mathrm{x}\mathrm{\right|}\mathit{\theta }\mathrm{\left)}{\ln }^{\mathrm{2}}\mathit{\rho }\mathrm{\right(}\mathrm{x}\mathrm{\left|}\mathit{\theta }\mathrm{\right)}\hspace{0.17em}{\mathrm{d}}^{\mathrm{d}}\mathrm{x}\mathit{.}$(15)This expression, together with Eq. (14)is useful to define acceptance boundaries for the likelihood value and therefore for model test.

Fig. 1 Maximum likelihood best-fit parameters for a set of 100 simulations. The dotted lines show the true, original values used: κ = 1.5 star pc-2, β = 1.8, A0 = 0.3 mag, and σ = 0.5 pc. Note how the points, representing the individual best-fits, cluster around the true values. The ellipses, showing the expected 3-σ errors as deduced from the Fischer information matrix, are in excellent agreement with the observed distribution of points.

5. Application: deriving the local Schmidt scaling law for a giant molecular cloud

In this section we apply the framework developed above to one relevant astrophysical problem, determining the local star formation scaling law for a nearby Giant Molecular Cloud, which in its simplest version can be written as in Eq. (5).

In nearby molecular clouds Σ_YSO is a directly measured quantity and can be related to Σ_SFR via Σ_SFR = Σ_YSO × ⟨ m ⟩ /τ, where ⟨ m ⟩ is the mean mass of a protostar and τ the age of the protostellar (class I) population in the cloud. We note here that we are interested in the local star formation scaling law within a GMC and this is physically different from the Kennicutt-Schmidt law which relates global quantities averaged over entire galaxies (e.g., Kennicutt 1998) or subregions of galaxies (e.g., Bigiel et al. 2008) which in either case do not resolve the individual molecular clouds. Indeed, the indexes (β) of the local and global relations are not necessarily physically related, nor expected to be the same, when comparing individual molecular clouds and normal galaxies (e.g., Calzetti et al. 2012).

Since we use infrared extinction measurements to trace the gas column density in the clouds we write Eq. (5)as ${\mathrm{\Sigma }}_{\mathrm{YSO}}\mathrm{\propto }{\mathit{A}}_{\mathit{K}}^{\mathit{\beta }}$ and note that Σ_gas ∝ A_K for a constant gas-to-dust ratio. Additionally we allow for the possibility of a star-formation threshold, i.e., a lower limit of the surface density for gas, or equivalently dust, below which no star formation takes place and no protostars are produced in situ (e.g., Lada et al. 2010, Heiderman et al. 2010). Finally, since it is physically reasonable to assume that star formation is not an instantaneous process, we allow for some diffusion of protostars from their birth sites. We model this diffusion process by smoothing the initial protostellar surface density, ${\mathrm{\Sigma }}_{\mathrm{YSO}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$, by a Gaussian spatial kernel. In summary, we have ${\mathrm{\Sigma }}_{\mathrm{YSO}}\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\kappa ,\beta ,}{\mathit{A}}_{\mathrm{0}}\mathit{,\sigma }\mathrm{)}\mathrm{=}\mathrm{\int }\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{e}}^{\mathrm{|}\mathit{x}\mathrm{-}{\mathit{x}}^{\mathrm{\prime }}{\mathrm{|}}^{\mathrm{2}}\mathit{/}\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}{\mathrm{\Sigma }}_{\mathrm{YSO}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{(}{\mathit{x}}^{\mathrm{\prime }}\mathrm{\left|}\mathit{\kappa ,\beta ,}{\mathit{A}}_{\mathrm{0}}\mathrm{\right)}{\mathrm{d}}^{\mathrm{2}}{\mathit{x}}^{\mathrm{\prime }}\mathit{,}$(16)where ${\mathrm{\Sigma }}_{\mathrm{YSO}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\kappa ,\beta ,}{\mathit{A}}_{\mathrm{0}}\mathrm{)}\mathrm{=}\mathit{\kappa }\mathrm{H}({\mathit{A}}_{\mathit{K}}\mathrm{(}\mathit{x}\mathrm{)}\mathrm{-}{\mathit{A}}_{\mathrm{0}}){\left(\frac{{\mathit{A}}_{\mathit{K}}}{\mathrm{1}\mathrm{mag}}\right)}^{\mathit{\beta }}\mathrm{(}\mathit{x}\mathrm{)}\mathit{.}$(17)In this equation H is the Heaviside step function (defined to be zero for a negative argument and one for positive one). The constants involved are the normalization κ (taken to be measured in units of star pc^-2 mag^− β, the star-formation threshold A₀ (in units of K-band extinction), the dimensionless exponent β, and the diffusion coefficient σ (measured in pc).

The data at our disposal will be the extinction map of a molecular cloud, i.e. A_K(x), and a catalog with the positions of protostars { x_n }. With these data we can fit the four parameters θ = { κ,β,A₀,σ } using both the approaches described in Sects. 3 and 4.

5.1. Simulations

Before using the techniques described in this paper on real data, we validated them on simulated data. The simulations were carried out by taking the extinction map of the Orion molecular cloud² from Lombardi et al. (2011), and by randomly generating protostars according to the law (16). Specifically, we fixed the parameters β, A₀, and σ to test values, and we set κ such that the expected number of protostars in the field was 300. For each simulation, we then drew the number of protostars from a Poisson distribution with the appropriate average, and we distributed the initial positions of the stars in the field following the density ${\mathrm{\Sigma }}_{\mathrm{YSO}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}$ of Eq. (17), i.e. without any diffusion represented by the σ parameters. Finally, we changed the initial positions of the protostars by drawing random offsets from a two dimensional normal distribution with variance σ², and by moving each protostar position according to the drawn offsets.

We then tried to recover the “unknown” parameters using both a maximum likelihood approach and a Bayesian inference. The input data provided to the code were the positions of the protostars, as generated using the procedure described above, and the extinction map of the cloud. Figure 1 reports the results obtained in a typical set of simulation with true parameters κ = 1.5 star pc^-2 mag^− β, β = 1.8, A₀ = 0.3 mag, and σ = 0.5 pc. To produce the plots, we performed 100 independent simulations, drawing each time a different set of protostars, and we report in the plot the locations of the best-fit estimates obtained from the maximization of the likelihood. We also plot the expected 3-σ error ellipse, as derived from the Fisher information matrix. As expected, and as shown by these plots (and analogous ones produced during our tests), the maximum likelihood estimate does not suffer any evident bias and is able to constrain all parameters in a very efficient way. Note that for the set of simulations shown in Fig. 1 we used on average 300 protostars, a number comparable with the number of class I objects known in Orion A (see below). During the tests we also verified that the value of the likelihood function at its maximum was compatible with the expectations of Eqs. (14)and (15).

Fig. 2 Differential distribution of ΣYSO as a function of extinction on a simulation carried out with input parameters κ = 1.5 star pc-2 mag− β, β = 1.8, A0 = 0.3 mag, and σ = 0.5 pc. The best-fit power law, shown as a line in this plot, is κ = (1.23 ± 0.12) star pc-2 mag− β and β = 2.71 ± 0.08.

Fig. 3 Posterior probability for the four density parameters obtained for the Orion A molecular cloud.

We also tested in a similar way the Bayesian approach. For these tests, we adopted both uniform and uninformative Jeffreys priors for the parameters, excluding for all of them negative values. In general, during our tests we found that original, true values of the parameters were always within the Bayesian 95% credible intervals.

As a comparison, Fig. 2 shows the projected density of protostars as a function of the projected density of the dust. To make this figure, we considered increasing contour levels of extinction, and we estimated the density of protostars within two consecutive levels of extinction by computing the ratio between the number of protostars and the area enclosed within the two contours; errors were estimated from simple Poisson statistics. The figure also shows the best-fit power law obtained from these data, which is not in agreement with the original input: the measured parameters where κ = (1.23 ± 0.12) star pc^-2 mag^− β and β = 2.71 ± 0.08, compared with the original input parameters of κ = 1.5 stars pc^-2 mag^− β and β = 1.8. This simple test shows that in presence of a threshold and/or of a diffusion of the protostars from their original loci of formation, the standard fit of the Σ_gas–Σ_YSO plot leads to unreliable results. Note that the poor performance of the simple fit of histograms is a combination of several factors: the use of a simple least-squares minimization algorithm does not take into account the Poisson statistics of the bins (see Laurence & Chromy 2010); the lack of modeling of the threshold and of the diffusion; the lost of information when binning the data. A partial removal of these factors occasionally makes the algorithm less biased (but only marginally so); in some cases (particularly when adding the parameter A₀ to the fit) the bias was actually larger, presumably as a result of a poor minimization of the least square algorithm. In general, we find that the fitting algorithm based on binned data is very unstable, leading to large differences in the results depending on the particular configuration used (i.e., on the specific location of the protostars { x_n } and on the particular choice of the bin sizes). Note instead that the use of a Poisson statistics with infinitesimal bins produces the same likelihood of Eq. (6)(see Appendix A.1), and therefore guarantees unbiased and robust results.

Additionally, we also checked the results of the Bayesian analysis when other observational effects were present in the data. Specifically, after generating the stars according to the local Schmidt law, we convolved the map with a Gaussian beam with FWHM = 6 arcmin to simulate the effect of a unresolved structures due to finite resolution. Moreover, we added artificial noise that takes into account the proper error propagation in the extinction map (see Lombardi & Alves 2001). Simulations showed that the Bayesian technique can cope well with these effects, and the results obtained are largely unaffected.

In summary, the simulations completely validated the approach described in this paper. We therefore applied our method to the best studied star-forming molecular cloud, the Orion complex.

5.2. Application to Orion A

In order to derive the local Schmidt scaling law in Orion-A, we used our 2MASS/Nicest extinction map (Lombardi et al. 2011). We preferred the Nicest (Lombardi 2009) algorithm over the Nicer (Lombardi & Alves 2001) one because the former produces maps that are less affected by systematic biases in the central regions of molecular clouds, where presumably most protostars are formed. Additionally, we used a catalog of protostars obtained from a survey carried out in the region with Spitzer Space Telescope (Megeath et al. 2012). From this catalog we selected only Class I sources (specifically, objects classified as “probable protostars”); moreover, we excluded 56 sources found by cross-correlation with a catalog of 624 foreground objects (identified as unreddened sources visible in front of the molecular cloud, see Alves & Bouy 2012). As a result, we were left with 329 objects enclosed within the polygonal area in Orion-A where the observations have been carried out (see Megeath et al. 2012 for details on the area selection).

We performed frequentist and a Bayesian analyses of these data using the techniques described in this paper. For the Bayesian analysis we fitted the model (16)with flat priors over all parameters (taken however to be positive) using the likelihood of Eq. (7). We explored the resulting posterior probability distribution with Markov Chain Monte Carlo, using a simple Metropolis-Hastings sampler. The resulting posterior probability, shown in Fig. 3, presents the following relevant results (all bounds are referred to a 95% confidence level interval):

• The exponent β is exquisitely close to the value 2: the result obtained is indeed β = 2.03 ± 0.15.

• We measure a star-formation coefficient κ = (1.65 ± 0.19) star pc^-2 mag^-2.

• The data show that there is no threshold for star formation: A₀ < 0.09 mag.

• Class I protostars seem to have undergone no detectable diffusion from their original positions beyond the resolution of our extinction maps.

The frequentist analysis provides consistent results. Here we just mention the fact that the best-fit log-likelihood is ; as a reference, the value predicted from Eq. (14)is , with an expected tolerance computed from the square root of Eq. (15)of ~52. We deduce therefore that the model proposed is in agreement with the observations.

Fig. 4 Extinction map of Orion A (taken from Lombardi et al. 2011), together with the location of the Class I protostars (from Megeath et al. 2012). The polygonal line shows the region where the protostar observations have been carried out.

5.3. Discussion

Using the analysis described above we can now write the fundamental star formation scaling relation for the Orion A molecular cloud: ${\mathrm{\Sigma }}_{\mathrm{YSO}}\mathrm{=}\mathrm{(}\mathrm{1.65}\mathrm{\pm }\mathrm{0.19}\mathrm{)}{\left(\frac{{\mathit{A}}_{\mathit{K}}}{\mathrm{mag}}\right)}^{\mathrm{2.03}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{0.15}}\hspace{0.17em}\mathrm{star}\mathrm{p}{\mathrm{c}}^{-2}\mathit{.}$(18)This result can be considered a complete description of the local Schmidt law characterizing this cloud because our analysis has enabled the derivation of robust values for both the power-law index, β, and the coefficient, κ, the constant of proportionality. While β determines how Σ_YSO varies with column density of the cloud, κ sets the overall scale or magnitude of the relation.

It is of interest to compare our results to those of a few recently published studies. Gutermuth et al. (2011) derived a value for β of 1.8 ± 0.01 for the Orion cloud from a least-squares fit to the Σ_YSO–Σ_gas relation, using the nearest neighbor method to derive Σ_YSO for each star and infrared extinction measurements to derive the corresponding Σ_gas. This result is in reasonable agreement with our findings, particularly given the large scatter in the Gutermuth et al. (2011) data. Steeper βs (~4), however, have been derived for other clouds using different methodologies (Heiderman et al. 2010; Harvey et al. 2013). Whether this represents evidence for variations in β between clouds is far from clear at this time. More robust statistics and application of a more consistent methodology to studies of other clouds would be needed to make a better comparison to our results for Orion and a better assessment of possible cloud-to-cloud variations in the star formation scaling law. We started already a follow-up analysis in this direction, and the results of the application of the statistical technique described in this paper to a set of clouds are presented in Lada et al. (2013), to which we refer the reader for a more the astrophysical conclusions of the analysis.

5.4. Concluding remarks

In this paper, we described the use of a likelihood analysis to derive a parametric density model that best describes the discrete spatial data. The likelihood can be used both in a frequentist and Bayesian analysis, and appears to be validated by its application to simulated data in which input parameters were accurately recovered. We applied the method to model the observed distribution of Class I protostars in the Orion A molecular cloud and derive the star formation scaling law for that cloud. Specifically we find: Σ_YSO = (1.65 ± 0.19) (A_K/mag)^{2.03 ± 0.15} stars pc^-2. Moreover, we found no evidence for an extinction threshold for star formation in the cloud and no evidence for any significant diffusion of the protostars from their birth sites.

References

Appendix A: Analytical derivations

Appendix A.1: Alternative derivation of Eq. (7)

In this section, we provide an alternative derivation of Eq. (7), that highlights the relationship between the likelihood and the fit of histograms. Suppose we that we know completely the density ρ(x) ≡ ρ(x | θ) responsible for a Poisson point process. One possibility to evaluate the probability to observe a given configuration { x_n } of points is to make a partition of the domain of ρ(x) into a small regions or bins of equal area a, chosen such that a ≪ 1/ρ(x) for any point x ∈ R^d. In this way, since the average number of stars for each bin is much less then unity, then each bin will contain at most one point. In this limit, each bin will have with probability P₀(x) = 1 − aρ(x) no stars, and with probability P₁(x) = aρ(x) one star. The joint probability for the whole configuration will be therefore the product of the individual probabilities for the individual bins, i.e. $\mathit{P}(\mathrm{\left\{}{\mathit{x}}_{\mathit{n}}\mathrm{\right\}}\mathrm{|}\mathit{\theta })\mathrm{=}\underset{\mathit{n}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}[\mathit{a\rho }\mathrm{(}{\mathit{x}}_{\mathit{n}}\mathrm{\left|}\mathit{\theta }\mathrm{\right)}]\underset{\mathit{m}}{\mathrm{􏽙}}[\mathrm{1}\mathrm{-}\mathit{a\rho }\mathrm{\left(}{\mathit{x}}_{\mathit{m}}\mathrm{\right|}\mathit{\theta }\mathrm{)}]\mathit{,}$(A.1)where the first product runs over the N bins with a single data point, and the second one over the other bins, i.e. all bins without any data point (we recall that since bins are taken to be small, no bin has more than a datapoint). This probability is proportional to a^N, and therefore as we take the limit a → 0 it is interesting to consider ℒ = P/a^N: that operation corresponds to use a probability density. If we switch to logarithms, the first product over n becomes a sum of lnρ at the locations of the points, while the second product, in the limit of a → 0, becomes an integral over the entire space (because the bins with points make a negligible contribution, see Fig. A.1): $\ln \mathrm{ℒ}(\mathrm{\left\{}{\mathit{x}}_{\mathit{n}}\mathrm{\right\}}\mathrm{|}\mathit{\theta })\mathrm{=}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}\ln \mathit{\rho }\mathrm{(}{\mathit{x}}_{\mathit{n}}\mathrm{\left|}\mathit{\theta }\mathrm{\right)}\mathrm{-}\mathrm{\int }\mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right|}\mathit{\theta }\mathrm{)}\hspace{0.17em}{\mathrm{d}}^{\mathit{d}}\mathit{x}\mathit{.}$(A.2)

Appendix A.2: Properties of the log-likelihood

In order to prove Eq. (13), we note that the average has to be carried out over all possible configurations of data points { x_n }, and therefore involves only the first term, i.e. the summation, of this equation; the integral is a constant term and the average is trivially itself. Let us call ℓ the first term of Eq. (13), and μ the second, so that lnℒ = ℓ + μ. By using the same reasoning of Eq. (6), we can write the average of ℓ by considering separately the probability to have exactly N points in the field, which follows a Poisson probability, and the probability for the distribution of each of these points. In summary, $\begin{array}{ccc}\mathrm{E}\mathrm{\left[}\mathit{\ell }\mathrm{\right]}& \mathrm{=}& \sum _{\mathit{N}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}{\mathrm{e}}^{\mathrm{-}\mathit{\mu }}\frac{{\mathit{\mu }}^{\mathit{N}}}{\mathit{N}\mathrm{!}}\mathrm{\int }\mathrm{d}{\mathit{x}}_{\mathrm{1}}\hspace{0.17em}\frac{\mathit{\rho }\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathrm{\right)}}{\mathit{\mu }}\mathrm{·}\mathrm{·}\mathrm{·}\mathrm{\int }\mathrm{d}{\mathit{x}}_{\mathit{N}}\hspace{0.17em}\frac{\mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\mathit{\mu }}\mathit{\ell }\\ & \mathrm{=}& \sum _{\mathit{N}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\frac{{\mathrm{e}}^{\mathrm{-}\mathit{\mu }}}{\mathit{N}\mathrm{!}}\mathrm{\int }\mathrm{d}{\mathit{x}}_{\mathrm{1}}\hspace{0.17em}\mathit{\rho }\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathrm{\right)}\mathrm{·}\mathrm{·}\mathrm{·}\mathrm{\int }\mathrm{d}{\mathit{x}}_{\mathit{N}}\hspace{0.17em}\mathit{\rho }\mathrm{\left(}{\mathit{x}}_{\mathit{N}}\mathrm{\right)}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}\ln \mathit{\rho }\mathrm{\left(}{\mathit{x}}_{\mathit{n}}\mathrm{\right)}\\ & \mathrm{=}& \sum _{\mathit{N}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}{\mathrm{e}}^{\mathrm{-}\mathit{\mu }}\frac{{\mathit{\mu }}^{\mathit{N}\mathrm{-}\mathrm{1}}}{\mathit{N}\mathrm{!}}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{\int }\mathrm{d}{\mathit{x}}_{\mathit{n}}\hspace{0.17em}\mathit{\rho }\mathrm{\left(}{\mathit{x}}_{\mathit{n}}\mathrm{\right)}\ln \mathit{\rho }\mathrm{\left(}{\mathit{x}}_{\mathit{n}}\mathrm{\right)}\\ & \mathrm{=}& \sum _{\mathit{N}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}{\mathrm{e}}^{\mathrm{-}\mathit{\mu }}\frac{{\mathit{\mu }}^{\mathit{N}\mathrm{-}\mathrm{1}}}{\mathrm{(}\mathit{N}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{!}}\mathrm{\int }\mathrm{d}\mathit{x}\hspace{0.17em}\mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right)}\ln \mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right)}\\ & \mathrm{=}& \mathrm{\int }\mathrm{d}\mathit{x}\hspace{0.17em}\mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right)}\ln \mathit{\rho }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{.}\end{array}$(A.3)Finally, putting back the term μ of lnℒ we obtain the desired result of Eq. (13).

Fig. A.1 Derivation of the likelihood of Eq. (7). Note that when we reduce the grid size from the black grid to the grey one, each grid cell eventually contains at most one point. Note also how the total area of the grid cells with points becomes negligible when the grid size is small.

The variance of lnℒ is identical to the variance of ℓ, which can be evaluated from E [ℓ²]. Following a calculation similar to Eq. (3)and using the independence of x_n from x_m for n ≠ m one obtains $\begin{array}{ccc}\mathrm{E}\mathrm{\left[}{\mathit{\ell }}^{\mathrm{2}}\mathrm{\right]}\mathrm{=}{\left[\mathrm{\int }\mathit{\rho }\mathrm{\left(}\mathrm{x}\mathrm{\right)}\ln \mathit{\rho }\mathrm{\left(}\mathrm{x}\mathrm{\right)}\hspace{0.17em}{\mathrm{d}}^{\mathrm{d}}\mathrm{x}\right]}^{\mathrm{2}}\mathrm{+}\mathrm{\int }\mathit{\rho }\mathrm{\left(}\mathrm{x}\mathrm{\right)}{\ln }^{\mathrm{2}}\mathit{\rho }\mathrm{\left(}\mathrm{x}\mathrm{\right)}\hspace{0.17em}{\mathrm{d}}^{\mathrm{d}}\mathrm{x}\mathit{.}& & \end{array}$(A.4)From this equation one immediately derives Var(ℒ) = Var(ℓ) = E(ℓ²) − E²(ℓ) as given in Eq. (15).

In a similar way one can derive Eq. (12). We start from the second equality of Eq. (10), and we note that the partial derivatives of lnℒ can be written as $\begin{array}{ccc}\mathrm{-}\frac{{\mathit{\partial }}^{\mathrm{2}}\ln \mathrm{ℒ}}{\mathit{\partial }{\mathit{\theta }}_{\mathit{i}}\hspace{0.17em}\mathit{\partial }{\mathit{\theta }}_{\mathit{j}}}& \mathrm{=}& \sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}\left[\frac{\mathrm{1}}{{\mathit{\rho }}^{\mathrm{2}}\mathrm{\left(}{\mathit{x}}_{\mathit{n}}\mathrm{\right)}}\frac{\mathit{\partial \rho }\mathrm{\left(}{\mathit{x}}_{\mathit{n}}\mathrm{\right)}}{\mathit{\partial }{\mathit{\theta }}_{\mathit{i}}}\frac{\mathit{\partial \rho }\mathrm{\left(}{\mathit{x}}_{\mathit{n}}\mathrm{\right)}}{\mathit{\partial }{\mathit{\theta }}_{\mathit{j}}}\mathrm{-}\frac{\mathrm{1}}{\mathit{\rho }\mathrm{\left(}{\mathit{x}}_{\mathit{n}}\mathrm{\right)}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\rho }\mathrm{\left(}{\mathit{x}}_{\mathit{n}}\mathrm{\right)}}{\mathit{\partial }{\mathit{\theta }}_{\mathit{i}}\hspace{0.17em}\mathit{\partial }{\mathit{\theta }}_{\mathit{j}}}\right]\\ & & \mathrm{-}\mathrm{\int }\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\rho }\mathrm{\left(}{\mathit{x}}_{\mathit{n}}\mathrm{\right)}}{\mathit{\partial }{\mathit{\theta }}_{\mathit{i}}\hspace{0.17em}\mathit{\partial }{\mathit{\theta }}_{\mathit{j}}}\hspace{0.17em}{\mathrm{d}}^{\mathit{d}}\mathit{x}\mathit{.}\end{array}$(A.5)When taking the average over the point positions { x_n } of this expression, the summation of the first line becomes an integral over ρ(x) d^dx, and as a result the last term of the first line cancels with the integral of the second line. We are therefore left out with Eq. (12).

All Figures

Fig. 1 Maximum likelihood best-fit parameters for a set of 100 simulations. The dotted lines show the true, original values used: κ = 1.5 star pc-2, β = 1.8, A0 = 0.3 mag, and σ = 0.5 pc. Note how the points, representing the individual best-fits, cluster around the true values. The ellipses, showing the expected 3-σ errors as deduced from the Fischer information matrix, are in excellent agreement with the observed distribution of points.

In the text

	Fig. 2 Differential distribution of ΣYSO as a function of extinction on a simulation carried out with input parameters κ = 1.5 star pc-2 mag− β, β = 1.8, A0 = 0.3 mag, and σ = 0.5 pc. The best-fit power law, shown as a line in this plot, is κ = (1.23 ± 0.12) star pc-2 mag− β and β = 2.71 ± 0.08.
In the text

	Fig. 3 Posterior probability for the four density parameters obtained for the Orion A molecular cloud.
In the text

	Fig. 4 Extinction map of Orion A (taken from Lombardi et al. 2011), together with the location of the Class I protostars (from Megeath et al. 2012). The polygonal line shows the region where the protostar observations have been carried out.
In the text

	Fig. A.1 Derivation of the likelihood of Eq. (7). Note that when we reduce the grid size from the black grid to the grey one, each grid cell eventually contains at most one point. Note also how the total area of the grid cells with points becomes negligible when the grid size is small.
In the text