aa r X i v : . [ a s t r o - ph ] J a n To Be Published in The Astrophysical Journal
The Distribution of Stellar Mass in the Pleiades
Joseph M. Converse and Steven W. Stahler
Astronomy Department, University of California, Berkeley, CA 94720 [email protected]
ABSTRACT
As part of an effort to understand the origin of open clusters, we present a statisticalanalysis of the currently observed Pleiades. Starting with a photometric catalog of thecluster, we employ a maximum likelihood technique to determine the mass distributionof its members, including single stars and both components of binary systems. Wefind that the overall binary fraction for unresolved pairs is 68%. Extrapolating toinclude resolved systems, this fraction climbs to about 76%, significantly higher thanthe accepted field-star result. Both figures are sensitive to the cluster age, for whichwe have used the currently favored value of 125 Myr. The primary and secondarymasses within binaries are correlated, in the sense that their ratios are closer to unitythan under the hypothesis of random pairing. We map out the spatial variation ofthe cluster’s projected and three-dimensional mass and number densities. Finally, werevisit the issue of mass segregation in the Pleiades. We find unambiguous evidence ofsegregation, and introduce a new method for quantifying it.
Subject headings: open clusters and associations: individual (Pleiades) — stars: massfunction, statistics — binaries: general
1. Introduction
Open clusters, with their dense central concentrations of stars, are relatively easy to identify.Over a thousand systems are known, and the census is thought to be complete out to 2 kpc (Brown2001; Dias et al. 2002). Because the clusters are no longer buried within interstellar gas and dust,their internal structure and dynamics is also more accessible than for younger groups.Despite these favorable circumstances, many basic questions remain unanswered. Most fun-damentally, how do open clusters form? All observed systems have undergone some degree ofdynamical relaxation. Thus, the present-day distribution of stellar mass differs from the one justafter disruption of the parent cloud. Recovering this initial configuration will clearly be of value inaddressing the formation issue. But such reconstruction presupposes, and indeed requires, that we 2 –gauge accurately the stellar content of present-day clusters. Even this simpler issue is a non-trivialone, as we show here.We consider one of the most intensively studied open clusters, the Pleiades. Again, our ultimategoal is to trace its evolution from the earliest, embedded phase to the present epoch. Focusing hereon the latter, we ask the following questions: What are the actual masses of the member stars?Do they follow the field-star initial mass function? How many of the members are single, and howmany are in binary pairs? Are the primary and secondary masses of binaries correlated? What isthe overall density distribution in the cluster? What is the evidence for mass segregation, and howcan this phenomenon be quantified?All of these questions have been addressed previously by others. Deacon & Hambly (2004)constructed a global mass function for the Pleiades. Their method was to assign masses basedon the observed distribution of R -magnitudes. A more accurate assessment should account forthe photometric influence of binaries. Several studies directed specifically at binaries have probedselected regions for spectroscopic pairs (Raboud & Mermilliod 1998; Bouvier et al. 1997). However,the overall binary fraction has not been carefully assessed, despite some preliminary attempts(Steele & Jameson 1995; Bouvier et al. 1997; Moraux et al. 2003).A fuller investigation of these and the other issues raised requires statistical methods; theseshould prove generally useful in characterizing stellar populations. Section 2 describes our approach,which employs a regularized maximum likelihood technique (e.g. Cowen 1998). A similar methodhas been applied to other astronomical problems, including the reconstruction of cloud shapes(Tassis 2007) and the investigation of binarity within globular clusters (Romani & Weinberg 1991).Our study is the first to apply this versatile tool to young stellar groups. In doing so, we also relaxmany of the restrictive assumptions adopted by previous researchers. Section 3 presents our derivedmass function for the Pleiades, along with our results for binarity. The density distribution, as wellour quantification of mass segregation, are the topics of Section 4. Finally, Section 5 summarizesour findings, critically reexamines the binarity issue, and indicates the next steps in this continuingstudy.
2. Method of Solution2.1. Stellar Mass Probability Function
The basic problem is how to assign stellar masses to all the point-like sources believed to becluster members. In many cases, the source is actually a spatially unresolved binary pair. Morerarely, it is a triple or higher-order system; for simplicity, we ignore this possibility. The availableobservational data consists of photometry in several wavebands. Given the inevitable, random errorassociated with each photometric measurement, it is not possible to identify a unique primary andsecondary mass for each source. Instead we adopt a statistical approach that finds the relative 3 –probability for each source to contain specific combinations of primary and secondary masses.We introduce a stellar mass probability density, to be denoted Φ( m p , m s ). This two-dimensionalfunction is defined so that Φ( m p , m s ) ∆ m p ∆ m s is the probability that a binary system exists withprimary mass (hereafter expressed in solar units) in the interval m p to m p + ∆ m p and secondarymass from m s to m s + ∆ m s . Single stars are viewed here as binaries with m s = 0. We normalizethe function over the full mass range: Z m max m min dm p Z m p dm s Φ( m p , m s ) = 1 . (1)Note that we integrate the secondary mass m s only to m p , its maximum value. Furthermore, weset the lower limit of m s to 0, in order to account for single stars. It is assumed that Φ = 0 for0 < m s < m min . Here, the global minimum mass m min is taken to be 0.08, the brown dwarf limit.We consider m max to be the highest mass still on the main sequence. The age of the Pleiades hasbeen established from lithium dating as 1 . × yr (Stauffer et al. 1998). This figure representsthe main-sequence lifetime for a star of 4 M ⊙ (Siess et al. 2000), which we adopt as m max .We examine separately the handful of stars that are ostensibly more massive than m max , andhence on post-main-sequence tracks. For these 11 sources, we assigned approximate masses fromthe observed spectral types, and obtained data on their known unresolved binary companions; thisinformation was taken from the Bright Star Catalogue (Hoffleit et al. 1991). These systems werethen added by hand to our mass functions. Finally, we ignore the brown dwarf population, thoughtto comprise from 10 to 16% of the total system mass (Pinfield et al. 1998; Schwartz & Becklin2005).The most direct method for evaluating Φ( m p , m s ) would be to guess its values over a discretegrid of m p - and m s -values. Each ( m p , m s ) pair makes a certain contribution to the received flux invarious wavebands. Thus, any guessed Φ( m p , m s ) yields a predicted distribution of fluxes, whichwill generally differ from that observed. One changes the guessed function until the observed fluxdistribution is most likely to be a statistical realization of the predicted one.Unfortunately, this straightforward approach is impractical. The basic difficulty is the mass-sensitivity of stellar luminosities. For secondary masses m s only modestly less than m p , the binaryis indistinguishable photometrically from a single star having the primary mass. A 0 . M ⊙ main-sequence star, for example, has a K -band flux of 10.81 mag at the 133 pc distance of the Pleiades(Soderblom et al. 2005). Pairing this star with a 0 . M ⊙ secondary (which is not yet on the mainsequence) only changes the flux to 10.68 mag. In summary, the function Φ( m p , m s ) evaluated inthis way is unconstrained throughout much of the m p − m s plane. 4 – Since binaries of even modestly large primary to secondary mass ratio are difficult to recog-nize observationally, we need to infer their contribution indirectly, within the context of a largertheoretical framework. The physical origin of binaries is far from settled (Zinnecker & Mathieu2001). There is a growing consensus, however, that most pairs form together, perhaps within asingle dense core. The accumulating observations of protobinaries, i.e., tight pairs of Class 0 orClass I infrared sources, have bolstered this view (e.g., Haisch et al. 2004).If binaries indeed form simultaneously, but independently, within a single dense core, there isno reason to expect a strong correlation between the component masses. (Such a correlation wouldbe expected if the secondaries formed within the primaries’ circumstellar disks, for example.) Acredible starting hypothesis, then, is that each component mass is selected randomly from itsown probability distribution. If the formation mechanism of each star is identical, then thesedistributions are also the same. That is, we postulate that the true binary contribution to Φ( m p , m s )is φ ( m p ) φ ( m s ), where the single-star probability density φ is properly normalized: Z m max m min φ ( m ) dm = 1 . (2)Of course, not all sources are unresolved binaries. Let b represent the fraction of sources thatare. We suppose that this fraction is independent of stellar mass, provided the mass in question canrepresent either the primary or secondary of a pair. While this hypothesis is reasonable for low-massstars, it surely fails for O- and early B-type stars, which have an especially high multiplicity (e.g.,Mason et al. 1998; Garcia & Mermilliod 2001). Such massive objects, however, are not present inclusters of the Pleiades age.Accepting the assumption of a global binary fraction, we have a tentative expression for thefull stellar mass probability:Φ( m p , m s ) = 2 b φ ( m p ) φ ( m s ) + (1 − b ) φ ( m p ) δ ( m s ) . (3)Here, the first term represents true binaries, and the second single stars of mass m p . The factorof 2 multiplying the first term is necessary because of the restricted range of integration for m s inequation (1). That is, this integration effectively covers only half of the m p - m s plane. On the otherhand, the normalization condition of equation (2) applies to both the primary and secondary star,and covers the full range of mass, from m min to m max , for each component.We shall see below that the strict random pairing hypothesis, as expressed in equation (3),does not yield the optimal match between the predicted and observed distribution of magnitudes.The match can be improved, in the statistical sense outlined previously, if one allows for a limiteddegree of correlation between the primary and secondary masses within binaries. In other words,there is an apparent tendency for more massive primaries to be accompanied by secondaries thathave a greater mass than would result from random selection. 5 –A simple way to quantify this effect is to consider the extreme cases. If there were perfect cor-relation between primary and secondary masses, then the contribution to Φ( m p , m s ) from binarieswould be b φ ( m p ) δ ( m p − m s ). With no correlation at all, Φ( m p , m s ) is given by equation (3). Weaccordingly define a correlation coefficient c , whose value lies between 0 and 1. Our final expressionfor Φ( m p , m s ) uses c to define a weighted average of the two extreme cases:Φ( m p , m s ) = 2 b (1 − c ) φ ( m p ) φ ( m s ) + b c φ ( m p ) δ ( m p − m s ) + (1 − b ) φ ( m p ) δ ( m s ) . (4)Note that the last righthand term, representing the probability of the source being a single star, isunaffected by the degree of mass correlation within binaries. Reconstructing the stellar mass probability Φ( m p , m s ) requires that we evaluate the constants b and c , as well as the single-star probability φ ( m ). To deal with this continuous function, wedivide the full mass range into discrete bins of width ∆ m i . Integrating over each bin, we find y i ,the probability of a star’s mass being in that narrow interval: y i ≡ Z m i +∆ m i m i φ ( m ) dm . (5)We symbolize the full array of y i -values by the vector y , and similarly denote other arrays below.Our task, then, is to find optimal values not only for b and c , but also for all but one element of y .The normalization of φ is now expressed by the constraint X i y i = 1 , (6)which sets the last y -value.For each choice of b , c , and y , equation (4) tells us the relative probability of binaries being atany location in the m p − m s plane. After dividing the plane into discrete bins, each labeled by anindex α , we define µ α as the predicted number of systems associated with a small bin centered onan ( m p , m s ) pair. If µ tot is the total number of systems, i.e., of unresolved sources in all magnitudebins, then our chosen b , c , and y i -values yield the relative fractions µ /µ tot .As an example, consider a bin α in which m p and m s have different values lying between m min and m max . Then the system is an unequal-mass binary, for which equation (4) gives µ α µ tot = 2 b (1 − c ) y p y s . (7)Here, y p is the element of y corresponding to the selected m p , while y s is similarly associated with m s . For a bin where m p = m s , the corresponding relation is µ α µ tot = b (1 − c ) y p + b c y p . (8) 6 –Note the additional term accounting for correlated binaries. Note also that a factor of 2 has beendropped from equation (4), since we are integrating only over that portion of the mass bin with m s < m p . Finally, if the system is a single star, so that m s = 0, we have µ α µ tot = (1 − b ) y p . (9)Our observational data consists of a catalog of n tot sources, each of which has an apparentmagnitude in at least two broadband filters. (In practice, these will be the I - and K -bands; seebelow.) As before, we divide this two-dimensional magnitude space into small bins. Our choice of b , c , and y leads not only to a predicted distribution in mass space, but also in magnitude space.Let ν β be the predicted number of sources in each magnitude bin, now labeled by the index β .Then we may write the transformation from the mass to the magnitude distribution as ν β = X α R βα µ α , (10)which may be recast in the abbreviated form ν = R µ . (11)Here, R is the response matrix , whose elements R αβ give the probability that a source in a massbin α is observed in a magnitude bin β . In detail, this probability utilizes a theoretical isochronein the color-magnitude diagram (see Section 3.1). We must also account for random errors in themeasured photometry. In other words, a given magnitude pair can have contributions from a rangeof mass pairs. It is for this reason that each element of ν involves a sum over all α -values.We previously showed how to obtain the relative mass distribution µ /µ tot , not the actual µ itself. However, it is the latter that we need for equation (11). To find µ tot , we sum equation (10)over all β -values, and demand that this sum be n tot , the total number of observed sources: n tot = X β ν β = X β X α R βα µ α , so that n tot = µ tot X β X α R βα (cid:18) µ α µ tot (cid:19) . (12)In summary, choosing b , c , and y gives us µ /µ tot through through equations (7), (8) and (9). Wethen solve equation (12) for µ tot . Supplied with knowledge of µ , we finally use equation (11) tocompute ν . Having chosen b , c , and y , how do we adjust these so that the predicted and observed magni-tude distributions best match? Our technical treatment here closely follows that in Cowen (1998, 7 –Chapter 11), but specialized to our particular application. Let n β be the number of sources actuallyobserved in each two-dimensional magnitude bin. We first seek the probability that the full array n is a statistical realization of the predicted ν . The supposition is that each element ν β representsthe average number of sources in the appropriate bin. This average would be attained after a largenumber of random realizations of the underlying distribution. If individual observed values followa Poisson distribution about the mean, then the probability of observing n β sources is P ( n β ) = ν βn β e − ν β n β ! . (13)This probability is highest when n β is close to ν β .The likelihood function L is the total probability of observing the full set of n β -values: L ≡ Y β P ( n β ) . (14)We will find it more convenient to deal with a sum rather than a product. Thus, we useln L = X β n β ln ν β − ν β − ln ( n β !) . (15)The strategy is then to find, for a given n , that ν which maximizes ln L . For this purpose, wemay neglect the third term in the sum, which does not depend on ν . We thus maximize a slightlymodified function: ln L ′ ≡ X β n β ln ν β − ν β . (16)Since, for a given n , each P ( n β ) peaks at ν β = n β , maximizing ln L ′ is equivalent to setting ν equal to n in equation (11), and then inverting the response matrix to obtain µ . Such a directinversion procedure typically yields a very noisy µ , including unphysical (negative) elements. Thesolution is to regularize our result by employing an entropy term S : S ≡ − X i y i ln ( y i ) . (17)The function S is largest when the elements of y are evenly spread out. Adding S to ln L ′ andmaximizing the total guarantees that the y -values are smoothly distributed, i.e., that φ ( m ) is alsoa smooth function.In practice, we also want to vary the relative weighting of S and ln L ′ . We do this by defininga regularization parameter λ , and then maximizing the function Γ, whereΓ ≡ λ ln L ′ + S . (18)For any given value of λ , maximizing Γ yields an acceptably smooth y that reproduces well theobserved data. For the optimal solution, we find that value of λ which gives the best balancebetween smoothness of the derived φ ( m ) and accuracy of fit. We do this by considering anotherstatistical measure, the bias. 8 – Our observational dataset, n , is an imperfect representation of the unknown probability density φ ( m ) in two senses. As already noted, n may be regarded as only one particular realization of theunderlying distribution. Even this single realization would directly reveal φ ( m ) (or, equivalently, y ) if the sample size were infinite, which of course it is not.Imagine that there were other realizations of φ ( m ). For each, we employ our maximum likeli-hood technique to obtain y . Averaging each y i over many trials yields its expectation value, E ( y i ).However, because of the finite sample size, E ( y i ) does not necessarily converge to the true value, y true i . Their difference is the bias, b i : b i ≡ E ( y i ) − y true i . (19)The values of the biases, collectively denoted b , reflect the sensitivity of the estimated y tochanges in n . Following Cowen (1998, Section 11.6), we define a matrix C with elements C iβ ≡ ∂y i ∂n β . (20)The bias is then given by b = C ( ν − n ) . (21)To evaluate the derivatives in C , we consider variations of the function Γ about its maximum. Inmatrix notation, C = −A − B , (22)where the matrix A has elements A ij ≡ ∂ Γ ∂y i ∂y j , (23)and B is given by B iβ ≡ ∂ Γ ∂y i ∂n β . (24)Since Γ is a known function of both the y and n , the derivatives appearing in both A and B may be evaluated analytically. Another matrix that will be useful shortly is D , whose elements D αi ≡ ∂µ α ∂y i (25)are also known analytically from equations (7)-(9).To determine the regularization parameter λ appearing in Γ, we seek to minimize the biases.In practice, we consider the weighted sum of their squared values: χ b ≡ X i b i W ii , (26) 9 –and vary λ to reduce this quantity. Here, the W ii are diagonal elements of W , the covariance ofthe biases. Recall that the elements of b are here considered to be random variables that changewith different realizations.We find the covariance matrix W by repeated application of the rule for error propagation(Cowen 1998, Section 1.6). We begin with V , the covariance of n . Since these values are assumedto be independently, Poisson-distributed variables, V has elements V αβ = ν β δ αβ . (27)Here we have used the fact that the variance of the Poisson distribution equals its mean, ν β . Inthis equation only, both α and β range over the ν -values, i.e., V is a square matrix.We next consider Y , the covariance of y . This is given by Y = C V C T . (28)Finally, we obtain the desired W by W = F Y F T . (29)The matrix F in this last equation has elements F ij ≡ ∂b i ∂y j . (30)Differentiating equation (21) and applying the chain rule, we find F to be F = C R D − I , (31)where I is the identity matrix. Thus far, we have focused on determining global properties of the cluster, especially the massfunction φ ( m ). We also want to investigate the spatial distribution of stellar masses. For thispurpose, we need not perform another maximum likelihood analysis. The reason is that we cantreat the mass distribution at each radius as a modification of the global result.We divide the (projected) cluster into circular annuli, each centered on a radius r . What is µ rα ,the number of sources in an annulus that are within mass bin α ? (As before, each bin is labeledby the masses of both binary components.) The quantity we seek is µ rα = X β Q αβ ν rβ . (32) In practice, we require that χ b be reduced to N , the number of free parameters in our fit. As noted by Cowen(1998, Section 11.7), the average b i -value at this point is about equal to its standard deviation, and so is statisticallyindistinguishable from zero.
10 –Here, Q αβ is the probability that a source observed within magnitude bin β has component masseswithin bin α . In principle, this probability depends on radius. For example, the source could havea high probability of having a certain mass if there exist many such stars in that annulus, evenstars whose real magnitude is far from the observed one. If we discount such extreme variations inthe underlying stellar population, then we may approximate Q by its global average.The factor ν rβ in equation (32) is the estimated number of sources at radius r in magnitude bin β . We only compute, via the maximum likelihood analysis, ν β , the estimated number of sourcesin the entire cluster. But we also know n rβ , the observed number of sources in the annulus. If thetotal number of observed sources, n β , is non-zero, then we take ν rβ ≡ n rβ n β ν β . (33)In case n β = 0, then n rβ vanishes at all radii. We then assume ν rβ ≡ = n r tot n tot ν β , (34)where n r tot is the observed source number of all magnitudes in the annulus. In the end, equation (32)attributes the radial mass variation to changes in the local magnitude distribution, rather than toimprobable observations of special objects.It is clear that the global Q must be closely related to the response matrix R , which is theprobability that a source with a given mass has a certain magnitude. The precise relation betweenthe two follows from Bayes Theorem: Q αβ = R βα µ α /µ tot ν β /ν tot . (35)The numerator of the fraction is the probability that a source at any radius lies within the massbin α , while the denominator is the probability of it lying within magnitude bin β . Note that µ tot and ν tot need not be identical. The first quantity is the estimated number of sources covering allpossible masses. The second is the observed number in the magnitude range under consideration.In practice, this range is extensive enough that the two are nearly the same. We thus write Q αβ = R βα µ α ν β . (36)Using this last equation, along with equation (33), equation (32) now becomes µ rα = µ α X β R βα n rβ n β . (37)For those terms where n β = 0, equation (34) tells us to replace the ratio n rβ /n β by n r tot /n tot . 11 –Summing µ rα over all α and dividing by the area of the annulus gives the projected surfacenumber density of sources as a function of radius. The total mass in each annulus is∆ m r = X α µ rα m α , (38)where m α is the sum of the masses of both binary components in the appropriate bin. Division of∆ m r by the annulus area gives the projected surface mass density. Under the assumption of spher-ical symmetry, the corresponding volume densities are then found by the standard transformationof the Abel integral (Binney & Merrifield 1998, Section 4.2.3).
3. Application to the Pleiades: Global Results3.1. The Response Matrix
We begin with the observational data. Figure 1 is a dereddened (
I, I − K ) color-magnitudediagram for the Pleiades, taken from the recent compilation by Stauffer et al. (2007). Shown areall sources which have high membership probability, as gauged by their colors, radial velocities,and proper motions (see, e.g., Deacon & Hambly (2004) for one such proper motion study.) Thelower open circles correspond to probable brown dwarfs; we exclude such objects from our study.Most brown dwarfs are too faint to be observed, and the population, in any case, is more sparselysampled. (The magnitude cutoffs corresponding to a 0 . M ⊙ object are M I = 12 and M K = 9.)After also exluding the 11 bright, post-main-sequence stars, shown here as large, filled circles, wehave a total sample size of n tot = 1245.The solid curve near the lower boundary of the stellar distribution is a combination of thetheoretical zero-age main sequence for m ∗ > The isochrone is that for the measured cluster ageof 125 Myr. Our basic assumption is that the observed scatter about this curve stems from twoeffects - binarity and intrinsic errors in the photometric measurements. We do not consider, fornow, possible uncertainty in the cluster’s age. (See Section 5 for the effect of this uncertainty.) Wealso ignore the finite duration of star formation. This duration is roughly 10 yr (Palla & Stahler2000), or about 10 % of the cluster age. Both theoretical results are presented in magnitudes. We have applied corrections to the theoretical K -bandmagnitudes to make them consistent with the 2MASS K s -band used in Stauffer’s catalog. See Cohen et al. (2003)for this transformation. We further ignore the effect of differential reddening across the cluster. Stauffer et al. (2007) adjusted individuallythe fluxes from sources in especially obscured regions, bringing their effective extinction to the observed average A V of 0.12. We therefore constructed Figure 1 by applying uniformly the corresponding A I - and A K -values of 0.06 and0.01, respectively.
12 –After doing a polynomial fit to the mass-magnitude relations found by Siess et al. (2000) andBaraffe et al (1998), we have analytic expressions for M ∗ I ( m p , m s ) and M ∗ K ( m p , m s ), the absolute I - and K -magnitudes for a binary consisting of a primary mass m p and secondary m s . Here,the superscripts indicate that the magnitudes are theoretically derived. Both M ∗ I and M ∗ K arecalculated by appropriately combining the individual absolute magnitudes for m p and m s .We do not actually observe M ∗ I or M ∗ K for any source. What we have are dereddened, apparentmagnitudes in these wavebands. Using the Pleiades distance of 133 pc, these apparent magnitudesare readily converted to absolute ones, M I and M K . The salient question is: Given a source withintrinsic magnitudes M ∗ I and M ∗ K (or, equivalently, with masses m p and m s ), what is the probabilitythat it is observed to have magnitudes M I and M K ?Here we confront the issue of photometric errors. We assume the errors in the two wavebandsto be normally distributed. Then the relevant probability density is S ( M I , M K ; m p , m s ) = 12 π σ I σ K exp " − ( M I − M ∗ I ) σ I − ( M K − M ∗ K ) σ K . (39)Here, S ∆ M I ∆ M K is the probability of observing a source in magnitude bin β , centered on( M I , M K ), and having widths ∆ M I and ∆ M K .The quantities σ I and σ K in equation (39) are the standard deviations of the photometricmeasurements. According to Stauffer et al. (2007), the average standard deviation in the I -bandis about 0.15. Figure 2, constructed from Table 2 of Stauffer et al. (2007), shows that σ K isgenerally lower, and rises steeply with M K for the dimmest sources. The two branches of thecurve presumably represent the results from two different observations. We do a polynomial fit tothe upper, majority, branch, and thus have an explicit expression for σ K ( M K ).Suppose now that m p and m s are centered within a mass bin α , which has widths ∆ m p and ∆ m s . Then the response matrix R αβ is obtained by integrating S ( M I , M K ; m p , m s ) over themagnitude bin, then averaging over the mass bin: R αβ = R m p +∆ m p m p dm p R m s +∆ m s m s dm s R M I +∆ M I M I dM I R M K +∆ M K M K dM K S ∆ m p ∆ m s . (40)The magnitude integrals can be expressed in terms of error functions if we reinterpret σ K as being afunction of M ∗ K rather than M K . The remaining numerical integrals over m p and m s are performedby finding, for each ( m p , m s ) pair, the M ∗ I - and M ∗ K -values from our polynomial fits to the mass-magnitude relations. This rise in σ K occurs because the observed K -magnitudes are approaching the sensitivity limit of the observa-tions. Many of the I -band measurements come from POSS II plates, for which the limit is 18.5 (Hambly et al. 1993).Another large source of data was the observations of Pinfield et al. (2000), whose limiting magnitude was 19.7. Ourlower cutoff for brown dwarfs corresponds to an apparent I -magnitude of 17.7, so the rise in our σ I should be modest.
13 –
With the response matrix in hand, we are ready to input the catalog of source magnitudes.Before turning to the Pleiades itself, we first employed a number of synthetic datasets, in order totest various aspects of the code. We shall describe these tests shortly. First, however, we summarizethe standard procedure we adopted for the analysis of any cluster, real or synthetic.The basic problem, we recall, is to guess the optimal values of b , c , and y that maximize thefunction Γ, as given in equation (18). The entropy part of Γ, labeled S , is directly a function of y (eq. (17)), while the modified likelihood function L ′ depends on the observed magnitude distribution n and the guessed one ν (see eq. (16). The guessed y does not yield ν itself, but the guessed massdistribution µ , through equations (7)-(9). It is in this transformation that the binary fraction b andcorrelation coefficient c appear. Finally, ν is obtained from µ via the response matrix (eq. (11)).We begin the maximization procedure by first setting the regularization parameter λ to zero.Since Γ = S in this case, the optimal set of y -values will be uniformly distributed, while subject tothe normalization constraint of equation (6). We guess b , c , and y , and vary them to maximize Γ.For the actual maximization, we employ a standard simplex algorithm (Press et al. 2002, Chapter10). The resulting best-fit parameters are then perturbed and the maximization rerun. This check,which may be redone several times, is done both to confirm convergence and to avoid becomingtrapped in small, local maxima of the function Γ. We record the relevant covariances and biases,to be used in estimating errors in predicted quantities and to set the optimal λ -value.The next step is to increase λ slightly. We maximize Γ in the same way as before, againrecording covariances and biases. We again increase λ , repeating the entire procedure until χ b , theweighted sum of the biases, starts to become acceptably small. At this point, the best-fit b , c , and y have been established.As a first test of the procedure, we introduced an artificial cluster whose single-star probability, φ ( m ), we selected beforehand. Sources were chosen randomly to have masses according to thisdistribution. In a certain fraction of the sources, our preset binary fraction b , a second star wasadded to the source. The mass of this object was also randomly chosen from φ ( m ). Thus, thecorrelation coefficient c was initially zero. Given both masses in a source, its intrinsic M ∗ I and M ∗ K are readily obtained. These magnitudes are smeared out into neighboring bins accordingto Gaussian errors with the appropriate standard deviations σ I and σ K . Thus, the “observed”magnitude distribution, n , is established.Figure 3 shows two representative examples. In the left panel, the chosen φ ( m ), is a powerlaw: φ ( m ) ∝ m − . . On the right, we used a log-normal distribution: φ ( m ) = Cm exp (cid:20) − (log m − log m ◦ ) σ m (cid:21) , (41)where C is the normalization constant. The central mass was chosen as m ◦ = 0 . σ m was 0.4. The binary fraction b was chosen to be 0.30 in the power-law example, 14 –and 0.68 for the log-normal distribution. The total source number n tot was 10,000 in both cases.The smooth curve in both panels is φ ( m ), while the data points are the best-fit values of y i / ∆ m i , where ∆ m i is the bin width. Shown also is the estimated error for each value. This wasderived from the covariance matrix Y , introduced in Section (2.1). Specifically, the plotted erroris √Y ii / ∆ m i . We divide each y i and its associated error by ∆ m i because y i is integrated over thebin (eq. (5)).It is evident that the code reproduces well the assumed φ ( m ) in these two examples. Notethat most of the scatter seen in both plots, especially in the left panel, was already present in theinput data, which were finite realizations of the analytic distributions. The derived (i.e., predicted)binary fractions, b = 0 . ± .
008 and b = 0 . ± . n tot to 1245, the actual number in our Pleiades sourcecatalog. In this smaller sample, the errors in our predicted mass function and binary fractionincreased, roughly as n tot − / .Figure 4, taken from a dataset with n tot = 1245, shows in more detail how the regularizationparameter λ was chosen. The figure also illustrates some of the subtlety involved in this procedure.Plotted here, as a function of λ , is χ b , defined in equation (26). As λ is gradually increased, χ b takes a sudden, sharp dip. After climbing back, χ b then more slowly declines, eventually fallingbelow N = 21, the number of tunable parameters in this maximization ( b , c , and 19 y -values).It is the second threshold ( λ = 0 .
027 in this case) that marks the true stopping point. Theearlier dip in χ b is due, not to a decrease in the biases, but to a sharp increase in the covariances W . This increase commonly occurs when the likelihood term ln L ′ starts to become comparable tothe entropy S in the full function Γ. At that point, y makes an abrupt shift away from its earlier,nearly uniform, distribution. With further increase in λ , y settles down gradually to its optimalform.Continuing our synthetic data tests, we next introduced a correlation between the primary andsecondary masses. First, we generated uncorrelated pairs, as above. Generalizing the prescriptionof de La Fuente Marcos (1996), we then altered the secondary mass in each source according to m s → m s (cid:18) m p m s (cid:19) γ . (42)Here, γ , a preset number between 0 and 1, represents our imposed degree of correlation. Thus,setting γ = 0 yields the previous, uncorrelated case, while, for γ = 1, every binary has equal-masscomponents. We ran our routine with a variety of input single-star mass functions, binary fractions,and degrees of correlation.Our general result was that the predicted y still reproduced well the synthetic φ ( m ). Thebinary fraction b was similarly accurate. Most significantly, the predicted c -value tracked theinput quantity γ . Figure 5 shows this relation. We conclude that our statistical model, whilecrudely accounting for correlation by inserting a fraction of equal-mass pairs, nevertheless mimics 15 –a smoother correlation, such as would be found naturally. The shaded patch in the figure is theprobable region occupied by the real Pleiades; Section 3.4 below justifies this assessment. One price we paid for our simplified account of correlation was that our matching of φ ( m )was less accurate than for randomly paired input binaries. Consider, for example, the log-normalfunction of equation (41). While our best-fit y still reproduced φ ( m ) reasonably well, the outputfunction peaked at too high a mass compared to m ◦ . The filled circles in Figure 6 shows thatthis shift, ∆ m ◦ , increased with the input γ -value. Concurrently, our output function was toonarrow compared to the input σ m . The (negative) difference, ∆ σ m , displayed as open circles inFigure 6, was also more pronounced at higher γ . These systematic errors need to be consideredwhen analyzing a real cluster. The two patched areas in the figure again represent the likely regimeof the Pleiades, as we explain shortly. We now present the results of applying our maximum likelihood analysis to the Pleiades itself,i.e., to the I - and K -magnitudes of 1245 sources from the catalog of Stauffer et al. (2007). Ourbest-fit binary fraction was b = 0 . ± .
02, while the correlation coefficient was c = 0 . ± . y i / ∆ m i . As in Figure 3, these pointsare a discrete representation of the single-star mass function φ ( m ). The large error bars on the twopoints at highest mass are due to the small number of sources gauged to be in the respective bins.The smooth, solid curve in Figure 7 is a log-normal mass function that best matches the empirical y . Referring again to equation (41), we find that m ◦ = 0 . ± .
04 and σ m = 0 . ± .
02. Thepresence of a finite binary correlation affects both estimates. Judging from Figure 6, our m ◦ isoverestimated by about 0.06, while σ m should be raised by 0.08.Each of our mass bins has contributions from both the primary and secondary components ofbinary pairs, as quantified by equation (4). Integrating the full stellar mass probablity Φ( m p , m s )over all secondary masses, we obtain φ p ( m p ), the probability distribution of primary masses: φ p ( m p ) = Z m p dm s Φ( m p , m s ) . (43)Note that this distribution includes the possibility that the star is single ( m s = 0). The prescription for mass correlation given in equation (42) is more realistic than introducing a subpopulationof identical-mass binaries (eq. (4)). We employed the latter device only for convenience. If we had parametrized thecorrelation through γ , equations (7) and (8) would have been numerical integrals, and the derivative matrix D inequation (25) would also have required numerical evaluation.
16 –The solid curve in Figure 8 is a log-normal fit to the empirical φ p ( m p ). Shown for comparisonas a dashed curve is the fit for φ ( m ) from Figure 7. Relative to the latter function, the primarydistribution falls off at lower masses. This falloff simply reflects the fact that less massive objectsare more likely to be part of a binary containing a higher-mass star, and thus to be labeled as“secondaries.” In any event, we now see why the peak of φ p ( m p ), m ◦ = 0 . ± .
02, is elevatedwith respect to the peak of φ ( m ). Similarly, the primary distribution is also slightly narrower thanthe single-star mass function, with σ m = 0 . ± . φ p ( m p ) may be compared to those ofMoraux et al. (2004). These authors fit the entire mass function. Since, however, they did notaccount for binarity, their results are more closely analogous to our primary distribution. Theirbest-fit m ◦ of 0.25 is close to ours, while their σ m of 0.52 is higher, mostly because of their inclusionof the highest-mass members. These parameters are also close to those given by Chabrier (2003)in his log-normal fit to the field-star initial mass function ( m ◦ = 0 . σ m = 0 . φ ( m ) to the field-star initial mass function(dashed curve). The latter, which has been raised in the figure for clarity, is taken from Kroupa(2001), who did correct for binarity. It is apparent that φ ( m ) itself veers away from the IMF forboth low- and high-mass objects. When these are added in, the resemblance improves. The opencircles in Figure 9 are Pleiades low-mass stars and brown dwarfs found by Bihain et al (2006). Wehave normalized their data, taken from a limited area of the cluster, so that their total number ofstars matches ours within the overlapping mass range. No such normalization was necessary for the11 B-type stars (filled circles), which are from the catalog of Stauffer et al. (2007) but not includedin our maximum likelihood analysis. Adding both these groups not only improves the match to theIMF, but also reveals a gap in the stellar distribution between about 2 and 5 M ⊙ . A similar gap isseen in the Pleiades mass function of Moraux et al. (2004, see their Figure 1).Our estimate for the total cluster mass, based solely on the 1245 catalog sources, is 738 M ⊙ ,with a 4% uncertainty. Adding in the brightest stars brings this total to 820 M ⊙ , with the samerelative error. Tests with synthetic data indicate that the systematic bias due to binary correlationraises this figure by roughly 50 M ⊙ , to 870 M ⊙ . Addition of the brown dwarfs would cause afurther, relatively small, increase. For comparison, Pinfield et al. (1998) found 735 M ⊙ in stars,and an upper limit of 131 M ⊙ for the brown dwarf contribution. Raboud & Mermilliod (1998) usedthe virial theorem to estimate a total mass of 720 M ⊙ , with a 28% uncertainty. Direct integrationof their mass function gave 950 M ⊙ , with an 18% fractional error. We may similarly calculate a secondary mass distribution φ s ( m s ) by integrating Φ( m p , m s ) over m p , from m s to m max . The function φ s has an excess of low-mass stars and drops very steeply at high masses, as most such objectsare primaries.
17 –
The global binary fraction, b = 0 .
68, obtained in our analysis represents most, but not all, ofthe full binary population. Omitted here are spatially resolved systems. For these, the primary andsecondary appear as separate sources in the catalog of Stauffer et al. (1998). Counting resolvedpairs raises the total fraction to about 76%, as we now show.The smallest angular separation between stars in the catalog is 10 ′′ . At the Pleiades distanceof 133 pc, the corresponding physical separation is 1400 AU. An edge-on circular binary of exactlythis orbital diameter will still be unresolved, since the components spend most of their time closertogether. The true minimum separation in this case is 2200 AU. Here, we have divided 1400 AUby 2 /π , which is the average of | sin θ | , for θ randomly distributed between 0 and 2 π .Of course, only a relatively small fraction of binaries have separations exceeding 2200 AU.The average total mass of our unresolved systems is 0 . M ⊙ . A binary of that total mass and a2200 AU diameter has a period of 1 . × yr. What fraction of binaries have even longer periods?Our average primary mass is 0 . M ⊙ , corresponding to a spectral type of M1. Fischer & Marcy(1992) studied the period distribution of binaries containing M-type primaries. They claimed thatthis distribution was indistinguishable from that found by Duquennoy & Mayor (1991) for G-typeprimaries. In this latter sample, 11% of the systems had periods greater than our limiting value.If the Pleiades periods are similarly distributed, then the total fraction of binaries - both resolvedand unresolved - becomes 0 . / (1 − .
11) = 0 . . × yr. Using the period distribution of Duquennoy & Mayor (1991)to extrapolate their observed binary fraction of 28% yields a total fraction of 60%. Mermilliod et al.(1992) observed spectroscopic pairs with periods under 3 yr. A similar exercise again yields 60%.We note, however, that this ostensible concurrence of results is based on very broad extrapolationsfrom limited data. (See Fig. 4 of Bouvier et al. (1997).)Our derived binary fraction also exceeds that found in the field-star population. Duquennoy & Mayor(1991) found that 57% of G stars are the primaries of binary or higher-order systems. Note thatour b represents the total probability that a star is in a binary, whether as the primary or secondarycomponent. Since G stars are rarely secondaries, the comparison with Duquennoy & Mayor (1991)is appropriate. On the other hand, M stars are frequently secondaries, so we would expect thefraction of binaries with M-type primaries to be reduced. Lada (2006) has found that only 25% ofM-stars are the primary components of binaries. Our own analysis yields a binary fraction of 45%for M-star primaries, still in excess of the field-star result.If our finding of a relatively high binary fraction proves robust, it may provide a clue to theprogenitor state of the Pleiades and other open clusters. A similar statement applies to the corre-lation between component masses within binaries. Our adopted method of gauging this correlation 18 –- inserting a fraction of equal-mass pairs in the mass function - is admittedly crude. Nevertheless,the strong result ( c = 0 . ± .
06) is significant. Referring back to Figure 5, we find that thePleiades correlation is equivalent to setting γ equal to about 0.65 in the alternative description ofequation (42). Whatever the origin of the Pleiades binaries, the primaries and secondaries werenot formed by completely independent processes.
4. Application to the Pleiades: Radial Distributions4.1. Number and Mass Profiles
We now employ the procedure outlined in Section 2.4 to investigate both the surface andvolumetric density as a function of the projected distance from the cluster center. The filled circlesin Figure 10, along with the associated error bars, represent the surface number density of sources,measured in pc − . The solid curve is a density profile using the empirical prescription of King(1962). Here, the core radius is 2.1 pc, while the tidal radius is 19 pc. For comparison, Adams et al.(2001) also fit the surface number density profile of their low-mass stars to a King model, with acore radius of 2.3-3.0 pc. Our profile is also at least roughly consistent with the cumulative numberdistribution displayed by Raboud & Mermilliod (1998). Our best-fit tidal radius is slightly largerthan the 17 pc cited by these authors.The surface mass density is plotted in an analogous fashion, again as a function of the projectedradius. We show both the data points (small open circles) and, as the dashed curve, the best-fitKing model. Here, the core radius is 1.3 pc, and the tidal radius is 18 pc. Note that the massdensity profile falls off more steeply than the number density. Thus, stars near the center areabnormally massive, a trend we shall explore more extensively below.Figure 11 displays the corresponding volumetric densities. As we indicated, the deconvolutionfrom surface profiles assumes spherical symmetry. In fact, the Pleiades is slightly asymmetric, witha projected axis ratio of 1.2:1 (Raboud & Mermilliod 1998). This ellipticity is thought to stemfrom the tidal component of the Galactic gravitational field (Wielen 1974). Under the sphericalassumption, the filled circles and solid, smooth curve show the number density. Here, the Kingmodel is the same used for the surface number density in Figure 10, but deprojected into three-dimensional space.Figure 11 also shows, as the small open circles and dashed curve, the mass density as a functionof spherical radius. Again, the King model here is the deprojected version of that from Figure 10.The relatively rapid falloff in the mass, as opposed to the number, density is another sign of thetendency for more massive stars to crowd toward the center.The information we used in obtaining these profiles also gives us the spatial variation of thebinary fraction b . That is, we first used equation (37) to obtain µ r , the predicted mass distributionin each radial annulus. Recall that the distribution refers to both primaries and secondaries, as 19 –well as single stars. The binary fraction can thus be computed locally. To within our uncertainty,about ± .
05 at each radial bin, we find no variation of b across the cluster. We have mentioned, in a qualitative manner, that more massive cluster members tend to residenearer the center. In Figure 12, we explicitly show this trend. Here, we plot h m p + m s i , the averagesystem mass (primary plus secondary) as a function of the projected cluster radius. It is apparentthat h m p + m s i monotonically falls out to about 4 pc. Beyond that point, the average mass isroughly constant.The pattern here is consistent with mass segregation, but is not a clear demonstration of thateffect. The problem is that Figure 12 gives no indication of the relative populations at differentannuli. If the outer ones are occupied by only a small fraction of the cluster, is mass segregationpresent? To gauge any variation in the mass distribution of stars, that distribution must becalculated over an adequate sample size.Previous authors have also claimed evidence of mass segregation, using various criteria. Adams et al.(2001) looked at the distribution of surface and volumetric number densities for a number of differ-ent mass bins. Raboud & Mermilliod (1998) divided the population by magnitude into relativelybright and faint stars. They calculated the cumulative number as a function of radius for bothgroups, and found the bright stars to be more centrally concentrated. Finally, Pinfield et al. (1998)fit King profiles to the surface density of various mass bins. As the average mass increases, the coreradius shrinks.Figure 13 gives a simpler and more clear-cut demonstration of the effect. Here, we consider f N , the number of sources enclused in a given projected radius, divided by the total number ofsources in the cluster. We also consider f M , the analogous fractional mass inside any projectedradius. The figure then plots f M versus f N . In the absence of mass segregation, f M would equal f N at each annulus. This hypothetical situation is illustrated by the dotted diagonal line.In reality, f M rises above f N , before they both reach unity at the cluster boundary (see uppersmooth curve, along with the data points). This rise, as already noted, indicates that the innermostportion of the cluster has an anomalously large average mass, i.e., that mass segregation is present.Moreover, the area between the solid and dotted curves is a direct measure of the effect. In the caseof “perfect” mass segregation, a few stars near the center would contain virtually all the clustermass. The solid curve would trace the upper rectangular boundary of the plot, and the enclosedarea would be 0.5. We thus define the Gini coefficient , G , as twice the area between the actual f M − f N curve and the central diagonal. For the Pleiades, we find that G = 0 . ± . The name derives from economics, where the coefficient is used to measure inequality in the distribution of wealth
20 –It is possible, at least in principle, that this effect is due entirely to a few exceptionally massivestars located near the center. In fact, this is not the case. We have artificially removed the 11brightest sources (all late-B stars) and recalculated f M versus f N . The result is shown by thedashed curve in Figure 13. While the rise above the diagonal is diminished, it is still present. Thatis, the intermediate-mass population exhibits segregation, as well.An interesting contrast is presented by another populous group, the Orion Nebula Cluster(ONC). The distribution of stellar masses in this far younger system was recently studied byHuff & Stahler (2006). Figure 2 in that paper compares the stellar populations in the inner andouter halves of the cluster. Apart from a few high-mass objects, the two populations are essentialidentical.We may also construct an f M − f N curve for the ONC, as shown here in Figure 14. The solidcurve again lies well above the fiducial diagonal, ostensibly indicating mass segregation. However,removal of just the four Trapezium stars gives a dramatically different result (dashed curve) thatis virtually indistinquishable from the diagonal. All stars except this tiny subset are similarlydistributed. The cluster is too young to have undergone true mass segregation, a conclusion drawnpreviously from N-body simulations (Bonnell & Davies 1998). The Trapezium represents a specialpopulation, one that probably formed just prior to cloud dispersal (Huff & Stahler 2007).
5. Discussion
In this paper, we have applied a versatile statistical tool, the maximum likelihood technique, toassess the distribution of stellar mass and the incidence of binaries in the Pleiades. We began witha near-infrared catalog of cluster members. Our basic assumption was that all cluster membersshare the same evolutionary age, and that any dispersion in the color-magnitude plane stems frombinarity and random photometric errors. We were then able to infer the most probable distributionof masses, both for the cluster as a whole, and as a function of distance from its center. Finally,we introduced a simple method for gauging the degree of mass segregation in the cluster.One of our surprising results is the relatively high fraction of binaries. We estimate that 68%of all systems in the cluster are unresolved binaries; this figure climbs to about 76% if resolvedpairs are included. These fractions are significantly higher than the accepted field-star result, sowe should scrutinize them carefully. Could they stem from an underestimate of the photometricerror at faint magnitudes? Since the error in I is greater than K , we artificially increased thedispersion σ I . We kept σ I at 0.15 until M I = 9 .
5, below which we increased it linearly, reaching σ I = 0 .
20 at M I = 12. After redoing the maximum likelihood analysis, the global binary fraction (Sen 1997). As we show in the Appendix, G is also half the mean mass difference of radial shells, where that massdifference is normalized to the average system mass in the cluster. Note that the axis labels in Figure 2 of Huff & Stahler (2006) were inadvertantly switched.
21 – b for unresolved pairs is unchanged.Another potential difficulty is our neglect of the physical thickness of the cluster. We haveassigned all members a distance of 133 pc, although there will naturally be some variation. However,this effect is also relatively small. From Section 4.1, the volumetric number density falls off withradius approximately as a King model with core and tidal radii of 2.1 and 19 pc, respectively.Consider the front half of a spherical cluster with such a density distribution. It may readily beshown that the mean distance from the plane of the sky of any cluster member is d = 2 . D , the induced magnitude spread is 5 log [( D + d ) /D ], which is 0.04 inour case. Although the actual spread in magnitudes is not Gaussian, we have added this figure inquadrature to both σ I and σ K , and rerun the analysis. Again, the binarity is unaffected.The errors due to both photometry and finite cluster thickness induce a symmetric spread instellar magnitudes. That is, they scatter as many sources below the fiducial isochrone as above it.Thus, they cannot reduce the estimated binarity, which stems from an excess of stars above theisochrone. One systematic error that would affect b is an overestimation of the cluster distance.If D were lowered, the absolute magnitudes of all sources would decrease equally, as would theinferred b -value. Quantitatively, the distance would have to decrease by about 15 pc to bring thebinary fraction down to the field-star result for G-dwarf primaries. An error of this size for theaverage distance is excluded by current observations, for which the estimated uncertainty is only1 pc (Soderblom et al. 2005).Since our method relies solely on photometry to assess binarity, we cannot distinguish betweenphysically linked pairs and chance alignments. As mentioned in Section 3.4, the resolution limit ofour data is 10 ′′ , or ∆ r ◦ = 1400 AU at the distance of the Pleiades. Consider a star at a radius r from the cluster center. Its average number of neighbors within ∆ r ◦ is π ∆ r ◦ n s ( r ),where n s ( r ) isthe projected surface number density of the cluster. Since each ring of width dr contains 2 π n s r dr stars, integration over all members yields the total number of chance alignments: N chance = 2 π ∆ r ◦ Z R n s ( r ) r dr . (44)where R is the cluster’s outer radius, Using n s ( r ) from Figure 10, we find N chance = 2 .
4. Thus,chance alignments have no quantitative impact.Yet another source of systematic error is the cluster age. We have adopted the lithium-basedfigure of 125 Myr from Stauffer et al. (1998). Earlier estimates, using the main-sequence turnoff,yielded a range of answers. For example, Meynet et al. (1993) found 100 Myr. Even this minorreduction affects our results, since it lifts the low-mass end of the isochrone toward higher luminosity.For an age of 100 Myr, our analysis gives b = 0 . ± .
02 and c = 0 . ± .
06. The binary fractionis augmented to 0.64 when we include resolved pairs. From Figure 3 of Stauffer et al. (1998), a100 Myr age corresponds to a lithium edge at M I = 11 .
7, or I = 17 . A. Interpreting the Gini Coefficient
In Section 4.2, we introduced the Gini coefficient geometrically, as twice the area between the f M − f N curve and the diagnonal line representing zero mass segregation. Alternatively, G may bedefined in terms of the mean mass difference between radial shells in the cluster. Here we describemore precisely, and prove the equivalence of, this second interpretation.Altering our previous notation, we now let m ( r ) be the average system mass (primary plussecondary) in a shell with outer radius r . If there are many shells, then we may define a systemnumber density n ( r ), such that the number of systems between r and r + dr is n ( r ) dr . The totalnumber of systems at all radii is N tot = Z ∞ n ( r ) dr . (A1)This total was called µ tot in the text. The average system mass throughout the entire cluster is¯ m = 1 N tot Z ∞ m ( r ) n ( r ) dr . (A2)We will be concerned with the relative mean difference in the mass of shells. This is¯∆ = 1¯ m N Z ∞ dr Z ∞ dr ′ | m ( r ) − m ( r ′ ) | n ( r ) n ( r ′ ) . (A3)We will prove that the Gini coefficient, as defined in Section 4.2, is also ¯∆ / G , we utilized the cumulative fractional number f N and cumu-lative fractional mass f M . These may be written in terms of the system number density: f N ( r ) = 1 N tot Z r n ( r ′ ) dr ′ , (A4) f M ( r ) = 1¯ m N tot Z r m ( r ′ ) n ( r ′ ) dr ′ . (A5)We will later need the differentials of f N and f M , which are df N = n ( r ) N tot dr , (A6) df M = m ( r ) n ( r )¯ m N tot dr . (A7)In terms of f N and f M , the geometric definition of G is then G = 2 Z ( f M − f N ) df N . (A8)If the mass is centrally concentrated, as expected in a real stellar cluster, then f M ≥ f N at allradii, and G ≥
0. Hypothetical clusters in which larger masses are preferentially located fartherfrom the center would have G ≤
0. 24 –We now return to equation (A3) and manipulate it to obtain G , as defined above. A centralassumption we will make is that m ( r ) declines monotonically with r . We may then expand therighthand side of equation (A3). We split the integral over r ′ into two parts, one for r ′ ≤ r and theother for r ′ > r . Under our assumption, m ( r ′ ) ≥ m ( r ) in the first integral, and m ( r ′ ) < m ( r ) inthe second. After pulling n ( r ) from the r ′ -integration, we further split the integrands to find¯∆ = 1¯ m N tot Z ∞ ( −I + I + I − I ) n ( r ) dr . (A9)The first two terms of the integrand are I ≡ m ( r ) N tot Z r n ( r ′ ) dr ′ = m ( r ) f N ( r ) , (A10) I ≡ N tot Z r m ( r ′ ) n ( r ′ ) dr ′ = ¯ m f M ( r ) . (A11)The third term is I ≡ m ( r ) N tot Z ∞ r n ( r ′ ) dr ′ (A12)= m ( r ) N tot Z ∞ n ( r ′ ) dr ′ − m ( r ) N tot Z r n ( r ′ ) dr ′ (A13)= m ( r ) − m ( r ) f N ( r ) , (A14)while the fourth is I ≡ N tot Z ∞ r m ( r ′ ) n ( r ′ ) dr ′ (A15)= 1 N tot Z ∞ m ( r ′ ) n ( r ′ ) dr ′ − N tot Z r m ( r ′ ) n ( r ′ ) dr ′ (A16)= ¯ m − ¯ m f M ( r ) . (A17)Putting equations (A10), (A11), (A14), and (A17) back into equation (A9) yields¯∆ = 1¯ m N tot Z ∞ [2 ¯ m f M ( r ) − m ( r ) f N ( r ) + m ( r ) − ¯ m ] n ( r ) dr (A18)= 2 (cid:18)Z f M df N − Z f N df M (cid:19) , (A19)where we have used both equations (A1) and (A2) to cancel the last two terms on the right side ofequation (A18), and equations (A6) and (A7) to transform the remaining two terms. The secondintegral in the last equation is Z f N df M = Z f N f M =1 f N f M =0 d ( f N f M ) − Z f M df N (A20)= 1 − Z f M df N , (A21) 25 –since f M and f N attain their upper and lower bounds simultaneously. Using this result, equa-tion (A19) becomes ¯∆ = 2 (cid:18) Z f M df N − (cid:19) . (A22)Finally, we note that 2 Z f N df N = 1 . (A23)Thus, equation (A22) becomes¯∆ = 4 (cid:18)Z f M df N − Z f N df N (cid:19) (A24)= 4 Z ( f M − f N ) df N (A25)= 2 G , (A26)as claimed.
REFERENCES
Adams, J. D., Stauffer, J. R., Monet, D. G., Skrutskie, M. F & Beichman, C. A. 2001,AJ, 121,2053Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, A&A, 337,403Bihain, G. et al. 2006, A&A, 458, 805Binney, J. & Merrifield, M. Galactic Astronomy, Princeton: Princeton U. PressBonnell, I. A. & Davies, M. B. 1998, MNRAS, 298, 93Bouvier, J., Rigaut, F. & Nadeau, D. 1997, A&A, 323, 139Brown, A. G. A. 2001, Rev. Mexicana Astron. Astrofis., 11, 89Chabrier, G. 2003, PASP, 115, 763Cohen, M., Wheaton, M. A., & Megeath, S. T. 2003, AJ, 126, 1090Cowen, G. 1998, Statistical Data Analysis, Oxford: Oxford U. PressDeacon, N. R. & Hambly, N. C. 2004, A&A, 416, 125de La Fuente Marcos, 1996, A&A, 314, 453xDias, W.S., Alessi, B.S., Moitinho, A., & Lepine, J.R.D. 2002, A&A, 389, 871 26 –Duquennoy, A. & Mayor, M. 1991, A&AS, 88, 241Fischer, D. A. & Marcy, G. W. 1992, ApJ, 396, 178Garcia, B. & Mermilliod, J. C. A&A, 368, 122Haisch, K. E., Greene, T. P., Barsony, M. & Stahler, S. W. 2004, AJ, 127, 1747Hambly, N. C., Hawkins, M. R. S. & Jameson, R. F. 1993, A&AS, 100, 607Hoffleit, D. et al. 1991, Bright Star Catalogue, New Haven: Yale University PressHuff, E. H. & Stahler, S. W. 2006, ApJ, 644, 355Huff, E. H. & Stahler, S. W. 2007, ApJ, 666, 281King, I. 1962, AJ, 67, 471Kroupa, P. 2001, MNRAS, 322, 231Lada, C. J. 2006, ApJ, 640, 63Mason, B. D., Henry, T. J., Hartkopf, W. I., Ten Brummelaar, T. & Soderblom, D. R. 1998, AJ,116, 2975Mermilliod, J.-C., Rosvick, J. M., Duquennoy, A. & Mayor, M. 1992, A&A, 265, 513Meynet, G., Mermilliod, J.-C., & Maeder, A. 1993, A&AS, 98, 477Moraux, E., Bouvier, J., Stauffer, J. R., & Cuillandre, J.-C. 2003, A&A, 400, 891Moraux, E., Kroupa, P., & Bouvier, J. 2004, A&A, 426, 75Palla, F. & Stahler, S. W. 2000, ApJ, 540, 255Pinfield, D. J., Jameson, R. F., & Hodgkin, S. T. 1998, MNRAS, 299, 955Pinfield, D. J., Hodgkin, S. T., Jameson, R. F., Cossburn, M. R., Hambly, N. C., & Devereux, N.2000, MNRAS, 313, 347Pinfield, D. J., Dobbie, P. D., Jameson, R. F., Steele, I. A., Jones, H. R. A., & Katsiyannis, A. C.2003, MNRAS, 342, 1241Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2002, Numerical Recipes inC: Second Edition, Cambridge: Cambridge U. PressRaboud, D. & Mermilliod, J.-C. 1998, A&A, 329, 101Romani, R. W. & Weinberg, M. D. 1991, ApJ, 372, 487 27 –Schwartz, M. J. & Becklin, E. E. 2004, AJ, 130, 2352Sen, A. 1997, On Economic Inequality, Oxford: Clarendon PressSiess, L., Dufour, E., & Forestini, M. 2000, A&A, 358, 593Soderblom, D., Nelan, E., Benedict, G., McArthur, B., Ramirez, I., Spiesman, W., & Jones, B.2005, AJ, 129, 1616Stauffer, J. R., Schultz, G, & Kirkpatrick, J. D. 1998, ApJ, 499, 199Stauffer, J. R. et al. 2007, ApJS, 172, 663Steele, I. A. & Jameson, R. F. 1995, MNRAS, 272, 630Tassis, K. 2007, MNRAS, 379, L50Wielen, R. 1974, in Stars and the Milky Way System, ed. L. N. Mavridis (Berlin: Springer-Verlag),p. 326Zinnecker, H. & Mathieu, R. D. (eds) 2001, The Formation of Binary Stars, San Francisco: ASP
This preprint was prepared with the AAS L A TEX macros v5.2.
28 –Fig. 1.— Near-infrared color-magnitude diagram for the Pleiades. Small dots represent the 1245stars in our sample. Open circles are the 41 likely sub-stellar objects which have been removed fromthe sample. Filled circles are the 11 brightest stars, which are likely post-main-sequence objects.The 125 Myr isochrone for stars with masses between 0.08 and 4.0 M ⊙ is shown as the smooth,solid curve. 29 –Fig. 2.— Observational error in the K -band measurements as a function of absolute magnitudefor all 1417 stars in the catalog of Stauffer et al. (2007). The smooth curve is the approximate fitused in our maximum likelihood analysis. 30 –Fig. 3.— Sample synthetic data results for the single star mass function φ ( m ). In both panels, thesmooth curve is the input function from which the synthetic data were drawn. Shown also are thebest-fit values, along with errors, for our discrete representation of the function. In the left panel,the input φ ( m ) is a power law with slope -2.8; in the right panel, it is a log-normal function withpeak m ◦ = 0 . σ m = 0 .
4. 31 –Fig. 4.— Weighted sum of the biases as a function of the regularization parameter λ . The syntheticinput here was a log-normal function with 1245 sources. The dashed, horizontal line is drawn at χ b = N , where N = 21 is the number of free parameters in our fitting. The short, vertical arrowindicates the final value of λ used for this synthetic dataset. 32 –Fig. 5.— Comparison of our fitted correlation coefficient c with γ , the imposed degree of correlationin the synthetic dataset. The gray area indicates the region in which the Pleiades most likely falls. 33 –Fig. 6.— Systematic errors in parameters of the log-normal fit to φ ( m ), plotted as a function ofthe synthetic binary correlation γ . Filled circles show ∆ m ◦ , the error in m ◦ . Open circles show∆ σ m , the error in σ m . The gray areas indicate the regions in which the Pleiades most likely falls. 34 –Fig. 7.— Best-fit single star probability density φ ( m ) for the Pleiades. Actual bin values y i / ∆ m i are shown with associated errors. The smooth curve is a log-normal approximation to the results. 35 –Fig. 8.— Comparison of log-normal fits to the primary probability density log φ p ( m p ) (solid curve)and the single star probability density log φ ( m ) (dashed curve). The primary function peaks atlarger mass and has a smaller width. Note that φ p ( m ) includes single stars. 36 –Fig. 9.— Comparison of the Pleiades single star probability density to the field-star initial massfunction of Kroupa (2001), where the latter has been shifted upward for clarity. Shown here is thelog-normal approximation to log φ ( m ) (solid curve), augmented with the data of Bihain et al (2006)for low-mass members and brown dwarfs (open circles) and our 11 brightest, post-main-sequencestars (filled circles). 37 –Fig. 10.— Surface density distribution in the Pleiades. The filled circles represent the surfacenumber density (pc − ), displayed on a logarithmic scale. Open circles are the mass density, in M ⊙ pc − . The solid and dashed smooth curves are King model fits. 38 –Fig. 11.— Volume density profiles. The filled circles represent the number density (pc − ), againdisplayed logarithmically. Open circles are the mass density, in M ⊙ pc −3