text:
# HST/WFPC2 proper motions in two bulge fields: kinematics and stellar population of the Galactic bulgeBased on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Insitute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-2655. These observations are associated with proposal GO-8250.

## 1 Introduction

The bulge of our Galaxy is interesting for a number of reasons. It is the nearest galactic bulge to the sun, and represents a unique place to study the stellar populations and stellar dynamics of such objects in detail. Such analysis can provide important information for our understanding of how bulges formed their stellar populations, what gravitational potential they sit in, and how they came to have the structure they do. In the case of our Galactic bulge, stellar motions are of great interest both in testing the hypothesis that our bulge is in fact a bar (e.g. Zhao, Rich, & Spergel 1996) and in the modeling of microlenses, the overwhelming majority of which are seen in the direction of the bulge (Alcock et al., 2000).

There are several complications in trying to build an understanding of our bulge. Close to the galactic plane, extinction by foreground dust is large and uneven (even on very small scales) and effectively limits optical work to a few “windows” of low and relatively uniform extinction. More troubling is the difficulty that the stellar population seen in these directions samples everything along the line of sight. Simple exponential models for the disk, for example, predict that in Baade’s Window at least half of the stars visible are actually disk stars, and not bulge stars. Employing photometry alone, it is not possible to effectively sort the populations; bulge and disk populations overlap in color, especially near the turnoff (Holtzman et al., 1998) greatly complicating the use of HST-derived CMDs for age determination. Blue stragglers extending brighter than the turnoff in an old population overlap with the main sequence locus of a young population.

Among the first stellar populations imaged in 1994 with the repaired WPFC2 on board HST were fields near two of the low extinction regions originally identified by Walter Baade: what is now known as Baade’s Window $`(l,b)=(1.13^{},3.77^{})`$, and Sgr I $`(l,b)=(1.25^{},2.65^{})`$. A discussion of the photometry and luminosity function of the Baade’s Window field is given in Holtzman et al. (1998).

We noticed that the fields were ripe for a revisit with the aim of measuring proper motions and we proposed successfully (GO-8250). Although HST has been used to measure proper motions of field stars in the rough vicinity of the Galactic Center, these two fields are of special interest because they lie well within the COBE bulge (Dwek et al., 1995) and, in the case of Baade’s Window, abundances are measured (McWilliam & Rich, 1994) and many other studies have been done. Combining proper motions with precision photometry might make it possible to separate the observed populations based on their kinematics. The only previous study of bulge proper motions in Baade’s Window was based on photographic plates (Spaenhauer, Jones, & Whitford, 1992) and reported motions only of stars thought to be candidate red giants in the bulge, by Arp (1965). By systematically excluding the bluer disk stars, and by only measuring a few hundred of the brightest giants, this study, while pioneering, leaves much of the problem ripe for inquiry. Hence our decision to obtain second epoch images of both bulge fields, using WFPC2 on board HST.

## 2 Observations and proper motion measurements

In August of 2000 we obtained HST/WFPC2 observations of two bulge fields for which archival images suitable as first-epoch observations existed. The coordinates and observations are listed in Table 1.

The archival data, though they consist of multiple exposures, are not dithered. This compromises our ability to centroid stars accurately, as the images on the wide-field camera CCDs are very much undersampled. We therefore chose to observe the second epoch data (henceforth referred to as Y2K data) in the F814W filter, which gives the widest, and hence least undersampled, point-spread function (PSF). We also decided to dither the Y2K exposures in a regular 2x2 pattern of offsets (with a step of 0.55 arcsec, or 5.5 pixels on the WF CCDs), ensuring that the uncertainties in the proper motions would not be dominated by the newer epoch.

It is easy to show that the 1-$`\sigma `$ uncertainty in the centroid of a stellar image of full-width at half maximum FWHM that is detected at significance (S/N) is given by

$$0.7\times \text{FWHM}/\text{(S/N)}$$

where the numerical coefficient has a constant value for Moffat-profile PSF’s with asymptotic power-law slopes ranging from 2 to infinity, the latter corresponding to a gaussian PSF (see Appendix). With HST, therefore (FWHM$`0.12`$arcsec in F814W), a star detected at 25-$`\sigma `$ could be centered to about 3.4mas (1/30th of a pixel). Over a 7-year baseline, proper motions can be measured for such stars to an accuracy of 0.7mas/yr. At the distance to the bulge of about 8kpc this corresponds to about 25km/s.

Anderson & King (2000) and Ibata & Lewis (1998) describe ways in which astrometry can be performed in dense star fields using multiple dither observations. Essentially, these methods simultaneously derive a PSF (convolved with the pixel shape) that is consistent with the multiple dither observations of the same star field, and centroids for all stars in the field.

In the case of undersampled data, such as our first-epoch data, it is necessary to have independent knowledge of the shape of the PSF in order to be able to determine an accurate centroid. In dense star fields this information can be extracted provided the undersampling is not too severe: the stars will be positioned randomly with respect to the pixel grid of the detector, and so the ensemble of stars does in fact sample the PSF at a wide range of phases. Our method consists in a modification of the Anderson & King (2000) approach. For each exposure we start from an analytic PSF model, and determine best-fit centers and intensities for the 100 brightest unsaturated stars in the frame that are not affected by neighbours. These centers and intensities are then used to build an ensemble-averaged PSF from these stars. An analytic fit of somewhat higher order is then made to this new sampling of the PSF, with a robust outlier rejection. This new PSF is then used to define new centers and intensities, etc. This process iterates to convergence quite rapidly. In practice we use a gaussian PSF multiplied by polynomials as our model.

Armed with a PSF for each observation, we determine the proper motions as follows. First we combined all undithered data from the 1995 epochs, filter by filter, and identify 20-$`\sigma `$ detections in the combined image using the DAOFIND package. This master coordinate list for each CCD is then transformed using bright reference stars to align it with the positions in each individual exposure. For each F814W exposure a PSF fit is then performed near each transformed master list position, and the best-fit center and intensity recorded. The result is a set of up to 30 (21) position measurements for each star in the Baade Window (Sgr-I) fields, referred to the pixel positions of the different exposures, and spread over three (two) epochs.

These position measurements contain the proper motion information we seek, but also plate solution mappings from one exposure to another, and residual systematic centroiding errors. An iterative procedure, consisting of the following loop, separates the various effects: (i) find an average transformation that maps the star positions of each exposure onto those of a reference exposure, and apply this transformation to the positions; (ii) determine proper motions for each star based on these transformed positions, and determine position residuals for each position measurement from the best-fit proper motion; (iii) determine the average position residual as a function of pixel phase of each stellar image, and subtract this residual from each measurement; (iv) determine the average position residual as a function of phase with respect to the ‘34th-row’ effect (Anderson & King, 1999) and subtract this. At each pass through the loop a new proper motion solution is made for each star, using a weighted linear least-squares fit with rejection of outliers. The weights are deduced from the estimated centroiding error of each image on the basis of the number of photons detected, including an estimated residual systematic error (see Figure 1).

Effects (iii) and (iv) are only detectable with dithered data, and so we have had to assume that the systematic rediduals with pixel phase and 34th-row phase (both on the order of 2–5% of a pixel) apply also to the earlier epochs. Any error we make in this assumption is not easily detectable in our data: any such positional error in the first epoch will simply be absorbed into the proper motion measurement. Two pieces of evidence give us confidence that our results are not seriously affected by such residual systematics:

1. The multiple dithers observed in our Y2K data provide a test of the procedure. With many dithered observations at the same epoch, residual errors show up as extra scatter about the best-fit position that cannot be absorbed into the proper-motion solution. Figure 1 shows that for the fainter stars our errors are safely dominated by photon statistics, and that systematic errors set in at the 0.03-pixel level on the undersampled WF detectors (corresponding to 20km/s over a 7-year baseline at the distance of the bulge).
2. Comparisons of the distributions of the proper motions between the PC and WF detectors shown below do not indicate that our WF proper motions contain larger errors than those derived from the PC.

Using this technique, we obtained proper motion measurements for some 20,000 stars in the Sgr-I field, and for 15,000 stars in Baade’s window. In addition, we used the 1994/1995 archival data for each field to derive magnitudes in the F814W and F555W filters, using the magnitude calibrations given in Holtzman et al. (1995).

Note that the proper motions derived are relative: we arbitrarily assign the mean proper motions of the complete sample of stars in each field to be zero in both components. Absolute (i.e., with respect to an inertial frame) proper motion measurements would require extragalactic sources such as qso’s to be identified in these fields. Perhaps future spectral or multi-color studies can provide this extragalactic reference frame.

## 3 Analysis

Proper motion distributions and color-magnitude diagrams are presented in Figures 2 and 3. The velocity dispersions per CCD are listed in table 2, and are consistent with each other—there is no indication of systematic differences between the velocity dispersions measured on the WF and PC CCD’s for the same fields, but there are field-to-field differences in the measured velocity dispersions. The proper motion distributions are skewed, with a longer tail towards positive $`\mu _l`$. We argue below that this is in large part due to foreground contamination.

At first inspection, the color-magnitude diagram shows a turnoff population with a blue extension and a red giant branch. A first cut at separating the kinematics of the various stellar populations can be attempted using those stars brighter than the old turnoff point, which we expect to show clear differences in kinematics (Figure 4). In fact, the blue main sequence stars have a mean positive longitudinal proper motion of some 5mas/yr with respect to the red giants. This is most simply understood as a foreground, thin disk population which rotates in front of the bulge stars.

The bright blue main sequence stars are the part of the younger disk population that stand out in the color-magnitude diagram (CMD), but of course disk stars are found among fainter stars as well. This is best illustrated in Figs. 5 and 6, which are binned color-magnitude diagrams, with each bin color-coded by various kinematic quantities. These figures clearly show the extent of disk contamination across the CMD, as well as a number of other kinematic features:

1. The disk kinematics of the bright blue main sequence stars (large positive $`\mu _l`$) are shared by redder, fainter stars located above the bulge main sequence. This is naturally explained as the main sequence of disk stars in front of the bulge. As expected, the effect is more pronounced in the Sgr-I field, which lies closer to the galactic plane.
2. There is a gradient in the mean $`\mu _l`$ of faint main sequence stars with magnitude: at a given color, the faintest main sequence stars drift towards negative $`l`$. This suggests that we are measuring stars clear through the bulge rotation field, even to the far side of the bulge, where stars are rotating behind the galactic minor axis.
3. The proper motion dispersions of the main sequence stars decrease with magnitude: at a given color the faintest main sequence stars have smaller dispersion in proper motion. We interpret this as a distance effect: the faintest stars are further away and hence a given velocity translates into a smaller proper motion.
4. There is a distinct suggestion in these diagrams of a kinematically homogeneous bulge population which follows a single isochrone, characteristic of an old stellar population.

In order to investigate the kinematics further, we make cuts of main sequence stars for which a crude distance modulus can be calculated. We find that the quantity

$$M^{}=m_{F814W}2\left(m_{F814W}m_{F555W}\right)$$

removes the slope of the main sequence in the CMD, so that it may be used as a simple relative distance indicator. As a function of this distance modulus we plot the proper motions in Figs. 7 and 8. In both Baade’s Window and Sgr I, the same picture emerges. We see a smooth gradient of rotation down the line of sight. This is either due to contamination of a non-rotating bulge by disk stars, or represents a true rotation of the bulge. The notable similarity of the slope of both ‘rotation curves’ in both fields (even though the Sgr-I field is expected to be more disk-dominated) argues for an intrinsic bulge rotation.

Radial velocity surveys find rotation in the bulge population; perhaps the best example is that of Izumiura et al. (1995), which shows the rotation curve based on a survey of SiO masers. The practicial difficulty with measuring the rotation curve of the pure bulge population is that as one moves off-axis, the density of the bulge drops rapidly and sample contamination (from the halo and disk) becomes a major problem. Our study gives one of the cleanest demonstrations of rotation in a bulge star dominated sample.

Some undulations in the vertical mean proper motion are also seen. There is the hint of a rise in the vertical proper motion for those stars which might correspond to the “middle” of the bulge. At this point it is not clear whether these trends in the vertical mean proper motion are real, or a manifestation of an unknown systematic residual error.

It is notable that the proper motion dispersions decrease for the faintest $`M^{}`$. This is as expected if these stars are indeed at larger distances: these stars would lie a factor two further away than stars with $`M^{}`$ ca. 1.5 magnitudes brighter.

Detailed modelling of these results is complicated by the uncertain scatter in these crude photometric parallaxes. However, we can estimate values for the velocity dispersions of bulge stars in each field by considering the $`M^{}`$ bin with zero mean $`\mu _l`$. These stars should lie as close to the Galactic minor axis as these lines of sight reach. The numbers are given in Table 3. We see evidence for a subtle anisotropy in the bulge, but this could be an effect of the velocity gradient through the $`M^{}`$ bin.

The tangential dispersions in Table 3 are similar to those for the full sample (Table 2). This coincidence is caused by two effects which act in opposite directions. The full sample contains many distant stars, whose proper motions are small and hence reduce the dispersion; on the other hand the line-of-sight gradient in tangential velocity tends to increase the width of the distribution. Only the first effect operates in the vertical direction, so $`\sigma (\mu _b)`$ values are larger in Table 3 than they are in the full sample.

## 4 The stellar population of the bulge

With our picture of the kinematics in the bulge region, we now study the stellar populations. We isolate a ”pure bulge” sample by keeping only stars with negative $`\mu _l`$ and moderate $`|\mu _b|`$. Given the rotation gradient of the bulge, this kinematic cut may also preferentially remove near-side bulge stars from the sample, so it may be slightly biased in distance. Far-side disk contamination should be very small as the line of sight rises away from the galactic plane behind the bulge.

The color-magnitude diagrams that these cuts produce in Baade’s Window and the Sgr I field are shown in Figure 9. The $`\mu _l,\mu _b`$ cuts cleanly remove the bright blue portion of the main sequence and foreground M dwarfs, by kinematic selection only. The color-magnitude diagram that remains strongly resembles that of the old bulge globular cluster NGC 6553 (cf. Zoccali et al. 2001). A handful of stars brighter than the old turnoff remain (field blue stragglers or outliers?) but the dominant blue extension is clearly demonstrated to be a foreground population. The present color-magnitude diagram suggests that blue stragglers are not common in the field population (for example, not nearly as frequent as found in the cluster NGC 6553 by Zoccali et al. 2001). However, it does not put interesting limits on their numbers. To do so would require radial velocities and perhaps line strength measurements. When additional analysis is finished, we will be able to derive the age of the bulge field relative to NGC 6553 and 6528; a cursory examination of Figure 9 shows the characteristic turnoff point of a clearly old stellar population.

The impression that two very different populations are being identified via the proper motions is amplified when one considers the luminosity functions for the bulge and disk populations. Recall that (especially at the bright end) incompleteness and photometric errors are identical, being only magnitude dependent. Figure 10 shows our plot for the Baade’s Window disk and bulge samples, and the disk’s extension (relative to the sharp rise of the main sequence turnoff) is evident. These plots give further strength to the visual impression from Figure 9.

We emphasize that our population separation is based on kinematics, not on any photometric or spectroscopic classification. While disk and bulge exist cospatially, they separate well kinematically, and this is therefore at least a complementary, and we would argue a more robust, technique for isolating a ‘pure bulge’ population.

## 5 Conclusions

With HST photometry and proper motions determined with high enough precision, it is possible to separate the disk and bulge populations by their kinematics alone. The long standing question regarding the nature of the blue main sequence extension in the bulge field population is settled: the great majority of those stars evidently belong to the foreground disk. When these stars are excluded, the old turnoff population in the bulge remains, and no measureable population of blue stragglers or intermediate age stars is present. This was first demonstrated for the Baade’s Window field by Ortolani et al. (1995). Feltzing & Gilmore (2000) came to the same conclusion based on counts of stars brighter and fainter than the turnoff point in their WFPC2 data. However, a proper accounting for the foreground disk has been a persistant issue, and our application of the kinematic data strengthens greatly the conclusion that the bulge is dominated by old stars.

We take our study a step further, and find direct evidence for the rotation of the bulge population. The observed proper motion anisotropy of bulge stars is largely caused by the line-of-sight gradient of the rotation of the bulge; when this is removed a nearly isotropic velocity distribution of the bulge stars results. The velocity dispersion declines from Sgr I to Baade’s Window.

When appropriate samples in Baade’s Window are compared, our proper motion dispersions agree with those found by Spaenhauer, Jones, & Whitford (1992) and by Feltzing & Johnson (2001). Our results compare well with predictions of the Zhao (1996) and Zhao, Rich, & Spergel (1996) bulge model, but since our proper motion samples are more than an order of magnitude larger than existing radial-velocity and proper-motion samples of bulge stars, and extend well below the turnoff, the time is ripe for more involved modelling. Extension of this work to further bulge fields (e.g. Zoccali et al. 2001), and combination of these results with spectroscopy (for radial velocities, metallicity, and improved distance estimates) should open the way for a new chapter in our study of the Galactic bulge.

Support to RMR for proposal GO-8250 was provided by NASA through a grant form the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.

## Appendix A Accuracy of PSF-fitting photometry and positions

Here we calculate the accuracy of object centroids derived by means of PSF fitting.

Let $`f_i`$ be the data: intensities on the pixels $`i`$ at positions $`x_i,y_i`$ on the image plane. We try to model these data as $`AP(x_i\mu _x,y_i\mu _y)\mathrm{\Delta }^2`$, where $`\mathrm{\Delta }`$ is the pixel width, and $`P`$ is the PSF normalized to total intensity 1. $`A`$ is the intensity of the star, and $`(\mu _x,\mu _y)`$ are the position of the star, to be fitted for. Write

$$\chi ^2=\underset{i}{}\left(f_iAP(x_i\mu _x,y_i\mu _y)\mathrm{\Delta }^2\right)^2/\sigma _i^2$$

where $`\sigma _i`$ is the error on the measured intensity in pixel $`i`$. Then the minimum of $`\chi ^2`$ gives the best-fit PSF, and the second partial derivatives of $`\chi ^2`$ with respect to any two parameters give twice the inverse covariance matrix. In what follows we ignore covariances between the different variables $`A`$, $`\mu _x`$ and $`\mu _y`$, as is appropriate for bisymmetric PSF’s.

At the best-fit value (we may assume without loss of generality that $`\mu _x=\mu _y=0`$) we obtain the inverse variance on the star’s intensity $`A`$ as

$$\text{Var}(A)^1=0.5^2\chi ^2/A^2=\underset{i}{}P(x_i,y_i)^2\mathrm{\Delta }^4/\sigma _i^2$$

which reduces to

$$=\left(P^2𝑑x𝑑y\right)\mathrm{\Delta }^2/\sigma ^2$$

if the PSF is fully sampled and $`\sigma _i`$ constant (i.e., background- or read noise-limited data).

The inverse variance of the best-fit position is similarly (removing terms which go to zero at the best fit)

$$\text{Var}(\mu _x)^1=0.5^2\chi ^2/\mu _x^2=\underset{i}{}A^2(P/x)^2\mathrm{\Delta }^4/\sigma ^2=\left((P/x)^2𝑑x𝑑y\right)A^2\mathrm{\Delta }^2/\sigma ^2$$

(and similarly for $`\mu _y`$). These two relations can be combined to give

$$\delta \mu _x=\frac{\delta A}{A}\times \sqrt{\frac{P^2𝑑x𝑑y}{(P/x)^2𝑑x𝑑y}}$$

where $`\delta A=\text{Var}(A)^{1/2}`$ is the 1-$`\sigma `$ error on $`A`$, etc.

The centroid error therefore depends on the significance $`A/\delta A`$ of the detection of the star, and on a geometric factor governed only by the shape of the PSF.

For a gaussian PSF, dispersion $`s`$, we find

$$\text{Var}(A)^1=\mathrm{\Delta }^2/(4\pi s^2\sigma ^2)$$

hence

$$\delta A=\sqrt{4\pi }s\sigma /\mathrm{\Delta }$$

and

$$\text{Var}(\mu _x)^1=A^2\mathrm{\Delta }^2/(8\pi s^4\sigma ^2)$$

hence

$$\delta \mu _x=\sqrt{8\pi }s^2\sigma /(\mathrm{\Delta }A)=\sqrt{2}s\frac{\delta A}{A}$$

For a Moffat function PSF of the form

$$P=\frac{\beta 1}{\pi a^2}\left(1+\frac{r^2}{a^2}\right)^\beta $$

we get

$$\delta \mu _x=a\frac{\delta A}{A}\sqrt{\frac{2\beta +1}{\beta (2\beta 1)}}$$

If we write $`a`$ as $`\frac{1}{2}\text{FWHM}/\sqrt{2^{1/\beta }1}`$, we find that for $`\beta >1.5`$, which covers PSF’s with tails as shallow as $`r^3`$

$$\delta \mu _x=\frac{\delta A}{A}\times 0.67\times \text{FWHM}\pm 10\%.$$

(The gaussian is the limit $`\beta \mathrm{}`$; in this case the coefficient is 0.6.)


format: text

text:
# Analysis of colour-magnitude diagrams of rich LMC clusters: NGC 1831

## 1 Introduction

The Large Magellanic Cloud (LMC) presents three essential characteristics that make it an excellent complementary laboratory for studying the formation and evolution of galaxies and stellar systems in general: a) it is close to the Galaxy; b) it has markedly different morphological, chemical and kinematical properties from our Milky Way; c) it presents a large variety of stellar clusters, displaying distinct physical characteristics among themselves and when compared to those in the Galaxy (Westerlund 1990). Due to the diversity in ages and metallicities, LMC clusters are found at distinct evolutionary stages (Westerlund 1990, Olszewski et al. 1991). The determination of a cluster’s present physical properties, such as density profile, shape, internal velocity distribution and its position dependent Present Day Mass Function (PDMF), provide us with essential links needed to assess the role of gravitational dynamics, including effects of mass segregation and stellar evaporation (Heggie & Aarseth 1992, Spurzem & Aarseth 1996, de Oliveira et al. 2000). Thus, once these present properties are known, modelling techniques like N-body simulations allow us to recover the initial conditions under which clusters formed (Goodwin 1997, Vesperini & Heggie 1997, Kroupa et al. 2001). In this context, the initial mass function (IMF) and its possible universality, are key pieces in the study of stellar contents of distant galaxies (Kroupa 2001).

Through the analysis of colour-magnitude diagrams (CMDs) one has a great variety of physical information about a cluster. Isochrone fits help constraining the system’s age, metallicity, distance and reddening. Furthermore, observational luminosity functions (LFs) have allowed derivation of the stellar mass function (Elson et al. 1995, De Marchi & Paresce 1995, Santiago et al. 1996, Piotto et al. 1997, de Grijs et al. 2002a). However, the transformation of luminosity into mass depends on age and metallicity, the uncertainties in these parameters being therefore incorporated into the inferred mass function. Besides, the theoretical uncertainties in the mass-luminosity relation itself, specially for the low-mass stars (Piotto et al. 1997, Baraffe et al. 1998), added to the effect of unresolved binarism, further hampers real mass function determination through observational luminosity functions.

From both observational and theoretical points of view, the advent of the Hubble Space Telescope (HST), associated with the constant improvement in theories of stellar interiors, atmospheres and evolution, require ever more sophisticated methods of CMD analysis. In this context, computational modelling techniques have allowed a wider use of CMDs as tools to constrain physical properties of stellar populations and systems. Model CMDs, along with statistical methods of comparing them to observed ones, have been useful means to investigate the star formation history within a galaxy (Gallart et al. 1996, Gallart et al. 1999, Hernandez et al. 1999, Holtzman et al. 1999, Hernandez et al. 2000) or to constrain structural parameters and the stellar luminosity function in the Milky Way (Kerber et al. 2001).

With these issues in mind, our work aims at extracting as much physical information as possible about rich LMC clusters, by means of objective comparison of their observed CMDs with artificial ones. The idea is to simultaneously infer PDMF slope, age, metallicity, distance, reddening and unresolved binary fraction for each system studied. The present work introduces the techniques developed for that purpose and shows the results for NGC 1831, one of the richest LMC clusters for which we have deep HST data.

Previous works are evidence of the large difficulties in extracting physical parameters for NGC 1831. Techniques based on CMDs from CCD photometry (Mateo 1987,1988; Chiosi 1989; Vallenari et al. 1992; Corsi et al. 1994), integrated spectroscopy or colours (Bica et al. 1986; Meurer et al. 1990; Cowley & Hartwick 1992; Girardi et al. 1995) and spectroscopy of individual giant stars (Olszewski et al. 1991) were employed with this aim, constraining the values of the main parameters: $`8.50\stackrel{<}{}\mathrm{log}(\tau )\stackrel{<}{}8.80`$ ($`300\stackrel{<}{}\tau \stackrel{<}{}650`$ Myr); $`0.002\stackrel{<}{}Z\stackrel{<}{}0.020`$ ($`1.00\stackrel{<}{}[\mathrm{Fe}/\mathrm{H}]\stackrel{<}{}0.01`$); $`0.00\stackrel{<}{}E(BV)\stackrel{<}{}0.07`$. For the distance to NGC 1831, there are not reliable determinations, the standard procedure being the adoption of typical values of the intrinsic distance modulus, $`(mM)_0`$, for the LMC centre. In this aspect, the most reliable estimate seems to be that of Panagia et al. 1991, $`(mM)_0=18.51`$, since it is based on purely geometrical arguments applied to high quality imaging and spectral data on supernova SN1987a. In terms of dynamical evolution for this system, Elson et al. (1987) estimated $`6.5\stackrel{<}{}\mathrm{log}(t_{\mathrm{cross}})\stackrel{<}{}7.0`$ and $`9.6\stackrel{<}{}\mathrm{log}(t_{\mathrm{rh}})\stackrel{<}{}10.0`$ for the crossing time and two-body relaxation time, respectively, the range quoted being due to different mass-luminosity relations. Comparing with its estimated age values, these results suggest that NGC 1831 is a system dynamically well mixed, but not totally relaxed through stellar encounters. Hence, NGC 1831 is sufficiently old to have suffered mass segregation, affecting the PDMF slope at different distances to its centre, but perhaps still young enough that the initial conditions could be preserved in its outer regions. Similarly, external effects may not have had enough time to affect the cluster dynamics either.

One of the main objectives of this paper is to verify the effect of mass segregation in NGC 1831, quantifying the variation in PDMF slope with projected distance to the cluster’s centre. This determination may yield strong links to IMF reconstruction efforts based on N-body models. The paper is outlined as follows: in Sect. 2 we present the data and the methods of accounting for sample incompleteness and field star contamination in the observed CMD; in Sect. 3 we present the algorithm used for CMD modelling and the statistical tools used for model vs. data comparisons; in Sect. 4 we discuss control experiments used for to verify the validity of the method; finally, in Sect. 5 we present the results for the NGC 1831 data, which are discussed in Sect. 6.

## 2 The observed CMD

We have data taken with the Wide Field and Planetary Camera 2 (WFPC2) on board HST for 8 rich LMC clusters and nearby fields. These data are part of the GO7307 project, entitled “Formation and Evolution of Rich Star Clusters in the LMC” (Beaulieu et al. 1999, Beaulieu et al. 2001, Johnson et al. 2001). For each cluster, images were obtained using the F555W (V) and F814W (I) filters. Most of the photometry had been previously carried out: cluster stellar LFs were built and analyzed by Santiago et al. (2001), de Grijs et al. (2002b); field stellar populations were studied by Castro et al. (2001). Exposure times, field coordinates, image reduction and photometry processes are described in detail by those authors.

Fig. 1 shows the observed WFPC2 CMD for stars in the direction of NGC 1831 in panel (a) (hereafter the on-cluster field) and for a nearby field (hereafter the off-cluster field) studied by Castro et al. (2001) in panel (b). The on-cluster sample shown here is the final composition of the CEN and HALF images described by Santiago et al. (2001). These have the Planetary Camera (PC) centred on the cluster’s centre and half-light radius, respectively. The off-cluster field is located at about 7.3 arcmin away from the cluster’s centre. A clear cluster main sequence (MS) is visible in the figure, stretching from $`m_{555}18.5`$ down to $`m_{555}25`$. The cluster MS turn-off is also clearly visible at the upper MS end. Notice, however, that saturation becomes a problem in the HALF field for $`m_{555}\stackrel{<}{}19`$ ($`m_{555}\stackrel{<}{}17.8`$ for the CEN field). Hence, all our subsequent analysis will be based on the CMD fainter than this limit. A branch of evolved stars is seen as well, especially in the range $`18\stackrel{<}{}m_{555}\stackrel{<}{}19`$, where the cluster red clump is located. The subgiant branch at fainter magnitudes is due to field contamination and is largely made up of older ($`\tau >3`$ Gyr) stars.

The on-cluster data suffer from two important effects: sample incompleteness and contamination by field stars. Our CMD modelling algorithm does not incorporate such effects. Therefore it is crucial to adequately correct the observed CMD for them in order to place models and data on equal footing. Quantifying systematic and random photometric uncertainties and either correcting for them or applying them to model CMDs is extremely important as well, as they are responsible for most of the observed CMD spread. These data corrections are the subject of the next subsections.

### 2.1 Random photometric uncertainties

The random photometric uncertainties are the main source of spread in our HST/WFPC2 CMDs. Therefore, a suitable assessment of these uncertainties in both filters is crucial for correctly incorporating this effect into the artificial CMDs.

For the on-cluster field, two independent photometric measurements were available for a fraction of the stars due to the overlap region imaged by both the HALF and CEN fields (Santiago et al. 2001). Thus, we used the stars belonging to this overlap region to estimate the typical photometric uncertainties in the data. For each filter and at each magnitude bin, we estimated the dispersion, $`\sigma ^{}`$, in the distribution of differences between the independent magnitude measurements. For simplicity we assume that $`\sigma ^{}`$ is the composition of two equally-sized realizations of photometric error, $`\sigma `$. Thus, we get $`\sigma ^2=2\sigma ^2`$, and therefore

$$\sigma =\frac{\sigma ^{}}{\sqrt{2}}.$$

We emphasize that the two filters were treated as absolutely independent. As a consequence, the uncertainty in colour $`m_{555}m_{814}`$, $`\sigma _{\mathrm{colour}}`$, will be the quadratic sum of the individual uncertainties:

$$\sigma _{\mathrm{colour}}^2=\sigma _{555}^2+\sigma _{814}^2.$$

Fig. 2 shows the derived uncertainties for MS stars in the two filters, $`\sigma _{555}`$ and $`\sigma _{814}`$, and for the colour, $`\sigma _{\mathrm{colour}}`$, using the expressions above.

For the off-clusters stars, we did not have two independent and overlapping WFPC2 images and, as a result, we could not apply the same method to quantify their photometric uncertainties. The solution found was to employ the same estimate as for the on-cluster stars. As the off-cluster image is deeper and sparser than the on-cluster ones, we can expect that this approximation yields an overestimate of the photometric uncertainties in the off-cluster data. However, the off-cluster CMD is used only for statistical subtraction of field contamination from the on-cluster CMD. We will see later that this particular correction technique is fairly independent of the estimated photometric uncertainties.

### 2.2 Systematic photometric uncertainties

Systematic effects in WFPC2 data have been found by several authors. Johnson et al. (2001) measured an exposure time effect varying from 0.01 to 0.06 for different WFPC2 chips and filters in images of NGC 1805 and NGC 1818, as part of this project. de Grijs et al. (2002b) finds similar trends, but of slightly larger amplitude for the same data. Previously, Casertano & Muchtler (1998) found an exposure time effect but in the opposite sense as Johnson et al. (2001). Colour shifts of $`0.04`$ have also been measured between different WFPC2 chips, possibly due to errors in CTE and aperture corrections or in zeropoints (see also Johnson et al. 1999).

Any photometric biases, as a function of either exposure time or chip, must be eliminated from our observed CMD, since the model CMDs which will be compared to it do not incorporate them.

As there does not seem to exist a consensus on the corrections to be applied, our approach was to empirically measure such biases and to apply appropriate shifts to the data when necessary. We searched for both exposure time and chip vs. chip effects. No systematic effects were found in the off cluster data. The main source of bias in the on-cluster data was found to be an offset between the PC and the Wide Field Camera (WFC) chips in the sense that MS stars with the same $`m_{555}`$ magnitude tend to be bluer by 0.05-0.10 mag when imaged with the PC than with the WFC. This applies to both HALF and CEN images. As the PC in the HALF image is centered on the cluster half-light radius, this effect is unlikely to be due to differences in crowding. In order to correct the data for this effect we first defined MS fiducial lines, taking the median $`m_{555}m_{814}`$ colour at different magnitude bins. This was done separately for each chip and each image (CEN or HALF). We noticed that the WFC MS lines were more stable, always occupying loci in the CMD very close to each other. Thus, we transposed the PC fiducial lines to the corresponding WFC locus. The corrections are shown in Fig. 3. Panel (a) shows the uncorrected PC and WFC (a mean locus of the 3 chips) lines for both the CEN and HALF images. Panel (b) shows the corrected fiducial lines.

### 2.3 Selection effects

#### 2.3.1 Sample incompleteness

Sample incompleteness occurs essentially due to two factors: overlapping of stellar profiles due to crowding and background noise. Therefore, in a cluster incompleteness will depend not only on apparent magnitude, but also on spatial position. Faint stars in dense stellar regions are the ones that suffer most from this effect. Completeness effects in the on-cluster data were previously discussed and measured by Santiago et al. (2001). These authors carried out experiments where artificial stars were added to the cluster image and subjected to the same sample selection as the real stars. The completeness c of each real star was estimated as the fraction of artificial stars of similar magnitude and location which were successfully recovered in the experiments. The estimated weight w for each star is then simply given by the inverse of its completeness ($`w=1/c`$). This weight corresponds to the total number of stars similar to the observed one which should be detected and measured in an ideal image.

Assigning a weight to each observed star is enough for the sake of luminosity and mass functions. However, correcting a CMD for incompleteness requires an extra step, namely to place an extra number of stars on the CMD according to the position of each observed star and its previously computed completeness weight. In order to fill the observed CMD with the missing stars, we first defined a fiducial line for the MS. As described in Sect. 2.2, this line was defined by taking the median $`m_{555}m_{814}`$ colour at different magnitude bins. Given the $`i^{th}`$ MS star, its position along the fiducial line is provided by its $`m_{555}`$ magnitude and the corresponding $`m_{555}m_{814}`$ colour. If its completeness weight is $`w_i`$, then $`w_i1`$ extra stars were spread out from its position along the MS line taking the measured photometric uncertainties (as described in Sect. 2.1) into account. We assumed a Gaussian distribution of uncertainties both in $`m_{555}`$ magnitude and $`m_{555}m_{814}`$ colour. Sample completeness falls to less than $`50\%`$ for $`m_{555}\stackrel{>}{}23.8`$. We completed the CMD down to $`m_{555}=24.5`$ and then cut it at $`m_{555}=23.5`$, therefore avoiding boundary effects. A complete sample of stars both in number and in position resulted from this method.

As for the off-cluster CMD, it is complete at least down to $`m_{555}=24.5`$ (Castro et al. 2001).

#### 2.3.2 Field star subtraction

Field star subtraction from the cluster sample is carried out in two separate steps. We first cut-off all stars in the on-cluster CMD which are much farther from the MS than expected by photometric errors. We therefore eliminate all evolved stars, as well as objects which are likely to be foreground stars or remaining non-stellar sources in the sample (distant galaxies, spurious image features, etc).

The second step involves the statistical removal of field stars located along the MS. We thus compare the distribution of stars in the on-cluster CMD to that of the off-cluster CMD. The comparison method is based on the hypothesis that the positions of the off-cluster stars represent the most likely positions for field stars on any similar CMD. We thus try to estimate the probability of each on-cluster star to be a field star and according to this probability we randomly remove stars from the on-cluster CMD.

We consider pairs of on-cluster/off-cluster stars. For each pair we compute the expected scatter in CMD position between the pair members under the assumption that they are independent photometric realizations of the off-cluster star. We will then have

$$\sigma _{\mathrm{c},555}=\sqrt{2}\sigma _{\mathrm{off},555}$$

and

$$\sigma _{\mathrm{c},\mathrm{colour}}=\sqrt{2}\sigma _{\mathrm{off},\mathrm{colour}}$$

respectively for the expected scatter in $`m_{555}`$ magnitude and $`m_{555}m_{814}`$ colour, where the “off” subscript in the expressions above refers to the off-cluster pair member.

For the $`i^{th}`$ off-cluster CMD star we then consider the $`N_i`$ on-cluster stars inside a $`3\sigma _{\mathrm{c},555}`$ x $`3\sigma _{\mathrm{c},\mathrm{colour}}`$ box centered on it. Using a Gaussian error distribution in magnitude and colour, we estimate the probability $`p_{i,j}`$ that the $`j^{th}`$ on-cluster CMD star, inside this box, is a second photometric measurement of the $`i^{th}`$ field star. Therefore,

$$p_{i,j}exp[\frac{(mag_imag_j)^2}{2(\sigma _{\mathrm{c},555})^2}]exp[\frac{(colour_icolour_j)^2}{2(\sigma _{\mathrm{c},\mathrm{colour}})^2}]$$

where the normalization of $`p_{i,j}`$ is such that

$$\underset{j=1}{\overset{N_i}{}}p_{i,j}=1$$

Doing the same for all $`N_{\mathrm{off}}`$ off-cluster CMD stars, we estimate the probability $`P_j`$ that the $`j^{th}`$ on-cluster CMD star is any one of the $`N_{\mathrm{off}}`$ field stars. Hence,

$$P_j=\underset{i=1}{\overset{N_{\mathrm{off}}}{}}p_{i,j},$$

where

$$\underset{j=1}{\overset{N_{\mathrm{on}}}{}}P_j=N_{\mathrm{off}}$$

and $`N_{\mathrm{on}}`$ is the total number of stars in the on-cluster CMD. In practice, not all off-cluster stars will have on-cluster stars within their $`3\sigma _{\mathrm{c},555}`$ x $`3\sigma _{\mathrm{c},\mathrm{colour}}`$ boxes. The actual off-cluster stars taken into account will therefore be $`N_{\mathrm{off}}^{}<N_{\mathrm{off}}`$.

Based on the $`P_j`$ probabilities, and scaling the number of field stars to the solid angle of the on-cluster field, we randomly remove

$$N_{\mathrm{field}}=N_{\mathrm{off}}^{}\frac{\mathrm{\Omega }_{\mathrm{on}}}{\mathrm{\Omega }_{\mathrm{off}}}$$

stars from the on-cluster CMD, where $`\mathrm{\Omega }_{\mathrm{on}}`$ and $`\mathrm{\Omega }_{\mathrm{off}}`$ are the solid angles covered by the on-cluster and off-cluster fields, respectively.

Fig. 4 shows the results of correcting an observed NGC 1831 CMD for incompleteness and field contamination. Panel (a) shows the original CMD obtained from photometry, excluded of non-MS stars and corrected only for the systematic photometric effects described in Sect. 2.2. Panel (b) presents the complete CMD, i.e., corrected for sample incompleteness and cut at $`m_{555}=23.5`$; the extra stars added by completeness correction represent about 23% of the total within this magnitude range. Panel (c) shows the 89 stars ($``$ 3%) in the on-cluster CMD that were considered as LMC field stars in the field star subtraction process; finally the clean and final cluster CMD is shown in panel (d).

We tested the field star subtraction algorithm for different box sizes and assumptions regarding the photometric scatter. The results are insensitive to the details in the algorithm.

## 3 CMD modelling and statistical tools

### 3.1 CMD Modelling

We model the MS of a cluster as a single stellar population, characterized by stars of the same age and chemical composition. The first step is the choice of an isochrone, which defines a sequence of magnitude and colour as a function of mass for this population. For the present work we used Padova isochrones (Girardi et al. 2000) because they present masses inside the observed MS mass range ($`M\stackrel{<}{}2.5`$ $`M_{}`$) and are expressed in the vegamag WFPC2 photometric system. Notice that Padova models assume convective overshooting in the stellar interiors. This assumption may influence age determinations, specially when based in the position of turn-off and He-burning stars of younger populations (Testa et al. 1999, Barmina et al. 2002). However, as NGC 1831 is at least several Myr old and our modelling makes use of the entire MS, we believe that these model uncertainties tend to be of smaller importance in our analysis.

The free model parameters are: metallicity ($`Z`$), age ($`\mathrm{log}(\tau )`$), intrinsic distance modulus $`(mM)_0`$, reddening (E(B-V)), PDMF slope ($`\alpha =\mathrm{d}\mathrm{log}\mathrm{\Phi }(M)/\mathrm{d}\mathrm{log}M`$) and unresolved binary fraction ($`f_{\mathrm{bin}}`$). The PDMF was considered as a power-law ($`\mathrm{\Phi }(M)=\mathrm{\Phi }_0M^\alpha =\mathrm{d}N/\mathrm{d}M`$), where the only free parameter is the slope $`\alpha `$. $`f_{\mathrm{bin}}`$ is defined as the systemic binary fraction, $`f_{\mathrm{bin}}=N_{\mathrm{bin}}/(N_{\mathrm{bin}}+N_{\mathrm{sing}})`$, where $`N_{\mathrm{bin}}`$ and $`N_{\mathrm{sing}}`$ are, respectively, the number of unresolved pairs and single stars.

The process of generating artificial stars works as follows:

(1) we fix $`Z`$ and $`\mathrm{log}(\tau )`$ for the stellar population by means of a chosen isochrone;

(2) we randomly draw a stellar mass according to the PDMF and get the absolute magnitudes in the two desired filters through the mass-luminosity relation given by the isochrone;

(3) for $`f_{\mathrm{bin}}`$ randomly chosen cases, we repeat step (2), representing a companion star in a binary system, and combine the two luminosities in both filters;

(4) we apply the intrinsic distance modulus $`(mM)_0`$ and reddening vector($`A_\mathrm{V},E(BV)`$) to the system, defining its theoretical CMD position. For this purpose, we use $`R_\mathrm{V}=A_\mathrm{V}/E(BV)=3.1`$ and the photometric transformation to the vegamag WFPC2 system according to Holtzman et al. 1995a;

(5) we introduce the photometric uncertainties by spreading the star with a Gaussian distribution of errors with $`\sigma _{555}`$ and $`\sigma _{814}`$ as empirically determined (see Sect. 2.1). This yields observational versions of the magnitude and colour;

(6) finally, we verify if the star is inside the MS observational ranges in magnitude and colour defined for the data and throw it away if it is not.

For each model realization we generate the same number of MS stars as observed in the real CMD, corrected for sample incompleteness and field contamination, and located inside the $`19.0m_{555}23.5`$ range. This range in apparent magnitude corresponds to $`0.5\stackrel{<}{}M_{555}\stackrel{<}{}5.0`$ and $`0.9\stackrel{<}{}M\stackrel{<}{}2.3`$ $`M_{}`$.

The best models are chosen by a direct comparison of the observed CMD with the artificial ones. The statistical tools used in this comparison are presented in the next section. The model vs. observation comparison strategy is as follows:

(1) the global parameters for the cluster, $`\mathrm{log}(\tau )`$, $`Z`$, $`E(BV)`$ and $`(mM)_0`$, are determined using the CMD of the central cluster region ($`R30`$ arcsec, where $`R`$ is the projected distance from the cluster’s centre), where field contamination and statistical noise are minimized (see Table 1);

(2) for the best combinations of the global parameters, the position dependent parameters $`\alpha `$ and $`f_{\mathrm{bin}}`$ are then derived using the CMDs in concentric rings of variable radii.

Table 1 shows important information about the cluster regions used in this modelling process. The first column gives the range in $`R`$, Col. (2) the original number of CMD stars in the $`19.0m_{555}23.5`$ range, Col. (3) the number of stars in the completeness corrected CMD (see Sect. 2.3.1), Col. (4) the number of assigned field stars (Sect. 2.3.2) and, finally, Col. (5) the number of stars in the CMD used in the modelling process. Notice that, not unexpectedly, field contamination becomes a serious issue for the outermost region, since statistical removal of field stars reduces the CMD numbers by about 40%. On the other hand, crowding in the central regions yields larger photometric incompleteness, as reflected by the clear increase in the number of stars between Cols. (2) and (3).

This modelling strategy allows an efficient and systematic search of best fit models in a 6-dimensional parameter space and makes use of the 2-dimensional information provided by the CMD data. Furthermore, this strategy naturally splits the parameters into those that define the position of the MS in the CMD plane (the global ones) and those that influence the way stars are distributed within the MS locus (the position dependent ones).

Fig. 5 shows four CMDs for the central cluster region. The one in the bottom right (panel d) is the real data, whereas the other three are realizations from different models generated by the modelling process described above. The input model parameters are shown in each panel. The three model CMDs do in fact look different, their MS ridge lines having different shapes and occupying different positions along the CMD plane.

### 3.2 Statistical tools

One of the main goals of this paper is to establish an objective comparison method between models and data. This required developing and applying statistical criteria that discriminate the model CMDs that most adequately reproduce the observed one. Ideally these comparison criteria should be both simple and easy to implement but yet make use of as much information provided by the bidimensional colour-magnitude plane as possible. We stress that these methods, in principle, are not restricted only to CMD analyses, but may be applied to the comparison of any two bidimensional distributions of points. In similar context as in this work, statistical techniques of CMDs comparison have been successfully applied by Kerber et al. (2001), Lastennet & Valls-Gabaud (1999), Saha (1998), Hernandez et al. (1999), allowing a model reconstruction of the main CMD features displayed by the component stars within a galaxy or a cluster.

Three statistics of CMD comparison were considered in this paper: $`S^2`$, $`PSS`$ and $`L`$. The first two are explained in more detail by Kerber et al. (2001). $`S^2`$ is essentially a dispersion between model and data points in the CMD plane; $`PSS`$ is proportional to the joint probability that the two CMDs being compared are drawn from the same population.

As for the $`L`$ statistics, it is an empirical version of the likelihood statistics described and used by Hernandez et al. (1999). For each model, 300 realizations with the same number of artificial stars as the real data (hereafter $`N_{\mathrm{obs}}`$) were generated. Dividing the CMD plane into $`N_\mathrm{b}`$ boxes, the model probability of one star, randomly chosen from any of these 300 realizations, to be inside the $`k^{th}`$ box is given by $`p_k=N_k/(300N_{\mathrm{obs}})`$. $`N_k`$ is the sum of stars from all realizations which fall in the $`k^{th}`$ box. Thus, $`p_k`$ may be interpreted as a probability function along the CMD plane.

The likelihood L of a given model is then defined as

$$L=\underset{i=1}{\overset{N_{\mathrm{obs}}}{}}p_{\mathrm{obs},i}$$

or

$$\mathrm{log}L=\underset{i=1}{\overset{N_{\mathrm{obs}}}{}}\mathrm{log}(p_{\mathrm{obs},i}),$$

where the product and sum are over the $`N_{\mathrm{obs}}`$ observed stars, and $`p_{\mathrm{obs},i}`$ is the model probability function evaluated at the CMD position of the $`i^{th}`$ observed star.

So, for each model we have three distinct statistical values. In order to establish the best models we build diagnostic diagrams (hereafter DDs), which are planes where we confront any two of these statistics. The best models will naturally have large $`L`$ and $`PSS`$, and small $`S^2`$ values. The method was tested by means of control experiments, that are shown in Sect. 4.

### 3.3 The model grids

The model input values for $`\mathrm{log}(\tau )`$, $`Z`$, $`(mM)_0`$, $`E(BV)`$ were chosen in order to bracket those found in the literature. In this respect, the web page www.ast.cam.ac.uk/STELLARPOPS/LMCdatabase by Richard de Grijs was very useful as it includes a very large compilation of parameter values and references on the 8 LMC rich clusters imaged by the GO7307 project.

In accordance with this database, we set the range of possible values for each physical parameter and created a regular model grid within this range. We expect this systematic grid to prevent biases in the parameter values determination.

Using the cluster central region ($`0R30`$ arcsec), we explored the following space defined by the global parameters:

$`Z=0.004,0.008,0.012`$

$`\mathrm{log}(\tau /\mathrm{yr})=8.70,8.75,8.80`$

$`(mM)_0=18.30,18.40,18.50,18.60,18.70`$

$`E(BV)=0.00,0.02,0.04,0.06,0.08,0.10`$

The position dependent parameters were kept fixed as $`\alpha =2.30`$ and $`f_{\mathrm{bin}}=0.30`$. Therefore this initial grid has 270 models.

Once the CMD comparison statistics defined in Sect. 3.2 are computed for each model, the DDs are built and the best models are identified, we investigate the dynamically affected, position dependent parameters by generating artificial CMDs to be compared to the observed CMDs within the concentric rings. In this second step we explored the following bidimensional parameter space:

$`\alpha =1.40,1.70,2.00,2.30,2.60`$

$`f_{\mathrm{bin}}=0.25,0.50,0.75`$

A total of 180 models were built for each ring, the only difference between one ring and another being the number of artificial stars.

## 4 Control Experiments

As mentioned in Sect. 3.2, we tested the validity of our statistical methods using control experiments. These experiments consist of drawing one realization of some specific model and calling it the “observed CMD”. We then verify if the DDs recover the generating model (hereafter input model) as the best one describing the “observed CMD”.

Figs. 6 and 7 show the results of a control experiment involving the 270 models to be latter applied to the NGC 1831 central region. All panels show DDs of $`\mathrm{log}S^2`$ vs. $`\mathrm{log}L`$, each point in the DD representing a particular model. The panels on the right are blow-ups of those on the left, showing in detail the region where the best models are located; this region corresponds to larger $`\mathrm{log}L`$ and smaller $`\mathrm{log}S^2`$ values. The different symbols in a single panel are coded according to the values of one of the four global model parameters ($`Z`$, $`\mathrm{log}(\tau )`$, $`(mM)_0`$ and $`E(BV)`$), therefore allowing the effect of varying each parameter to be separately assessed. Notice that the figures do not show the entire model grid in order to avoid cluttering. The grid regions discarded from the DDs are those of systematically high $`\mathrm{log}S^2`$ and low $`\mathrm{log}L`$ values. The “observed CMD” was taken to be a realization of the model with $`Z=0.008`$, $`\mathrm{log}(\tau )=8.75`$, $`(mM)_0=18.50`$ and $`E(BV)=0.06`$. This input model is shown as the large star in the blowup panels.

A tight correlation between the two statistics is clearly seen in all panels, adding reliability and stability to the choice of the best fitting models. The control experiments also show that one is in fact capable of recovering the input model based on the values of the statistics used: it is the model with largest $`\mathrm{log}L`$ and smallest $`\mathrm{log}S^2`$, as ideally we would expect.

Another important result from Figs. 6 and 7 is the combining and/or canceling effect of some global parameters, yielding models of comparable quality. As an example, the effects of metallicity $`Z`$ and reddening $`E(BV)`$ tend to cancel each other. Some models with $`Z`$ ($`E(BV)`$) lower (higher) than the input value, along with some high $`Z`$ and low $`E(BV)`$ ones, rank among the best models in the DDs. This degeneracy in the DDs is not surprising since the effect of increasing $`Z`$ is to make stars redder and fainter at a given mass, roughly opposite to the effect of decreasing $`E(BV)`$.

Fig. 8 presents similar DDs as in Figs. 6 and 7, but shows the model grid to be applied to the concentric regions of NGC 1831 (in a total of 180 models). The symbols now indicate different values of the PDMF slope $`\alpha `$ (panels (a) and (b)) and binary fraction $`f_{\mathrm{bin}}`$ (panels (c) and (d)). As before, the panels on the right are blow-ups of the ones on the left, showing the best models only. The input model (large star), in this case, is the one with $`\alpha =2.00`$, $`f_{\mathrm{bin}}=0.50`$ (and $`Z=0.012`$, $`\mathrm{log}(\tau )=8.75`$, $`(mM)_0=18.60`$, $`E(BV)=0.02`$). It is again located at the extreme upper left corner of the DD, confirming the applicability of our statistical approach.

However, the panels on the right show that the three best models have the same $`\alpha `$ ($`=2.00`$) but different $`f_{\mathrm{bin}}`$ values, revealing a larger difficulty in determining the latter than the former. This occurs because the effect caused by binaries is of the same order as or smaller than the photometric uncertainties in our WFPC2 CMDs. Consequently, as will be discussed in Sect. 5.2, the $`f_{\mathrm{bin}}`$ determination by means of our CMDs becomes a difficult task.

Notice that, in comparison with the DDs for the central region, the DDs in Fig. 8 present a larger dispersion. This is caused by the much smaller number of stars used in the set of models with varying $`\alpha `$ and $`f_{\mathrm{bin}}`$; in this case, all model CMDs (including the “observed” one) have 1154 stars, therefore mimicking the situation of the second concentric region to be studied in NGC 1831 (see Table 1).

Fig. 9 shows DDs involving the $`PSS`$ statistics. The upper panels show $`PSS`$ vs. $`\mathrm{log}L`$ and $`PSS`$ vs. $`\mathrm{log}S^2`$ plots for the central region. The lower panels show the same plots for one of the concentric rings. The same models as in Figs. 6, 7, and 8 are depicted but without symbol coding as a function of parameter values. The correlation between $`PSS`$ and the other statistics is again quite tight. In fact, the results based on the DDs are insensitive to the particular choice of statistics to be plotted. This is a very important result, since it further enhances the reliability of our statistical CMD modelling techniques. For the sake of simplicity, we hereafter restrict ourselves to $`\mathrm{log}L`$ vs. $`\mathrm{log}S^2`$ DDs only.

## 5 Results

### 5.1 Central region

In this section we model the CMD of stars belonging to the region within $`R30`$ arcsec from the centre of NGC 1831. In accordance with Santiago et al. (2001), this corresponds roughly to $`R\stackrel{<}{}7`$ pc or within 2 half-light radii. As indicated in Table 1, there are 2221 stars in this region in the magnitude range $`19.0m_{555}23.5`$. As also mentioned earlier, this central cluster region suffers from insignificant contamination by LMC field stars. On the other hand, due the high stellar density, incompleteness effects become more important. The chosen faint magnitude cut-off represents the magnitude at which completeness falls at $`50\%`$, in an attempt to reduce the relevance of this effect on our results.

The DDs for the model vs. real CMD comparison are presented in Figs. 10 and 11. As in the control experiments, the correlation between $`\mathrm{log}L`$ and $`\mathrm{log}S^2`$ is quite noticeable, the best models being again at the upper left region in the panels. These figures follow the same conventions and notations as Figs. 6 and 7, therefore allowing us to assess the effect of varying each global cluster parameter separately. For example, panels 10(a) and 10(b) clearly show that the best models have $`Z=0.012`$. The effects of the other parameters are not as striking as in the case of metallicity. Yet, panels 10(c) and 10(d) favour an age in the range $`8.75\mathrm{log}(\tau )8.80`$. Likewise, the three lowest values of intrinsic distance modulus (two of which are not even shown in the figure) can essentially be ruled out (panels 11(a) and 11(b)), placing NGC 1831 near or beyond the distance to the LMC centre. Finally, $`E(BV)0.02`$ is favoured by our modelling approach (panels 11(c) and 11(d)). Notice that this latter parameter is again anti-correlated with $`Z`$, in the sense that the models in the blowup panels with $`Z=0.008`$ (circles in panel 10(b)) have larger $`E(BV)`$ values (circles or crosses in panel 11(d)).

These best choices of the global parameters, as discussed in Sect. 3.1, besides being useful constraints by themselves, serve as model inputs to the position dependent parameters $`\alpha `$ and $`f_{\mathrm{bin}}`$, whose modelling is based on the concentric regions.

### 5.2 Concentric regions

We now analyze the NGC 1831 CMDs in the concentric regions listed in Table 1. As just mentioned, the free model parameters in this case are the PDMF slope $`\alpha `$ and the unresolved binary fraction $`f_{\mathrm{bin}}`$. For all regions we corrected the observed CMD to the effects of sample incompleteness and field star contamination. Fig. 12 shows their final CMDs. As mentioned before, Table 1 shows the number of points in all the steps along the data correction process.

In the modelling, we apply the same model grid for every region, only changing the number of artificial stars generated. The results from the previous section were used to restrict the possible values of the cluster s global parameters, allowing us to investigate the dynamically dependent parameters for the system.

Figs. 13 and 14 present the DDs resulting from model vs. data CMD comparison in the four regions. The different symbols in each panel correspond to different values of $`\alpha `$. As usual, the panels on the right show the best models in detail. The results for the two innermost regions are presented in Fig. 13. Panels (a) and (b) indicate that the best models for innermost region ($`R15`$ arcsec) have $`1.40\alpha 1.70`$. In the ring $`15<R30`$ arcsec, shown in panels (c) and (d) of the same figure, the best models present similar $`\alpha `$ values: $`1.40\alpha \stackrel{<}{}1.70`$.

In Fig. 14 we have the results for the two outermost rings. Panels (a) and (b) show the DDs for the stars with $`30<R60`$ arcsec. Now, there is evidence for a sharp change in the PDMF s slope: the best models have $`2.30\alpha 2.60`$. Finally, the DD for the most peripheric ring has a large dispersion (panels 14(c) and 14(d) ). This spread is also seen in the best PDMF slope: values in the range $`2.00\stackrel{<}{}\alpha 2.60`$ are seen in the upper left corner of the DD. This may reflect large uncertainties in the field star subtraction, since it leads to a large reduction of CMD points in this ring.

For each ring, we assign a representative $`\alpha `$ and associated uncertainty using the 10 models with the largest values of likelihood $`L`$. The final slope is taken to be the weighted average value among these best models, and its associated uncertainty is the dispersion around the average. The weight assigned to each model was the inverse of the difference in $`\mathrm{log}L`$ between the observed CMD and a CMD from a typical model realization. This difference is a measure of the discrepancy between model and observed CMDs. Table 2 lists the final $`\alpha `$ values and uncertainties for each ring (Cols. 3 and 4). The first 2 columns in the table show the regions in arcsec and parsecs (assuming $`(mM)_0=18.61`$ for NGC 1831 as adopted in Sect. 6), respectively.

As for the fraction of unresolved binaries, $`f_{\mathrm{bin}}`$, the results, as anticipated, are not conclusive in any ring, given the CMD spread. A more precise treatment of photometric errors may help constrain this parameter.

## 6 Summary and Conclusions

In this work we analyzed a deep CMD of NGC 1831, a rich LMC cluster, obtained from HST/WFPC2 images in the F555W and F814W filters. We inferred physical parameters such as metallicity, age, intrinsic distance modulus, reddening and PDMF slopes by comparing the cluster CMD with artificial ones. We presented in detail the techniques used to build the model CMDs, which use these parameters as model input and take into account observational uncertainties as in the real data. The parameter space explored by our regular model grids bracketed the values found in the literature.

The model vs. data comparison required correcting the observed CMD for selection effects caused by photometric incompleteness and CMD contamination by stars belonging to the LMC field. We also presented the statistical techniques used to compare the real CMD, corrected for the aforementioned selection effects, to the artificial ones. These statistical tools allowed us to discriminate the models that best reproduce the data. They are based on simple and objective statistical quantities and make use of the full bidimensional distribution of points in the CMD.

The best parameter values inferred for NGC 1831 are in the ranges $`8.75\mathrm{log}(\tau )8.80`$, $`18.50(mM)_018.70`$, $`0.00E(BV)0.02`$. Using the weighted average value among the 10 best models, as described in Sect. 5.2, we have $`\tau =588\pm 43`$ Myr, $`(mM)_0=18.61\pm 0.07`$, $`E(BV)=0.013\pm 0.015`$. As for the metallicity, all 10 best models have $`Z=0.012`$.

As discussed in Sect. 4, there is a coupling between reddening and metallicity, and therefore the determination of the former limits the values of the latter. As a consequence, the derived high metallicity for NGC 1831 brings about a low reddening determination.

The combined values of the global parameters obtained in this work for NGC 1831 suggest that this cluster is: a) metal-richer and older than in most previous estimates; b) placed near or beyond the distance to the LMC centre; c) found in a low reddening region.

Our estimated age and metallicity for NGC 1831 are in perfect accordance with the age-metallicity relation for LMC clusters obtained by Olszewski et al. (1991).

The best models for the global parameters were then used to build artificial CMDs of regions separated according to distance from the cluster centre and with varying values of the PDMF slope and fraction of unresolved binaries. For the PDMF slope $`\alpha `$ our statistical modelling shows that significant mass segregation exists in NGC 1831: for the central regions (out to 30 arcsec $`7.6`$ pc), we derive $`\alpha 1.70`$, whereas $`\alpha \stackrel{>}{}2.20`$ for the outer regions. As for binarism, our results were not as conclusive. One explanation is that the spread in the CMD caused by unresolved binaries is of similar or smaller amplitude than the empirically derived photometric uncertainties in the data. Therefore, this parameter is highly sensitive to the adopted photometric errors, rendering its realistic estimate a task for yet deeper CMDs or for data for which photometric uncertainties may be better estimated and modelled.

In a preliminary analysis of these data, Santiago et al. (2001) observed the effect of mass segregation in NGC 1831 with a steepening of the LF slope as a function of distance from the centre. However, no PDMF was derived. A global PDMF slope for NGC 1831 is presented by Mateo (1988), through the conversion of the luminosity function down to $`V23`$ ($`M\stackrel{>}{}M_{}`$) into a PDMF. The resulting slope is $`\alpha 4.0`$, therefore considerably steeper than our position dependent ones. Global slope values of other clusters from Mateo (1988) are in the range $`2.5\stackrel{<}{}\alpha \stackrel{<}{}4.6`$. In contrast, Elson et al. (1989) find much shallower slopes $`0.8\stackrel{<}{}\alpha \stackrel{<}{}1.8`$ for another sample of rich LMC clusters. These earlier results are based on ground-based observations, for which the effects of crowding are much more serious than in the present work.

Our determined position dependent PDMFs constrain the $`\alpha `$ values within $`0.2`$. These results can be used along with N-body simulations in order to recover the initial conditions, in particular the cluster IMF. We are applying the techniques shown in this paper to the others LMC clusters imaged with HST/WFPC2 as part of the GO7307 project in order to infer the same physical information as for NGC 1831. These future results can be very useful in investigations on the IMF universality.

###### Acknowledgements.

We acknowledge CNPq and PRONEX/FINEP 76.97.1003.00 for partially supporting this work. We are grateful to R. de Grijs, S. Beaulieu, R. Johnson, G. Gilmore for useful discussions. BXS is, as always, in debt with the late Becky Elson for all he learned from her.