A new method to unveil embedded stellar clusters
Marco Lombardi Charles J. Lada Joao Alves
astro-ph.IM astro-ph.GA



 marco.lombardi@unimi.it 
  University of Milan, Department of Physics, via Celoria 16, I-20133
  Milan, Italy Harvard-Smithsonian Center for Astrophysics, Mail
  Stop 72, 60 Garden Street, Cambridge, MA 02138 University of
  Vienna, Türkenschanzstrasse 17, 1180 Vienna, Austria



  In this paper we present a novel method to identify and characterize
  stellar clusters deeply embedded in a dark molecular cloud.  The
  method is based on measuring stellar surface density in wide-field
  infrared images using star counting techniques. It takes advantage
  of the differing H-band luminosity functions (HLFs) of field stars
  and young stellar populations and is able to statistically associate
  each star in an image as a member of either the background stellar
  population or a young stellar population projected on or near the
  cloud.  Moreover, the technique corrects for the effects of
  differential extinction toward each individual star. We have tested
  this method against simulations as well as observations. In
  particular, we have applied the method to 2MASS point sources
  observed in the Orion A and B complexes, and the results obtained
  compare very well with those obtained from deep Spitzer and
  Chandra observations where presence of infrared excess or
  X-ray emission directly determines membership status for every star.
  Additionally, our method also identifies unobscured clusters and a
  low resolution version of the Orion stellar surface density map
  shows clearly the relatively unobscured and diffuse OB 1a and 1b
  sub-groups and provides useful insights on their spatial
  distribution.  
A new method to unveil embedded stellar clusters
    Marco Lombardi1, Charles J. Lada2, and João
  Alves3
    Received ***date***; Accepted ***date***
========================================================













§ INTRODUCTION



Embedded star clusters offer one of the best opportunities to
understand star formation.  Within these astrophysical laboratories
hundreds of stars are formed in volumes below 1 ^3.
Overall, it is estimated that 80%–90% of young stellar objects
(YSOs) form in embedded clusters (;
although the exact fraction is somewhat sensitive to the definition of
a cluster, see ).  It is now established
that the fate of these objects is directly linked to the evolution of
the molecular gas, which is responsible for most of the gravitational
potential that binds the stars in the cluster.  Therefore, it is
important to study these objects in their early phases, before they
are subject to infant mortality <cit.>.

Operationally, clusters are often defined and identified as groups of
stars whose stellar surface density exceeds that of field stars of the
same physical type (see, e.g., ).[In this paper
  we focus on the optical identification of a cluster, and we
  explicitly ignore issues such as the gravitational stability of
  overdensities of stars.  Therefore, we will call “cluster” any
  overdensity, irrespective of its boundness.] In this respect,
embedded clusters pose special problems because they are buried in
dust and gas, and therefore are often not even visible in optical
images; moreover, their shape is often elongated and clumpy,
reflecting the initial structure of the dense molecular gas
<cit.>.  However, even with
infrared observations, discovering deeply embedded clusters and
measuring their basic parameters (such as surface density and size)
can still be a challenge since the associated dust extinguishes both
the cluster members and the field stars behind the cloud (see, e.g.,
). In fact, typical optical or
near-infrared observations stellar fields around molecular clouds show
underdensities at the location of the clouds because of the effects of
extinction on the density of background stars. In such a situation the
observed cluster stellar surface density can be comparable or even
less than the unobscured field stellar surface density (see
Fig. <ref>).



Different authors have used different techniques to take into account
the effects of extinction in the identification and characterization
of star clusters. <cit.> built density maps in
the direction of nearby molecular clouds using the 2MASS Point Source
Catalog, and corrected the effect of dust extinction using publicly
available CO maps (converted into infrared extinction with a constant
X-factor). In his technique, the extinction correction is only applied
to field stars: the cluster stellar density is obtained by 
subtracting an extinction-corrected model of the field stellar density
to the observed stellar density. Moreover, the use of CO maps and the
largely uncertain X-factor represent a primary source of error.

An opposite approach has been adopted by <cit.>,
who studied star clusters embedded in the Rosette nebula using 2MASS
data.  In this case, the local extinction was determined directly
from the stellar near-infrared colors, and a correction for the
effects of extinction was applied to each star. In practice, as noted
by the authors themselves, the correction applied would only be
correct for background stars, while it is used for all stars
(including the foreground, which should not be corrected, and embedded
ones, which should be partially corrected).  Things are further
complicated because molecular clouds are know to have steep gradients
in their column densities, and these, if undetected (because of
resolution effects, which is surely the case of the Rosette cloud at
the 2MASS resolution), would introduce further biases in the
extinction map and in the correction to the cluster richness.

<cit.> use yet another technique (similar to
) to study the properties of IC 348 in
Perseus. The basic idea is that two independent measurements of
extinction in molecular clouds are available, one from the star color
excess (here applied to H - K color), and one from the star number
counts. In absence of contaminants, both measurements should
agree. The presence of a cluster, however, will affect both the color
excess measurement (in a way that depends on the intrinsic color of
the cluster members) and the star count method (in a way that depends
on the location of the cluster within the cloud). The difference of
the two extinction estimates is a proxy for the cluster density. The
various assumptions made when applying this technique (in particular,
the degree of embedding of the cluster within the cloud) needs to be
resolved by calibrating it using independent measurements of cluster
richness: this clearly limits its application.

In this paper we present a new methodology to identify and
characterize clusters in wide-field or all sky, multi-band, infrared
imaging surveys.  The method is based on the production of
extinction-corrected maps of stellar surface density, and takes
advantage of the different H-band luminosity functions of embedded
clusters with respect to that of the background population.
Additionally, in contrast to the methods described above, it is able
to correct for cluster members unidentified because of extinction in a
way that takes into account the position of the cluster within the
cloud along the line of sight.  The technique is based on a rigorous
mathematical framework; this provides a clear advantage, in that we
can perform a detailed statistical analysis of the method. But the
detailed mathematical derivation of this method might not be of as
much interest to astronomers as its implementation.  However, we
stress that those readers not interested in the detailed mathematical
aspects of the derivation can still benefit from the method proposed
here, because its application is simple and straightforward:
Eq(<ref>), with optionally a second expression for the noise
estimate, Eq. (<ref>).  We also provide a pseudo-code for it in
the appendix.

This paper is organized as follows. In Sect. 2 we present the general
framework and discuss the standard technique generally employed to
identify clusters.  In Sect. 3 we present our new method and we
discuss its practical implementation in Sect. 4.  Simple numerical
tests used to validate the method are described in Sect. 5, while the
results of an application to the Orion molecular complex are discussed
in Sect. 6.  Finally, we summarize the results obtained in this paper
in Sect. 7.



§ THE STANDARD TECHNIQUE




















Traditionally, star clusters are identified as overdensities in the
distribution of stars.  Although many different options are available
to estimate the local angular density of stars σ
(see, e.g., ), we consider in this paper the simple “moving
average” estimate[Throughout this paper we will use “hats”
  to denote measured quantities.  Hence, σ(x⃗) is the
  true star density at the angular position x⃗, while
  σ̂(x⃗) is the measured star density at the
  same position.]

    σ̂(x⃗) = ∑_n=1^N W(x⃗ - x⃗_n
      )   .

Here n runs on the individual stars, {x⃗_n } are the
locations of the various stars, and W is a window function,
normalized to unity:

    ∫ W(x⃗')  ^2 x' = 1   .

So, by construction, W has a unit of one over area.  As a simple
example, suppose that W is taken to be a top-hat function:

    W(x⃗') =
      
        1 / (π r^2)    if |x⃗'| < r   ,
    
        0                    otherwise  .

In this case, the sum of Eq. (<ref>) runs just over the N_r stars
within an angular distance r from the point of the estimate, and the
observed estimates reduces to

    σ̂= N_r/π r^2  .

For a constant density of stars, this quantity is unbiased, in the
sense that the mean (ensemble average) E of σ̂ is
just the true density of stars:

    [ σ̂] ≡⟨σ̂⟩ =
      [N_r]/π r^2 = σ  ,

and the associated variance is

    [σ̂] = [N_r]/(π r^2)^2 =
      σ/π r^2  .

This equation shows that the error associated with the measured star
density decreases as N_r^-1/2: therefore, if it
is known that the density of stars is constant, to determine its value
it is sensible to use relatively large window functions W.

More generally, if the star density is variable, the average measured
density can be shown to be a convolution of the true density with the
window function W:

    [ σ̂] (x⃗) = ∫
      W(x⃗ - x⃗') σ(x⃗')  ^2 x'   ,

and the associated two point correlation function is

    [σ̂] (x⃗, x⃗') 
       ≡[ ( σ̂(x⃗) - 
       [σ̂] (x⃗) ) (
      σ̂(x⃗') -  
       [ σ̂] (x⃗') ) ]
      
        = ∫ W(x⃗ - x⃗”) W(x⃗' - x⃗”) 
      σ(x⃗”)  ^2 x”  .

Equation (<ref>) shows that W sets the scale for the resolution
of the density map; similarly, Eq. (<ref>) shows that points
close enough in the density map will have correlated noise.
Therefore, if one aims at finding changes in the star density (such as
a star cluster), the window function should have a scale that is
comparable or smaller than the typical size of the variations of the
star density.  However, this is in tension with the noise properties
of Eq. (<ref>), because a small window function implies a large
noise.

In order to be able to detect a star cluster, the measured density of
stars at the location of the cluster must differ from the average
density much more than the standard deviation of σ̂.  Hence,
the quantity [σ̂](x⃗, x⃗) sets the sensitivity
of the cluster finding algorithm.  In general the true density σ
= σ_field + σ_cluster is the sum of the
field density and of the cluster density, and in some cases
σ_cluster≪σ_field.  In these
conditions the error is dominated by the shot-noise due to the field
star population.  In reality, many other effects can prevent the
discovery and characterization of star clusters:


  * Extinction by dark nebulæ, which reduces the surface density
  of background sources;

  * The galactic structure, which induces smooth large-scale
  variations;

  * Differences in the sensitivity across a mosaic image due to
  changes in the observational conditions;

  * Other systematical observational effects, such halos produced by
  bright objects within the image;




 §.§ Extinction correction



Among the effects listed above, the first one is particularly
important for young clusters, since these objects tend to be deeply
embedded and thus can escape a detection; additionally, detected
clusters are plagued by large uncertainties in their physical
properties (number of stars, mass, and size). For this reason,
<cit.> developed a technique to perform a simple
correction of the results produced by Eq. (<ref>). They noted
that the density of background stars observed through a molecular
cloud decreases by a factor 10^-α A, where α is the
exponential slope of the number counts and A is the extinction, both
in the band used for the observation.  Therefore, to account for the
undetected stars one can just multiply the local density estimate
σ̂(x⃗) by 10^α A(x⃗), where A(x⃗) is
the local estimate extinction (i.e., the extinction derived from an
extinction map at the location of the x⃗).

The problem of this approach is that it uses the same correction
factor for foreground stars, embedded objects, and background stars,
which generally results in an overestimate of the local corrected star
density. Additionally, the same correction is applied to young stellar
objects (YSOs) and field stars, which however have number count
distributions very different. This leaves a large uncertainty on the
corrected σ̂(x⃗) and, ultimately, on the
characterization of each cluster.




§ MAXIMUM-LIKELIHOOD APPROACH





 §.§ Constant extinction, no weighting



As discussed in the previous section, the error associated to
σ̂ is often dominated by shot noise due to the field star
population.  However, as mentioned, YSOs have photometric properties
that differ significantly from these of field stars: the former have a
H-band luminosity function (hereafter HLF) that can be approximated
by a Gaussian distribution <cit.>, while the
latter have a HLF that follows, up to H ∼18mag, an
exponential with slope α∼ 0.33 (see, e.g.,
).  This suggests that we might use
suitable cuts in the H-band to reduce the number of field stars,
while essentially retaining all YSOs, thus gaining in terms of noise.
However, this naive procedure is difficult to apply, since both YSOs
and field stars are also affected by extinction, which would change
the shape and the peaks of the respective HLFs.  Therefore, we are
forced to use a more systematic approach.

Let us first model the observed H-band luminosity function for a set
of angularly close stars (i.e., in a patch of the sky).  We suppose
that the true, unextinguished HLF can be described as a mixture of L
different HLFs, corresponding to different stellar populations (in the
situation considered later on in this paper we will use just two
components, one for the field stars and one for the YSOs, but
deeper observations might require the inclusion of a third component
corresponding to galaxies and following exponential number counts with
a slope of 0.6; see ). The observed HLF
differs from the true one because of extinction and because of
incompleteness in the faint end; additionally, we will also have
photometric errors, but generally these will be small compared to the
typical width of the various HLF components, and for the sake of
simplicity are therefore neglected here.[Photometric errors
  can be easily included in Eq. (<ref>) by replacing p_i(m - A)
  there with the convolution of this function with the photometric
  uncertainty.  In presence of different uncertainties for each star,
  one will need to specialize Eq. (<ref>) to each star.  A
  similar technique can be used to include the effects of errors on
  the extinction measurements.] We can thus write finally the
observed HLF, measured in units of stars per area per magnitude bin, as

    σ(m) = c(m) ∑_i=1^L σ_i p_i(m - A)   ,

where A is the H-band extinction,[In this section we will
  take the extinction to be identical for a set of angularly close
  stars; we will then relax this assumption and use the individual
  extinction in the direction of each star.] c(m) is the
completeness function (i.e. the probability to detect a star with
apparent magnitude m), p_i(m) is the probability distribution for
the i-th component of the HLF, and σ_i is a coefficient that
indicates the predominance of the i-th component in the HLF.  In
order to give σ_i a simple interpretation, it is useful to take
p_i suitably normalized.  In particular, we assume that

    ∫ c(m) p_i(m)   m = 1   .

With this choice, σ_i can be identified as the observed
angular density of stars for the i-th component in absence of
extinction.

Our goal is to infer the parameters {σ_i } from a set of
observed H-band magnitudes { m_n } in a patch of the sky of area
S; we will then repeat this operation in different directions in
order to build maps of densities for the various components.  This
will allow us to identify overdensities in each component, such as in
the YSO one.  To this purpose, we use Eq. (<ref>) to write
the log-likelihood following the prescriptions of
<cit.>.  The entire analysis is done in the
magnitude space, by using S σ as the probability distribution
density for the magnitudes:

    lnℒ =    ∑_n=1^N ln S σ(m_n) - ∫
      S σ(m)   m 
     =    ∑_n=1^N [ ln S c(m_n) ∑_i=1^L
      σ_i p_i(m_n - A) ] 
        - ∫ S c(m) ∑_i=1^L σ_i p_i(m - A)   m
        .

This likelihood can be used in Bayes' theorem to infer a posterior
probability for the densities {σ_i }.  More simply, we can
just find the most likely parameters {σ_i } using a maximum
likelihood approach.  To this purpose, we consider

    ∂lnℒ/∂σ_i = ∑_n=1^N
      c(m_n) p_i(m_n - A)/σ(m_n) - S ∫ c(m) p_i(m-A)   m   .

The maximum likelihood solution is given by the set of densities
σ_i that maximize ℒ or, equivalently, lnℒ, i.e. by the zero of Eq. (<ref>); that is, we need
to solve

    ∑_n=1^N c(m_n) p_i(m_n - A)/σ(m_n) = S ∫
      c(m) p_i(m-A)   m   .

Unfortunately, the solution of this equation for L > 1 cannot be
provided in analytic form, and must be obtained by numerical methods.
We will discuss below in Sect. <ref> simple ways to
find it.

We can obtain an estimate of the errors associated to the measurements
σ_i from the Fisher information matrix (see, e.g.,
), that we recall is defined as

    I_ij = [ ∂lnℒ/∂σ_i∂lnℒ/∂σ_j] = - [ ∂^2 lnℒ/∂σ_i  ∂σ_j]   .

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:

    [σ] ≥ I^-1  .

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.

The Fisher information matrix can be readily evaluated from our data
using Eq. (12) of <cit.>:

    I_ij = S ∫c^2(m) p_i(m-A) p_j(m-A)/σ(m)  m   .

This relation is interesting from several points of view.  First, note
that the Fisher matrix contains elements outside the diagonal, unless
the probability distributions for the various components do not
overlap, i.e. have non-intersecting support: this would basically
mean that we could immediately tell to which component belongs a star
from its unextinguished magnitude.  Second, note that all elements of
the matrix are non-negative and therefore, in case of L=2
components, [σ] ≃ I^-1 will have non-positive
elements off-diagonal: in other words, the measurements of the two
densities σ_1 and σ_2 will be affected by a negative
correlation.  This is expected and corresponds to a classification
error: if p_1 overlaps with p_2, we cannot uniquely associate
stars to each different population and in general an overestimate of
one population is done at the expenses of an underestimate of the
other population.

It is instructive to consider the form of the information matrix in
the special case where L=1.  When just one component is
used, then I is a scalar and it reduces to

    I = S/σ∫ c(m) p(m-A)   m =
      S^2/[N]  ,

where [N] = S σ is the average number of stars observed in
the area S.  Its inverse, I^-1, is therefore [N] /
σ^2, as expected from a simple Poisson statistics.

With L ≥ 2, in principle we can encounter cases where the Fisher
information matrix is singular.  Given the analytic form of I, this
happens in particular when the two components have the exact same
probability distributions within the support of c(m), i.e. c(m)
p_i(m) = c(m) p_j(m): in this case, the corresponding rows (or
columns) of I are identical.  In such a situation, clearly, it
is virtually impossible to classify objects as one of the two
degenerate components, and therefore the uncertainty on the respective
densities σ_i and σ_j are infinite.

For completeness, we also report the expected maximum value of the
log-likelihood

    [ℒ] = L/2 + S ∫σ(m) [ ln
      S σ(m) - 1 ]   m   ,

and the associated variance

    [ℒ] = S ∫σ(m) ln^2 S σ(m)   m   .

These equations can be used to verify that the chosen model
(<ref>) can account well for the data.




 §.§ Differential extinction and spatial weight



So far, we have assumed that all stars are subject to the same
extinction A; moreover, we have not weighted stars depending on
their angular position as done in Eq. (<ref>).  In this section
we intend to remove these two limitations and consider the full
problem. 

The simpler approach to include the effects of differential extinction
is to consider the joint density in magnitudes and positions.
In this framework, we can rewrite Eq. (<ref>) for log-likelihood
as

    lnℒ =    ∑_n=0^N lnσ(m_n, x⃗_n)
      - ∫ m ∫^2 x'  σ(m, x⃗')   .

In this equation the quantity σ(m, x⃗') represents the
predicted density of stars with magnitude m at the location x'.
Similarly to Eq. (<ref>), we write this quantity as a mixture of 
different densities, corresponding to different stellar populations:

    σ(m, x⃗') = c(m) ∑_i=1^L σ_i(x⃗')
      p_i ( m - A(x⃗') )   ,

where σ_i(x⃗') represents the density of stars of class i
at the sky position x⃗'.  In order to proceed, we need to model
these densities.  A simple and consistent way of doing this, is to
suppose that lnσ_i(x⃗') can be written as the weighted
sum of two terms: one, lnσ(x⃗), associated to the density
at the point of interest x⃗ (the point where we intend to
evaluate the local densities of stellar populations); and one, lnτ_i(x⃗'), which parametrizes the local changes of the densities.
As a result, we write

    σ_i(x⃗') = ( σ_i(x⃗)
      )^ω(x⃗ - x⃗')( τ_i(x⃗')
      )^1 - ω(x⃗ - x⃗')  .

The function ω describes the spatial correlation between
densities at different positions and plays a central role in
identifying which stars contribute to the density estimate at x⃗.

We can now insert Eqs. (<ref>) and (<ref>) in
Eq. (<ref>) and find the maximum likelihood solution over
σ_i(x⃗).  Calling A_n ≡ A(x⃗_n) the extinction
suffered by the n-th star, we find

    ∂lnℒ/∂σ_i(x⃗) =
      ∑_n=1^N ω(x⃗ - x⃗_n) p_i (m_n - A_n)
        σ_i(x⃗_n) / σ_i(x⃗)/∑_j σ_j(x⃗_n) p_j(m_n - A_n)
          - ∫ m   c(m) ∫^2 x'  ω(x⃗')
      p_i( m - A(x⃗') ) σ_i(x⃗') / σ_i(x⃗)
        .

We now make an important assumption.  Because of the form of the
parametrization (<ref>), a solution for the maximum likelihood
can only be found if we specify the functional form of the functions
τ_i(x⃗').  However, these functions are truly unknown, since
they parametrize the local changes of the various densities, which in
turn depend on the local density map.  Using a statistical approach,
however, it is natural to assume that these functions have the same
distribution of σ_i(x⃗), the quantity we are interested
in.  As a result, if we take an average of Eq. (<ref>), terms
such as σ_i(x⃗') / σ_i(x⃗) cancel out:

    ∂lnℒ/∂σ_i(x⃗) =
      ∑_n=1^N ω(x⃗ - x⃗_n) p_i (m_n - A_n)/∑_j
        σ_j(x⃗) p_j(m_n - A_n)
          - ∫ m   c(m) ∫^2 x'  ω(x⃗')
      p_i( m - A(x⃗') )   .


Before proceeding, it is useful to consider the solution of the
maximum likelihood approach in the simple case where there is a single
population of stars and where the extinction vanishes, A(x⃗') =
0.  We find in this case

    ∂lnℒ/∂σ(x⃗) =
      ∑_n=1^N ω(x⃗ - x⃗_n)/σ(x⃗) -
      ∫^2 x'  ω(x⃗')   ,

where we have used the normalization property (<ref>).  The
solution of this equation is immediately found as

    σ(x⃗) = ∑_n=1^N W(x⃗ - x⃗_n)   ,

where

    W(x⃗) = ω(x⃗)/∫ω(x⃗')  ^2
        x'

We have therefore recovered Eq. (<ref>), with the correct
normalization (<ref>) for the weight W.

In the general case, the maximum-likelihood solution of
Eq. (<ref>) must be obtained numerically.  The errors associated
to the solutions can be estimated using the Fisher matrix.  In our
case, we can obtain the Fisher matrix from

    I_ij(x⃗) =    ∬ m  ^2 x'  1/σ(m,
        x⃗')∂σ(m, x⃗')/∂σ_i(x⃗)∂σ(m, x⃗')/∂σ_j(x⃗)
     =    ∬ m  ^2 x'   c^2(m) ω^2(x⃗ -
      x⃗') ×
           p_i ( m - A(x⃗') ) p_j ( m -
        A(x⃗') )/σ(m,
        x⃗')σ_i(x⃗') σ_j(x⃗')/σ_i(x⃗) σ_j(x⃗)  .

As before, we can replace x⃗ to x⃗' in all arguments of
σ_i, with the justification that the statistical properties of
σ are invariant upon translation.  We will discuss in the next
section useful expressions to evaluate this quantity in practical
cases.





 §.§ Implementation



The algorithm proposed in this paper is essentially the search of the
solutions of Eq. (<ref>) as a function of the star densities {σ_i } corresponding to the various components or star
populations.  The same procedure must be applied to different patches
of the sky, so that maps of the star densities can be obtained.
These, in turn, will allow us to identify and characterize star
clusters, and in particular embedded ones.

Although the practical implementation of the algorithm follows this
schema, a number of technical and theoretical aspects must be
correctly addressed in order to optimize the detection and make the
technique efficient.

First, we note that a simple way to obtain the (positive) solutions of
Eq. (<ref>) is through the use of a recursive formula.  To this
purpose, multiply both members of this equation by σ_i(x⃗),
and solve for this quantity, thus obtaining the expression

    σ_i(x⃗) ←∑_n=1^N ω(x⃗ - x⃗_n) σ_i(x⃗)
          p_i(m_n - A_n)/∑_j=1^L σ_j(x⃗) 
          p_j(m_n - A_n)/∫ m   c(m) ∫^2 x'  ω(x⃗')
        p_i( m - A(x⃗') )  .

Unfortunately, we are unable to use this equation because we only know
the extinction at the locations of the stars A_n ≡ A(x⃗_n):
this prevents us from evaluating the integral over x^2 in the
denominator.  We can, however, move the denominator inside the sum,
and evaluate the second integral by replacing A(x⃗) with A(x⃗_n), the extinction at the direction of each star.  Additionally,
using the same argument that has been employed in Eq. (<ref>),
that is the similarity between σ and τ, it is convenient
to replace σ_i(x⃗) with σ_i(x⃗_n) in the
numerator.  This procedure leads to the iteration

    σ_i(x⃗) ←∑_n=1^N 
        W(x⃗ - x⃗_n) σ_i(x⃗_n)
          p_i(m_n - A_n)/∑_j=1^L σ_j(x⃗_n)
          p_j(m_n - A_n)/∫ m   c(m) p_i(m - A_n)  .

Equation (<ref>) is the solution proposed in this paper to
estimate the local density of stars.  As indicated by the left arrow
symbol, we can obtain the set of values {σ_i } by starting
with some arbitrary (positive) values for these quantities, and then
by calculating updated values of σ_i by applying
Eq. (<ref>).  The convergence is usually reached within a few
tens of iterations.

Note that Eq. (<ref>) has a simple interpretation.  Let us
ignore for a moment the weight W, i.e. let us assume that all stars
have the same weight.  The sum in Eq. (<ref>) is carried out
over all stars in the patch of the sky, but each star is counted only
partially (i.e., contributes with a term between zero and unity in the
sum): precisely, each star contributes by the computed probability
that the star be associated with the i-th component.  The way this
probability is computed is actually a simple application of Bayes'
theorem, where p_i(m_n - A_n) plays the role of the likelihood,
σ_i(x⃗_n) is proportional to the prior that the star is of
class i, and the sum over j in the denominator is proportional to
the evidence.  The result of the sum of these terms is divided by the
by the result of the integral: this is a correcting factor that takes
into account the fact that, because of extinction and incompleteness,
we will miss a fraction of stars.  Note also that Eq. (<ref>)
can be also considered as a special case of a K-means soft
clustering algorithm where the only unknown quantities are the classes
σ_i <cit.>.

Before proceeding, it is useful to recall the hypotheses of this
algorithm and its strengths.  First, we assume that we have some
knowledge of the H-band luminosity function for the various
populations of stars that are likely to be present in the field. In
practice, we will generally use two probabilities, one for the field
stars, and one for the YSOs. Second, we assume that we have measured
the extinction A_n of each star: note that this is not the
average extinction at the location of the star, which might be very
different because of perspective effects: for example, a foreground
star in front of a cloud would have A_n ≃ 0, while the average
extinction would be significant. This way, the algorithm can directly
account for foreground contamination: foreground stars will not be
corrected in their counts, since the integral in the denominator of
Eq. (<ref>) will evaluate to unity. Similarly, stars within
molecular clouds will be corrected only for the amount of material
that is really in front of them.

Finally, we stress that the iterative procedure proposed here only
find positive solutions for the values σ_i.  Although
reasonable, nevertheless this choice inevitably introduces biases in
the results: for example, in a region where no YSO is present,
because of errors we will still measure small positive values for the
density of YSOs.  However, numerical tests have shown that the
bias amount is limited; moreover, a reduction of the bias is
associated to a large increase in the scatter.  Therefore, we will
force the σ_i to be positive and use Eq. (<ref>) for the
solution.

The uncertainties associated to Eq. (<ref>) can be computed from
the Fisher matrix.  For practical applications, it is convenient to
rewrite Eq. (<ref>) by replacing the integrals over m
and ^2 x' with a sum over the observed objects.  This leads to
the approximated Fisher matrix expression

    I_ij = ∑_n=1^N ω^2(x⃗ - x⃗_n) c^2(m_n)
        p_i(m_n - A_n) p_j( m_n - A_n)/σ^2(m_n, x⃗_n)  .

In this equation, we take the spatial function ω to be
normalized such that

    ∫ω(x⃗')  ^2 x' = ∫ω^2(x⃗')  ^2 x'   ;

that is, in terms of W,

    ω(x⃗) = W(x⃗)/∫ W^2(x⃗')  ^2 x'  .

Table <ref> reports a few common choices for the spatial
function ω(x⃗) and the corresponding weight function
W(x⃗), both correctly normalized.  As usual, the covariance
matrix associated with the measurements of the densities σ_i
can be computed from the inverse of the Fisher matrix, I^-1.




§ SIMULATIONS




















In order to test our method and verify its robustness we have
performed a set of simulations.  We have considered a small patch of
the sky with two different stellar populations: field stars, with
exponentially increasing number counts p_1(m) ∝ 10^α m,
with α = 0.33mag^-1; and YSOs, with Gaussian
number counts p_2(m) ∝exp(-(m-m_0)^2 / 2 s^2 ),
with m_0 = 12mag and s = 1.65mag.  We have
distributed both populations randomly in the small patch of the sky
considered and we have added to each star a random extinction drawn
from the probability distribution p_A(A) ∝ A^-2, in the
range A ∈ [0.1, 2.0] mag.  This choice is meant to simulate
the effects of differential extinction for objects within molecular
clouds.  Finally, we have imagined that our observations are complete
up to 15mag, and that no stars can be observed beyond that
magnitude: in other words, we have modeled the completeness function
as a Heaviside function c(m) = H(15 - m).  This way, our final
catalog contains, for each star, the angular position in the sky, the
H-band magnitude, and the measured extinction.  Note that the
parameters used here are chosen to simulate a real situation
corresponding to the sample application of
Sect. <ref>, that is the Orion molecular cloud
complex observed with 2MASS.

We have used these data in our algorithm, represented by
Eqs. (<ref>) and (<ref>).  As weight function W(x⃗)
we have used a Gaussian, and we have chosen the angular units so that

    ∫ W^2(x⃗')  ^2 x' = [ ∫ W(x⃗')  ^2
      x' ]^2   .

This choice guarantees that the effective area of the weight function
is unity, i.e. the effective number of stars used for the analysis,
in absence of extinction, would be just σ_1 + σ_2.   In
reality, the presence of the extinction reduces this number by a
factor that depends on the relative ratio between σ_1 and
σ_2 (typically, by ∼ 20%).

We have then performed different simulations, with various choices for
σ_1 and σ_2, to verify the ability of the algorithm to
recover the input densities.  Figures <ref>–<ref>
show the observed distributions of σ̂_1 and σ̂_2,
together with the predicted ones (Gaussian distributions centered on
the true values of σ_1 and σ_2, with the variances
predicted from the Fisher matrix I).  In general, we can note a few
points:


  * When one of the input densities is below ∼ 7, there is a
  clear excess of the corresponding quantity for small measured
  values.  This is a consequence of the fact that the proposed
  algorithm only returns positive values for the densities.

  * Except for the point above, the predicted distributions
  reproduce very well the measured data.  The agreement is
  particularly good when the input densities are large.

  * Overall, the total density σ_1 + σ_2 is better
  constrained than the individual densities σ_1,2.

Figures <ref> and <ref> show the biases and the errors
on the measured densities σ̂_1,2 for σ_1 = 10 and
σ_2 varying in the interval [0, 20].  We can see that there
is a measurable bias for σ_2 < 5, while the results are
essentially unbiased above this value.  Correspondingly, we observe in
the same range σ_2 ∈ [0, 5] a measured scatter in measured
densities that is significantly smaller than what predicted from the
Fisher matrix.  For larger values of the input density the error
estimate gets closer to the measured errors, but still remains
slightly above.  This is actually good, because implies a small
overestimate of the error which will make the entire method more
robust for cluster detections (that is, it will decrease the number of
false positive).  In summary, all these simulations confirm that the
method works and that the error estimate is quite accurate.




§ SAMPLE APPLICATION: ORION
















  
  




  
  




We have applied the method proposed in this paper to 2MASS data of the
Orion A and B molecular cloud complexes.  The regions selected ensure
that we can verify the reliability of the algorithm proposed here using
some of the best studied objects in the sky.  In particular, the
populations of embedded clusters for both clouds have been the
targets of extensive observational campaigns using ground-based, 
near-infrared <cit.>) and millimeter-wave
<cit.> surveys as well as space-based, mid-infrared
Spitzer Space Telescope <cit.>, and Chandra X-ray 
<cit.> surveys.  Additionally, the distance
of these clouds ensures that the 2MASS H-band data are, in absence of
extinction, complete for YSOs: that is, the cluster HLF at the distance of
Orion is essentially entirely within the 2MASS H-band limiting
magnitude (∼15mag).



Figures <ref> and <ref> show the density of fiducial
YSOs measured in Orion A and B.  These maps have been produced by our
pipeline together with Nicer <cit.> and
Nicest <cit.> extinction maps from the
2MASS Point Source Catalogue (see  for
details on the data selection and extinction map construction).  The
algorithm has been run in conditions similar to the simulations
described above: that is, we have used two different stellar
populations, one associated with the background and characterized by
exponential number counts, and one associated with the YSOs and
characterized by Gaussian number counts (with parameters consistent
with the H-band luminosity function of
).  Using an angularly close control
field we also measured the distribution of intrinsic colors of stars
and the shape of the completeness function: the latter has been
modeled using an complementary error function erfc, as described in
Appendix <ref>. This choice makes it possible
to use the entire 2MASS catalogue without any magnitude cut (which
would increase the noise in the final data products).  The maps
produced have a pixel size of 1.5arcmin and a weight function W
∝ω in turn proportional to a Gaussian with 𝐹𝑊𝐻𝑀
= 3arcmin.  We have used a relatively large beam in these maps
to increase the sensitivity of our algorithm and to minimize the
effects of the biases shown in the simulations described in
Sect. <ref>, while still retaining in most situations
the ability to distinguish different clusters (i.e., avoid confusion
effects at the distance of Orion).

Since we have at our disposal the covariance map of these
measurements, we have assessed the reliability of each density peak in
these figures.  The red contours in the figures show the areas (larger
than 2 pixels) where the local YSO density exceeds
1.5stars/pixel, corresponding approximately to a
signal-to-noise ratio larger than 3.5: note how some regions in black
in Fig. <ref> do not reach the threshold because of the large
error associated with them (mostly due to the high extinction values
there).

Table <ref> shows the YSO clusters identified in the Orion A and
B areas using our simple prescription, together with the most relevant
parameters. In some cases we clearly see that angularly close clusters
appear as a single contour in our maps: the simple procedure used here
to define clusters, the relatively coarse resolution used, and the
cluster morphology itself prevent us from deblending some close
objects. An extreme case of this situation might be the ISF (the
Integral Shaped Filament) cluster, where the limitations due to
angular resolution would make it difficult to resolve and separate
smaller clusters if they exist in such a very populous region. We note
that the ISF cluster encompasses M42, the Trapezium and ONC clusters
as well as an extended population of YSOs along the ISF.

The radius R reported in the table corresponds to the radius of a
circle that would occupy the same area as the identified cluster,
i.e. to the connected region of the sky where the inferred
density of YSOs exceeds the background by 3σ.  At the estimated
distance of Orion, 413pc <cit.>, 1'
corresponds to 0.12pc: therefore, the clusters identified have
radii spanning from ∼2.4pc to ∼0.15pc.

The well known clusters in these clouds are correctly identified by
our procedure.  It is interesting to compare Table <ref> with
clusters identified independently using much more secure data.  Among
the ones at our disposal, the recent catalog of embedded YSOs
obtained by <cit.> using the Spitzer
  Space Telescope and the Chandra observatory is probably
the most secure and complete: we will therefore focus on this catalog.
Since our definition of a cluster is based on an somewhat arbitrary
parameters (signal-to-noise threshold, minimum number of pixels, no
correction for the stars missed at the boundaries), and since
different, more-or-less, arbitrary parameters are also used by
<cit.>, we find it more appropriate and fair to
make a comparison after we first homogenize the data. Specifically, we
take <cit.> YSO list and we make out of it
a density map using a Gaussian kernel of the same size of the one used
for our map. Figures <ref> and <ref> show the results
obtained for Orion A and B, which clearly compares very well with our
own maps, derived purely from the 2MASS point source catalog. The most
relevant structures are present in both maps and have very similar
shapes; the only differences are the noise present in our maps
(however at a relatively low level), and the limited coverange of
the Spitzer derived density maps.

The qualitative similarity of these maps can be quantified if compare
clusters identified in both maps using the same
criteria. Table <ref> shows a list of clusters identified in the
smoothed <cit.> maps using a fix density
threshold (set to 1.5stars/pixel). In this table we compare the
number of Spitzer YSOs with the number of YSOs predicted from
the integral of the σ_YSOs over the area of each
cluster as defined from Megeath et al. density map, together with the
computed 1-σ error. It is clear that in almost all cases we
find an excellent agreement, although in many cases our estimates are
slightly larger than the ones by Megeath et al. We can speculate that
this is due to the presence of class III YSO, which likely would
be missed by Spitzer. Indeed, a comparison of the two panels
of Fig. <ref> shows that the bottom panel, corresponding to
our density map, has spatially more extended clusters than the top
panel, corresponding to Megeath et al. density map.

As discussed earlier on, our algorithm is a statistical one and works
best when it is applied to a sizeable number of stars. However, we can
also push it and associate to each single star a probability of being
a YSO: to this purpose, for the n-th star we can compute

    P_i = σ_i(x⃗_n) p_i(m_n - A_n)/∑_j=1^L
        σ_j(x⃗_n) p_j(m_n - A_n)  .

Note how this quantity resembles the term within the outer sum of
Eq. (<ref>).  Figure <ref> shows the distribution of
P_YSO (that is, the distribution in the probabilities
assigned to each object to be a YSO) for the
<cit.> YSO candidates and for all the other
objects. It is clear how all the other objects have P_YSO
that is concentrated around zero, while the YSO candidates have
a much broader distribution that extends to unity.  For these latter
objects the distribution, in addition to a substal peak at
P_YSO = 1, shows a long tail up to small values of
P_YSO: this is not unexpected, since our identification
is only statistical (and therefore we cannot identify unambigously
YSOs). Note also how the relatively low values of
P_YSO for some genuine YSOs in our algorithm are
compensated by the small tail in the distribution of field stars (this
works because there are many more field stars than YSOs, a fact
that is accounted for in the algorithm).



 §.§ Sensitivity to the distributed population
 

Recently, , have identified a rich and
well-defined stellar population of about 2 000 objects, mostly M
stars without extinction or infrared-excesses, as the low-mass
counterpart to the Orion OB Ib subgroup (the Orion belt
population). This low-mass population is not obviously clustered but
instead appears to be distributed across ∼ 10 square degrees and
the authors speculate that it could represent the evolved counterpart
of a Orion nebula-like cluster.  While more data is needed to test
this scenario, it is clear that much can be learned about the origin
of stellar OB associations and the dispersal of clusters into the
Galactic field if one is able to trace in a robust manner the
distribution of the slightly older and more expanded populations
surrounding star-forming molecular clouds.

We now investigate the extent to which the technique proposed here is
suitable for detection of looser, more expanded distributions of young
stars, in particular the low-mass counterpart to the Orion OB
association presented in .  For this
purpose, we have built a lower resolution map of the region, employing
a FWHM of 30arcmin.

Figure <ref> shows that, surprisingly, we are well able to
recover the stellar over-density of the Ori Ib population, and for the
first time, the stellar over-density of the Ori Ia population. These
over-densities are likely to be created by low-mass stars as 2MASS is
still sensitive to the peak of the IMF for the putative distance and
age of these subgroups. An analysis of the substructure seen in the
distributed population visible in Figure <ref> above the noise
pattern is beyond the scope of this paper, but will best addressed
once Gaia parallaxes are generally available. Of relevance for this
paper is that the ability of the method to trace the dispersed
population from existing all sky data opens a new window on the
unsolved problem of the origins of OB association and cluster
dispersal in to the field.










§ CONCLUSIONS



The following items summarize the main results presented in this
paper:


  * We have developed a new method to discover and characterize
  deeply embedded star clusters.

  * The method is able to statistically classify objects as field
  stars or YSOs and corrects for the effects of differential
  extinction.

  * We have provided expressions for the covariance of the inferred
  densities and we have validated both the method and the analytic
  expression for the covariance with a set of simple but realistic
  simulations.

  * We have applied the new method to 2MASS point sources observed
  in the Orion A and B and we have shown that we can identify and
  characterize well protostellar clusters in these regions, as well as
  detect much looser associations such as the OB 1a and 1b subgroups.


Finally, we note that the method proposed here can be easily extended
to multi-band observations by using suitable probability distributions
p_i(m⃗) for various populations as a function of the magnitude
vector m⃗.  Its implementation would be essentially identical,
with minor modifications to take into account the different effects of
extinction on different bands.  The natural use of such an extension
would be in the context of techniques such as the one proposed by
<cit.> which are able to recover the extinction
from a complex analysis of multi-band data.


  This research has made use of the 2MASS archive, provided by
  NASA/IPAC Infrared Science Archive, which is operated by the Jet
  Propulsion Laboratory, California Institute of Technology, under
  contract with the National Aeronautics and Space Administration.e
  Additionally, this research has made use of the SIMBAD database,
  operated at CDS, Strasbourg, France.






§ IMPLEMENTATION



In this section we consider a few details and analytical expressions
useful to evaluate some of the expressions needed to implement the
method proposed here.



 §.§ Algorithm



The algorithm proposed in this paper is essentially an iterative
application of Eq. (<ref>) to the set of densities {σ_i
(x⃗) }. A few notes, however, are necessary to better implement
the algorithm:


  * The method assumes that the extinction values in the H band are
  available for each star used. In our implementation these values are
  obtained through the use of the Nicer algorithm. Other
  techniques can be used, as long as the extinction is refereed to the
  single star and is not an average extinction at the
  location of the star. In our implementation, this requirement limits
  the method to stars for which, in addition to the H-band
  measurements, have at least another band photometry.

  * Since the various densities are spatially variable, one needs to
  repeat the iteration of Eq. (<ref>) to each point in the map.

  * Equation (<ref>) itself needs to be iterated a few times
  (typically around ten, in some cases a little more) before reaching
  a good convergence.  This generally suggests the use of simple
  analytical models for the various functions involved in this
  equation (the probability distributions p_i(m) and the
  completeness function c(m), discussed below).





 §.§ H-band probability distributions and completeness function





In this paper we have considered two useful models for the H-band
luminosity function: the exponential distribution and the normal
(Gaussian) one.

We parametrize the exponential distribution as

    p(m) = k β^β m∝ 10^α m  ,

where k is a normalization constant that will be found lather on and
β = αln 10.

The normal distribution is parametrized as

    p(m) = k/√(2 π s^2)^-(m - m_0)^2 / 2 s^2  ,

where again k is a normalization constant to be investigated later.

Finally, we model the completeness function c(m) in terms of the
complementary error function . This can justified from an
empirical point of view since the error function has a shape that
resemble the completeness function.  It is also reasonable from a
statistical point of view if the photometric errors are Gaussian: in
this case, one can suppose that the probability that an object be
detected is an integral over a relevant part of the
Gaussian. Specifically, we assume for c(m) the following functional
form:

    c(m) = 1/2( m -
          m_c/√(2 σ^2_c))   ,

where m_c is the 50% completeness limit and
σ_c sets the width of the completeness
function. Figure <ref> shows that our simple model for the
control field HLF, the product of p(m) c(m) with p(m) following an
exponential distribution and c(m) described in terms of an erfc,
reproduces well the data.




 §.§ Statistical errors and completeness function



As noted above, statistical errors can be included in our algorithm by
convolving the probability distributions p_i with appropriate
kernels.  We expect to have two main sources for statistical errors:
photometric errors in the H-band magnitude measurements and errors
in the extinction measurements.  While the formers are relatively
simple to characterize (typically, they will be provided in the star
catalogue), the latter are a combination of different sources of
errors: the scatter in the intrinsic color of sources and the
photometric errors in each band used.  If extinctions are derived
using the Nicer algorithm, as assumed here, then the variance
associated to each extinction measurement can be computed from the
expression

    (A) ≡σ^2_A = 1/k⃗^T C^-1k⃗  ,

where, following the notation of <cit.>, we have
called k⃗ is the reddening vector and C the combined
covariance matrix of the observed colors (including both the intrinsic
scatter in the color of unextinguished stars and the individual
photometric errors).

Looking again at Eq. (<ref>), i.e. the main equation
representing our method, we see that in reality the H-band
magnitudes and the associated extinction measurements typically enter
the problem through the combination m - A; moreover, while in the
numerator the combination is m_n - A_n and involves therefore
measurements of both m and A for individual stars, in the
denominator the combination is m - A_n, with m integration
variable. Therefore, in presence of errors these two parts of the
expression must be computed separately.

Let us consider first the denominator. Since there the only argument
of p_i with associated errors is A_n, the measured extinction for
the n-th star, we just need to convolve the magnitude distribution
with a Gaussian kernel with variance provided by Eq. (<ref>). 

In the case of the numerator the situation is slightly more
complicated, because there we find the combination m_n - A_n.  Since
A_n is computed from the observed magnitude of each star (including
of course the H-band magnitude), m_n and A_n are correlated.  A
simple calculation shows that the variance associated to the
combination m_n - A_n is

    (m_n - A_n) = k_H^2 σ^2_A + σ^2_H [ 1 - 2  
      (b_J-H - b_H-K)   k_H ]   .

In this expression the last term is related to the correlation between
A and m and is equal to σ^2_H, the square of the
photometric error on the H-band magnitude. The factors b_J-H and
b_H-K are quantities that can be computed from Eq. (12) of
<cit.>, while k_H is the reddening vector for
the H band, i.e. the ratio A_H / A.



 §.§ Convolution integrals



The convolution of the luminosity functions (<ref>) and
(<ref>) with Gaussian kernels representing the statistical
errors take simple forms.  Calling σ^2_p the variance of the
convolution kernel, given depending on the term by Eq. (<ref>)
or by Eq. (<ref>), and using the tilde to represent the
convolved H-band luminosity functions, we find for the exponential
distribution

    p̃(m) = k β^β m + σ^2_p β^2 / 2  .

Similarly, for the Gaussian model we find

    p̃(m) = k/√(2 π (s^2+σ^2_p))^-(m -
        m_0)^2 / 2 (s^2 + σ^2_p)  ,




 §.§ Completeness integrals



Finally, we need to compute the integral in the denominator of
Eq. (<ref>), involving the H-band luminosity function p(m)
and the completeness function c(m). A simple change of variable
casts this integral (similar to a convolution) into a more convenient
form:

    ∫ c(m) p(m-A)   m = ∫ c(m+A) p(m)   m ≡
      K(A)   .

In this equation we have called K(A) the result of the integral,
retaining its explicit dependence on A. We now consider the
derivative K' of K:

    K'(A) = ∫ c'(m+A) p(m)   m   ,

and note that, since c is the error function, c' is a Gaussian
function. Since Eq. (<ref>) is also (essentially) a convolution,
this allows us to compute K'(A) using formulae similar to the ones
of the previous section: as a result, we see that K'(A) remains
either an exponential or a Gaussian, depending on the function form of
p(m).  

Finally, we need to integrate back K'(A) to obtain K. The result
we obtain from this procedure is, in case of an exponential
distribution,

    K(A) = k ^β (m_c - A) + σ^2_cβ^2 / 2  .

For the Gaussian distribution the result contains again the error
function:

    K(A) = k/2(  m_0 - m_c +
          A/√(2(σ_c^2 + s^2)))   .

If statistical errors have to be taken into account, these expressions
simply change into

    K(A) = k ^β (m_c - A) + (σ^2_c + σ^2_p) β^2 / 2  ,
  
and

    K(A) = k/2(  m_0 - m_c +
          A/√(2(σ_c^2 + s^2 + [b]σ^2_p)))   .


  
aa 


