text:
# Constraints on Thick Disc and Halo parameters from HST photometry of field stars in the Galaxy

## 1 Introduction

Our Milky Way galaxy is known to have a large disk with spiral structure, typical of luminous galaxies of Sb type. It has also been known for a long time to contain a population of old ($`\tau `$ $`\stackrel{>}{}`$$`10`$ Gyrs) stars distributed in an spheroidal component. More recent models of its structure have included at least a third component, of intermediate properties, often called the thick disc  (Gilmore & Reid 1983). Our understanding of the spatial, kinematical and chemical structure of the Galaxy has continuously improved in more recent years as a result from the quest for a more quantitative and detailed picture of how the Galaxy (and galaxies in general) formed and evolved. However, despite the enlarged and improved stellar samples with accurate photometric, parallax and velocity measurements, several issues concerning the structure, origin and evolution of the main components and populations in the Galaxy still remain to be settled (Gilmore et al. 1995, Majewski et al. 1996, Norris 1999).

The very existence of the thick disc as a discrete component, with stars kinematically and chemically distinct from those of the thin disk and halo, is subject to controversy (Norris 1987, Carney et al. 1989, Reid 1998). There is little doubt, however, that models with one planar component, with a single density profile, do not provide a suitable description of the structure of the Galaxy, since such models do not successfully fit recent star count data (Santiago et al. 1996a, Buser et al. 1999).

The structure of the thick disc is generally described as a double exponential, with horizontal and vertical scales in the ranges $`r_02.54.0`$ kpc and $`600`$ pc $`\stackrel{<}{}`$$`z_0`$ $`\stackrel{<}{}`$$`1600`$ pc, respectively (Chen 1996, Robin et al. 1996, Buser et al. 1999, Norris 1999). The wide range of scale height values reflects in part the difficult task of defining and separating thick disc stars from those of other structural components. The value of $`z_0`$ for the thick disc is anti-correlated with the local density of intermediate population stars, $`n_{0D}`$, which normalizes the assumed density profile: models with larger $`z_0`$ usually have smaller $`n_{0D}`$ values and vice versa.

The spatial structure, stellar luminosity distribution and shape of the outer regions of the spheroid, the stellar halo, have also been under closer scrutiny in recent years (Wetterer & McGraw 1996, Gould et al. 1998, Stetson et al. 1999, Moore et al. 1999, Samurovic et al. 1999). Largely due to Hubble Space Telescope (HST) observations, it is now known that the stellar luminosity function (LF) both of globular clusters and field disk stars decreases at the low-luminosity end (De Marchi & Paresce 1995a,b; Elson et al. 1995; Santiago et al. 1996a,b; Gould et al. 1996,1997, Piotto et al. 1996, 1997, Mendez & Guzman 1998). However, there is no consensus yet about the shape of the stellar LF at the faint end for other components (Dahn et al. 1995, Gould et al. 1998). Also, the flatenning of the stellar halo and the density profile function that best describes it are subject to uncertainty.

Improved observational data will more efficiently constrain the models if use is made of objective and statistically sound means of comparing model predictions to observations. Recent efforts in the direction of efficiently modeling observed colour-magnitude diagrams (CMDs) or other N-dimensional spaces of observables have been developed (Saha 1998, Hernandez et al. 1999, Lastennet & Valls-Gabaud 1999, Stetson et al. 1999). These methods usually try to make use of non-parametric statistics and a minimum amount of initial assumptions.

In this paper we analyze a composite CMD of 32 stellar fields towards various directions in the Galaxy, all imaged with HST’s WFPC2. We apply simple and objective statistical tools to compare the observed CMD with those generated using structural models of the Galaxy. In Sect. 2 we discuss the data used, whereas in Sect. 3 we present the statistical methods and models. In Sect. 4 we discuss the models found to best reproduce the data and the consistency among the different statistics used. Finally, in Sect. 5 we present our final conclusions and future perspectives.

## 2 The WFPC2 data

Our WFPC2 fields are part of the Medium Deep Survey (MDS) database. The data were extracted from 32 deep fields, 17 of which, mostly at high Galactic latitudes, were analyzed by Santiago et al. (1996a). The new 15 fields are listed in table 1. Their MDS id is given, along with their Galactic longitude and latitude ($`l_{II}`$ and $`b_{II}`$, respectively), total number of exposures and exposure times and the number of stars found in each. All fields have been observed with the $`F606W`$ and $`F814W`$ filters and have at least 2 exposures per filter. Several of the 15 new fields have $`|b_{II}|`$ values in the range $`[15^{},30^{}]`$, thus containing a larger fraction of thick disc stars than those studied by Santiago et al. (1996a).

Sample definition, photometry, instrumental and completeness corrections, and calibration to the standard system for the new fields were carried out in a similar fashion as in Santiago et al. (1996a). In brief, the DAOPHOT package within IRAF was used to automatically make an object list (whose peak intensities were $`5\sigma `$ above sky background) and to measure magnitudes both using aperture photometry and point spread function (psf) modeling. A single psf template was built using a compilation of very high signal-to-noise isolated stars from the different fields. Star/galaxy separation proceeded by applying cuts in the space of parameters output by the psf fitting task ALLSTAR. Fig. 1 shows the cuts applied to the $`\delta mvsI_{814}`$ (panel a) and $`\chi ^2vsI_{814}`$ (panel b) relations in the field uo50. From these panels we also determined cut-off magnitudes both in the bright end (due to saturation, where $`\chi ^2`$ sharply increases) and in the faint end (where disentangling stars from other sources becomes impossible). For the fields with too few stars for that to be done, we used the cut-off magnitudes from the more crowded fields, scaled up or down according to exposure time.

Photometric corrections for limited aperture and charge-transfer effect (CTE), as prescribed by Holtzman et al. (1995a,b) were applied to the data. The magnitudes and colours were converted to the standard system, also using the transformations listed by those authors. Finally, completeness levels were obtained as a function of the $`F814W`$ magnitudes and used to compensate for the loss of stars in the faintest magnitude bins. Typically, completeness started to drop abruptly from 100% at $`I_{814W}24`$

In Fig. 2 we show the joint CMD for all 32 WFPC2 fields. It contains over 1000 stars distributed within a large range of magnitudes and colours. The magnitudes and colours were obtained by aperture photometry, are in the standard photometric system and corrected for the several effects mentioned above. Notice, however, that no correction to redenning was applied to the data. We preferred to incorporate redenning effects to our model predictions, described in the next section.

## 3 Models for the Galaxy and Statistical Techniques

### 3.1 Modeling the Galaxy

Our main goal is to compare the observed CMD of HST/WFPC2 stars with theoretical CMDs based on models for the structure of the Galaxy. The model number counts as a function of colour and magnitude can be obtained from numerical integration of the Fundamental Equation of Stellar Statistics. The main ingredients entering this equation are the stellar luminosity function, $`\mathrm{\Phi }(M)`$ (LF), and the density distribution $`\rho (𝐫)`$ of stars for each structural component in the model. We consider models with 3 major components: a thin disc and a thick disc, both with double exponential density profiles, and a spheroidal or halo component, with a deprojected de Vaucouleurs profile. The main parameters governing these profiles are the horizontal and vertical exponential scales, $`r_0`$ and $`z_0`$ respectively, and the effective de Vaucouleurs radius, $`R_e`$, for the halo. Also relevant are the normalization factors, $`n_0`$, for the different components, which are usually expressed as percentage values of the local spatial density of stars.

The LF used in our modeling is the one originally from Wielen et al. (1983) up to $`M_V=12`$. For lower luminosities we treat it as a power law in luminosity, $`log\mathrm{\Phi }(M_V)\alpha M_V`$ with slope $`\alpha `$ considered as a free parameter. In our data $`vs`$ model comparisons (see below), we also allow for variations in $`z_0`$, $`r_0`$, $`R_e`$ and in the normalization of the thick disc ($`n_{0D}`$) and spheroidal components ($`n_{0H}`$). For more on the models, profiles and LFs we refer to Reid & Majewski (1993), Majewski et al. (1996) and Santiago et al. (1996a).

The numerical integration of stellar statistics was carried out for all 32 fields, taking into account their variable available magnitude ranges (set by saturation and detection) and completeness functions. Galactic extinction was also incorporated to the models using the Burstein & Heiles (1982) E(B-V) maps. Conversion of extinction vectors and model predictions from BV to VI was carried out in the same way as in Santiago et al. (1996a). The output from these integrals are tables with $`(VI)`$ colour and $`I`$ magnitude counts for each field which can then be added together to yield a composite model CMD.

We explored an extensive grid of models, covering large regions of the parameter spaces formed by the thick disc and halo structural parameters and LF slopes. Thick disc scale height and normalization values varied within the ranges $`500z_02000`$ pc and $`1\%n_{0D}8\%`$, respectively. But notice that these two parameters were paired so that their anti-correlation was respected. A total of 32 thick disc models resulted. For each one of them, we considered 4 choices of the horizontal thin disk e-folding length: $`r_0=2500,3000,3500,4000`$ pc, and 2 choices for the LF slope at its low-luminosity end ($`M_V12`$): $`\alpha =0.2`$ and $`\alpha =0`$, totalling 256 thick disc models.

As for the halo we used $`R_e=1500,3000,4500`$ pc for its effective radius and $`n_{0H}=0.1,0.125`$ and $`0.2\%`$ for the local normalization. Another halo parameter varied was its axial ratio, with values $`c/a=0.5,0.65,0.8,1.0`$, yielding a total of 36 stellar halo models.

This extensive parameters grid was defined so as to include the variety of parameter values often quoted in the literature (Norris 1999) and to accommodate most thick disc and halo parameters whose values are still uncertain.

### 3.2 Statistical Tools and Techniques

Fig. 3 shows a comparison of a theoretical CMD, whose main parameters are indicated, with the observed one. The simulated CMD was built from the expected model number counts as a function of colour and magnitude using Monte-Carlo techniques: we randomly threw points on the CMD plane with probabilities proportional to the expected number in each CMD position. The total number of points in the simulated CMD was a Poisson deviate of the total model counts.

Our main challenge is to find an objective and efficient statistical tool to quantify how different or similar these two two-dimensional distributions of points are. A visual comparison reveals regions in the CMD, such as $`0VI0.8`$ and $`18I22`$, where there is a clear excess of model points. There are other regions (like $`0VI0.8`$ and $`23.5I24.5`$, for instance) with more observed than model stars. The different models in the grid show a smooth and continuous variation in their CMDs, making it hard to visually select best fitting models.

We here propose two statistics which are objective and simple tools for comparing two distributions of points on a plane. The first is simply a measure of the dispersion between number counts within different data bins. Let $`m_i`$ and $`o_i`$ be the number of model and observed stars in the $`i^{th}`$ of $`N`$ bins regularly spaced on the CMD plane; we then define the dispersion as:

$$s^2=\underset{i=1}{\overset{N}{}}(m_io_i)^2$$

The second statistics was used by Santiago & Strauss (1992) in a different context, to study spatial segregation of galaxies. It also compares the number of CMD points found in different bins. Let $`p_i`$ be the percentile position of $`|o_im_i|`$ within the distribution of $`|m_im_{ij}|`$ ($`j=1,N_{real}`$) values, where $`m_{ij}`$ are Poisson realizations of the expected model counts $`m_i`$ in the $`i^{th}`$ bin. We then define

$$pss=\underset{i=1}{\overset{N}{}}log(100p_i(\%))$$

We typically compute $`pss`$ from a distribution of $`N_{real}=300`$ Poisson realizations of each model. Even though both statistics depend on data binning, we have observed that both are quite insensitive to the binning scheme. They are also applicable to datasets with varying number of points. This is in contrast with what we found for the W statistics (Saha 1998), whose results seem to be sensitive to the total number of points in both distributions being compared.

Clearly, the better the agreement between model and observed CMDs, the smaller $`s^2`$ and the larger $`pss`$ will be. We tested both statistics by means of controlled experiments, in which we generated a Monte Carlo realization from a particular model (as described in the beginning of this section) and compared it to all the models. Fig. 4 shows the result of a typical such experiment. It shows a diagnostic diagram of $`pssvss^2`$ values obtained by comparing the simulated “data” with the theoretical CMD for all 256 models in the thick disc grid described in the previous section. There is a clear anti-correlation between the 2 statistics, as expected.

The large dot in Fig. 4 corresponds to the particular model used to generate the “data” points. It is clearly among the extreme values of both statistics, again as expected. However, the fact that it is not the most extreme point in the upper left of Fig. 4 is a measure of the resolution provided by the diagnostic diagram used.

## 4 Results

### 4.1 Thick Disc

In this section we apply the $`pss`$ and $`s^2`$ statistics to the comparison between our observed CMD and the theoretical CMDs from the grid of models presented earlier. Our main tool is the $`pssvss^2`$ diagram. In computing $`pss`$ we ran $`N_{real}=300`$ Poisson realizations for each model in order to build the distribution of $`m_{ij},j=1,300`$ model number counts in each bin. We adopted a 14x10 binning covering the I vs (V-I) plane in the ranges $`18I25`$ and $`0(VI)4`$.

Fig. 5 shows the position of all 256 models we considered by varying the thick disc structural parameters $`z_0`$, $`n_{0D}`$ and $`r_0`$ and its stellar LF slope. There is again a clear correlation among the two statistics. Furthermore, the position within this diagram is correlated with the parameter values. The diagram splits into two branches, the upper one corresponding to models with $`z_01200`$ pc, which usually feature a larger than observed number of stars. The lower branch (small pss and $`s^2`$ values) contains models with a deficiency of stars relative to the observed sample. The best models in the upper left have $`800z_01200`$ pc and $`4n_{0D}8\%`$, with little or no dependence on the horizontal exponential scale, $`r_0`$.

These results are further confirmed by Fig. 6, where the $`pss`$ values (assuming a decreasing LF at the faint end) are shown for different locations within the $`z_0vsn_{0D}`$ plane for fixed $`r_0`$. The diagrams do not differ significantly, with the best $`pss`$ values located within the range $`800z_01200`$ pc and $`4n_{0D}8\%`$ for all panels. Note that the anti-correlation between these two structural parameters persists, a $`z_0800`$ pc with $`n_{0D}8\%`$ being essentially as good a model as one with $`z_01200`$ pc and $`n_{0D}4\%`$.

There is also a clear trend for models with decreasing faint end LF slope to be favoured by our observed CMD. In Fig. 7 we show the difference in $`pss`$ values between models with decreasing and flat LF slopes at the faint end, all the other parameters being the same. The large majority of the 128 such differences are positive, indicating that the data are best described by a decreasing LF. In fact, of the 35 models located closer to the upper left boundary of the $`pssvss^2`$ diagram, only 1 had a flat thick disc LF slope. A similar result of a decreasing faint end for the LF has also been found for the halo and/or disk stars by Bahcall et al. (1994), Santiago et al. (1996a), Gould et al. (1997) and Mendez & Guzman (1998) using different data sets.

The consistency of the results just shown was tested by adopting both coarser (70 bins) and narrower (280 bins) binning schemes, as well as different color ranges than the one quoted early in this section. There is no significant change in the $`pssvss^2`$ diagram or in the position within it of models with different parameters.

### 4.2 Halo

In order to increase the fraction of halo stars in our observed CMD, we have considered only the 18 WFPC2 fields with $`|b|40^{}`$ and which point away from the centre of the Galaxy. This reduced, higher latitude, sample contains over 200 stars. Their CMD was then compared to those generated by the 36 halo models discussed at the end of Sect. 3.1.

The $`pssvss^2`$ diagram for these 36 models is shown in Fig. 8. We again notice the anti-correlation between the two statistics. And again there is a correlation between position within this diagram and model parameters. The 2 models above the main branch have $`n_{0H}=0.2\%,c/a=1.0`$ and large values of $`R_e`$ ($`R_e=3000`$ or $`4500`$ bf pc). Most importantly, models with $`c/a=0.5`$ are concentrated towards the lower part of the diagram, therefore being ruled out. Due to the smallness of the sample, however, it was not possible to single out any other region of parameter space favoured by the data.

## 5 Summary and Conclusions

We carried out V and I photometry for over 1000 field stars in the Galaxy down to faint ($`I25`$) magnitudes and distributed over different directions on the sky imaged with deep HST/WFPC2 exposures. The CMD for the data, corrected for instrumental effects and calibrated, was compared with theoretical CMDs expected from a grid of models for the structure of the Milky-Way. The main grid parameters are the thick disc scale height $`z_0`$, its local density normalization, $`n_{0D}`$, the stellar halo effective radius $`R_e`$, normalization $`n_{0H}`$ and axial ratio $`c/a`$, as well as the slope of the luminosity function at its low-luminosity end ($`M_V12`$).

We used two simple statistical tools to compare model and observed CMDs, in an attempt to restrict the range of acceptable values for the structural and LF parameters of the halo and thick disc. Our analysis indicates acceptable models in the range $`800z_01200`$ pc, $`4n_{0D}8\%`$ with these two parameters being anti-correlated. These ranges in both parameters are substantially narrower than those quoted in a recent review by Norris (1999). In particular, our results do not accommodate the very thick disc of $`M>0.3M_{}`$ stars suggested by Gyuk & Gates (1999) in their reinterpretation of microlensing results. The $`z_0`$ range inferred from our work is also in disagreement with earlier works using ground-based star counts or samples of RR Lyrae stars, which have yielded $`z_01500pc`$ (Hartwick 1987, Rodgers 1991, Reid & Majewski 1993). On the other hand, as pointed out by Majewski (1993), introduction of metallicity gradients in our models could yield larger $`z_0`$ values. Another possibility, along similar lines, is that the thick disc may have a complex structure with possible subcomponents, such as the metal-weak thick disc (Chiba & Yoshii 1998, Norris 1999). Finally, our results are only marginally consistent with that of Robin et al. (1996), who find $`z_0=760\pm 50`$ pc and $`n_{0D}=5.6\pm 1\%`$. A better agreement is found with Buser et al. (1999), whose best model for the thick disc has $`n_{0D}=5.9\pm 3\%`$ and $`z_0=910\pm 300`$ pc.

As already established for disk stars by Santiago et al. (1996a), Gould et al. (1996,1997) and Mendez & Guzman (1998) (from studies of HST/WFPC2 fields), we find the thick disc LF to be decreasing beyond $`M_V12`$ with large confidence. This is also in agreement with most studies of globular cluster luminosity functions.

We also attempted to restrict the structural parameters of the halo using a smaller sample with $`|b|>40^{}`$. We were able to rule out flattened oblate halo models, with $`c/a`$ $`\stackrel{<}{}`$$`0.5`$, or models with both high normalization and high effective radius. The smallness of the sample, however, prevents further conclusions about halo structure.

A way forward will be to devise more and better statistical tools to be applied to larger stellar samples. Particularly useful in constraining the structure of the main components of the Galaxy will be the inclusion of future kinematic data to be obtained from new orbital astrometric missions such as SIM and GAIA (Majewski 2000).

###### Acknowledgements.

We are very grateful to David Valls-Gabaud for useful discussions and suggestions. We also thank Gerry Gilmore and Jim Lewis for the original software used for modeling the structure of the Galaxy. We also acknowledge CNPq and PRONEX/FINEP 76.97.1003.00 for partially supporting this work.