NICEST, a near-infrared color excess method tailored for small-scale structures
aa r X i v : . [ a s t r o - ph ] S e p Astronomy&Astrophysicsmanuscript no. art3 c (cid:13)
ESO 2018October 16, 2018 N icest , a near-infrared color excess method tailored for small-scalestructures. Marco Lombardi , European Southern Observatory, Karl-Schwarzschild-Straße 2, D-85748 Garching bei M¨unchen, Germany University of Milan, Department of Physics, via Celoria 16, I-20133 Milan, ItalyReceived ***date***; Accepted ***date***
ABSTRACT
Observational data and theoretical calculations show that significant small-scale substructures are present in dark molecular clouds.These inhomogeneities can provide precious hints on the physical conditions inside the clouds, but can also severely bias extinctionmeasurements. We present N icest , a novel method to account and correct for inhomogeneities in molecular cloud extinction studies.The method, tested against numerical simulations, removes almost completely the biases introduced by sub-pixel structures and by thecontamination of foreground stars. We applied N icest to 2MASS data of the Pipe molecular complex. The map thereby obtained showssignificantly higher (up to 0 .
41 mag in A K ) extinction peaks than the standard N icer (Lombardi et al. 2001) map. This first applicationconfirms that substructures in nearby molecular clouds, if not accounted for, can significantly bias extinction measurements in regionswith A K > ff ect, moreover, is expected to increase in more distant molecular cloud, because of the poorer physicalresolution achievable. Key words. dust, extinction – methods: statistical – ISM: clouds – infrared: ISM – ISM: structure – ISM: individual objects: Pipemolecular complex
1. Introduction
Although the details of star and planet formation are very lit-tle known, there is a large consensus that these objects arecreated inside dark molecular clouds from the contraction andfragmentation of dense, cold cores. As a result, in the lastdecades molecular clouds have been studied in detail using manydi ff erent techniques, from optical number counts (Wolf 1923;Bok 1937), to radio observations of carbon monoxide (CO)molecules (Wilson et al. 1970; Frerking et al. 1982), to near-infrared (NIR) extinction measurements (Lada et al. 1994).Molecular hydrogen and helium represent approximately99% of the mass of a cloud, but the absence of a dipole momentin these molecules makes them virtually undetectable at the lowtemperatures ( T ∼
10 K) present in these objects. Hence, the(projected) mass distribution of dark clouds is usually inferredfrom the distribution of relatively rare tracers (such as CO ordust grains) under the assumption that these are uniformly dis-tributed in the cloud.As pointed out by Lada et al. (1994), the reddening of back-ground stars in NIR bands is a simple and reliable method tostudy the dust distribution and thus the hydrogen column den-sity. Compared to optical bands, NIR bands are less a ff ected byextinction and are less sensitive to the physical properties of thedust grains (Mathis 1990), and thus their color excess can bedirectly translated into a hydrogen column density. In a seriesof papers we reconsidered and improved the NIR color excess(N ice ) technique, by generalizing it to use three or more bands(N icer , see Lombardi & Alves 2001) and by taking a maximum-likelihood approach (Lombardi 2005). Send o ff print requests to : M. Lombardi Correspondence to : [email protected]
Although the N ice and N icer techniques have been verysuccessful in studying molecular clouds (see, e.g., Alves et al.2001), they are plagued by two complications: small-scale inho-mogeneities and contamination by foreground stars. NIR stud-ies typically assume a uniform extinction over (small) patchesof the sky; however, in presence of substructures such as steepgradients, unresolved filaments, or turbulence-induced inhomo-geneities (Lada et al. 1994; Larson 1981; Heyer & Brunt 2004),the background stars used to estimate the cloud extinction wouldnot be uniformly distributed in the patch, but would be preferen-tially observed in low density regions. Similarly, in presence offoreground stars, their relative fraction increases with the dustcolumn density. Unfortunately, both e ff ects will bias all colorexcess estimators toward low column densities; moreover, thebias will be especially important in the dense regions of molec-ular clouds, which is where the star formation takes place.In this paper we describe in details the repercussions of smallscale inhomogeneities and foreground stars in extinction mea-surements, and propose a simple method to obtain estimates ofthe column density in presence of both e ff ects that, under sim-ple working hypotheses, is asymptotically unbiased. Our methoddoes not rely on any (usually poorly verified) assumption regard-ing the small scale structure of the cloud; moreover, it can be ap-plied to a general class of NIR extinction measurements (markedspatial point processes).The paper is organized as follows. In Sect. 2 we introducea general smoothing technique virtually used in all cases to cre-ate smooth, continuous maps from the discrete pencil-beam ex-tinction measurements carried out for each background star; wealso show that in two typical applications (moving average andnearest neighbours smoothing) a bias is expected. In Sect. 3 wepropose a new method that can correct this bias without makingany assumption on the underlying form of the cloud substruc- Marco Lombardi: N icest ture. This new method has been tested with numerical simula-tions, as described in Sect. 4, and with a real-case application,presented in Sect. 5. We discuss and comment the results ob-tained in Sect. 6, and finally we briefly present our conclusionsin Sect. 7.
2. Cloud substructure
Measurements of the reddening of stars observed through amolecular cloud provide estimates of the cloud column densi-ties on pencil beams characterized by a diameter of the order ofa fraction of milliarcsec. However, these high resolution mea-surements are strongly undersampled and are a ff ected by largeuncertainties due to the photometric errors and to the intrinsicscatter of star colors in the NIR. Smooth and continuous extinc-tion maps are normally obtained by interpolating and binningthe various pencil beams. Di ff erent authors use di ff erent recipesto smooth the data (see, e.g., Lombardi 2002 for a description ofseveral interpolation techniques), but in most cases the interpola-tion follows a standard scheme. Consider N K -band extinction measurements (cid:8) ˆ A n (cid:9) obtained, for example, from the color excessof N background (see Lada et al. 1994 or Lombardi & Alves2001). The extinction ˆ A at any location in the sky is evaluatedusing a weighted averageˆ A = P n w n ˆ A n P n w n . (1)The weights { w n } are usually chosen to be significantly di ff er-ent from zero only for stars angularly close to the given locationof the map (for example, Cambr´esy et al. 2002 assign a unityweight to the N nearest neighbors). Clearly, Eq. (1) does notinclude slightly more complex situations where, for example,one takes the median of the extinctions measured from objectsnearby a given position (e.g. Dobashi et al. 2008); however, mostof the discussion carried out for the weighted average actuallyapplies also to median or related estimators.The binning in Eq. (1) washes out the cloud substructure onscales smaller than the typical size where the { w n } are signifi-cantly di ff erent from zero, but this is needed in order to havesmooth maps and to increase the signal-to-noise ratio. However,Eq. (1) also introduces a significant bias on the estimated col-umn density. Suppose that, in the region of the sky that we areinvestigating (i.e., in the area covered by the N stars used to es-timate A ) the column density has significant variations. Becauseof this, the local density ρ ( x ) of stars is not homogeneous thoughthe cloud, but rather follows the scheme ρ ( x ) = ρ − α k λ A ( x ) , (2)where ρ is the density of stars where no extinction is present, α is the slope of the number counts, and k λ = A λ / A is the extinctionlaw in the band λ considered. This e ff ect, which is at the origin ofthe historical number count method to measure column densitiesin molecular clouds (Wolf 1923; Bok 1937), also induces a biasin Eq. (1), since regions with smaller extinctions and thus higherdensity of background stars will, on average, contribute more tothe sum in Eq. (1) than regions with large extinctions. Note thatsince the color excess method requires measurements in at leasttwo di ff erent bands, one should replace Eq. (2) with a version For simplicity, in this paper we drop the subscript K normally usedin literature to denote the K -band extinction A K ; it is also obvious thatthe same method applies to extinction measurements referred to anyband (e.g., the visual extinction A V ). that provides the expected density of stars observed in all bandsrequired for the application of the method (e.g., H and K ). Forequally deep observations on all bands, the result has the sameform of Eq. (2), where k λ is the extinction law of the bluer bandconsidered (e.g., H ).The main issue considered in this paper is the bias intro-duced by unresolved substructures in Eq. (1). In general, the biasshould be defined here as the di ff erence between the expected,mean value of h ˆ A i and the true column density A at the same po-sition. However, clearly every smoothing technique introduces abias because of the smoothing itself, and this bias is normally ac-ceptable if it does not introduces a systematic e ff ect on the total column density. In other words, a required property is Z A ( x ) d x = *Z ˆ A ( x ) d x + = Z (cid:10) ˆ A ( x ) (cid:11) d x . (3)Clearly, in order for Eq. (3) to hold, the di ff erence h ˆ A i − A mustbe of alternating sign. For our purposes, thus, it is more usefulto define the bias as the di ff erence between the expected meanvalue h ˆ A i and the same quantity that we would theoretically ex-pect in presence of a uniform density of stars , i.e. by ignoringthe dependence of the density of stars on the local extinction de-scribed by Eq. (2) [see also below Eq. (7)]. This definition isuseful, since it isolates the standard e ff ects of a smoothing fromthe e ff ects introduced by the non-uniform sampling of the cloudstructure. In addition, a column density estimator ˆ A ( x ) that is un-biased according to this definition will also be unbiased accord-ing to Eq. (3) (provided the weights w n are spatially invariant).For, when the density of background stars is uniform, (cid:10) ˆ A ( x ) (cid:11) canalways be written as a convolution of the true column densitymap A ( x ) with some kernel K normalized to unity (Lombardi2002).In the following, we will calculate explicitly this bias in twocommon smoothing schemes. In this section we will consider a simple weighting scheme inEq. (1) where each weight w n = w ( x n ; x ) is a simple function ofthe location x n of the n -th star (typically, w ( x n ; x n ) will dependonly on the modulus | x n − x | ).As discussed above, our aim is to compare the expected col-umn density estimate with the one that we would obtain in ab-sence of any selection e ff ect in the number of background stars.We need thus to evaluate two averages: (i) (cid:10) ˆ A ( x ) (cid:11) , where ˆ A is cal-culated according to Eq. (1) and the ensemble average is takenover all positions of stars with the density given by Eq. (2) andover all individual column density measurements ˆ A n ; and (ii)the average ¯ A ( x ), which is basically the same quantity evaluatedwith a constant density of background stars ρ .In principle, both averages can be calculated analytically us-ing the method described by Lombardi & Schneider (2001; seebelow Sect. 3); in practice, for the purposes of this section asimple approximation can be used provided that in the weightedaverage of Eq. (1) a relatively large number of column densitiesare used with weights significantly di ff erent from 0. In this casewe find for the first average (cid:10) ˆ A ( x ) (cid:11) = R A ( x ′ ) w ( x ′ ; x ) ρ ( x ′ ) d x ′ R w ( x ′ ; x ) ρ ( x ′ ) d x ′ . (4)As shown by Eq. (2), ρ ( x ′ ) is a simple function of A ( x ′ ). Thissuggests that the equation can be recast in a simpler form by arco Lombardi: N icest defining p A ( A ; x ), the (weighted) probability that a point withextinction A is used to evaluate ˆ A ( x ), the column density at x : p A ( A ; x ) = R w ( x ′ ; x ) δ (cid:0) A − A ( x ′ ) (cid:1) d x ′ R w ( x ′ ; x ) d x ′ . (5)By using Eq. (5) in Eq. (4) we find then (cid:10) ˆ A ( x ) (cid:11) = R Ap A ( A ; x ) ρ ( A ) d A R p A ( A ; x ) ρ ( A ) d A = R Ap A ( A ; x )10 − α k λ A d A R p A ( A ; x )10 − α k λ A d A . (6)In presence of a uniform distribution of stars, ρ ( x ) = ρ , wewould instead obtain simply¯ A ( x ) = R A ( x ′ ) w ( x ′ ; x ) d x ′ R w ( x ′ ; x ) d x ′ = Z Ap A ( A ; x ) d A . (7)The two values presented in Eqs. (6) and (7) do not di ff er signifi-cantly provided the scatter of A in the patch considered around x is much smaller than ( α k λ ln 10) − . For example, if we considerthe H band of a typical region close to the Galactic plane witha standard reddening law (Rieke & Lebofsky 1985), we have α ≃ .
34 mag − and k H ≃ .
55 (Indebetouw et al. 2005), so thatthe maximum scatter allowed in A K to have a negligible biasis ≪ .
82 mag. Hence, clearly we need to be concerned by thisbias in regions from moderate to large extinction.It is also interesting to evaluate the bias (cid:10) ˆ A ( x ) (cid:11) − ¯ A ( x ) in pres-ence of small variations of A within the region considered. In thiscase, the probability distribution p A ( A ; x ) is peaked around themean value of Eq. (7), and we can expand to second order thetwo exponential functions present in the numerator and denomi-nator of Eq. (6). After a few manipulations we obtain then (cid:10) ˆ A ( x ) (cid:11) − ¯ A ( x ) ≃ − α k λ ln 10 Z p A ( A ; x ) (cid:2) A − ¯ A ( x ) (cid:3) d A = − β R w ( x ′ ; x ) ∆ ( x ′ ; x ) d x ′ R w ( x ′ ; x ) d x ′ , (8)where β ≡ α k λ ln 10 and ∆ ( x ′ ; x ) ≡ A ( x ′ ) − ¯ A ( x ). Hence, thedi ff erence between the two values is proportional to a weightedaverage of the scatter of A around its mean value ¯ A ( x ) at x , aquantity that, as we will see below, can be directly estimatedfrom the data. Note that the bias is, to first order, quadratic on ∆ ( x ′ ; x ) and thus will be particularly severe when steep gradientsare present in the underlying column density A ( x ). Cambr´esy et al. (2002) suggest to use a di ff erent prescription tomake smooth extinction map from the individual column densitymeasurements. For each point in the map, they take the averageextinction of the N angularly closest stars (nearest neighboursinterpolation). As argued by Cambr´esy et al. (2005, 2006), thismethod can potentially alleviate the bias introduced by the vary-ing background density of stars described by Eq. (2), becausemeasurements from low density regions will mostly appear iso-lated and will thus be used for relatively large patches of the sky.Interestingly, using the technique described by Lombardi(2002) it is possible to obtain an analytical estimate for the aver-age of ˆ A in the smoothing scheme proposed: (cid:10) ˆ A ( x ) (cid:11) = Z A ( x ′ ) K ( x ′ ; x ) ρ ( x ′ ) d x ′ , (9) -4 -2 0 2 4 6 x 〈 A ( x ) 〉 N = 1 N = 3 N = 7 N = 15 Fig. 1.
The bias of the nearest neighbors interpolation on a toymodel where the true extinction is A ( x ) = . x < . x ≥
0. The solid lines show the average measuredextinction of the nearest neighbours method, calculated usingEq. (9) for various number of neighbours N . The dashed linesshow the extinction that one would measure if the density ofstars were uniform on the whole field. Note that the solid linesare always below the dashed lines except at the extremes of thisplot (where all estimators agree with the true value of A there).For this plot we set α = . k = .
55, and ρ = K ( x ′ ; x ) is given by K ( x ′ ; x ) = e − µ ( x ′ ; x ) N N − X n = (cid:2) µ ( x ′ ; x ) (cid:3) n n ! . (10)In this equation, N is the number of stars used in the nearestneighbors interpolation and µ ( x ′ ; x ) is the average number ofstars observed in B ( x ′ ; x ), the disk centered in x and of angularradius | x − x ′ | : µ ( x ′ ; x ) = Z B ( x ′ ; x ) ρ ( x ′′ )d x ′′ . (11)As shown by Lombardi (2002), the kernel K ( x ′ ; x ) is normal-ized, i.e. Z K ( x ′ ; x ) ρ ( x ′ ) d x ′ = , (12)an obvious result if one consider the measurement of the extinc-tion on a uniform cloud where A ( x ) = const. In our case, thisproperty can also be proved directly from Eqs. (10) and (11).First, note that µ ( x ′ ; x ) = µ ( r ; x ) depends only on x and on r = | x − x ′ | , and thus the same applies to K ( x ′ ; x ) = K ( r ; x ).This suggests that we can recast the integral of Eq. (12) as Z ∞ K ( r ; x ) ∂µ ( r ; x ) ∂ r d r = N − X n = n ! N Z ∞ e − µ ( r ; x ) (cid:2) µ ( r ; x ) (cid:3) n ∂µ ( r ; x ) ∂ r d r = . (13)The fact that K ( x ′ ; x ) only depends on x and r = | x − x ′ | (andnot on the direction of the vector x − x ′ ) also implies that thenearest neighbors interpolation su ff ers from the bias discussedin this paper. Indeed, in the integral of Eq. (9) the kernel K gives Marco Lombardi: N icest -2 0 2 4 x 〈 A ( x ) 〉 N = 1 N = 3 N = 7 N = 15 Fig. 2.
Similarly to Fig. 1, but for the median-nearest neighborsinterpolation. Again, the solid lines, representing the expectedmean measurements, are always below the dashed lines, repre-senting the same measurements in the ideal case where the den-sity of background stars is uniform within the whole field. Notethat because of the properties of the median, for large values of N the “transition” from A = . A = . x , irrespec-tive of the specific value of ρ on the various points of the circles(what matters here is only the average value of ρ on a circle,and not the specific value on the various points). Hence, we donot expect that this technique can solve the bias introduced bythe correlation between the extinction A and the density of back-ground stars ρ expressed in Eq. (2), especially if steep gradientsare present in the intrinsic cloud column density. We also ex-pect that the bias present in the nearest neighbours interpolationincreases with the number of neighbours N . Indeed, when N in-creases, the “size” of the kernel K ( r ; x ), i.e. the range in r wherethe kernel is significantly di ff erent from 0, also increases becauseof the e ff ects of the polynomial terms in Eq. (10). This, as shownby Eq. (9), implies that the kernel averages out more distant re-gions in the sky, where di ff erences among the local densities canbe significant.In order to better understand the points mentioned in the pre-vious paragraph, it is useful to look at simple example. Considera cloud characterized by a Heaviside column density: A ( x ) = ( A if x < , A if x ≥ . (14)In this case we can basically carry out all the calculations analyt-ically, and obtain for each value of x the expected average mea-sured extinction h ˆ A i . We do not report here the relevant equa-tions, but rather show the results obtained in a typical case inFig. 1. As shown by this plot, the nearest neighbours are clearlybiased: for example, the various curves are not symmetric around x =
0, and in addition h ˆ A i = A α k λ A + A α k λ A α k λ A + α k λ A < A + A . (15)In order to better understand the bias, we also plot in Fig. 1 theexpected results in case of a uniform density, obtained assuming α = x = h ˆ A i = ( A + A ) /
2. A com-parison of the solid and dashed lines also confirms that curves -4 -2 0 2 4 6 x 〈 A ( x ) 〉 σ = 0.5σ = 1.0σ = 2.0 Fig. 3.
Similarly to Fig. 1, but for the moving average smooth-ing with a Gaussian window function with di ff erent dispersionparameters σ . As for Figs. 1 and 2, the solid lines, represent-ing the expected mean measurements, are always below thedashed lines, representing the same measurements in the idealcase where the density of background stars is uniform within thewhole field.with large N , being less steep, are biased on a larger intervalaround x =
0. This is expected, since as shown by Eq. (10) theintrinsic size of the smoothing kernel K increases with N .One might wanders if the bias decreases by using a medianestimator instead of a simple mean, as done by Cambr´esy et al.(2005, 2006). These authors still identify, for each sky position x , the N nearest neighbour stars, and then evaluate a simple me-dian of the relative extinction measurements. Luckily, it is possi-ble to evaluate the exact statistical properties of the median esti-mator (see Lombardi 2002, Appendix A). For this purpose, it isuseful to evaluate the probability distribution p N ( A ; x ) that oneof the N neighbours around the position x has a column density A : p N ( A ; x ) = Z d µ e − µ µ N − N ! Z B ( µ ; x ) d x ′ ρ ( x ′ ) δ (cid:0) A − A ( x ) (cid:1) , (16)where B ( µ ; x ) is the disk centered at x and with radius r definedby µ = µ ( r ; x ) (note that since µ ( r ; x ), for x fixed, is a mono-tonic function of x ′ this definition is unique). The cumulativedistribution associated with p N ( A ; x ) is given by P N ( A ; x ) = Z d µ e − µ µ N − N ! ν ( µ, A ; x ) . (17)In this equation, ν ( µ, A ; x ) is the integral of ρ ( x ) carried out overall points of B ( µ ; x ) where A ( x ′ ) < A : ν ( µ, A ; x ) ≡ Z B ( µ ; x ) ∩{ x ′ | A ( x ′ ) < A } ρ ( x ′ ) d x ′ . (18)Note that ν ( µ, A ; x ) → µ as A → ∞ . Finally, the median of the N = k − p m ( A ) = k Nk ! p N ( A ; x ) (cid:2) P N ( A ; x ) (cid:3) k − (cid:2) − P N ( A ; x ) (cid:3) k − . (19)Although it is di ffi cult to obtain a general expression for thebias introduced by the median-nearest neighbours estimator, itis clear that to first approximation the same arguments discussedabove for the simple mean apply (note also that, trivially, when arco Lombardi: N icest N = N . For comparison, in Fig. 3 we also show the results ob-tained from the simple moving average smoothing discussed inSect. 2.1: interestingly, the plot is very similar to the one for a thesimple average nearest neighbors interpolation (Fig. 1). This, inreality, is not too much surprising because the expected averagevalues (cid:10) A ( x ) (cid:11) measured with these two techniques can be de-scribed in terms of a simple convolution with some appropriatekernels [cf. Eq. (4) and (10); see also Lombardi 2002].In summary, a careful analysis of the nearest neighborsshows that this interpolation technique, for the specific case ofextinction measurements, is still a ff ected by a significant bias.At least for the case of a simple mean estimator, the origin of theproblem can be traced to the di ff erent behaviour of the kernel K ( x ; x ′ ) with respect to the density ρ ( x ′ ) in Eq. (9): in particu-lar, the kernel K does not respond directly to variations of ρ , andinstead depends only µ ( x ′ ; x ), i.e. on the total “mass” within adisk of radius r = h x − x ′ i .
3. N icest : over-weighting high column-densitymeasurements
As shown by Eq. (8), the bias discussed here strongly dependson the small-scale structure of the molecular cloud through theprobability distribution p A ( A ). Although on large scales manymodels predict a log-normal distribution for p A , clearly weshould expect significant deviations from the theoretical ex-pectations on the small scales often investigated in NIR stud-ies. In addition, even at large scales the superposition of dif-ferent cloud complexes, or a significant “thickness” of a cloudalong the line of sight, can produce significant deviations from alog-normal distribution (V´azquez-Semadeni & Garc´ıa 2001; seealso Lombardi et al. 2006 for a case where the log-normal distri-bution is not a good approximation). Hence, we need to be ableto correct for the substructure bias in a model independent way . In order to better understand, and then fix, the bias present inEq. (1) in the case of the moving average smoothing, we use thetheory developed in Lombardi & Schneider (2001, 2002, 2003).The key equations to obtain the ensemble average of ˆ A ( x ) aresummarized below: Q ( s ; x ) = Z (cid:2) e − sw ( x ′ ; x ) − (cid:3) ρ ( x ′ ) d x ′ (20) Y ( s ; x ) = exp (cid:2) Q ( s ; x ) (cid:3) , (21) C ( w ; x ) = − P ( x ) Z ∞ e − ws Y ( s ; x ) d s , (22) (cid:10) ˆ A ( x ) (cid:11) = Z A ( x ′ ) w ( x ′ ; x ) C (cid:0) w ( x ′ ; x ) , x (cid:1) ρ ( x ′ ) d x ′ . (23)Let us briefly comment this set of equations: – The equations are invariant upon the transformation w ( x ′ ; x ) qw ( x ′ ; x ), with q positive constant. This is ex-pected, since an overall multiplicative constant in the weight function does not e ff ect Eq. (1). For simplicity, in the fol-lowing we will assume, without loss of generality, that theweight function is normalized according to Z w ( x ′ ; x ) ρ ( x ′ ) d x ′ = . (24) – The quantity Q ( s ; x ) can be shown to be related to theLaplace transform of p w ( w ; x ), the probability distributionfor the values of the weight function w ( x ′ ; x ) with x fixed. – P ( x ) is the probability that no single star is present within thesupport π w ( x ) of w ( x ′ ; x ), defined as π w ( x ) ≡ { x ′ | w ( x ′ ; x ) > } . Hence, P ( x ) can be simply evaluated from P ( x ) = exp " − Z π w ( x ) ρ ( x ′ )d x ′ . (25)Note that P ( x ) = – C ( w ; x ) is a simple Laplace transform of Y ( s ; x ). This is thekey quantity that enters the final result (23) together withthe original smoothing function w ( x ′ ; x ) and the star density ρ ( x ′ ). – A key parameter in the calculation of C ( w ), and thus on thefinal result h ˆ A ( x ) i , is the so-called weight number of objects N ( x ), defined as N ( x ) ≡ "(cid:2) − P ( x ) (cid:3) Z (cid:2) w ( x ′ ; x ) (cid:3) ρ ( x ′ ) d x ′ − , (26)where as usual w ( x ′ ; x ) is taken to be normalized accordingto Eq. (24). Informally, N ( x ) counts the number of stars thatcontribute (with a weight significantly di ff erent from zero)to the average ˆ A ( x ). – In the limit N ( x ) ≫ C ( w ; x ), which for P ( x ) = C ( w ; x ) ≃ + w + (cid:2) N ( x ) (cid:3) − (1 + w ) . (27)This expression shows that to lowest order C ≃ − w + N − ,where both terms w and N − introduce a relative correctionof order N − in the final result (cid:10) ˆ A ( x ) (cid:11) .The results summarized above rigorously prove the state-ments of Sect. 2.1. In particular, in the limit of a large weightnumber of objects N ( x ), we have C ( w ) ≃
1, and thus Eq. (23)together with the normalization (24) reduces to Eq. (4). In prac-tice, numerical simulations show that Eq. (4) is already accuratefor relatively small weight numbers (of the order of
N ∼
10; seeFig. 3 of Lombardi & Schneider 2001).
The bias of Eq. (1) is essentially due to a change in the densityof extinction measurements as a function of A [Eq. (2)]. In theframework of the moving average smoothing, it is thus reason-able to try to “compensate” this e ff ect by including in the weightsof the estimator (1) a factor 10 α k λ A : w n = w ( x n ; x )10 α k λ ˆ A n . (28)With this modification, we expect ˆ A to be unbiased provided that N is large (so that we can take C ≃
1) and that measurement er-rors on ˆ A n can be neglected. Instead, in presence of significanterrors on the individual column density estimations, we expected Marco Lombardi: N icest -4 -2 0 2 4 6 x 〈 A ( x ) 〉 σ = 0.5σ = 1.0σ = 2.0 Fig. 4.
Similarly to Fig. 3, but for the N icest moving averagesmoothing (see Sect. 3). In this case a significant bias is observedonly for σ = .
5, which however correspond to a very smallweight number
N ≃ A to be biased toward large extinctions. This happens because ofthe non-linearity of the introduced factor in ˆ A n : for example, ˆ A n is symmetrically distributed around A n , 10 α k λ ˆ A n ˆ A n will be skewedagainst values larger than 10 α k λ A n A n . In summary, the modifiedweight of Eq. (28) removes (to a large degree) the bias due to thecloud substructure, but introduces a new bias related to the scat-ter of the individual extinction measurements. While the formeris basically unpredictable, the latter can be estimated accuratelyprovided we know the expected errors on ˆ A n .In order to better understand the properties of the estimator(1) with the new weight, let us evaluate the average (cid:10) ˆ A (cid:11) . Forsimplicity, initially we will ignore measurement errors and focuson the e ff ects introduced by a non-uniform density. We first notethat Eq. (28) is equivalent to the use of a weight w n = w ′ ( x n ; x ),where w ′ ( x ′ ; x ) ≡ w ( x ′ ; x ) /ρ ( x ′ ) . (29)We require the new weight w ′ ( x ′ ; x ) to be normalized accordingto Eq. (24), which implies Z w ( x ′ ; x ) d x ′ = . (30)Interestingly, the normalization condition (30) for w does not de-pend on the density ρ anymore, and it is thus possible to choosefor w a spatially invariant function w ( x ′ ; x ) = w ( x ′ − x ) (a typicalchoice would be a Gaussian, w ∝ e −| x ′ − x | / σ ). We have then (cid:10) ˆ A ( x ) (cid:11) = Z A ( x ′ ) w ′ ( x ′ ; x ) (cid:2) − w ′ ( x ′ ; x ) + N − ( x ) (cid:3) ρ ( x ′ ) d x ′ = Z A ( x ′ ) w ( x ′ ; x ) (cid:2) − w ( x ′ ; x ) /ρ ( x ′ ) + N − ( x ) (cid:3) d x ′ . (31)After a few manipulations we can rewrite this expression in theform (cid:10) ˆ A ( x ) (cid:11) = ¯ A ( x ) − ρ Z w ( x ′ ; x )10 α k λ A ( x ′ ) (cid:2) A ( x ′ ) − ¯ A ( x ) (cid:3) d x ′ ≃ ¯ A ( x ) − α k λ ¯ A ( x ) ρ Z w ( x ′ ; x ) (cid:2) ∆ ( x ′ ; x ) + β ∆ ( x ′ ; x ) (cid:3) d x , (32) where in the second equality we have expanded the exponen-tial term to first order using the definitions following Eq. (8). Insummary, the bias (cid:10) ˆ A ( x ) (cid:11) − ¯ A ( x ) of the proposed estimator iscomposed of two terms: – A weighted average of ∆ ( x ′ ; x ). Since the weighted averageuses w ( x ′ ; x ) as weight, this term is not expected to vanishidentically unless the weight is a top-hat function [cf. Eq. (8),where instead the weighted average involves w ( x ′ ; x ) andtherefore no linear contribution in ∆ is present]. Still, we donot expect any systematic bias related to this term, becauseon average ∆ ( x ′ ; x ) by construction has alternating sign. – A weighted average of ∆ ( x ′ ; x ). This term is manifestly neg-ative and therefore introduces a bias in the result. The bias isvery similar to the original bias discussed in Eq. (8): it isproportional to β ≡ α k λ ln 10, and depends on the scatter ∆ ( x ′ ; x ) ≡ A ( x ′ ) − ¯ A ( x ) of the local column density withrespect to its average value. However, there is a significantdi ff erence: the bias is inversely proportional to ρ − α k λ ¯ A ( x ) ,i.e. to the local average density within the weight function;moreover, the average is taken over w instead of w . Thesetwo di ff erences, together, indicate that the bias of Eq. (32),compared to the on of Eq. (8), is smaller by a factor ∼ N ( x ).In other words, the modification of the weight operated byEq. (28) reduces the bias present in the original simple mov-ing weight average by a factor N ( x ), which in typical cases issignificantly larger than unity. On the other hand, the fact thatthe method is ine ff ective then N is low is evident from the def-inition of w n : in particular, as ρ → w n = w ( x n ; x )10 α k λ ˆ A n significantly di ff erent from zero, but thentrivially ˆ A = ˆ A n and no correction is e ff ectively performed. In Fig. 4 we show the e ff ects of the proposed smoothing tech-nique on the toy-model discussed in Sect. 2.2. A comparisonwith Fig. 3 shows that the N icer interpolation is very e ff ectivein reducing the bias of the simple moving average interpolation,especially for relatively large values of the smoothing factor σ (and thus of the weight number of objects N ).We now consider the e ff ects of measurement errors in themodified estimator. For simplicity, we consider the limit of alarge number of stars ( N ≫ C ( w ) ≃ ffi cient when N isof order of unity). In this case, a simple calculation shows thatmeasurements errors introduce a bias B err of the order B err ≃ β P n w n σ n P n w n , (33)where { σ n } are the variances on { ˆ A n } . In other words, the methodproposed in this paper introduces a small bias toward large ex-tinction values. In order to estimate the order of magnitude of B err , we note that the median error in K band extinctions onN icer column density estimates from the 2MASS catalog istypically σ ≃ .
13 mag, and that less than 1% of stars have σ > . σ n = .
13 mag for all stars in Eq. (33), In reality, a detailed calculation shows that the key parameter hereis N e ff ( x ), which is defined analogously to N ( x ) but with w replacedby w e ff ≡ wC ( w ). As shown in Lombardi (2002), N e ff ≥
1, so the biasnever increases. In this example
N → N e ff →
1, which agrees with state-ment of the previous footnote. The expected error on ˆ A n is calculated using the standard N icer techniques, and thus includes both the typical 2MASS photometric un-certainties and the star color scatters as measured in control fields wherearco Lombardi: N icest we find a bias h ˆ A i − A ≃ .
02 mag (always in the K band).Interestingly, although very small this bias, in contrast to the onedescribed by Eq. (32), can be corrected to first order because itsexpression involves only known quantities.In summary, we propose to replace the estimator (1) withˆ A = P n w n ˆ A n P n w n − β P n w n σ n P n w n , (34)where w n ≡ w ( x ′ ; x )10 α k λ ˆ A n . This new estimator, that we nameN icest , significantly reduces the bias due small-scale inhomo-geneities in the extinction map and is (to first order) unbiasedwith respect to uncertainties on ˆ A n .If one is ready to accept a small bias on ˆ A for extinction mea-surement uncertainties, it is also possible to use Eq. (34) withoutthe last term. In this respect, we also note that the correctingfactor present in Eq. (34) is independent of the local extinction,and only depends on the individual errors on the measurements(these typically are a combination of the photometric uncertain-ties and of the intrinsic scatter of colors of the stellar populationconsidered). As a result, it is not unrealistic to expect that thisfactor is approximately constant within the field of observation,and that it is equally present in the control field used to fix thezero of the extinction.Finally, we note that in N icest a central role is played bythe slope of the number counts α , operationally defined throughEq. (2). In reality, however, it clear that Eq. (2) cannot strictlyhold, because there is an upper limit on the (un-reddened) mag-nitude of stars. Hence, depending on the limiting magnitude ofthe observations, there is an upper limit to the measurable extinc-tion of stars. As a result, in regions where the extinction is largerthan this limit, the density of background stars vanishes, with theresult that in these cases even N icest cannot be unbiased. The ex-act value of the upper measurable extinction depends on manyfactors, including the waveband λ and limiting magnitude m max of the observations, and the distance d of the molecular cloud.An approximate relation for the maximum extinction A max is A max ∼ k λ " m max − M max − d
10 pc ! , (35)where M max is the maximum absolute magnitude of stars in theband considered. The maximum extinction A max should be eval-uated in each specific case; for example, for typical 2MASS ob-servations, using m max ≃
15 mag, k H = .
55 (Indebetouw et al.2005), and M max ≃ − A max ≃
10 mag for a cloudlocated at 100 pc.
Foreground stars do not contribute to the extinction signal, butdo contribute to the noise of the estimators (1) and (34), andthus whenever possible they should be excluded from the analy-sis. Usually, foreground stars are easily identified in high ( A K > f of foreground objects.In principle, if this fraction is known, we could correct the col-umn density estimate ˆ A into ˆ A / (1 − f ): in this way, we would the extinction is negligible. Hence, this error includes the e ff ects of dif-ferent stellar populations on the scatter in the intrinsic star colors, aslong as the various stellar populations are represented in the controlfield. obtain an unbiased estimator of the column density at the priceof an increased noise (due to foreground objects). In practice,the fraction of foreground stars is itself an increasing function ofthe local column density, i.e. f = f ( A ), because of the decreaseddensity of background stars in highly extincted regions. Hence,the simple scheme proposed above is not easily implemented.Interestingly, N icest perfectly adapts to the presence of fore-ground stars. Suppose that in regions with negligible extinction a fraction f ≡ f ( A = < p A ( A = f ( A )of foreground stars increases but, because of the correcting fac-tor in w n , we still expect to measure (1 − f ) A : in other words,the estimator ˆ A / (1 − f ), with ˆ A given by Eq. (34) is expected tobe unbiased both for small-scale inhomogeneities and for fore-ground contamination. Note that correcting factor (1 − f ) − isusually very close to 1 for nearby molecular clouds, and can beeasily evaluated by comparing the density of foreground stars(easily measured in highly extincted regions) with the total den-sity of stars in regions free from extinction.
4. Simulations
In order to test the e ff ects of substructures and the reliability ofthe method proposed here in a real case scenario we performeda series of simulations. We started from a 2MASS / N icer mapof the Pipe nebula (Lombardi et al. 2006), a nearby complexof molecular clouds seen in projection to the Galactic bulge,an ideal case for extinction studies (see also Alves et al. 2008;Onishi et al. 1999). We took this map as the true dust columndensity of a cloud complex, and simulated a set of backgroundstars. In order to show the e ff ects of substructures, we used a verysmall density of background stars (25 stars deg − instead of theoriginal ∼ . × stars deg − used to build the 2MASS / N icer map). This, e ff ectively, correspond to a linear downsize of thestructures of the Pipe nebula by a factor ∼ α = .
34 in the K band and by intrinsic color and color scatters similar to the onesobserved in the 2MASS catalog (see Table 2 of Lombardi 2005for a complete list of the parameters used). We added to each starmagnitude the local value of the of the 2MASS / N icer extinctionand also some random measurement errors:ˆ J n = J n + k K A K ( x n ) + ε ( J ) n , (36)ˆ H n = H n + k H A K ( x n ) + ε ( H ) n , (37)ˆ K n = K + A k ( x n ) + ε ( K ) n , (38)where ε ( J , H , K ) n denote random variables used to model the pho-tometric errors of the stars, and where as usual we used thehat to denote measured quantities. For simplicity, here, we tookthe errors as normal distributed random variables with constantvariance 0 .
05 mag in all bands (the typical median variance of2MASS magnitudes); moreover, we assumed a hard complete-ness limit at 15 mag in all bands (i.e., if a magnitude exceeded 15we took the star as not detected in the corresponding band). Wethen applied the N icer algorithm to these data, thus obtaining,for each star, its measured extinction and related error. Finally,we produced smooth extinction maps by interpolating the indi-vidual extinction measurements with three di ff erent techniques: Marco Lombardi: N icest ° ° ° ° ° ° ° ° ° A K ° ° ° ° ° ° ° ° ° A K ° ° ° ° ° ° ° ° ° -0.10-0.08-0.06-0.04-0.020.00 ∆ A K ° ° ° ° ° ° ° ° ° -0.10-0.08-0.06-0.04-0.020.00 ∆ A K Fig. 5. Top left.
The original 2MASS / N icer map of the Pipe nebula used in the simulations. Top right.
The mean reconstructedmap that one would obtain if the density of background star were uniform. This is basically a smoothed version of the originalmap in the left. Levels are spaced at 0 . Bottom left.
The di ff erence between the average reconstructed maps, calculated frombackground stars following the density of Eq. (2), and the average map at top right. Levels are spaced at 0 .
02 mag.
Bottom right.
Same map as bottom left, but using the improved method presented in this paper. Levels are spaced at 0 . ° ° ° ° ° ° ° ° ° A K ° ° ° ° ° ° ° ° ° -0.10-0.08-0.06-0.04-0.020.00 ∆ A K Fig. 6. Left.
Average reconstructed map using a 1-star ( N =
1) Voronoi method, assuming a uniform distribution of stars. Levels arespaced at 0 . Right. Di ff erence between the average reconstructed Voronoi map, obtained from background stars followingthe density of Eq. (2), and the average true map to the left. Levels are spaced at 0 .
02 mag(1) simple moving average using Eq. (1); (2) modified movingaverage using N icest [Eq. (34)]; (3) nearest neighbour interpo-lation with a single star ( N = ff erent set of random background stars, and took the average maps ˆ A ( x ) obtained with the three techniques. We thencompared these average maps with similar maps obtained fromuniformly distributed stars (in order to produce such maps,we applied the completeness limit to un-extincted magnitudes).Figures 5 and 6 show the results obtained, which are in clear arco Lombardi: N icest ° ° ° ° ° ° ° ° ° -0.12-0.10-0.08-0.06-0.04-0.020.00 ∆ A K ° ° ° ° ° ° ° ° ° -0.12-0.10-0.08-0.06-0.04-0.020.00 ∆ A K Fig. 7.
Same as Fig. 5 (bottom), but with a fraction of f = .
05 of foreground stars. The increases of relative number of foregroundstars in the denser regions of the cloud makes the bias more pronounced. Still, N icest performs very well. Contours levels are spacedat 0 .
02 mag in the N icer plot (left) and at 0 .
005 mag in the N icest plot (right); contours at − . icest . In particular, both the simple moving aver-age and the nearest neighbor su ff er from significant biases, up to A K ≃ . ∼
10, a factor comparable to the weightnumber of stars N (which in our simulations is N ≃
10 in thehigher column density regions).The simulation performed here allows us to evaluate alsothe average quadratic di ff erence between the extinction map ob-tained for the various methods ˆ A ( x ) and the true (smoothed) map¯ A ( x ), i.e. the quantity (cid:10)(cid:2) ˆ A ( x ) − ¯ A ( x ) (cid:3) (cid:11) = Var (cid:2) ˆ A ( x ) (cid:3) + (cid:2)(cid:10) ˆ A ( x ) − ¯ A ( x ) x (cid:11)(cid:3) . (39)As shown by the above equation, the quantity considered herecan be written as the sum of the variance and the square of thebias of the extinction map. Our simulations show that, althoughN icest , as expected, has a slightly larger variance then N icer ,there is still a gain by a factor at least two in the squared di ff er-ence of Eq. (39). In other words, the significant reduction in thebias performed by N icest largely compensate the increase in thescatter introduced by this method. A naive comparison betweenthe results obtained from N icest and from the Voronoi methodwould provide even larger di ff erences in favour of N icest , butwe stress that since the Voronoi method interpolates the extinc-tions over a fixed number of stars (typically smaller than the oneemployed by N icest in our simulations) a direct comparison isnot possible.Figure 7 shows the result of simulation carried out with thepresence of foreground stars. In particular, we generated starsfollowing the same prescriptions described in the previous para-graph, but allowing for a fraction f = .
05 of foreground ob-jects. Because of the e ff ects of extinction, the e ff ective fractionof foreground objects increases in the most dense regions of thecloud, which usually are also the ones with the larger substruc-tures: the two biases thus add up and can be particularly severe.As shown by Fig. 7, N icest is e ff ectively able to cope with rel-atively large fractions of foreground stars while still providing avirtually unbiased column density estimate. ° ′ ° ′ ′ ′ ° ′ ° ′ ′ ′ ° ′ ° ′ ′ ′ ° ′ ° ′ ′ ′ ∆ A K Fig. 9.
The di ff erence between the modified extinction estimates[Eq. (28)] and the standard ones [Eq. (1)] around Barnard 59.The contour levels are at A = { . , . , . , . } mag of the mapof Fig. 8.
5. A sample application
The method presented in this paper was finally applied to thewhole 2MASS point source catalog for the Pipe nebula region.We used the 4.5 million stars of the 2MASS catalog located inthe window − ◦ < l < ◦ , + ◦ < b < + ◦ . (40)The analysis was carried out following the prescriptions ofLombardi & Alves (2001), but using the modified estimator (34)to evaluate ˆ A . In particular, we generated the final extinctionmap, shown in Fig. 8, on a grid of about 1000 ×
750 points, withscale 30 arcsec per pixel, and with Gaussian smoothing charac-terized by FWHM = α = . ± .
02 in the H band.The final, e ff ective density of stars of about 8 stars per pixelguarantees that the approximation used to derive the unbiased es-timator (34) is valid, and that a significant improvement over the icest ° ° ° ° ° ° ° ° ° Galactic Longitude 2 ° ° ° ° ° ° ° G a l ac ti c L a tit ud e ° ° ° ° ° ° ° ° ° ° ° ° ° ° ° ° A K Fig. 8.
The N icest extinction map of the Pipe nebula, using the modified estimator (28). ° ′ ′ ′ ′ ° ′ ′ ′ ′ ° ′ ° ′ ′ ′ ′ ° ′ ′ ′ ′ ° ′ ∆ A K Fig. 10.
The di ff erence between the modified extinction esti-mates [Eq. (28)] and the standard ones [Eq. (1)] around thepeak ID 3 of Lombardi et al. (2006). The contour levels are at A = { . , . , . , . } mag of the map of Fig. 8. standard N icer method can be expected. The largest extinctionwas measured close to Barnard 59, where ˆ A ≃ .
68 mag (a valuethat is 0 .
41 mag larger than what obtained in (Lombardi et al.2006)).Figures 9 and 10 show a comparison between dense regionsmapped using N icer and N icest . We note that, as expected, thetwo methods are equivalent in low-density regions, while thenew one consistently estimates larger column densities as A in-creases. The same e ff ect can be appreciated more quantitativelyfrom Fig. 11, where we plot the relationship between the averageestimates ˆ A obtained in Lombardi et al. (2006) and here.Finally, we argue that the plot of Fig. 11 is strongly related toFig. 9 of Lombardi et al. (2006) where we show the increase inthe scatter of ˆ A n for the di ff erent stars used in each point of theextinction map as a function of the average ˆ A in the point. This A - σ relation is most likely due to unresolved substructures in thecloud, that are expected to be more prominent in high-densityregions. Recently, we reconsidered the A - σ relation and defineda quantity called ∆ ( x ) (Lombardi et al. 2008). This quantity issimply related to the local scatter of measured extinctions, butalso properly takes into account the contribution of measurementerrors to the observed scatter. Interestingly, this quantity can bedefined in terms of simple observables, and can be shown to bedirectly related to the local scatter of column densities: ∆ ( x ) = P n w n ∆ n P n w n , (41) arco Lombardi: N icest A K ∆ A K Fig. 11.
The average di ff erence ∆ A between the estimators (28)and (1) as a function of the estimated column density A . Thelower, solid line is relative to the Pipe nebula map at FWHM = = . = . ∆ n ≡ A ( x n ) − ¯ A ( x ) is the di ff erence between the columnextinction at the position of the n -th star and the local averagecolumn extinction. A comparison of Eq. (41) with Eq. (8) clearlyshows that, except for a numerical factor β , the bias expected inthe standard N icer method is simply proportional to the ∆ ( x ) of(Lombardi et al. 2008).
6. Discussion
The numerical simulations and a first sample application ofN icest have shown that the method presented in this paper cansignificantly alleviate the bias introduced by small-scale struc-tures and by foreground stars in extinction studies. Of course,N icest too has some limitations, which however are largely un-avoidable (and thus inherent to any extinction-based method).First, we note that in order for N icest to work e ff ectively,the weight number of background stars N must be significantlylarger than unity: as shown in Sect. 3.2, N is directly relatedto the reduction in the bias provided by N icest with respect toN icer , and thus having N ∼
N ∼ f of foreground stars is large (e.g., because we are in a particu-larly dense core), on average we do not expect any backgroundstar, and we will thus consistently measure a vanishing extinc-tion. Hence, N icest , like any other extinction based method, onlyworks if there are a sizable number of background stars that canbe used for reliable extinction measurements.The above point might be interpreted as an exceedingly strin-gent requirement for the smoothing window used in N icest . Infact, if a weight function w ( x ′ ; x ) = w ( x ′ − x ) invariant upontranslation is used, this function should be taken broad enoughto guarantee that the weight number of background stars N ( x )is significantly larger than unit everywhere. This, in turn, wouldimply that N ≫ icest technique is still valid for weight functions which are notspatially invariant. Hence, one does not need to use a fixed win- dow size for w ( x ′ ; x ), but rather this function could be taken tochange shape for di ff erent locations x . For example, a simplescheme could be the use of a Gaussian shape for w ( x ′ ; x ) withthe typical scale chosen according to a local estimate of the den-sity of background stars, in a way such that N ( x ) ∼ const ∼ icer technique toreduce the inhomogeneity bias.Another point to keep in mind is that N icest gives a largeweight to red sources. This has two potential unwanted conse-quences. First, it introduces a small bias, that has been correctedto first order in Eq. (34). Second, as shown also by the numericalsimulations described in Sect. 5, it slightly increases the noise inthe final map. Again, since the noise in final map ˆ A ( x ) is propor-tional to 1 / N , the solution is to make sure that N is su ffi cientlylarge. We believe that a small increase in the noise is a fair priceto pay for the significant reduction in the bias provided by N icest over N icer (see also Eadie et al. 1971).A potentially more severe problem can be young stellar ob-jects (YSOs) with infrared excess. These sources can be presentin the most dense star-forming cores of molecular clouds, wherethey can severely bias our method. We note, however, that sincethese sources have peculiar colors and are not usually presentin the control field used for the calibration, they represent any-way a problem for all extinction based methods, including N icer and the nearest neighbours ones. A median estimator is in prin-ciple safer to use in these cases, as long as the YSOs present inthe core are a minority of the total number of background starsobserved in the region; unfortunately, the tendency of YSOs toappear in clusters does not help here (note also that a median isusually more noisy than the simple average). In any case, the ob-vious solution is thus to exclude as much as possible YSOs fromthe list of sources used in the extinction maps.Finally, let us briefly mention alternative possibilities to re-duce the bias considered in this paper. Since the substructurebias is due to a relationship between the local density ρ ( x ) andthe extinction map A ( x ) [Eq. (2)], one could try to use a tracer ofthe local density to perform the correction. We discussed alreadya method based on this concept, the nearest neighbours interpo-lation. However, as shown in Sect. 4, this method fails (see alsoSect. 2.2). In general, a problem with methods based on the localestimate of the density of background stars is that one needs tobe able to obtain ρ ( x ) on scales smaller than the ones that char-acterize the weight w ( x ′ ; x ), or otherwise it is not possible togive di ff erent weights to the di ff erent stars that contribute signif-icantly to the weighted average of Eq. (1). However, if the den-sity of background stars is large enough to have a good handle onthe local density ρ ( x ), one basically does not need to care aboutsubstructures, because it is always possible to make extinctionmaps at the same resolution as the density map. Hence, the onlye ff ective way to deal with substructures involves a use of the lo-cal extinction measurements, as done in N icest . Note also thatany density-based correction will be completely ine ff ective withforeground stars, in contrast to the method presented here.
7. Conclusions
The main results obtained in this paper can be summarized in thefollowing points: – We discussed the e ff ects of small scale inhomogeneities inthe NIR extinction maps based on color excess methods,and showed that large inhomogeneities can significantly biasstandard extinction maps toward low column densities. icest – We proposed a new estimator for ˆ A , N icest , and we showedthat it is (i) equivalent to the usual estimator in low-densityregions, (ii) unbiased and robust for any value of A , i.e.insensitive to the actual form and amount of substructurepresent in the cloud. – We showed that the new estimator is also suitable in presencecontamination by foreground stars. – We tested N icest against numerical simulations, and showedthat it e ff ectively reduces by a large factor the biases dueto substructures and to foreground stars. We also tested analternative approach, the nearest neighbor, and showed thatthe results obtained from this interpolation are still severelybiased. – We applied N icest to the Pipe nebula and showed a few pre-liminary properties of the resulting extinction map. We alsonoted a direct connection between the bias of the N icer andthe ∆ map defined by (Lombardi et al. 2008). Acknowledgements.
We thank J. Alves, C.J. Lada, and G. Bertin for stimulatingdiscussions and suggestions, and the referee, L. Cambresy, for helping us im-prove the paper. This research has made use of the 2MASS archive, provided byNASA / IPAC Infrared Science Archive, which is operated by the Jet PropulsionLaboratory, California Institute of Technology, under contract with the NationalAeronautics and Space Administration.
References
Alves, J., Lada, C. J., & Lada, E. A. 2001, Nature, 409, 159Alves, J., Lombardi, M., & Lada, C. 2008, in Handbook of low mass star forma-tion regions, ed. B. Reipurt, ASP Conference SeriesBok, B. J. 1937, The distribution of the stars in space (Chicago: University ofChicago Press, 1937)Cambr´esy, L., Beichman, C. A., Jarrett, T. H., & Cutri, R. M. 2002, AJ, 123,2559Cambr´esy, L., Jarrett, T. H., & Beichman, C. A. 2005, A&A, 435, 131Cambr´esy, L., Petropoulou, V., Kontizas, M., & Kontizas, E. 2006, A&A, 445,999Dobashi, K., Bernard, J.-P., Hughes, A., et al. 2008, A&A, 484, 205Eadie, W., Drijard, D., James, F., Roos, M., & Sadoulet, B. 1971, StatisticalMethods in Experimental Physics (Amsterdam New-York Oxford: North-Holland Publishing Company)Frerking, M. A., Langer, W. D., & Wilson, R. W. 1982, ApJ, 262, 590Heyer, M. H. & Brunt, C. M. 2004, ApJ, 615, L45Indebetouw, R., Mathis, J. S., Babler, B. L., et al. 2005, ApJ, 619, 931Lada, C. J., Lada, E. A., Clemens, D. P., & Bally, J. 1994, ApJ, 429, 694Larson, R. B. 1981, MNRAS, 194, 809Lombardi, M. 2002, A&A, 395, 733Lombardi, M. 2005, A&A, 438, 169Lombardi, M. & Alves, J. 2001, A&A, 377, 1023Lombardi, M., Alves, J., & Lada, C. J. 2006, A&A, 454, 781Lombardi, M., Alves, J., & Lada, C. J. 2008, A&A in pressLombardi, M. & Schneider, P. 2001, A&A, 373, 359Lombardi, M. & Schneider, P. 2002, A&A, 392, 1153Lombardi, M. & Schneider, P. 2003, A&A, 407, 385Mathis, J. S. 1990, ARA&A, 28, 37Onishi, T., Kawamura, A., Abe, R., et al. 1999, PASJ, 51, 871Rieke, G. H. & Lebofsky, M. J. 1985, ApJ, 288, 618V´azquez-Semadeni, E. & Garc´ıa, N. 2001, ApJ, 557, 727Wilson, R. W., Je ffff