A MONTE-CARLO METHOD FOR ESTIMATING STELLAR PHOTOMETRIC METALLICITY DISTRIBUTIONS

The metallicity distribution function (MDF) for stars in the Milky Way is of great importance for revealing the chemical structure of Galactic halo and disk systems, and provides essential clues to the assembly and enrichment history of the Galaxy (Carollo et al. 2007, 2010; Peng et al. 2012, 2013; An et al. 2013, 2015). Many multi-fiber spectroscopic surveys, including the Sloan Digital Sky Survey (SDSS; York et al. 2000), the Radial Velocity Experiment (RAVE; Steinmetz et al. 2006), and the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST; Deng et al. 2012; Zhao et al. 2012; Liu et al. 2014), have obtained metallicity estimates and other stellar parameters for millions of stars from low- and medium-resolution spectra. However, compared with photometric data, the number of spectra within the limiting magnitudes of the surveys is still too small to provide a detailed chemical map of the Galaxy, even in the relatively local region. The use of photometric data, which is available for far more stars than spectroscopic data, breaks through this limitation.

Considering that the exhaustion of metals in a stellar atmosphere has a detectable effect on the emergent flux (Schwarzschild et al. 1955), particularly in the blue region where the density of metal absorption is highest, the combination of spectroscopic data and photometric data can be used to derive estimates of [Fe/H] (Allende Prieto et al. 2006, 2008; Lee et al. 2008a, 2008b). For example, Siegel et al. (2009) derived approximate stellar metallicities through the measurement of the ultraviolet excess based on UBV data in SA 141. Ivezić et al. (2008) employed SDSS data to derive a metallicity estimator from u − g and g − r colors, and successfully mapped the metallicity distribution of millions of F/G main-sequence stars within a distance of ∼8 kpc from the Sun. Peng et al. (2012) also used BATC survey data to estimate the stellar photometric metallicity distribution. Gu et al. (2015) obtained a metallicity estimator using SCUSS (Zou et al. 2015, 2016) data, which can be applied up to fainter magnitudes due to the use of more accurate SCUSS u-band measurements, and used this estimator to explore the metallicity of the Sagittarius stream in the South Galactic cap. Using a minimum χ² technique, Yuan et al. (2015) estimated photometric metallicities simultaneously using the dereddened colors u − g, g − r, r − i, and i − z from the SDSS and metallicity-dependent stellar loci. An et al. (2013) calibrated stellar isochrones to derive metallicities of individual stars with SDSS ugriz photometry. An et al. (2015) applied this method to recently improved co-adds of ugriz photometry for Stripe 82 from SDSS, including a factor of two more stars than in their previous effort. The new analysis revealed an MDF for halo stars between 5 and 10 kpc from the Sun with peak metallicities at [Fe/H] ∼ −1.4 and [Fe/H] ∼ −1.9, which the authors associated with the inner-halo and outer-halo populations of the Milky Way, respectively.

The photometric metallicity calibrations developed by previous works are characterized by their assignment of a star-by-star metallicity estimate based on its color indexes. This inevitably introduces error because a single star's metallicity is actually uncertain even when its color indexes are fixed, varying around the metallicity estimate to form a distribution. In addition, the calibration methods used by Ivezić et al. (2008) and Gu et al. (2015) yield poor results for very metal-rich or very metal-poor stars. However, in order to investigate the chemical structure of the Galactic stellar populations, we only require knowledge of the MDF for a large statistical sample of stars. Here, we develop a Monte-Carlo method to estimate the photometric metallicity distribution of a large number of stars with available SDSS photometry.

This paper is organized as follows. In Section 2, we provide a brief overview of the SDSS and its photometric data. Details of our Monte-Carlo-based photometric metallicity calibration are presented in Section 3. Section 4 presents a comparison between this method and the more traditional polynomial-fitting method. As an example, in Section 5, we apply this method to derive the MDF for Galactic stars and consider its variation with height above the plane. A brief summary is given in Section 6.

The SDSS is a large international collaboration project, and it has obtained deep, multi-color images covering more than one-quarter of the celestial sphere in the North Galactic cap, as well as a small (∼300 deg²) but much deeper survey in the South Galactic hemisphere (York et al. 2000). The SDSS used a dedicated 2.5 m telescope at Apache Point Observatory, New Mexico (Gunn et al. 2006). The flux densities are measured in five bands (u, g, r, i, z) with effective wavelengths of 3551, 4686, 6165, 7481, and 8931 Å, respectively. The 95% completeness limits of the images are 22.0, 22.2, 22.2, 21.3, and 20.5 for u, g, r, i, and z, respectively (Abazajian et al. 2004). The relative photometric calibration accuracy for u, g, r, i, and z are 2%, 1%, 1%, 1%, and 1%, respectively (Padmanabhan et al. 2008). Figure 1 shows the error of the u-, g-, and r-band magnitudes as a function of the g-band magnitudes. Other technical details about SDSS can be found on the SDSS website http://www.sdss3.org/, which also provides an interface for the public data access.

Zoom In Zoom Out Reset image size

The Sloan Extension for Galactic Understanding and Exploration (SEGUE-1) obtained spectra of nearly 230,000 unique stars over a range of spectral types to investigate Galactic structure. Building on this success, SEGUE-2 spectroscopically observed around 119,000 unique stars, focusing on the stellar halo of the Galaxy, including stars with distance from 10 to 60 kpc from the Sun. We employ the complete set of derived stellar parameters for SEGUE-1 and SEGUE-2 from the newest version of the SEGUE Stellar Parameter Pipeline (SSPP; Beers et al. 2006; Allende Prieto et al. 2008; Lee et al. 2008a, 2008b, 2011; Smolinski et al. 2011), including effective temperature, surface gravity, and metallicity (parameterized as [Fe/H]). A thorough overview of the SEGUE effort can be found in Yanny et al. (2009).

We use the adopted stellar atmospheric parameters from the SSPP listed in the sppParams table. After excluding the repeated stars surveyed on different plates, we obtain a sample of 366,676 stars with SDSS u, g, r, i, and z magnitudes, as well as stellar parameters. Most stars in the sample have metallicities in the range −2.5 ≤ [Fe/H] ≤ 0.0. In this study, in order to generate the [Fe/H] probability distribution from the colors, we select main-sequence stars, adopting the following selection criteria:

• 1.  
  14 < g < 19.5;
• 2.  
  $0.2\lt g-r\lt 0.8;$
• 3.  
  $0.6\lt u-g\lt 2.2;$
• 4.  
  main-sequence stars are selected by rejecting those objects at distances larger than 0.15 mag from the stellar locus described by the following equation (Jurić et al. 2008):

  $$\begin{eqnarray*}g-r & = & 1.39\{1-\mathrm{exp}[-4.9{(r-i)}^{3}-2.45{(r-i)}^{2}\\ & & -1.68(r-i)-0.05]\};\end{eqnarray*}$$

  and
• 5.  
  we further refine the selection of main-sequence stars by rejecting those objects at distances larger than 0.3 mag from the stellar locus described by the following equation (Jia et al. 2014):

  $$\begin{eqnarray*}&&u-g=\mathrm{exp}[-{(g-r)}^{2}+2.8(g-r)-1].\end{eqnarray*}$$

As an illustration, Figure 2 shows the two-color diagrams for r − i versus g − r and u − g versus g − r. Our final sample includes 268,029 sample stars. Throughout this paper, it is understood that magnitude and color have been corrected for extinction and reddening following Schlegel et al. (1998).

Zoom In Zoom Out Reset image size

Figure 3 shows the spectrum-determined [Fe/H] distribution of stars versus u − g color. These stars are selected within small color ranges around $g-r=0.3$, $g-r=0.4$, and $g-r=0.5$. The solid lines represent the photometric metallicity estimator derived by Ivezić et al. (2008) with g − r color specified by 0.3, 0.4, and 0.5, respectively. From inspection of this figure, these lines can only roughly describe the relation between [Fe/H] and u − g color. Those stars with specified u − g and g − r colors exhibit a metallicity distribution which may result from several other factors. We build the metallicity probability distribution from the scattered points in each specified color bin and then use a Monte-Carlo method to infer the MDF for large numbers of stars. That is to say, what we obtain from our photometric calibration is not a one-to-one function, but is instead a probability-governed one-to-many function.

Zoom In Zoom Out Reset image size

The Monte-Carlo method relies on repeated random sampling to obtain numerical results. As mentioned in the selection criteria, the two colors we employ are confined to $0.2\lt g-r\lt 0.8$ and $0.6\lt u-g\lt 2.2$. We divided the u − g versus g − r panel into 0.05 × 0.05 mag² bins and designate each square unit by an index computed in the following manner:

$$\begin{eqnarray*}\mathrm{index} & = & \mathrm{int}((u-g-0.6)/0.05)\ast 12\\ & & +\;\mathrm{int}((g-r-0.2)/0.05),\end{eqnarray*}$$

where the symbol int stands for the integer portion. We obtain 384 square units. For a given sample of stars, we also equally divided the value of [Fe/H] from −3.5 to 0.5 into 80 bins with a bin width of 0.05 dex. In this manner, we obtained an array of index ×[Fe/H]. Each array element records the number of stars whose colors and [Fe/H] match their corresponding positions based on the above-selected sample of 268,029 stars. The resulting array, with its array elements holding the pertinent information, is called the "seed" array. Each element of the seed array is denoted by value[index][i], where i ranges from 0 to 79. The maximum of value[index][i] for i in the range [0, 79] can easily be found and is denoted by max[index]. From this seed array, we can evaluate the MDF for the photometrically surveyed stars.

For a certain index (the u − g and g − r interval), we use the Monte-Carlo method to generate a random number sequence whose probability distribution is determined by the [Fe/H] distribution recorded in the seed array. Suppose X and Y are two stochastic variables which can be assigned a value by a random number generating function rand_X() and rand_Y(), respectively. The function rand_X() is modulated to generate a uniform-probability distributed random integer number from 0 to 79, and rand_Y() from 0 to max[index]. In each trial, we obtain a random number pair (X = rand_X(), Y = rand_Y()). When Y ≤ value[index][X], we record X as a useful value, but otherwise we discard it. Through numerous trials, we obtain a sequence of random numbers {X₁, X₂, X₃, ⋯} that follows the same probability distribution as recorded in the seed array. Then, we can transform the obtained random number sequence into metallicities, [Fe/H], by [Fe/H] = 0.05 ∗ index − 3.5 + 0.025. Here, because value[index][i] can be equal to zero for some i values, we can discard them and record the non-zero elements and their positions in a new array. Using this trick, we greatly improve the sampling efficiency.

For a large number of photometrically surveyed stars, we first count the number of stars corresponding to each array element, and then, following the above procedure, obtain the metallicities of the counted number for all array elements. The metallicities naturally form a distribution which we discuss below.

To test the feasibility and show the superiority of the method described above, we perform a comparison between this Monte-Carlo-based method and the third-order polynomial-fitting method presented in Ivezić et al. (2008), shown below:

$$\begin{eqnarray*}[\mathrm{Fe}/{\rm{H}}] & = & -4.37-8.56x+15.5y+23.5{x}^{2}\\ & & -39.0{xy}+20.5{y}^{2}-10.1{x}^{3}\\ & & +12.1{x}^{2}y+7.33{{xy}}^{2}-21.4{y}^{3},\end{eqnarray*}$$

where x = u − g for $(g-r)\lt 0.4$ and $x=(u-g)-2(g-r)+0.8$ for $(g-r)\gt 0.4$, $y=g-r$.

We select main-sequence stars covering the color range $0.2\lt g-r\lt 0.6$ from the spectroscopically surveyed stars and evaluate their photometric metallicity distribution. Figure 4 shows the derived MDFs from the three different methods. The top three panels show the MDFs for stars with $0.7\lt u-g\lt 1.0$ and the bottom three panels for stars with $1.3\lt u-g\lt 1.6$. The left panels present the spectrum-based MDF, which can be regarded as "ground truth". The middle panels show the photometric MDF determined by our Monte-Carlo method, and the right panels show the photometric MDF from Ivezić et al.'s (2008) model. The peak values in all of the histograms are normalized to one, with the actual peak values labeled. As shown in Figure 4, the photometric MDFs based on our Monte-Carlo method much more closely resemble the spectroscopic MDFs than those obtained from the Ivezić et al. (2008) approach. Not only are the peak values of the photometric MDFs determined by our Monte-Carlo method approximately equal to those of the spectrum-based MDF, but the wings of the two distributions have almost the same profile. By contrast, there exists some clear discrepancies between the photometric MDFs from Ivezić et al.'s (2008) model and the spectrum-based MDFs, particularly at the very metal-rich and very metal-poor ends. We thus believe that polynomial-based fitting methods are clearly inferior and should no longer be generally used. For application to a large number of stars, there is great advantage in using the Monte-Carlo method for the derivation of photometric MDFs. Note that if the number of stars used to evaluate the photometric MDFs is small, we can easily increase the number of desired random numbers X using a multiplicative factor in the Monte-Carlo method.

Zoom In Zoom Out Reset image size

Figure 5 shows the distribution of metallicity residuals between the spectrum-based metallicity and that determined by the Monte-Carlo method for individual stars. We determine the metallicity of a single star from the peak value of the metallicity distribution in each specific color bin. The stars are selected covering the color range $0.2\lt g-r\lt 0.6$. The distribution is roughly a Gaussian profile with a mean zero-point offset of −0.058 dex and a dispersion of 0.368 dex.

Zoom In Zoom Out Reset image size

In the conventional photometric metallicity calibration through polynomial fitting, the error mainly arises from two sources. One source is from the fitting method itself (illustrated in Figure 3), and the other is from the errors in the color indexes. In SDSS, the photometric error of the u-band magnitude is relatively large (shown in Figure 1), limiting the application range of the photometric metallicity estimator. In the Monte-Carlo calibration method, we suppose that the spectrum-based metallicity distribution in a given u − g, g − r color bin is fairly reliable, and thus by reproducing its distribution, we effectively eliminate the errors arising from the fitting method. The error introduced by fluctuations in the Monte-Carlo method can be eliminated by increasing the number of desired random numbers. In addition, the third-order polynomial model is determined by only 10 coefficients, while the model (seed array) construction involves many more variables. This, to some extent, guarantees the accuracy of the photometric metallicity distribution determined by the Monte-Carlo method. The deviation of the photometric MDF determined by the Monte-Carlo method mainly arises from the errors in the color indexes, especially when estimating the MDFs for faint stars. As shown in Figure 4, the calibration based on the Monte-Carlo method can be used to derive photometric MDFs with sufficient accuracy for most purposes.

As an example, we use the method introduced in this study to derive the photometric MDF for stars as a function of distance from the Galactic plane. The main-sequence stars with $0.2\lt g-r\lt 0.6$, 6 kpc < R < 9 kpc, and in different $| Z|$ intervals are selected (R, Z represent cylindrical Galactocentric coordinates, with the Sun's coordinate (R, Z, ϕ) = (8 kpc, 0, 0)). As shown in Figure 6, the photometric MDFs in the top two panels can be well-fit by two Gaussians with peaks at about −0.6 and −1.2, respectively: one associated with the disk system and the other with the halo (and metal-weak disk). The MDFs in the two bottom panels are found to be better-fit by three Gaussians with peaks at about −0.6, −1.4, and −1.9, respectively. The two lower peaks may be associated with inner-halo and outer-halo populations, respectively (Carollo et al. 2007, 2010; An et al. 2013, 2015). These four histograms clearly show that the number ratio between disk stars and halo stars decreases with vertical distance from the Galactic plane. In the metal-rich and metal-poor ends, the number of stars decreases gradually. The number ratios as a function of $| Z|$ between the disk and halo above the Galactic plane could be recalculated from the photometric metallicity distribution, and the result can be compared with that from the method of star counting. This will be presented in our future papers.

Zoom In Zoom Out Reset image size

This paper presents a new method to estimate the photometric MDF for main-sequence stars. The method is tested using a spectroscopic sample of stars from SDSS. Compared with the method from Ivezić et al. (2008), the current method is more accurate, particularly for very metal-rich and very metal-poor stars. This method is sufficiently accurate to be used to investigate the distribution and chemical structure of the Galactic stellar populations. At the same time, as an example, we also apply the method to main-sequence stars with $0.2\lt g-r\lt 0.6$, 6 kpc < R < 9 kpc, and different $| Z|$ intervals. The metallicity distribution of sample stars near the Galactic plane can be well fit by two Gaussians with peaks at about −0.6 and −1.2, respectively: one associated with the disk system and the other with the halo/metal-weak thick disk. However, the metallicity distribution of the sample stars far from the Galactic plane can be well fit by three Gaussians with peaks at −0.6, −1.4, and −1.9, which supports the existence of two components in the halo: the inner halo and the outer halo. The number ratio between disk stars and halo stars varies with vertical distance from the Galactic plane. In the metal-rich and metal-poor ends, the number of stars decreases gradually. With the advantage of the method introduced in this paper, we can better study the photometric MDF for different stellar populations in the Galaxy and provide detailed constraints on the Galactic chemical evolution, which we will consider in future papers.

We especially thank the referee for insightful comments and suggestions which have improved the paper significantly. This work was supported by a joint fund of Astronomy of the National Natural Science Foundation of China and the Chinese Academy of Science, under grant U1231113. This work was also by supported by the Special funds of cooperation between the Institute and the University of the Chinese Academy of Sciences. In addition, this work was supported by the National Natural Foundation of China (NSFC, No.11373033, No.11373035), and by the National Basic Research Program of China (973 Program) (No. 2014CB845702, No.2014CB845704, No.2013CB834902).

Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III website is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Institute de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.