Photometric survey databases are useful resources for studying stellar populations and structures of the Milky Way. The size of spectroscopic samples has grown rapidly in recent years, but photometric surveys still cover a significantly larger volume of space and thereby can provide the least biased sample of Galactic stars. Multiband observations are particularly useful, since they can be used to constrain fundamental stellar parameters, such as effective temperature (T eff ) and metallicity ([Fe/H]), with sufficient accuracy. Specifically, the overall shape of a spectral energy distribution as traced by multiband photometry depends on T eff , and ultraviolet excess provides information about a starʼs metallicity. When these data are combined with all-sky, high-precision astrometry from Gaia, they can provide rich information on chemical and kinematical properties of stars, as we demonstrated in a series of papers (An & Beers 2020, 2021a, 2021b, hereafter Papers I, II, and III, respectively).
To obtain a clear view of Galactic stellar populations, one needs to establish an accurate relationship between photometry and fundamental stellar parameters on an empirical basis or by using theoretical predictions. To take advantage of each method, we adopted a hybrid approach in previous papers in this series to calibrate theoretical isochrones of the main sequence using observations of well-studied Galactic globular and open clusters. The models were taken from YREC (Sills et al. 2000) and were combined with MARCS (Gustafsson et al. 2008) synthetic spectra, in order to convert T eff and luminosities into photometric colors and magnitudes. Differences of the models from observations typically amount to a few hundredths of a magnitude for warm (T eff > 5000 K) stars, but they become as large as a few tenths of a magnitude for cooler stars. The model offsets are also a function of metallicity, in that more metal-rich cluster sequences tend to exhibit larger color deviations.
The observed offsets from the models are systematic in nature and cannot be simply reconciled by adjusting input cluster parameters, which implies that they may originate from errors in the input physics and physical parameters. To overcome these difficulties, we took the observed offsets as empirical correction functions that one needs to apply to our specific choice of theoretical stellar models. When models are used with empirical corrections, we obtain distances from Sloan Digital Sky Survey (SDSS) photometry that are consistent with Gaia parallaxes, and our photometrically derived metallicities ([Fe/H]) are also in overall agreement with spectroscopic measurements in SDSS (Paper I).
While we developed a method of deriving photometric metallicities from SDSS, this set of color-T eff corrections is only valid for observations taken in the SDSS ugriz filter set. Other photometric surveys also adopt filter sets similar to that of SDSS, but their transmission curves are not exactly the same as each other, leading to nonnegligible differences in magnitudes. In this sense, direct calibration of synthetic spectra can serve as an alternative way of establishing such relations in various filter passbands. Importantly, it enables us to combine various photometric survey databases and produce chemokinematical phase-space maps over the entire celestial sphere in an internally consistent manner.
One of the goals in this study is to generalize our empirical correction procedure and construct a set of corrected synthetic spectra in order to generate magnitudes in any given filter set with high confidence. An essential requirement to achieve this goal is to finely sample flux for a set of calibration stars over a wide range of wavelength and stellar parameters using multiband photometry, which has become practical in the era of massive photometric surveys. The basic idea of calibrating model fluxes of theoretical stellar spectra, as opposed to making corrections on individual color indices in the models, was introduced by Lejeune et al. (1997Lejeune et al. ( , 1998)). However, the current work is based on a significantly larger set of photometric, spectroscopic, and astrometric data, which were unavailable then.
The other goal of this work is to probe a multidimensional data cube of Galactic stars, constructed based on the revised metallicity estimates. In addition to chemical information from photometry, we exploit kinematic data from Gaia, as in the previous papers of this series. The majority of main-sequence stars in our sample are too faint to have radial velocity measurements. However, along the great circle perpendicular to the direction of disk rotation (l = 0°and 180°), rotational velocities (v f ) in the cylindrical coordinate in the rest frame of the Galaxy do not depend on radial velocity. This enables us to derive v f by utilizing the proper motions of stars near the Galactic prime meridian and construct phase-space diagrams in v f and [Fe/H], which can subsequently be used to characterize kinematical and chemical properties of individual populations.
Because our method relies on calibration of theoretical models for main-sequence stars, our current approach is effectively limited to a local volume (d < 6 kpc). Nonetheless, in our previous papers we demonstrated the usefulness of such data by providing an unbiased, global perspective on local stellar populations, including Gaia Sausage/Enceladus (GSE; Belokurov et al. 2018;Helmi et al. 2018) and the Splash (Belokurov et al. 2020). In particular, we used Gaussian functions in Paper III to isolate individual stellar populations, and we evaluated their fractional contributions as a function of distance from the Galactic plane (Z) and Galactocentric distance (R GC ). In this work, we present evidence for yet another coherent structure of stars, which appears to have been formed during a starburst episode in the early history of the Galaxy, possibly driven by the GSE merger. This paper is organized into two parts, which describe each of the above two subjects: the revised calibration of theoretical models (Section 2) and its application to photometric databases (Section 3). We summarize our results in Section 4.
Stellar metallicities presented in this study are computed using a set of theoretical stellar isochrones with revised empirical corrections. As described below, our new color-T eff -[Fe/H] relations, which convert theoretically predicted quantities (T eff and luminosity) into observables (colors and magnitudes), hinge on both stellar sequences and field stars with spectroscopic metallicity estimates (see An et al. 2009An et al. , 2013; Papers I and II, and references therein, for more information on our previous model corrections). The same set of base theoretical models (Sills et al. 2000;Gustafsson et al. 2008) with identical model parameters, including an agemetallicity relation and α-element abundance mixtures, are adopted in this work. In contrast to the models adopted in the previous papers of this series, this newer version of calibration utilizes photometry in various filter passbands over a wide range of wavelength and thereby enables fine-tuning of synthetic stellar spectra.
As summarized in Table 1, we adopt both cluster sequences and a set of individual stars with spectroscopic metallicity estimates for the calibration of models. As shown in the first two columns, a total of 19 passbands from various photometric surveys are utilized in this study: ugriz photometry from the SDSS DR14, uvgriz from the SkyMapper Sky Survey (SMSS) DR2 (Onken et al. 2019), grizy from the Pan-STARRS1 surveys (PS1; Chambers et al. 2016), and BV from the AAVSO Photometric All-Sky Survey (APASS) DR10 (Henden et al. 2018). The gri photometry in APASS is not used owing to large photometric zero-point offsets (see also Tonry et al. 2018). We add to the list standard-star photometry in BVI C constructed by P. Stetson (see Stetson 2000). 6 Some of these databases use the same notation for their filter passbands (ugriz), but their response functions are not identical. Below we make a distinction between these filters by specifying the survey names.
Table 2 lists the stellar sequences adopted in this study. As in our previous exercise, we employ a set of well-studied Galactic globular and open clusters (M15, M92, M13, M3, M5, M67, and NGC 6791) over a wide range of metallicity (-2.4 [Fe/ H] 0.4). We use fiducial sequences from An et al. (2008), which were derived from the SDSS imaging data. Zero-point corrections (An et al. 2013) are applied to tie An et al. cluster In addition to clusters, we employ the Gaia double sequence, which appears on a color-magnitude diagram from stars with large proper motions (Gaia Collaboration et al. 2018). As demonstrated in Paper II, each of the sequences represents two dominant populations in the local halo-GSE and Splash-and has [Fe/H] ≈ -1.3 and -0.4, respectively (see also Sahlholdt et al. 2019). This dynamically defined group of stars provides a powerful constraint on the shape of a sequence, bridging the gap between globular and open clusters at intermediate metallicities. We follow the procedure developed in Paper II to extract individual sequences from the Gaia double sequence (see Appendix A). In short, this technique relies on metallicitysensitive u-band photometry to separate the two chemically distinct populations. For SDSS and SMSS, we use their u-band data. For PS1 and APASS, we use SDSS u. Additional cuts on kinematics further help to isolate each of the populations. As in our previous work, we compute v f using Gaia's proper motions and parallaxes, but without radial velocity measurements, within ±30°along the Galactic prime meridian. To construct a clean sequence, we use objects with good astrometry, having less than 20% uncertainty in parallax and 30% uncertainty in proper motion. We also apply cuts on E(B -V ) < 0.1 and |b| > 20°. As in Paper II, we impose -150 km s -1 < v f -50 km s -1 and 120 km s -1 < v f 150 km s -1 on the sample to derive the blue (metal-poor) and red (metal-rich) main sequences, respectively, for which we compute weighted median colors in bins of 0.5 mag in M r (ΔM r = 0.2 mag in ur CMDs).
For individual calibration stars with spectroscopic metallicities, we utilize the Sloan Extension for Galactic Understanding and Exploration (SEGUE; Yanny et al. 2009;Rockosi et al. 2022) and the Galactic Archaeology with HERMES (GALAH) survey (Buder et al. 2021), as they are among the largest and most uniform spectroscopic data sets in the northern and southern hemispheres, respectively. Specifically, we adopt metallicity estimates from a rerun of the updated SEGUE Stellar Parameter Pipeline (SSPP; Lee et al. 2008aLee et al. , 2008b)), performed by one of the coauthors (Y. S. Lee). We apply cuts based on a signal-to-noise ratio (S/N) of the spectra (>30), uncertainty in T eff < 200 K, and > g log 3.5 to select mainsequence stars with high-quality parameter estimates. For the GALAH sample, we also require that a stellar parameter quality flag (flag_sp) and an overall iron abundance quality flag (flag_fe_h) are not set.
The advantage of the SEGUE and GALAH samples is corroborated by the availability of photometric data in various passbands (such as SMSS uv; see Table 1). All SEGUE stars are covered by SDSS and PS1 imaging surveys but can only be matched to objects in SMSS near the equatorial region. Likewise, the majority of GALAH stars have good matches to SMSS photometry, but only a small fraction of its survey area overlaps with SDSS. Primary stellar sources (type = 6) in the SDSS are kept, with a set of minimal quality flags in the r-band measurements to ensure that sources do not have issues such as deblending, interpolation, and saturation. Similarly, primary detections in PS1 from its stacked imaging catalog are taken. We select point-like sources by imposing a maximum 0.05 mag difference in i-band photometry between a point-spread function and Kron magnitudes. For SMSS, we apply cuts on a number of photometric quality flags to only retain good photometric measurements: class_star > 0.9, flags < 3, nch_max = 1, prox > 7.5, ngood_min > 1, and nimaflags = 0 in each passband.
Photometric catalog objects are matched with Gaia Early Data Release 3 (EDR3; Gaia Collaboration et al. 2021) and Data Release 3 (DR3; Gaia Collaboration et al. 2023) using a 1″ search radius. Zero-point corrections on parallax by the Gaia team (Lindegren et al. 2021) are adopted. In the following calibration, objects with good parallaxes (σ π /π < 0.2) are used; about 68% of them have parallaxes within 2% from those inferred based on Bayes's theorem (Bailer-Jones et al. 2021). The foreground reddening values in Schlegel et al. (1998) are adopted, except for clusters, along with extinction coefficients at R V = 3.1 in Schlafly & Finkbeiner (2011) for SDSS, PS1, and Johnson-Cousins (those listed as "Landolt") bands. For SMSS, we adopt values in Wolf et al. (2018). These coefficients assume 14% reduction in the original E(B -V ) in Schlegel et al. (1998).
In this study, we aim to provide a set of isochrones with empirically calibrated synthetic spectra over a wide range of T eff and [Fe/H], which in turn can be employed to derive such quantities in other stars. Since we restrict our analysis to mainsequence stars, T eff can be directly mapped onto mass, luminosity, and surface gravity ( g log ) of a star in the isochrones; we take T eff as an independent variable in the following model comparisons. We compute theoretical flux ratios in various filter passbands as a function of wavelength, T eff , and [Fe/H]; then we attribute any deviation from the calibration samples to systematic errors in the models.
For this purpose, we employ YREC isochrones (Sills et al. 2000) and synthetic spectra generated using MARCS model atmospheres (Gustafsson et al. 2008); see An et al. (2009) for more information on the construction of the MARCS model library. We adopt the same age-metallicity and [Fe/H]-[α/Fe] relations for Galactic stars as in our previous papers of this series (see also An et al. 2013): ([Fe/H], [α/Fe]) = {(-3.0, +0.4), (-2.0, +0.3), (-1.0, +0.3), (-0.5, +0.2), (-0.3, 0.0), (+0.4, 0.0)} and ([Fe/H], age) = {(-3.0, 13 Gyr), (-1.2, 13 Gyr), (-0.3, 4 Gyr), (+0.4, 4 Gyr)}, with a linear interpolation in this metallicity grid. Inhomogeneous α-element abundance ratios in the Milky Way have a net effect of changing the overall metallicity of a star, but its impact is only mild; Δ[Fe/H] ∼ ± 0.2 for Δ[α/Fe] ∼ ± 0.2 (e.g., Kim et al. 2002). The effect of age is minimized in this study by restricting our sample to low-mass main-sequence stars.
We derive synthetic stellar colors in SDSS ugriz from the MARCS library using filter-response curves on the project webpage. 7 When deriving a Vega magnitude, we use a Vega model from the Hubble Space Telescope CALSPEC library (Bohlin et al. 2014), with the suggested flux rescaling in Riello et al. (2021). Effective wavelengths (λ eff ) of the filter passbands are computed for each model, as they have a mild dependence on the underlying stellar spectrum (primarily on T eff ). All flux ratios are referenced to the SDSS or SMSS r band among different filter passbands because of the large amount of flux collected in this passband and its relatively weak metallicity sensitivity in the T eff regime considered in this study (4000 K < T eff < 7000 K). Its bolometric corrections are also less prone to systematic errors.
For a given [Fe/H], differences in colors between observational data and models are computed as a function of absolute magnitude in the r band (M r ). We then use an M r -T eff relation of the isochrones to infer T eff from M r . For all of the calibration samples, we determine T eff photometrically using stellar isochrones, even if there exist spectroscopic T eff for the SEGUE and GALAH samples. This is because individual spectroscopic T eff measurements are not typically available for stars in clusters, and it is well known that there exists a T eff scale difference of a few hundred kelvins between photometric and spectroscopic approaches (e.g., Pinsonneault et al. 2004). This significantly reduces systematic differences in color-T eff relations from our heterogeneous data sets and makes our correction procedures internally more consistent. The infrared flux method (IRFM) is another useful way of deriving T eff from photometry, but our approach has an advantage of making a specific prediction on T eff and g log of a star, both of which are necessary for generating precise synthetic model colors.
Figure 1 displays a T eff -[Fe/H] space covered by the stellar sequences and spectroscopic targets in our sample, where the T eff values are determined from our isochrones. Stellar sequences in SDSS are shown by orange boxes with a width of ±0.1 dex in [Fe/H]. The dotted lines represent the same sequences, but without valid u-band measurements, demonstrating the necessity of deeper u-band photometry. Stellar sequences have discrete metallicities, while the spectroscopic sample (the green 2D histogram for SEGUE and contours for GALAH targets) fills up the remaining space. Our spectroscopic sample is heavily biased toward more metal-rich stars, with a significantly lower number in the metal-poor regime ([Fe/H] < -1). The upper right corner is not covered by both samples, due to the increased metal content and relatively old ages of the stars.
For each star in the spectroscopic sample, an isochrone is generated by interpolating the model grid at the star's metallicity ([Fe/H]), and differences in flux are computed in various filter passbands, as done for the stellar sequences. However, to evenly sample cluster sequences and the spectroscopic targets, we bin each of the SEGUE and GALAH samples in [Fe/H] and compute mean flux offsets as a function of T eff (see Appendix A). In accordance with the metallicity of the observed stellar sequences (Table 2), the central metallicities are set to [Fe/H] = {-2.4, -1.6, -1.3, -0.4, 0.0, +0.3}. To match this binning and keep the average metallicity of a subset of stars as close as possible to the central metallicity values, the spectroscopic sample is divided into [Fe/ H] = {(-2.9, -2.1), (-1.8, -1.4), (-1.5, -1.1), (-0.5, -0.3), (-0.1, +0.1), (+0.2, +0.4)}. A large width (0.8 dex) is set in the lowest metallicity bin to compensate for the small number of stars in the sample. Because the flux difference from our model changes mildly with metallicity, our adopted bin sizes have little impact on the following calibration. Nonetheless, mean flux offsets from GALAH are taken only at [Fe/ H] = {-0.4, 0.0, +0.3} owing to the lack of metal-poor stars in the sample.
Figures 2 and3 show magnitude differences between models and observations, as a function of λ eff , at two selected T eff (5800 and 4500 K, respectively). In each panel, differences are shown for the sequences and an ensemble of spectroscopic targets with open and filled symbols, respectively; different symbols are used to indicate references for photometry. All flux differences are registered to the SDSS r band, or the SMSS r band for the GALAH sample, due to a relatively small number of cross-matches with SDSS objects in the southern sky. The differences from the model at these two r passbands are nearly the same and are defined to be zero (i.e., their bolometric corrections are assumed to be correct). 9Error bars for the stellar sequences in Figures 2 and 3 represent propagated uncertainties from photometry and input parameters ([Fe/H], (m -M) 0 , E(B -V ), and age; Table 2). For the fiducial sequences in SDSS and PS1, constant uncertainties of 0.02 mag are assumed in the color indices. Similarly, a 2% error is adopted for the mean colors of Stetson's BVI C cluster sequences as a conservative limit. Photometry of Gaia's double sequence is collected across a large area on the sky, and therefore an observed scatter is taken as uncertainties in the mean colors, unless propagated photometric uncertainties are larger, since it represents a sum of random and systematic zero-point errors. A comparison of APASS photometry to Stetson's standard photometry for the sample clusters (except NGC 6791) reveals an rms dispersion of 0.04 mag in BV. Thus, it is added in quadrature to a total error budget for the APASS-based double sequence.
For the spectroscopic sample, the error bars in Figures 2 and 3 also indicate the quadratic sum of random and systematic uncertainties. The random component includes uncertainties in photometry, spectroscopic [Fe/H], and Gaia parallax. For the systematic uncertainty, 0.1 dex in [Fe/H] is assumed to take into account a scale difference between our models and the spectroscopic determinations. A 20% uncertainty in age is adopted for all stars. Flux differences from these systematic uncertainties are typically less than 0.01 mag in gizyBV but are as high as 0.05 mag in uv at high metallicities. In APASS, 0.04 mag uncertainty in photometry is further incorporated into the final uncertainty (see above).
In Figures 2 and3, there are systematic differences between the two classes of samples. The stellar sequences tend to show smaller flux deviations from the models than the spectroscopic sample; for example, at [Fe/H] -0.4, our result indicates that the spectroscopic sample is fainter at λ < 4000 Å but brighter at λ > 7000 Å than the sequences. Indeed, the cluster sequences in PS1 grizy (Bernard et al. 2014) exhibit the smallest differences overall. The observed discrepancy between stellar sequences and spectroscopic samples can be caused by inconsistent metallicity scales. However, other sources of errors, such as adopted ages, may also contribute to the observed offset, although the absolute model deviation changes monotonically with age, without modifying the observed wavelength-dependent offsets.
In Figures 2 and3, the red line indicates an average model deviation as a function of λ eff . Mean magnitude differences between models and observational data are computed by linearly interpolating values at three adjacent filters in λ eff . They are smoothed by applying a boxcar average with a width of 1000 Å, which is comparable to the FWHM of a broad filter passband. Average differences indicate that our models greatly overestimate flux below 5000 Å, by up to 20%, while they underestimate flux at longer wavelengths, by 10% at the most. The BV photometry from APASS (open and filled downwardpointing triangles) exhibits consistently larger fluxes than the other calibration sample by 0.05-0.1 mag, but the differences from the mean line are within our estimated 1σ uncertainties. We suspect all-sky photometric zero-point errors (at roughly 4% levels) as a likely source of the systematics (see above). Each data point represents one of the metallicity groups in this work. As shown by the blue lines, a third-order polynomial function is used to depict the observed trend. Beyond the metallicity range covered by the sample, a constant offset is assumed in the model deviation. Figure 5 displays slices of these mean model deviations at some selected wavelengths. Red colors indicate that models overpredict the flux, while the blue colors show regions with underpredicted flux. In this way, we construct a three-dimensional data cube of model deviations as a function of T eff , [Fe/H], and wavelength.
To first order, the model deviations change monotonically with wavelength: models are brighter than observations at shorter wavelengths, while the sign reverses at longer wavelengths. This may suggest an offset in the T eff scale of the isochrones as a major source of the systematic mismatch. However, the required amount of offset must be very large, by about 400 K, even for warm stars (Figure 6); it is even larger for cooler stars. Even if there are systematic differences between different approaches of determining T eff , such as the IRFM and spectroscopic determination from excitation/ionization balance, this is beyond the accepted range of errors in the models. Therefore, it seems that the observed offsets originate from a combination of various sources of systematic errors, such as incorrect input physics or underestimated line absorption in the models. Boundary conditions in stellarinterior models may also be incompatible with the atmosphere models in this study. On the observational side, an inconsistent metallicity scale, incorrect assumptions on elemental abundance ratios, or errors in the assumed age could be responsible for the systematic T eff offsets.
The model differences are highly nonlinear in the [Fe/H] versus T eff space, and no simple function can be adopted to remedy the problem. For this reason, we take the observed flux offsets in Figure 5 as a correction matrix for our choice of stellar isochrones and synthetic spectra. More specifically, we employ a semiempirical approach to correct synthetic spectra based on observations, while keeping stellar-interior models intact. Unlike in our previous work, the data cube in Figure 5 provides a continuous function of magnitude correction in wavelength. Thus, the corrected synthetic spectra can be applied to any filter sets in the wavelength range covered by our calibration sample.
In comparison to purely theoretical models, the net result of our empirical correction is redder colors or a higher photometric T eff , due to overestimation of the model flux at shorter wavelengths and underestimation at longer wavelengths. Apart from this fundamental change in the models, the revised calibration also differs from our earlier versions of the empirical corrections. The biggest change is the inclusion of individual spectroscopic targets in the sample, which inevitably modifies the metallicity scale of the models. As shown earlier in Paper I, isochrones calibrated using fiducial clusters produce photometric metallicities that are in agreement with spectroscopic estimates in SEGUE within Δ[Fe/H] ∼ 0.1 dex at [Fe/ H] > -1.5, but the difference amounts to ∼0.4 dex at [Fe/H] = -2, in the sense of a lower metallicity from our cluster-based approach. Such a difference is a direct consequence of a systematic offset in the metallicity scale between the cluster sequences and the SEGUE stars.
Figure 7 compares photometric metallicities from Paper III with those in this work, based on the revised calibration. The comparisons are shown for the solutions based on Gaia parallaxes using SDSS and PS1 photometry. Only stars having parallax uncertainties less than 10% are included. Other constraints are the same as for the main sample, as described in the next section, except in the top panel, where the metallicity difference is displayed over the full range of M r .
The large deviations for bright stars are evident, which are caused by the lack of hot stars in the calibration sample (Figure 1). In our subsequent analysis, including the bottom panel of Figure 7, we adopt 4.5 < M r < 7.5 to avoid regions with potentially large calibration errors.
The weighted median difference in Figure 7 indicates that the two calibration versions agree at high metallicity ([Fe/H] > -1) but that photometric metallicities from this work become larger for metal-poor stars, amounting to 0.35 dex at [Fe/H] = -2. Because the GALAH sample is confined to metal-rich stars in our calibration, the systematic trend highlights an inconsistent metallicity scale between the cluster sequences and SEGUE stars. The SSPP estimates have been checked thoroughly using clusters and high-resolution spectroscopic abundance determinations, and the overall agreement is impressive (Rockosi et al. 2022). Nonetheless, such comparisons were performed mostly using giants and main-sequence turnoff stars in the metal-poor regime, due to the lack of metalpoor main-sequence dwarfs in the SEGUE sample with a sufficiently high S/N. Therefore, the difference could originate from an internally inconsistent metallicity scale between dwarfs and giants in the SEGUE sample.
In summary, because neither metallicity scale is preferred over the other, the revised isochrones obtained in the current experiment should be taken as an alternative to our earlier cluster-based calibration. More precisely, a metallicity distribution of metal-poor stars in this study is hinged on an intermediate metallicity scale between SEGUE and the clusterbased work, set by their relative weights to the final calibration sample. The sense is that metallicity estimates in this study are systematically higher than those in our previous work at [Fe/ H] < -1.
Figure 8 shows a comparison between our distance estimates based on SDSS (combined with PS1) and Gaia parallaxes. In this case, we use distances that are determined jointly with metallicity, without relying on Gaia parallaxes. The same set of data as in the bottom panel of Figure 7 is used, but a more stringent test is performed by restricting a comparison to those having a reduced χ 2 fit of the models less than unity. We note that the majority of unresolved photometric binaries should have been rejected in Figure 8 (and the bottom panel of Figure 7), as we display photometric metallicity estimates at Gaia parallaxes on the abscissa. This is because unresolved binaries are systematically brighter than single stars and therefore show poor model fits to the observed fluxes, if a distance is fixed at a true value in the metallicity estimates (see also An et al. 2013). Our models produce a local distance scale that differs by at most 5% from Gaia EDR3 parallaxes. Error bars represent a standard deviation of the differences in bins of [Fe/H], but the errors in the "mean" differences are very small (<1%). In fact, the above good agreement with Gaia distances is not unexpected, since our calibration relies on Gaia parallaxes for nearby spectroscopic samples. We can expect a similar level of agreement to that of Bailer-Jones et al. (2021), as their Bayesian distance estimates closely align with those in the Gaia catalog when the parallax accuracy is less than 10%.
Figure 8 shows that the weighted standard deviation of the differences in distance is approximately σ[(m -M) 0 ] = 0.10-0.14 mag across the range of metallicities displayed. However, the quadrature sum of uncertainties from both methods is estimated to be in the range of 0.2-0.3 mag. This suggests that our distance measurement uncertainties may be overestimated by a factor of approximately three. The comparison with distance estimates based on SMSS (combined with PS1) also shows a similar level of discrepancy. One possible explanation for this discrepancy is that our model fitting does not account for correlations between photometric measurements in different passbands.
In addition, we compare our photometric metallicity estimates with spectroscopic measurements in GALAH. For stars with [Fe/H] > -1, the weighted standard deviation of the difference in metallicity is σ([Fe/H]) ≈ 0.15 dex, when photometric metallicities are estimated using SDSS or SMSS photometry without Gaia parallax priors. On the other hand, the expected value from propagation of uncertainty measurements is nearly 0.25 dex, indicating that our estimated uncertainties of photometric metallicities are overestimated by a factor of about two. However, when Gaia parallax priors are used in the computation of photometric metallicity, the difference between the estimated standard deviation and the propagated value is marginal; our measurement uncertainties are overestimated by only up to approximately 25%, depending on the stellar metallicity.
In this section, we apply our newly calibrated set of models to large photometric catalogs and provide new insights for Galactic stellar populations in the local volume. In addition to chemical information from photometry, we exploit kinematic data from Gaia to generate phase-space maps at various distances from the Galactic plane and Galactocentric distances. As in our previous papers, we restrict our analysis to a strip within ±30°from the Galactic prime meridian (l = 0°and 180°), where a conversion from transverse motions into v f is reliable (see Paper III for more details). Below, we first inspect metallicity distributions and phase-space diagrams to validate our new calibration (Section 3.1) and present distributions of scale heights and lengths for the stars in each section of v f and [Fe/H] (Section 3.2). Based on phase-space diagrams of high proper-motion stars, we demonstrate that our data reveal yet another stellar population formed during a period of Galactic starburst activity (Section 3.3).
In the following applications, we use Gaia EDR3 as a main source catalog, since v f in the rest frame of the Galaxy is a major ingredient of our phase-space diagrams. We combine Gaia EDR3 astrometric data with photometry in SDSS, SMSS, and PS1 using a 1″ match radius. As there is a little overlap between SDSS and SMSS, two sets of catalogs-SDSS ∩ Gaia and SMSS ∩ Gaia, respectively-are created, in addition to a master photometric catalog from all survey data ((SDSS ∪ SMSS) ∩ Gaia). PS1 grizy photometry is added to each data set, when SDSS or SMSS photometry exists, to better constrain the stellar parameters.
To estimate metallicities for individual stars in each catalog set, we employ the calibrated models and conduct a grid search. Gaia parallaxes degrade rapidly beyond ∼2.5 kpc from the Sun, whereasuncertainties in distance exhibit a more gradual increase when distances are derived photometrically (see Appendix C). Therefore, we calculate photometric metallicities under two conditions: with or without priors from the Gaia parallax. The former approach enables us to obtain photometric metallicities based on the best available parallax data but is restricted to nearby stars. In contrast, the purely photometric approach determines both distance and metallicity simultaneously and covers a larger volume of space. To assess the uncertainty in metallicity, we determine the Δχ 2 = 2.3 boundary for 2 degrees of freedom when Gaia parallaxes are used in the parameter estimation (Δχ 2 = 3.53 when using a purely photometric approach). In cases where we utilize Gaia priors, we determine the difference in the metallicity estimation from the ±1σ uncertainty in parallax and add it to the uncertainty in quadrature.
We have established certain criteria in our analysis to ensure the quality and reliability of our results. First, we only consider sources that have been detected in at least five photometric passbands in SDSS or SMSS, which guarantees that sources are observed in at least u or v. Additionally, we require that solutions have a reasonable fit to the model, as indicated by a reduced χ 2 of the best-fit model being less than 3, provided that they are within the range of 4.5 < M r < 7.5. To ensure accuracy in our photometric estimates, we exclude low-latitude regions with |b| < 20°and some areas where cumulative extinction exceeds E(B -V ) = 0.1. Moreover, we set a minimum safeguard by setting a maximum allowable uncertainty in metallicity of less than 1.5 dex. By implementing these selection criteria, we aim to minimize the potential for systematic errors and ensure that our results are of high quality and accuracy.
Figure 9 shows mean metallicity maps of stars in the local volume (1.2 kpc < d < 3 kpc) in the Galactic coordinate system (Mollweide projection), based on our metallicity estimates for individual stars in each of the three combined catalogs. Panel (a) displays a map from SDSS ∩ Gaia (3.2 million stars), which mostly covers the northern Galactic hemisphere, while panel (b) shows the southern hemisphere from SMSS ∩ Gaia (0.6 million stars). Panel (c) displays a full coverage map from SDSS, SMSS, and PS1 (3.7 million stars). Here we compute the mean metallicity in each HEALPix cell by convolving the metallicity of each star using a normalized Gaussian function, with a standard deviation set to its uncertainty. All three maps are smoothed using a median filter with a 2°radius.
There is a limited overlap between SDSS and SMSS, mainly along the celestial equator (≈80,000 stars). In these overlapping areas, photometric [Fe/H] estimates derived from individual catalogs (SDSS ∩ Gaia or SMSS ∩ Gaia) agree with those based on a combined catalog (SDSS ∩ SMSS ∩ Gaia) within 0.1 dex. Metallicity differences for individual stars also do not exhibit a systematic trend with metallicity, which provides a confirmation of the internal consistency in our models. Nonetheless, metallicity distribution functions from these subsets are not identical, due to unequal depths and qualities of these photometric surveys, which result in mild systematic differences in the mean metallicities in Figure 9.
Reassuringly, Figure 9 reveals that more metal-rich stars are found near the Galactic plane, from which a mean metallicity gradient is evident from the low to high Galactic latitudes, as expected from a simple population gradient. This exercise proves not only that our technique can be used to determine metallicities of stars precisely from multiband photometry but also that we can use corrected synthetic spectra to combine data in various filter sets to generate an internally consistent all-sky metallicity map. Our calibration procedure is currently valid for main-sequence stars and has lower precision for giant stars. Consequently, the above mapping based on main-sequence stars probes a local volume out to ∼6 kpc from the Sun. Giant stars are excluded in our sample using color-magnitude relations based on Gaia parallaxes, although a purely photometric approach can also be employed to tag such stars, as demonstrated in Paper III.
Our photometric technique is a sensitive probe of photometric zero-point errors. It is particularly useful for large photometric surveys because it is nontrivial to have an internally consistent photometric zero-point across large areas on the sky. In Appendix B, we demonstrate the existence of spatially correlated photometric zero-point errors in SDSS u and SMSS uv based on our corrected models. The size of photometric zero-point offsets is a few hundredths of magnitude level but as high as 0.1 mag in some areas. By inverting the problem, zero-point offsets in photometry can be derived to make a uniform mean metallicity of nearby stars on the sky. This backward design on photometric zero-point corrections improves the quality of the metallicity mapping and somewhat narrows the gap in model deviations between the cluster sequences and the SEGUE sample. For this reason, we iterate the calibration procedure (Section 2) using zero-pointcorrected photometry in SDSS u and SMSS uv. In the following analysis, including Figure 9, all input photometry is corrected for the spatially correlated zero-point offsets. These steps closely parallel similar exercises in zero-point corrections for the SMSS uv in Huang et al. (2021Huang et al. ( , 2022)).
More quantitative comparisons with previous studies can be made using phase-space diagrams such as shown in Figure 10 In panel(b) of Figure 10, a group of metal-poor stars is seen at [Fe/H] ≈ -2.2, which exhibits a slow net prograde rotation (〈v f 〉 ≈ 70 km s -1 ). Together with a more metal-rich ([Fe/H] ≈ -1.5) counterpart, which is most clearly seen in the inner Galactic region (panel (a)), it was considered one of the primary constituents of the Galactic halo in our series of papers. In Paper I, we referred to these components as "inner and outer halos" within the dual halo paradigm (Carollo et al. 2007(Carollo et al. , 2010;;Beers et al. 2012), whereas we labeled them as "metal-poor and metal-rich halos" in Paper II. On the other hand, we made a presumption in Paper III that both of them constitute the main body of GSE; this interpretation was primarily driven by the absence of other analogous structures known in the same phase space.
However, recent evidence suggests that the slowly rotating metal-poor stars are likely a separate entity from GSE. Belokurov & Kravtsov (2022) used aluminum abundances to separate in situ stars in the halo from accreted stars and found a Although they span an extreme range in v f (from ∼ -150 to 300 km s -1 ) and their spectroscopic sample is limited to [Fe/H] > -1.5, their approximate mean v f ∼ 100 km s -1 and low metallicity ([Fe/H] < -1) suggest that these stars (dubbed "Aurora") are a part of the prograde, metal-poor components seen in our previous work. Given its lower metallicity than the Splash, it is likely an old (primordial) in situ halo that formed before the GSE merger at z = 1-2.
In support of this view, we employ a simple proper-motion cut in the sample to separate stars in GSE from the metal-poor halo distribution, as demonstrated in panels (c) and (d) of Figure 10. They show the same phase-space diagrams of stars as in panel (b), while having different ranges of proper motion, <10 mas yr -1 and >10 mas yr -1 , respectively. In panel (c), both the metal-poor in situ halo and thick disk are seen, connected by a narrow band of stars, which we assigned to the metal-weak thick disk (MWTD)in Paper III. On the other hand, an elongated structure along the v f = 0 km s -1 line stands out from the high proper-motion sample in panel (d), which encompasses a wide range of metallicity (-3 < [Fe/ H] < -0.5). The chemical and kinematical properties of these high proper-motion stars are analogous to GSE in the original works of Belokurov et al. (2018) and Helmi et al. (2018).
Although our v f measurement is strongly correlated with proper motions, the above separation can be understood by the highly radial orbits of GSE stars, which contrast with a nearly isotropic velocity distribution of metal-poor in situ halo stars. We also note that the separation does not change appreciably even if the heliocentric distances of the sample are further narrowed down, indicating that smaller (negative) v f is not merely caused by systematically shorter distances. In summary, all of these chemical and kinematical properties of individual populations in our phase-space diagrams are consistent with those found from previous (mostly spectroscopic) studies, which essentially validates our photometric [Fe/H] and v f estimates.
To examine structural properties of each group of stars, we compute a scale length (L) and height (H) in bins of v f and [Fe/ H], with bin sizes set to 25 km s -1 and 0.2 dex, respectively. The survey data volume has a cone-shaped geometry, so we only use stars within ±1 kpc from the solar radius at 1 < |Z| < 4 kpc in each Galactic hemisphere to compute the scale height. Similarly, the scale length is determined using stars at 7 kpc < R < 12 kpc and 1 kpc < |Z| < 3 kpc. In each v f -[Fe/H] bin, we perform a linear least-squares fitting of the logarithmic number density of stars as functions of Z or R, where the bin size is fixed at 500 pc. Although the density profile of bulk halo stars has been known to follow a power law, we utilize an exponential function as a proxy in the local volume near the Galactic plane to obtain a relative comparison between different populations. From the best-fitting model, we calculate the standard deviation in a logarithmic number density and assume that all data points have the same uncertainty as this value. We then estimate the uncertainties in L or H by using the best-fitting slope and its estimated uncertainty in a successive regression (see Appendix C).
Figure 11 displays a distribution of scale height in each Galactic hemisphere from the SDSS ∩ Gaia sample (distances taken from Gaia). Only those pixels with a fractional uncertainty less than 25% are included. A 0.2 dex × 25 km s -1 pixel is subdivided into 4 × 4 subpixels, which are smoothed using a five-point boxcar average, in order to have a smoothed, global look at changes in the structural properties. Likewise, Figure 12 shows a distribution of scale length in both hemispheres, binned and smoothed in the same manner as in Figure 11. A requirement that stars lie in the stripe along the Galactic prime meridian severely limits the number of stars available in each data set, resulting in ≈4 × 10 6 and 2 × 10 6 stars in the northern and southern Galactic hemispheres, respectively.
In Figures 11(a) and 12(a), three notable features are seen in the northern Galactic hemisphere, which are marked by the circled numbers. The intermediate scale height and scale length valley is particularly evident, extending from a region mainly occupied by GSE stars ("①") to a region populated by the Splash ("②"). This valley appears even more dramatic, as it contrasts with a scale height "highland" at [Fe/H] ∼ -2.2 and v f ∼ +80 km s -1 ("③"). For reference, the circled numbers are also marked in the following figures, including Figures 11(b) and 12(b).
The intermediate scale height valley directly shows that GSE stars are distributed farther from the Galactic plane than disk stars, but not as far as metal-poor in situ halo stars. This observation can be understood by a low-inclination, highly radial orbit of the GSE progenitor, which penetrated deep into the primordial Galaxy. In addition, a comparison between the northern and southern hemispheres in Figures 11 and12 clearly demonstrates a larger amount of debris in the northern hemisphere, due to a pileup of stars at the last apocenter of a highly eccentric orbit of the GSE progenitor (Naidu et al. 2021).
The asymmetric stellar distribution in the halo between the northern and southern Galactic hemispheres can be checked using SMSS ∩ Gaia (based on Gaia parallaxes). SDSS mainly covers the northern Galactic hemisphere but has a limited coverage in the south. On the other hand, SMSS covers almost the opposite side of the celestial hemisphere (see Figure 9), with a total of ≈6 × 10 5 stars in the northern Galactic hemisphere and ≈10 6 stars in the south in our sample. Figure 13 displays scale height and scale length distributions from SMSS ∩ Gaia in the southern Galactic hemisphere, generated following the same steps used for Figures 11 and12. The coherent structure in the parameter space covered by GSE is not seen as clearly as in the northern Galactic hemisphere from SDSS ∩ Gaia, but it resembles the distributions in the south, highlighting the highly eccentric orbit of the GSE merger. Nonetheless, we note that, due to the limited overlap along the Galactic prime meridian, the SDSS ∩ Gaia sample is biased toward the Galactic anticenter direction in the southern hemisphere, while the sample from SMSS ∩ Gaia is more populated toward the Galactic center. Therefore, a trace of the intermediate scale height valley in panel (a) could be real and may capture the GSE debris in the direction toward the Galactic center. Deeper SMSS photometry in future data releases would be useful to explore this volume in more detail.
In Figure 11(a), a vertical trough with intermediate scale heights lies at [Fe/H] ∼ -0.8, extending from v f ∼ 0 to 100 km s -1 . Originally, Belokurov et al. (2020) defined a region occupied by the Splash as -0.7 < [Fe/H] < -0.2 and -150 km s -1 < v f < +100 km s -1 , which overlaps with this trough. The fact that the intermediate scale height valley from GSE stretches out to a region populated by Splash stars supports a previous claim that the GSE merger has dynamically heated stars in the primordial disk of the Milky Way (Bonaca et al. 2017;Belokurov et al. 2020). In other words, because a nearly in-plane collision with a dwarf galaxy would leave behind heated stars confined to a disk plane, this apparent coincidence supports a view on a causal connection between the Splash and GSE. In addition, since such stars originated from the primordial disk of the Galaxy, which was smaller in size than the current stellar disk, this also naturally explains the short scale length of these stars. Reassuringly, the lowland in panel (a) of Figure12 stretches out to the low v f region at the metallicity covered by Splash stars (-0.7 < [Fe/H] < -0.2), although such a feature is not clearly seen in panel (b).
However, our data clearly indicate that dynamical heating took place over a wider range of metallicity than previously considered by Belokurov et al. (2020). If stars were born in the inner region and then were displaced by mergers, the lowland with small scale lengths in Figure 12 shows an approximate extent where the dynamical heating took place. If the low scale length region is fenced by a contour line at L = 1-1.2 kpc, panel (a) indicates that even lower-metallicity stars down to [Fe/H]∼ -2 originated from similar excitation mechanisms. This conclusion is further supported by the near coincidence between the horizontal trough along v f ∼ + 150 km s -1 at -2 [Fe/H] -1.2 and the MWTD, which is considered to be the relic of the primordial disk. Therefore, even if the GSE progenitor was a major source of the dynamical heating of a primordial disk, it seems unlikely that the orbital properties of stars in this region are altered by a single massive merger event. Instead, our result indicates that such a dynamical heating process on the primordial disk was in operation, even before the GSE merger, driven by more numerous minor mergers, as predicted by numerical simulations of the early universe (e.g., Grand et al. 2020).
The other interesting feature is the "highland" in the scale height distribution ("③" in Figure 11), which coincides with the metal-poor in situ halo (see Figure 10). Its scale height reaches H > 1.1 kpc, which is significantly higher than those of heated stars (∼0.8 kpc); therefore, it is unlikely that it formed through the same dynamical heating process. Moreover, the highland exhibits a gradient in the scale length, as shown in Figure 12(a), in the sense that more metal-rich stars are more strongly concentrated in the inner region of the Galaxy. This implies that more active star formation took place in the deeper potential well, while the Milky Way has grown by chaotic coalescence of numerous small gas-rich dwarf galaxies in the early universe. It is unclear, however, how the Milky Way attained the net angular momentum in the same direction as that of the Galactic disk at this stage. Nonetheless, its small prograde net rotation (〈v f 〉 ∼ +80 km s -1 ) suggests that it had maintained a puffy disk-like structure owing to turbulent nature of the gas-rich minor mergers.
Our phase-space diagrams also reveal a long and narrow sequence of stars when a simple cut on proper motion is made. All four panels in Figure 14 are drawn from the same SDSS ∩ Gaia sample, but with different cuts on a minimum value of proper motions (20, 40, 60, and 70 mas yr -1 , respectively). Moreover, we convolve each count of stars using a normalized Gaussian function, with standard deviations set based on the measurement uncertainties in both axes. See Appendix D for the version generated using a Monte Carlo simulation.
Figure 14(a) is dominated by disk populations (mostly thickdisk stars), but a low-v f , low-[Fe/H] tail begins to show up when a mild cut on proper motion (>25 mas yr -1 ) is imposed. At >40 mas yr -1 , a striking elbow-like feature emerges from these diagrams, where halo and disk stars form a narrow, continuous sequence over a wide range of [Fe/H] and v f . This feature is characterized by two joint, orthogonal branches; the horizontal arm is nearly parallel to the v f = 0 km s -1 line over a wide range of metallicity (-2 [Fe/H] -0.6), while the vertical arm has a narrow metallicity range (-0.6 [Fe/ H] -0.4). The sequence passes through GSE and the Splash and is eventually connected to the disk. Although our v f estimates depend on distance, the elbow-like feature becomes stronger with higher proper-motion cuts but changes little with distance from the Sun.
As shown in Figure 15, the sequence persists even when SMSS ∩ Gaia is used, indicating that it is present in both hemispheres. Aside from our photometric metallicities, the same sequence can also be seen in Figure 16, based on spectroscopic metallicities in Gaia DR3 using the General Stellar Parameteriser-Spectroscopy (GSP-Spec) module (Recio-Blanco et al. 2022) from the Radial Velocity Spectrometer (RVS) spectra (λ/Δλ = 11,500). Because these spectroscopic observations are available for bright stars (G < 14 mag), the sample is limited to relatively nearby stars (1 kpc < |Z| < 2 kpc); the lower |Z| cut is made to exclude numerous disk stars. The proper-motion cut is lowered to 25 mas yr -1 to retain as many stars as possible along the sequence. On the other hand, the constraint on Galactic latitudes is lifted, as the spectroscopic metallicities are only weakly dependent on foreground reddening. While Figure 16(a) is based on the projected v f as for our photometric samples along the prime meridian, panel (b) is based on full three-dimensional space motions from both radial velocities and proper motions in Gaia. The sequence appears almost identical in both panels, which not only supports the existence of this coherent structure but also validates our approach for obtaining v f from proper motions only.
Our phase-space diagram presents a continuous chain of stars, suggesting that they formed through successive metal enrichment along this pathway. As seen in panels (b)-(d) of Figure 14, the sequence runs along ([Fe/H], v f ) ≈ {(-2, 0), (-0.6, 0), (-0.4, 180 km s -1 )}. Although the complete sequence covers nearly 2 dex in [Fe/H], the vertical arm of the sequence has a remarkably narrow range of metallicities. This implies that these stars formed quickly, with insufficient time for successive metal enrichment (except α-element enhancement from core-collapse supernovae), while starforming clouds were collapsing or reorienting their angular momentum vectors rapidly from nearly zero to ∼180 km s -1 in v f . In other words, our result indicates that the young Milky Way went through a phase of starburst activity. For this reason, we refer to this structure as the Galactic starburst sequence (GSS). The GSS traverses through known stellar populations and structures, including the GSE at the lower metallicity range and disk stars at the metal-rich end, implying a chronological order of formation of these Galactic components that can be traced to a common origin.
To better understand the properties of the GSS, we put the above results together in Figure 17. The red solid contour shows a region with a scale height H = 1 kpc from Figure 11(a), while the blue dashed contour delineates a scale length L = 1 kpc from Figure 12(a). They are overlaid on top of contours of stars with high proper motions (>60 mas yr -1 ) in Figure 14. First, the width of the vertical arm of the GSS is narrower than the metallicity range of the heated population, the extent of which can be delineated by the blue dotted enclosure with a small scale length. This indicates that the Splash, which can be defined as a group of metal-rich stars with halo-like kinematics, is in fact composed of two distinct groups of stars. The first is a group of dynamically heated stars, and the other is the starburst population. Heated stars formed in the inner region of the primordial Galaxy and then were displaced to the current location in the halo, while starburst stars formed in a top-down fashion, as indicated by a weak positive correlation between [Fe/H] and v f in our phase-space diagrams.
Second, both GSE and the GSS lie along the intermediate scale height valley, as shown by the red solid line in Figure 17. The simplest explanation is that the GSE merger has triggered a starburst in the mixture of gas from the primordial Milky Way and gas donated by the GSE merger, since mergers were likely gas-rich at high redshifts. According to this scenario, the radially biased orbit of GSE should be responsible for the relatively high transverse motions of stars in the GSS, which contrast with a nearly isotropic velocity distribution of metalpoor in situ halo stars, because these stars formed out of metalenriched gas having similar orbital properties with the GSE progenitor. Therefore, they are expected to show up in the high proper-motion sample, even though their orbits have evolved significantly over time from halo-like to disk-like orbits.
Notably, the full GSS is not seen in panel (d) of Figure 10, although it is also made using high proper-motion stars. This difference is manifested by a high central concentration of Splash(-like) stars in the Galaxy, as shown in Paper III, and more clearly demonstrated by the short scale lengths in Figure 12. The GSS is not a dominant structure at large Galactocentric distances.
Thanks to Gaia and large-area photometric and spectroscopic databases, it is now possible to perform accurate comparisons of theoretical models to extensive, high-quality data. In this work, we use models generated using YREC and MARCS and quantify model deviations as a function of T eff , metallicity, and wavelength. The current approach relies on a comparison with multiwavelength ultraviolet, optical, and near-infrared color-T eff relations derived from Galactic cluster sequences, Gaia's double sequence, and a sample of spectroscopic data in SEGUE and GALAH, using photometric data in broadband filters from SDSS, SMSS, and APASS and standard-star photometry in the literature. Mean flux deviations are derived as a function of wavelength, which amount to up to ∼20% at <4000 Å, but significantly less in longer wavelengths. We find that no single factor can remove the observed offsets, but it is more likely a problem arising from a combination of various sources of errors in the models and/or observational data. Subsequently, we define the model offset as an empirical correction function for our specific choice of models.
By combining our technique with proper-motion measurements in Gaia, we construct phase-space (v f vs. [Fe/H]) diagrams of stars to provide a global perspective on the stellar populations in the Milky Way. In this way, we identify a long and narrow sequence of stars in a phase-space diagram, which we call the GSS. The GSS is not a representation of a single stellar population, but rather consists of several previously known Galactic stellar populations or components, arranged like pearls on a string. In particular, it overlaps with GSE in a valley with intermediate scale heights, suggesting that GSE has likely triggered successive formation of these stars. It also passes through the Splash, showing rapid evolution of v f of star-forming clouds within a narrow metallicity range, which testifies to a starburst event in the young Milky Way. The wide metallicity range of dynamically heated stars as traced by small scale length regions in a phase-space diagram indicates that the Splash is likely composed of two stellar populations with distinct origins-dynamically heated and starburst populations.
The red sequence of stars with high transverse motions (>200 km s -1 ) in Gaia (Gaia Collaboration et al. 2018) is possibly another manifestation of the GSS. As shown in Paper II, the blue and red sequences are separated by a few tenths of a magnitude in griz colors, and more strongly when the u band is included, indicating a metallicity offset by ∼1 dex. Based on their corresponding metallicity ranges, stars in the blue sequence belong to GSE, while stars in the red sequence constitute the Splash (e.g., Gallart et al. 2019). According to our work, a fair fraction of stars in the red sequence should constitute the vertical arm of the GSS.
Our new perspective into the phase-space diagram of stars in the local halo reveals two consecutive modes of chemokinematical evolution: a rapid chemical enrichment along the orbit of GSE, followed by starburst in rapidly evolving orbits of gas clouds. Our finding is in line with recent numerical simulations in the literature, which explicitly predict the existence of a merger-driven starburst in the Milky Way-size galaxies (Cooper et al. 2015;Bignone et al. 2019;Grand et al. 2020;Renaud et al. 2021). For instance, in Grand et al. (2020), a merger with a gas-rich GSE-like progenitor triggers starburst activity owing to the increased compression of gas clouds in both galaxies. The star formation rate suddenly increases by a factor of two, which lasts less than 1 Gyr. Strikingly, their simulations reproduce some of the observed key features in our data, in that both starburst and heated populations produced by the merger event are present in a narrow metallicity bin but span a wide range of v f . Furthermore, the heated population shows a correlation between [Fe/H] and v f in their simulation, while there is essentially no such dependence for the starburst population. Considering an excessively large scatter in v f at a given metallicity, this may well be explained by a narrow range in [Fe/H] of the vertical arm of the GSS. Another intriguing aspect is that the starburst population in the simulation exhibits a radially concentrated, rotationally supported disk. This prediction is also in agreement with the fact that the GSS lies along the short scale length, intermediate scale height valley.
If the vertical and horizontal arms of the GSS are attributed to stars formed during a merger-driven star formation and those accreted from the GSE progenitor, respectively, we also find a good agreement in the fraction of such stars with numerical simulations. We count each group of stars in our phase-space diagram in Figures 14 and 15 using a simple box criterion: -50 km s -1 < v f 200 km s -1 and -0.9 < [Fe/H] 0 for the starburst, and -60 km s -1 < v f + 60 km s -1 and -2 < [Fe/H] -0.9 for the accreted population, respectively. The observed ratios between these two populations are 1.1-3.9 from SDSS ∩ Gaia (for proper-motion cuts at 60 and 20 mas yr -1 , respectively) and 1.7-3.1 from SMSS ∩ Gaia (for 60 and 40 mas yr -1 cuts, respectively), in that the total mass of stars formed during a merger-driven star formation in the primordial Milky Way exceeds the stellar mass of accreted stars. They are quantitatively in agreement with a stellar mass ratio of ∼2-5 found in Milky Way-like simulations with radially anisotropic stellar halos in Grand et al. (2020, except the case of a merger with the lowest stellar mass); see also Orkney et al. (2022).
On observational grounds, Myeong et al. (2022) argued for a starburst event triggered by the GSE merger by analyzing spectroscopic databases from APOGEE and GALAH (including the α-elements, Al, and Ce) and dynamical information from Gaia (orbital energy). Based on unsupervised Gaussian mixture models, they showed that the local halo populations could be described by four Gaussian components, of which three were previously known-the GSE, the Splash, and the in situ halo ("Aurora")-while the former was further divided into metal-poor and metal-rich parts. The remaining component was found to reside between the GSE and the low-α (thin) disk in the chemical space. The authors argued that the stars belonging to this component named "Eos" were formed from a starburst. The metallicity range of ∼170 possible members is bounded by a 2σ uncertainty to approximately -1.0 < [Fe/ H] < -0.3, with a mean of -0.7, where we found a rapid change of dynamical properties of GSS stars. These results provide supporting evidence that Eos could represent a subcomponent of the GSS, assuming that it formed during the same starburst event. However, there is currently a lack of available dynamical information on Eos members, and further investigations are needed to confirm the relationship between Eos and the GSS.
Using the framework of the GSS, additional pieces of the puzzle from other recent observational studies can be put together in a coherent manner. For instance, Lee et al. (2023) delved into a chemical space originally occupied by the Splash in the spectroscopic database from the SEGUE and the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST; Cui et al. 2012). They found that the sample in this narrow metallicity bin can be split into the low-and high-[α/Fe] groups, given the large dispersion in radial velocity (v R ) of the former in the Galactocentric coordinate system. They used systematic changes in kinematics (orbital inclinations and eccentricities) and v R -[Fe/H] relations to argue that about half of the low-[α/Fe] population and the majority of the high-[α/ Fe] population are from the GSE progenitor and the dynamically heated disk population, respectively, while the rest are of different origin, most likely from starburst activity.
In addition, Ciucă et al. (2023) utilized precise stellar-age estimates based on the asteroseismic measurements in APOKASC-2 (Pinsonneault et al. 2018) and discovered a rapid chemical enrichment of stars from [Fe/H] ≈ -0.5 to ≈+0.2 at a look-back time of 10-12 Gyr, which they dubbed the "Blob." It is also accompanied by an increase of [Mg/Fe] and a small decrease in [Fe/H] ("Dip") prior to this period, which qualitatively agrees with mixing of fresh gas by a gasrich merger in numerical simulations. This "Blob" feature is most clearly seen in the inner Galactic region, which is in line with our finding that the GSS has a small scale length.
In other studies, age-metallicity relations also reveal a chain of metal-rich stars, establishing a link between the halo and the disk (Haywood et al. 2013;Nissen et al. 2020;Xiang & Rix 2022). Most recently, Xiang & Rix (2022) used α-elementrich stars with low orbital angular momenta and demonstrated the existence of a narrow and continuous age-metallicity relation in the Milky Way. According to their analysis, the Milky Way achieved a high-metallicity floor ([Fe/H] ∼ -1) about 13 Gyr ago, with successive metal enrichment over the following ∼5 Gyr. The majority of stars are found in a narrow metallicity range (-0.7 < [Fe/H] < -0.2), which leads to its possible connection to the vertical arm of the GSS. Importantly, its narrow and continuous channel of stars in Xiang & Rix (2022) and Ciucă et al. (2023) reinforces the physical nature of the GSS.
At large v f , our data show that the GSS is connected to the disk, providing direct evidence that disk stars were formed in part from gas clouds left behind after the starburst episode.
Interestingly, according to the dichotomy of disk stars into high-and low-[α/Fe] sequences (e.g., Hayden et al. 2015, and references therein), the tip of the GSS has a metallicity ([Fe/H] ≈ -0.4) similar to that of the low-metallicity end of the low-α disk. This implies that gas clouds in the protogalactic disk were diluted by fresh, low-metallicity materials accreted by the GSE merger, as has often been invoked in numerical simulations (Brook et al. 2007;Buck 2020); see also Chiappini et al. (1997). Inflow from the circumgalactic medium (Grand et al. 2020) or a gaseous outer disk (Renaud et al. 2021) after the merger may also be plausible (see also discussions in Myeong et al. 2022). Intermediate [α/Fe] is then a consequence of corecollapse supernovae during the starburst. The inflow of metalenriched gas onto the thin disk can naturally explain the G-dwarf problem as well-the apparent excess of metal-rich stars in the local disk compared to the prediction from a simple closed-box chemical model (e.g., Greener et al. 2021).
In this context, we conjecture that the "Nyx" stream identified by Necib et al. (2020) is closely related to the vertical arm of the GSS. In their study, a group of ∼100 stars in the solar vicinity shows coherent radial (134 km s -1 ) and azimuthal (130 km s -1 ) motions. Their metallicities peak at [Fe/H] = -0.55, with a dispersion of 0.13 dex, and are mostly confined to a plane with a maximum vertical distance of Z = 1.7 kpc. The latter value is consistent with a small scale height of stars having similar [Fe/H] and v f to Nyx stars (H ≈ 0.6 kpc; Figure 11). Nonetheless, most of the Nyx stars have 100 km s -1 v f < 250 km s -1 , with a few extreme cases, and therefore trace only the upper half of the vertical arm of the GSS. Necib et al. (2020) postulated that the Nyx stream is a remnant of a disrupted dwarf galaxy. However, a more recent analysis in Zucker et al. (2021) argues against their extragalactic origin, based on the fact that a subset of these stars are indistinguishable from thick-disk stars in the elemental abundance space, having larger α-element abundances than those of accreted stars at a given metallicity.
The extensibility of calibrated synthetic spectra presented in this work will enable accurate prediction of stellar magnitudes for filter passbands in various photometric surveys. Our metallicity mapping technique based on the empirically calibrated isochrones will also serve as a useful resource for studying the demographics of stellar populations that are yet to be discovered in the local universe from the upcoming surveys such as the Legacy Survey of Space and Time (LSST; Ivezić et al. 2019). Appendix B Photometric Properties of SDSS, SMSS, and PS1
Because photometric data in various filter passbands are combined to derive a set of stellar parameters in this study, it is necessary to adopt accurate photometric uncertainties in the χ 2 statistics. Here we compare SMSS and PS1 photometry with SDSS and compute a mean magnitude difference and a dispersion to infer the size of the true uncertainties. The comparison with SMSS is limited to narrow regions, since the imaging stripes of SDSS overlap only a little with SMSS footprints along the celestial equator. The overlap with PS1 is more extensive, as both SDSS and PS1 cover the northern hemisphere.
Figure 24 shows statistical properties of the comparison of SDSS photometry with SMSS (left) and PS1 (right), respectively. The top panels display a distribution of a magnitude difference of relatively bright (14.5 mag < r < 18 mag) stars in each passband. To take into account nonnegligible color terms between different filter passbands (i.e., color transformations), the magnitude differences in each 1°-wide strip in R.A. are fit using a third-order polynomial as a function of gr in 0.3 < gr < 1.1, and its mean trend is removed. The SMSS v is compared to SDSS u, and PS1 y to SDSS z. A weighted median difference is computed in each (Δα, Δδ) = (1°, 1°) region, and an ensemble of these differences are fit using a Gaussian function, as shown by solid lines in the top panels of Figure 24. Standard deviations of the best-fitting Gaussian functions provide a measure of the spatial variation of the photometric zero-points across the sky. They are mmag, respectively, which are comparable to the quoted systematic uncertainties in these surveys.
The bottom panels of Figure 24 show ratios between the standard deviation of a magnitude difference and a propagated uncertainty for each filter passband in each (Δα, Δδ) = (1°, 1°) patch. Again, SMSS v is compared to SDSS u, and PS1 y to SDSS z. Some of the brightest objects are rejected in this comparison, due to unrealistically small uncertainties in the SDSS PSF i magnitude (<0.005 mag). As shown in the bottom left panel, the ratios for the SMSS and SDSS grz are near unity, indicating that photometric uncertainties are comparable to the observed scatter, while photometric uncertainties in uv are likely underestimated. As shown in the bottom right panel, the differences between propagated uncertainties and observed dispersions are even larger for PS1 passbands.
Based on the above comparisons, an uncertainty "floor" (σ f ) is computed in each passband in order to make a median of a standard deviation equal to a median of propagated mmag for SMSS and PS1, respectively. In all cases, zero-point uncertainties (σ zp ) are overwhelmed by the uncertainty floors (σ f ). We add both uncertainties in quadrature to the original photometric uncertainties in SMSS and PS1 and use them throughout this work.
By design, the zero-point of SMSS DR2 photometry was set based on synthetic griz photometry in the PS1 system, calculated from all-sky Gaia photometry (Onken et al. 2019). The wavelengths of SMSS griz passbands overlap with Gaia and PS1 passbands, so in principle one can tie them together without losing information on the properties of stellar spectra. On the contrary, photometric zero-points in u band and v band remain largely unconstrained owing to the lack of shortwavelength passbands in Gaia and PS1.
Figure 25(a) shows the mean metallicity distribution of stars (1.2 kpc < d < 3 kpc) from the SMSS in the Galactic coordinate system, which demonstrates the necessity for a second-order photometric zero-point correction. Here we use metallicities from a fully photometric solution (0.6 million stars), as it is more prone to photometric errors than the case based on Gaia parallaxes. The mean metallicity in each pixel is derived from a generalized histogram of photometric metallicities, which accounts for an uncertainty in metallicity by taking it as a standard deviation of a normal probability distribution. We implement HEALPix (Górski et al. 2005) in the Galactic coordinate system, for which we set a resolution parameter N side to 32, corresponding to a constant pixel size of 3.3 deg 2 . Although the bright survey limit restricts the sample to relatively nearby stars, nearly uniform metallicities of stars are contrary to what is expected in the local volume.
To obtain zero-point corrections on SMSS uv photometry, we assume that nearby stars (0.5 kpc < d < 1 kpc) have the same metallicity in every direction, as they are mostly thin-disk stars. Any deviation in the mean metallicity is attributed entirely to a zero-point error in the u and v bands because of their larger zero-point uncertainties (Figure 24) and stronger sensitivities on metallicity than other passbands. We also assume the same amount of offset in both passbands. The sensitivity of our metallicity estimate on zero-point is estimated using high-latitude stars (|b| > 60°) by comparing to a case assuming an arbitrary 0.06 mag offset in the u and v bands, from which we find Δ[Fe/H]/Δmag = -2.3. Using this, the metallicity map is forced to match a reference metallicity ([Fe/ H] = -0.28), which is taken from the average metallicity in the above volume. As shown in Figure 25(b), we adopt N side = 16 (a pixel area of 13.4 deg 2 ) for the zero-point correction map because a higher spatial resolution results in the loss of some pixels with small numbers of stars, while information on spatial dependence is lost on a lower-resolution map. In a new metallicity map based on zero-point corrections (panel (c)), a global change of the mean metallicity from low-to highlatitude regions is evident, as expected from a simple population gradient from the disk to the halo. Bottom: a distribution of the ratio between a standard deviation of a magnitude difference (σ(Δmag)) and a propagated uncertainty (σ(phot)) in each (1°, 1°) tile. The mean ratio in each filter passband is shown in the inset. For comparison,Figure 25(d) shows the zero-point corrections in Huang et al. (2021Huang et al. ( , 2022)), who used spectroscopic estimates from the GALAH survey and Gaia parallaxes to derive photometric offsets in each of the SMSS filters. It shows u-band corrections, but similar patterns and amplitudes are seen in the v band. The median difference from our map is negligible (0.004 and 0.001 mag in u and v band, respectively, in the sense of our study minus their values), and a standard deviation amounts to 0.025 mag in both bands. However, high-latitude regions (|b| > -50°) are only sparsely populated by GALAH targets, and therefore their correction functions are weakly constrained in the Galactic pole region. In this respect, our correction map provides more complete information for our chemo-kinematic sample along the Galactic prime meridian.
In the same manner as for the SMSS u and v photometry, SDSS u-band photometry is recalibrated for a small but significant offset in our metallicity map. Figure 26(a) shows a mean metallicity distribution from stars at 1.2 kpc < d < 3 kpc (5.1 million stars). There are strips with distinctly lower or higher metallicities than surrounding areas (e.g., the stripe along l = 30°; see also Figure 1 in Paper I). These strips are parallel to the scanning direction of the SDSS imaging footprints (each 2°. 5 wide and ∼120°long), suggesting that the spatially correlated offsets in metallicity are induced by photometric zero-point errors in the metallicity-sensitive u band.
The zero-point correction map in panel (b) is derived in the same way as for the SMSS u and v passbands. As there are more stars available for the construction of a metallicity distribution function in SDSS, we adopt a finer pixel size for HEALPix, N side = 32. The metallicity sensitivity Δ[Fe/H]/ Δmag = -2.6 is obtained from a case with a 0.06 mag offset in u and is used to make a uniform metallicity distribution of nearby stars (0.5 kpc < d < 1 kpc) at the ensemble average ([Fe/H] = -0.25). In this way, only the relative zero-point offsets are rectified, while the global mean remains intact. The revised metallicity map shown in panel (c) is significantly more smooth than in panel (a) and no longer shows artificial structures. Huang et al. (2022). In all panels, the Galactic center is at the center, and the north Galactic pole is to the top. Areas at low Galactic latitudes (|b| < 20°) and large cumulative extinction (E(B -V ) 0.1) are excluded. Note that the above metallicity maps (panels (a) or (c)) are not exactly the same as in Figure 9, due to different ways of estimating metallicities without and with Gaia parallaxes, respectively.
Figure 27 shows the fractional uncertainties of distance and projected v f of stars used in Section 3. The top panels display Gaia parallaxes and v f measurements computed from Gaia's proper motions based on Gaia parallaxes. In the middle and bottom panels, the same cases are shown from photometric approaches based purely on photometry from SDSS (middle panels) and SMSS (bottom panels) for all stars used in this work. Whenever possible, PS1 photometry is used, as in our main analysis. However, as found in the comparison with Gaia parallax in Figure 8, photometric distance uncertainties are likely overestimated by a factor of about three when photometry is used to estimate both distance and metallicity simultaneously. As v f estimates are linearly dependent on distance, our original uncertainties in v f are also overestimated by the same factor. The distribution in Figure 27 shows rescaled uncertainties based on these factors.
As we move away from the Sun, both geometric and photometric distance measurements become less precise. However, Gaia parallaxes are affected more significantly, with uncertainties deteriorating to 15% beyond a distance of 2.5 kpc. In contrast, photometric distances based on SDSS exhibit a more gradual increase in uncertainties with distance for most stars. SMSS, with its limited survey depth, can provide useful distance estimates only for nearby stars. While there is a second clump of stars with larger distance uncertainties in both SDSS and SMSS samples, due to large errors near the main-sequence turnoff, its impact on the overall sample is negligible. Following the ridge line of the majority of stars in the top and middle left panels, one can see that the photometric approach based on SDSS provides better distances than Gaia beyond 1 kpc.
To ensure the accuracy of our analysis in Section 3, we select a sample of stars with good Gaia parallax measurements and corresponding v f measurements. However, this sample is limited to a nearby volume within a distance of approximately 2 kpc. We use this sample to study high proper-motion stars in Figures 141516, where accurate v f measurements and precise sample cuts based on proper-motion measurements are necessary. Conversely, for a more extensive sample that requires a larger volume coverage, such as in estimating scale height and length in Figures 111213, we rely solely on photometric estimates of distance and v f .
Figure 28 illustrates the distribution of fractional uncertainties in our measurements of scale height and length based on SDSS photometry (Figures 1112). The uncertainties in the northern Galactic hemisphere are relatively small, due to the larger number of stars observed in SDSS. Additionally, we mark the approximate positions of the three major features observed in this work with circled numbers, indicating that such features are not affected by accidental inaccuracies in certain pixels on the phase-space diagram. Uncertainties associated with measurements of the distance (left panels) and projected v f (right panels) of the sample stars used in this analysis. The top to bottom panels display measurements obtained from Gaia EDR3, photometric distance with its associated v f derived from SDSS photometry, and measurements obtained from SMSS photometry, respectively. In the last two cases, PS1 photometry is used in the parameter estimation, when available. The original uncertainty estimates of photometric distance and its associated v f are rescaled by a factor of three (see text). The contour levels on each panel correspond to percentiles that encompass 60%, 70%, 80%, 90%, and 95% of the stars. rescale our original uncertainties by a factor of 4 and 2 in the left and right panels, respectively. This exercise demonstrates that adopting original uncertainties in the analysis fails to produce a clear separation between the thick disk and the halo. This offers additional evidence that our initial uncertainty estimates on [Fe/H] and distance may be overestimated, potentially due to correlations among photometric measurements in various passbands.
Deokkeun An https:/ /orcid.org/0000-0002-8072-7511 Timothy C. Beers https:/ /orcid.org/0000-0003-4573-6233 Young Sun Lee https:/ /orcid.org/0000-0001-5297-4518 Thomas Masseron https:/ /orcid.org/0000-0002-6939-0831
https://www.canfar.net/storage/list/STETSON/Standards.
The Astrophysical Journal, 952:66 (33pp), 2023 July 20 An et al.
https://www.sdss.org/instruments/camera/#Filters
http://svo2.cab.inta-csic.es/svo/theory/fps3
Similarly, because cluster fiducial sequences in PS1(Bernard et al. 2014) are not directly tied to SDSS photometry, color indices are registered using the PS1 r band, instead of the fiducial SDSS r band. Nonetheless, the magnitude offsets from models for the SDSS r band and PS1 r band are nearly the same (see Appendix A), so switching between the two passbands has little impact on the model comparison.