Data-driven Stellar Models

Much in astronomy relies on accurate determination and knowledge of stellar properties. The most accurate methods for determining stellar properties from observables are spectroscopic modeling and astroseismology. However, spectra and precise time-series photometry are unavailable for most objects. Broadband photometry, while providing coarser-grained information about stellar properties, is available in much greater abundance. To good approximation, stars can be characterized by . To constrain stellar parameters with broadband photometry m requires a forward model,

that maps stellar parameters to absolute magnitude M in the observed passbands.

Theoretical stellar models provide one way of mapping from stellar parameters to absolute magnitudes. Beginning with initial mass, age, and abundances, one can use a stellar evolutionary model to predict the bolometric luminosity and radius (which map into the photospheric T_eff and ), along with the atmospheric abundances. Synthetic stellar atmospheric models then map these parameters to broadband photometric magnitudes. These predictions can then be used to infer stellar parameters, as done for example by the StarHorse code (Santiago et al. 2016; Queiroz et al. 2018). A number of systematic effects can, and do, affect these theoretically derived absolute magnitudes, including poorly modeled molecular lines and microturbulence in low-temperature stars, inaccurate instrumental and atmospheric transmission curves, and errors in photometric zero-point calibrations (Casagrande & VandenBerg 2014, 2018). One method of dealing with these inaccuracies in theoretical models is to apply empirical corrections based on the observed photometry of stars of known type.

In this work, we cut out the middleman, and learn the mapping from stellar parameters θ to absolute magnitudes directly from the data, without explicit reference to theoretical models. This method should be simpler, more rigorous and quicker to adapt to new photometric data sets than the more traditional approach of tweaking theoretical models to better match the data.

We derive an empirical mapping from stellar atmospheric parameters to absolute magnitudes, using stars with well-measured parallaxes and determined by spectroscopy (hereafter referred to as “spectroscopic features”). This method can be applied to derive models of any combination of photometric passbands, and can be expanded to include dependence on more detailed element abundances (such as ). Our method represents the function f : θ ↦ M using a simple neural network architecture. At the same time, we must and do learn a simple model of dust extinction in the chosen photometric passbands. Because our mapping from spectroscopic features to photometry is learned directly from observations, the resulting stellar colors are more accurate—by construction—than those produced by theoretical models. This should enable more reliable determination of stellar parameters from broadband photometry, allowing the community to better leverage large photometric surveys.

We demonstrate the performance of this approach with a data set consisting of stars with Gaia parallaxes, optical, and near-infrared photometry from Gaia, the Panoramic Survey Telescope and Rapid Response System 1 (Pan-STARRS 1, hereafter PS1), the Two Micron All Sky Survey (2MASS), and the Wide-field Infrared Survey Explorer (WISE), and spectroscopic type determinations from the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST), the Apache Point Observatory Galactic Evolution Experiment (APOGEE), and Galactic Archaeology with HERMES (GALAH).

This paper is organized as follows. In Section 2, we lay out our machine-learning method in a manner that is agnostic to the specific photometric, parallax, and spectroscopic data sets used. In Section 3, we discuss the composition of our input data set, including spectroscopic features, parallaxes, photometry, and reddening estimates. In Section 4, we discuss the training procedure for our model. In Section 5, we present our trained model, which we validate in Section 6. Section 7 discusses possible uses and extensions of our approach.

Our goal is to learn a mapping from spectroscopic features to absolute magnitude, as well as a vector R that describes how dust affects stellar magnitudes in each passband. We assume that extinction is given by A = E R, where E is dependent on the amount of dust in front of the star. For an infinitesimally narrow passband, the vector R would not depend on spectroscopic features, and would be purely a property of the dust. For realistic passbands, however, R depends weakly on the stellar source spectrum (at the ∼10% level for the Gaia G-band, and less for other bands used in this work—see Appendix B). We will therefore model , with strong regularization to ensure that the dependence on the spectroscopic features is small. With these components, our forward model for observed stellar photometry is given by

where μ is the distance modulus. The causal structure of our model is represented graphically in Figure 1.

Zoom In Zoom Out Reset image size

It is worthwhile to reflect for a moment on why it would be possible to learn such a model. We know from stellar astrophysics that most of the variation in stellar properties is determined by three parameters: initial mass, initial metallicity, and age. There are higher-order parameters that affect stellar properties to a smaller extent (such as rotation and individual elemental abundances). Consequently, stars lie—to rough approximation—in a three-dimensional space. As initial mass and age are not easily measurable, they are not available for a large enough sample of stars to train a data-driven model of stellar photometry. It is, however, possible to measure a different set of parameters, describing the instantaneous present-day properties of the stellar atmosphere, from spectra: effective temperature, surface gravity, and metallicity. While these spectroscopic features determine stellar colors, there is no guarantee that they will map to a unique luminosity, as the latter depends on stellar radius. Fortunately, the mapping from to radius turns out to be nondegenerate over most of the domain, so that absolute magnitudes are nearly uniquely determined by these three features. In addition, the effect of dust extinction along the line of sight is to move stars along a single axis in multiband magnitude space. Thus, in the absence of observational uncertainties, we would expect stars to lie along a four-dimensional manifold in magnitude space. Though we take a data-driven approach in this work, the structure of our model—with unextinguished stellar photometry lying along a three-dimensional manifold in magnitude space, and extinction moving stars along a nearly type-independent axis in magnitude space—is informed by the underlying physics.

We learn both the mapping M(θ) and extinction vector R(θ) from a training data set of stars with spectroscopically determined type, (for example, T_eff, , and , though our method does not depend on the precise choice of parameterization of stellar type), observed apparent magnitudes , measured reddening , and for a subset of the stars, observed parallax . We defined our model above in Equation (2), assuming that θ, μ, and E were exactly known. As this is not true for realistic observations, we Taylor expand our model to first order in the observational errors:

where hatted quantities are noisy measurements of the corresponding un-hatted quantities, and δm, δθ, δϖ, and δE are the errors (measurement truth) in the corresponding measured quantities. This first-order expansion is valid when δθ is small and δϖ ≪ ϖ.

Assume that the observational errors, δE, δϖ, δθ, and δm, have Gaussian distributions and are independent of one another, with (co)variances of , , C_θ, and , respectively. In the first-order Taylor approximation, the observed apparent magnitudes are then also Gaussian distributed, with covariance given by

where i and j refer to passbands i and j, and δ_ij is the Kronecker delta. Using this covariance matrix, we can calculate a Gaussian likelihood for any observed star in our training data set. Figure 2 shows the contributions of each term to the covariance matrix for an example star. Note that this covariance matrix depends on the gradients of our model for M(θ) and R(θ), as well as on the value of R(θ). As we adjust the mappings M(θ) and R(θ), the covariance matrix of must also change.

Zoom In Zoom Out Reset image size

This method works well if σ_ϖ ≪ ϖ. As parallax uncertainty grows, the uncertainty in distance modulus μ becomes increasingly non-Gaussian, and we can no longer assume the uncertainty in apparent magnitude to be Gaussian. However, this non-Gaussianity does not propagate into stellar colors, as they do not depend on distance modulus. Instead of modeling m, we model one apparent magnitude and colors. We define c ≡ B m, where

This makes the first entry in c the apparent magnitude of the first passband, while the rest of the entries are colors: c_i = m_i − m₀. The covariance matrix describing the uncertainty in c is given by

After this transformation, the non-Gaussianity stemming from the parallax-to-distance-modulus conversion is confined to the first element of c, c₀, representing apparent magnitude. For stars with large uncertainties in ϖ, we only model stellar colors, ignoring c₀. Thus, we learn color information from all stars, but only learn about the absolute magnitudes of stars from those with well-measured parallaxes.

This now leaves us to find the functions and that maximize the likelihood of our stellar data:

where . Many classes of fitting functions could be used to model M(θ), such as polynomials or Gaussian processes. Here, we employ a shallow feed-forward neural network. In this context, the neural network should simply be thought of as a highly flexible, regularized fitting function. Our model for is the exponential of a linear function of the spectroscopic features, with strong regularization of the linear weights. In the limit that the weights go to zero, R becomes a constant vector, independent of the spectroscopic features. We apply an exponential activation function to the layer, which forces the entries of R to be positive. Our stellar model is thus a highly flexible fitting function, while our extinction model is highly regularized and expected to be approximately represented by a constant vector. We use the negative logarithm of the likelihood (Equation (7)) as the loss function, with additional regularization terms described below.

As noted earlier, the covariance matrix (Equation (4)) requires knowledge of the gradients of the model functions and . With a trained neural network representing the model, these gradients are easy to calculate for any . However, when we begin training the model, we do not know the values of these gradients. We therefore iterate between training the model while holding the covariance matrices fixed, and updating the covariance matrices while holding the model fixed. When calculating the covariance matrices before the initial iteration, we assume that the gradients are zero and that R = 0, meaning that only uncertainties in observed photometry and parallax contribute to C_m. In iteration n, we use the covariance matrices determined using the model resulting from iteration n − 1.

After each iteration, we additionally update our estimates of the reddening of each star. Given a model and , we can calculate the predicted zero-extinction magnitudes of a star, m_pred, by setting in Equation (3). The covariance matrix of this prediction, C_pred, is given by Equation (4), with . We can then infer the updated stellar reddening using the difference between the observed and predicted zero-extinction magnitudes, taking into account the covariance matrix C_pred. If we put a Gaussian prior on , with mean and variance , then the resulting posterior on is Gaussian, with mean and variance given by

We use the initial estimate of reddening to set the mean and variance of our prior. Keeping the mean of the prior constant across iterations is important in order to avoid a degeneracy in the model between individual stellar reddenings and the length of the vector R. If we were to update the prior after each iteration, the reddenings of all the stars in the training data set could drift either up or down, compensated exactly by an inverse scaling of R. We clip the inferred reddenings, such that . We additionally impose a floor of on the inferred reddening variances.

2.1. Neural Network Structure

The structure of our neural network is illustrated in Figure 3. The vector θ contains the spectroscopic features, with T_eff, , and each shifted and scaled so that the entire input data set has a median of zero and a standard deviation of unity. The sector of the neural network takes in the vector , passing it through two hidden layers to obtain the absolute magnitude in the first passband and the colors relative to the first passband: B M. The sector of the neural network feeds through a single dense layer with exponential activation (in order to ensure positivity), to obtain the vector R. This vector is multiplied by the input to obtain extinction, and then transformed to a mixture of magnitude- and color-space using the matrix B. The network adds the output of the two sectors to obtain , which by Equation (3) should be equal to , up to an error term with covariance C_c, given by Equation (6).

Zoom In Zoom Out Reset image size

Let L be the Cholesky decomposition of :

We operate with L^T on the output of the network, obtaining . We compare this output with our expected output, , using the mean square error (MSE) loss function. This is equal to the χ² error corresponding to the stellar likelihood, divided by the number of passbands:

Minimizing the MSE loss between the output of the network and the expected output is thus equivalent to maximizing the stellar likelihood. During training, the network varies both the stellar model and the extinction vector in order to minimize the sum of χ² over all the stars (equivalent to maximizing the product of all the stellar likelihoods).

The loss function contains regularizing terms in order to prevent the model from overfitting the data. In the dense layers that make up the stellar model (see Figure 3), we apply an L2 weight regularization of 10⁻⁴ to each dense layer. In the dense layer of the extinction model, we instead apply an L1 weight regularization (of 10⁻²), which penalizes even small variations in the reddening vector R, encouraging it to be independent of the spectroscopic features. Specifically, the per-star loss function is

where n_bands = 13 is the number of photometric passbands, W_s is the set of all weights in the stellar model, W_r is the set of all weights in the extinction model, λ_s = 10⁻⁴, and λ_r = 10⁻². During training, we minimize the value of this loss function, averaged over the stars in the training data set.

3.1. Spectroscopic Features

3.2. Photometry and Parallaxes

Our method can be applied to any arbitrary set of photometric passbands. We demonstrate our method using 13 photometric passbands, from Gaia (G, BP, RP), PS1 (g_P1, r_P1, i_P1, z_P1, y_P1), 2MASS (J, H, K_s) and unWISE (W₁, W₂), spanning wavelengths of ∼330 nm to ∼5000 nm. Figure 4 shows how these passbands cover the observed wavelength range.

Zoom In Zoom Out Reset image size

Gaia is a space telescope that is conducting a large-scale astrometric stellar survey, measuring the parallaxes and proper motions of billions of stars to unprecedented precision (Gaia Collaboration et al. 2016). In addition, Gaia provides photometry in a wide band, G, and spectrophotometry in two bands that essentially cover the two halves of the G passband, BP and RP (Liu et al. 2012). Gaia DR2 contains five-parameter astrometric determinations for 1.3 billion sources (Gaia Collaboration et al. 2018; Lindegren et al. 2018). We make use of all three photometric passbands from Gaia, as well as Gaia’s parallax measurement.

We choose Gaia G as the first passband in our vector c, meaning that c₀ represents the G-band apparent magnitude, while the remaining entries in c represent colors with respect to G. Because we learn absolute magnitudes in the first passband, we would like stars with well-measured parallaxes to also be measured in this passband, making G the natural choice. We additionally learn colors relative to the first passband, restricting our sample to stars that are detected in this passband. The range of apparent magnitudes probed by Gaia G spans the range of apparent magnitudes of stars targeted by LAMOST, APOGEE, and GALAH, so that stars observed by these spectroscopic surveys are likely to have Gaia G-band magnitudes. The ordering of the remaining passbands in c makes no difference to our method.

We make the following quality cuts on Gaia sources, as recommended by Arenou et al. (2018):

1. 1.  
   visibility_periods_used > 8,
2. 2.  
   

where m_G is the Gaia G-band magnitude. The first cut removes sources with an insufficient number of Gaia observations to reliably constrain the astrometric model, while the second cut removes sources which are poorly fit by the astrometric model. We additionally require three or more observations in the G-band: g_n_obs > 2.

PS1 is a 1.8 m optical and near-infrared telescope located on Mount Haleakala, Hawaii (Chambers et al. 2016), equipped with Gigapixel Camera #1, consisting of an array of 60 CCD detectors, each with 4800 pixels on each side (Tonry et al. 2006; Onaka et al. 2008; Stubbs et al. 2010; Chambers et al. 2016). From 2010–2014, PS1 conducted a multi-epoch, five-band (g_P1, r_P1, i_P1, z_P1, y_P1) survey of the sky north of decl. δ = −30° (Chambers et al. 2016). We use the mean of the photometric single-epoch PS1 photometry used to create DR1 (Flewelling et al. 2016), excluding PS1 photometry flagged as unphotometric, saturated, or containing bad pixels.

2MASS is an all-sky near-infrared (J, H, and K_s) survey (Skrutskie et al. 2006), which made use of two 1.3 m telescopes at Mount Hopkins, Arizona, and Cerro Tololo, Chile. The focal plane of each telescope was equipped with three 256 × 256 pixel arrays, with a 2′ pixel scale. The entire sky was covered six times with dual 51 ms and 1.3 s exposures, in order to cover a wide range in apparent magnitudes. The survey achieved a 10σ point-source depth of ∼15.8, ∼15.1, and ∼14.3 mag (Vega) in J, H, and K_s, respectively. 2MASS photometry was calibrated to 0.02 mag accuracy, with per-source photometric uncertainties for bright sources below 0.03 mag (Skrutskie et al. 2006). We make use of the 2MASS “high-reliability catalog.”⁴ When 2MASS reports that a source is extended, we remove it entirely from our training catalog.

WISE (Wright et al. 2010) is a 40 cm aperture space telescope, which in 2010 surveyed the full sky in four passbands, centered on 3.4, 4.6, 12, and 22 μm (and named W₁, W₂, W₃, and W₄, respectively), until the depletion of its cryogenic coolant. The spacecraft conducted several months of additional observations in W₁ and W₂, which do not require cryogenic cooling, in order to search for near-earth objects (Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE); Mainzer et al. 2011), and was subsequently put into hibernation. The spacecraft was reactivated in 2013 in order to continue the NEOWISE mission (Mainzer et al. 2014). Since then, WISE has gathered extensive photometry in the W₁ and W₂ bands across much of the sky. The “unWISE” Catalog (Schlafly et al. 2019) is based on a reprocessing of all publicly available WISE W₁ and W₂ photometry. This reprocessing is both deeper than previous WISE catalogs, due to the co-addition of a larger number of observations, and performs better in the Galactic plane, due to its use of a photometric pipeline optimized for crowded fields (Schlafly et al. 2018). We use unWISE photometry for which flags==0 and fracflux > 0.5, and reject photometric bands for which χ²/dof > 5. We convert unWISE magnitudes from Vega to the AB system by adding 2.699 and 3.339 mag to the reported W₁ and W₂ magnitudes, respectively (Schlafly et al. 2019).

In order to remove unrealistically small photometric uncertainties, we transform the reported uncertainties in the passbands from all surveys by adding 0.02 mag in quadrature. We handle unobserved passbands (and those that fail our quality cuts) by setting the photometric uncertainty to infinity. We additionally treat any passband with uncertainty greater than 0.2 mag as if it were unobserved. In order to limit the non-Gaussian error effects described in Section 2 and Appendix A, which occur for sources with weakly constrained parallaxes, we disregard parallax measurements for which . For these sources, we set σ_ϖ to infinity. This renders the uncertainty in absolute magnitude infinite, but does not affect the uncertainties in the colors.

3.3. Reddenings

3.4. Input Catalog

We construct our input catalog by matching each spectroscopic survey separately to the Gaia DR2 catalog, using a 0.2″ matching radius. We require Gaia G-band photometry, and beyond that, opportunistically use whichever Gaia BP/RP, PS1, 2MASS, and unWISE photometry meets our quality cuts for each star. We perform the catalog matching with the Large Survey Database (Juric 2012), using a radius of 0.2″ to match PS1 to Gaia, and a radius of 0.4″ to match 2MASS and unWISE to Gaia. When available, we include high-quality Gaia parallax measurements in the catalog (as described in Section 3.2). For each source, we include the reddening estimates described in Section 3.3. We thus obtain one input catalog for each spectroscopic survey, which we stack to create one catalog. This results in a small number of repeated stars (39,580, or 1.4% of the total), which we expect to have a negligible impact on our results. Our resulting data set has 2,888,361 stars. The number of stars with spectroscopy from each survey and photometry in each passband is given in Table 1. Figure 5 shows the distribution of training data in stellar parameter and reddening space.

Zoom In Zoom Out Reset image size

Table 1.  The Number of Stars with Spectroscopy from the Given Source or Photometry in the Given Passband

	Source	N	% of Total
	LAMOST	2,586,577	90%
Spectroscopy	GALAH	122,411	4%
	APOGEE	179,373	6%
Gaia parallaxes	ϖ	2,442,856	85%
	G	2,888,361	100%
Gaia	BP	2,884,921	99.9%
	RP	2,885,425	99.9%
	gP1	2,360,242	82%
	rP1	2,032,495	70%
PS1	iP1	1,853,962	64%
	zP1	2,121,541	74%
	yP1	2,428,790	84%
	J	2,742,702	95%
2MASS	H	2,687,149	93%
	Ks	2,576,599	89%
unWISE	W1	839,532	29%
	W2	1,857,203	64%

Download table as:  ASCIITypeset image

We split our input catalog into three parts by randomly assigning sources to the training set (70% of sources), validation set (20% of sources), or test set (10% of sources). We minimize the χ² of our model on our training set, using the validation set to track the progress of our training and to diagnose possible overfitting. All results, goodness-of-fit measures and validation plots in the following use the test set.

We implement our model in Tensorflow 2 (Abadi et al. 2015) with Keras (Chollet et al. 2015), using the Adam optimizer (Kingma & Ba 2014). We train the model in 20 iterations, with each iteration consisting of 25 training epochs. In the first iteration, we employ a learning rate of 0.001, and assume that ∇_θM = 0 and R = 0 when calculating the covariance matrix C_m (see Equation (4)). This means that in the initial iteration, uncertainties in the spectroscopic features and reddening do not propagate into the covariance matrix of the photometry. In order to minimize the impact of the missing reddening uncertainty term in the covariance, we exclude stars with σ_E > 0.2 from the initial training iteration. After each iteration, we multiply the learning rate by a factor of e^−1/5. We use the trained network to update our estimate of the reddening of each source (see Equations (8) and (9)) and to recalculate the terms in the covariance matrix that depend on M, R, and . Because and are represented by neural networks, obtaining the gradient terms in the covariance matrix, ∇_θM and ∇_θR, is trivial. Finally, after each iteration, we calculate the χ²/dof of each star. Because we re-estimate the reddening of each star, the number of degrees of freedom is one less than the number of observed passbands. We exclude stars that fail a cut from the next iteration. The threshold of the cut is reduced each iteration, beginning at χ²/dof > 100 after the first iteration, and decreasing logarithmically each iteration, such that it reaches χ²/dof > 5 by the 15th iteration, after which the threshold is held constant. In the final iteration, 2.7% of stars are excluded based on this cut.

We check for overfitting by comparing the loss on the training and validation data sets. A model is overfit if it has learned peculiarities of the training data set, rather than only the features that generalize to new data that the model was not trained on. Thus, over-fitting generally causes the loss obtained on the training data set to be significantly lower than the loss obtained on the validation data set. In our final iteration, the losses are 0.501 and 0.502 on the training and validation data sets, respectively, indicating that our model does not overfit the training data. The contribution of the regularization terms (see Equation (15)) to this loss is small, so that the loss is approximately equal to 〈χ²/n_bands〉 (where the average is taken over stars). Because photometric bands that are not observed do not contribute to χ², but are counted in n_bands, the loss function is significantly less than unity. Taking into account only observed passbands, and disregarding stars with χ²/dof > 5, we find that 〈χ²/dof〉 = 0.87 on the test data set. This is slightly less than the expected value of unity, and likely indicates that some of the error sources (photometric, spectroscopic, parallax, or reddening uncertainties) in our input catalog are systematically overestimated.

Our trained model encodes both the mapping f: θ ↦ M, from the spectroscopic features (T_eff, , and ) to absolute magnitudes in 13 passbands, and the extinction vector R (including a weak dependence on the spectroscopic features). The left panel of Figure 6 shows a Gaia color–magnitude diagram of the observed test data set photometry. We can map the measured spectral types of the stars in the test data set, , to unextinguished absolute magnitudes, to show what the color–magnitude diagram would look like in the absence of dust and photometric errors. The right panel of Figure 6 shows this ideal color–magnitude diagram, as well as the “average” extinction vector (obtained using the mean spectroscopic features in our data set).

Zoom In Zoom Out Reset image size

In Figure 7, we visualize the mapping from spectroscopic features to color–magnitude space. We predict the zero-extinction absolute magnitudes of each star in the test data set, and plot the mean spectroscopic features of stars falling in each region of color–magnitude space. These predictions have reduced aleatoric errors compared to the observations. In particular, they do not contain the final three error terms in Equation (3): photometric errors, error from the reddening estimate, and parallax errors. The only aleatoric error term which is present is the error propagated from the spectroscopic features to the absolute magnitudes:

By construction, however, this error term does not move a star’s photometry off of the stellar locus. It moves each star in a random direction along the stellar locus, but an ensemble of stars plotted in this manner will still fall along the modeled stellar locus. The epistemic errors deriving from any inaccuracies in our model itself remain in these reduced-error predictions.

Zoom In Zoom Out Reset image size

By leveraging observations of hundreds of thousands of stars, our model is able to learn the shape of the de-noised stellar locus with far greater precision than the typical photometric uncertainties of the input data. This can be illustrated by comparing observed, noisy photometry with predicted zero-extinction, de-noised photometry. The left panel of Figure 8 shows observed test-set photometry, colored by . In the right panel, we use the observed spectroscopic features of each star to predict its photometry. The projected reddening vector is plotted for each combination of colors. Figure 9 shows the same results, colored by . We recover color trends with that are not immediately apparent in the noisy observational data, with more metal-rich stars being slightly redder in g_P1 − r_P1 than metal-poor stars. This color trend with metallicity can also be seen in the right panel of Figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 10, we show the dependence of our inferred extinction vector, R on the spectroscopic features for a subset of photometric passbands used in this work. On the y-axis, we plot the fractional variation of the components of the reddening vector,

where R₀ is the value of R calculated at the mean spectroscopic features in our input catalog. The x-axis shows T_eff. Each pixel in the plot is colored by the average of stars falling in the pixel. We overplot the theoretical result obtained using a blackbody spectrum, measured transmission curves for each passband (Cohen et al. 2003; Wright et al. 2010; Rodrigo et al. 2012; Tonry et al. 2012; Rodrigo & Solano 2013; Maíz Apellániz & Weiler 2018), and the Fitzpatrick (1999, “F99”) extinction relation. As expected, the variation in the reddening vector is primarily dependent on T_eff, as temperature has the largest effect on the overall shape of the stellar spectrum. The reddening vector varies the most (by ±15%) in Gaia G, due to the width of the passband (see Appendix B). The variation in the reddening vector is similar to theoretical expectations in the bluest passbands, where the variation is strongest: G, BP, RP, and g_P1. However, our model learns a slight negative trend in ΔR/R with temperature in i_P1, z_P1, and y_P1, which is unphysical, and likely due to unmodeled systematics in the training data. Our model also learns a slight difference in the reddening vector for dwarfs and giants, which manifests itself as a dependence of R on , and which may be unphysical. This effect is also greatest in Gaia G, where R differs by ∼7% for 4000 K dwarfs and giants. The unphysical relations learned by the model can be suppressed with stronger L1 weight regularization in the extinction model. By contrast, the strong trends in ΔR/R with temperature that our model learns for G and BP—which roughly match theoretical expectations—are robust even under far stronger regularization.

Zoom In Zoom Out Reset image size

We can also compare our inferred extinction vector to previous results. In Figure 11, we compare our extinction vector to the family of vectors derived in Schlafly et al. (2016, “S16”). S16 works in color-space, rather than with absolute magnitudes, and therefore derives the reddening—as opposed to extinction—vector. We therefore compare the slope of our inferred reddening curve, defined as the ratio of reddening in subsequent colors, such as . In most color combinations, the slope of our reddening vector falls close to that of the mean (R_V ≈ 3.3) S16 reddening vector. However, we do see a somewhat different slope in two color combinations, and .

Zoom In Zoom Out Reset image size

Our trained neural network model and a tutorial on how to use it are available at doi:10.5281/zenodo.3902382.

In order to validate the performance of our model in explaining the data, we investigate the behavior of residuals in predicted stellar colors in the test data set, as a function of the spectroscopic features and reddening. For each color, X − Y, we consider the score of the prediction:

where is determined from the covariance matrix C_m (see Equation (4)), and the predicted colors are determined from the observed spectroscopic features and estimated reddening using the trained model and learned extinction vector: . We use the best-fit reddening estimates for this comparison, determined according to Equations (8) and (9).

Figure 12 shows the distribution of these scores for four colors, as a function of the spectroscopic features and estimated reddening. For colors that do not involve the Gaia G passband, the median residuals are close to zero, and there are no large trends as a function of the spectroscopic features or reddening. However, we observe trends in the median residuals of the r_P1 − G color with reddening, E. Similar trends are observed in all colors involving Gaia G. These trends are presumably due to the width of the Gaia G passband. In particular, our assumption that extinction is described by a vector R(θ) breaks down for extremely wide passbands, as extinction depends nonlinearly on the column density of dust. This effect is described in more detail in Appendix B.

Zoom In Zoom Out Reset image size

For most colors, the 16th and 84th percentiles of the scores approximately span the range ±1, as expected. In some colors, such as g_P1 − r_P1, the scores are closer to zero than expected, indicating that our uncertainties are overestimated. This could be due to the 0.02 mag floor we place on photometric uncertainties, to our floor on reddening uncertainty, or to the fact that we use diagonal covariance matrices for the spectroscopic features.

At low temperatures (T_eff ≲ 4500 K) and surface gravities (), the g_P1 − r_P1 residuals become larger, possibly due to decreased reliability in the spectroscopic features in this regime. In particular, line blanketing by molecular species at low temperatures renders accurate determination of spectroscopic features more difficult.

Our model is valid only in regions of θ-space in which we have training data. In particular, the spectroscopic pipelines from which we draw our features are tailored to stars with 4000 K ≲ T_eff ≲ 7500 K. The region covered by our training data also differs for each photometric passband, and depends on the overlap between photometric and spectroscopic surveys. We quantify the region of θ-space where our photometric predictions should be accurate, by measuring the density of stars in the input catalog. For each photometric passband, we construct a kernel density estimate of the input stars in θ-space, using a bandwidth of (50 K, 0.05 dex, 0.05 dex) in . For each passband, we normalize the peak density to unity. Figure 13 shows the dependence of photometric residuals on this input density for a subset of passbands. As expected, photometric residuals are smallest in the regions of spectroscopic-feature space with the highest density of training data. We find that our color predictions degrade more gracefully with increasing input density than our predictions of absolute G-band magnitude. One may set a density threshold below which one considers our model predictions unreliable for a given passband, based on the acceptable level of scatter in the photometric residuals. At an input density (for parallax observations) of 10⁻², . Residuals in X − G colors reach this level at input densities (in band X) of approximately 10⁻³. As a general rule of thumb, we therefore recommend to apply an input density threshold of 10⁻² when absolute magnitudes are critical, and of 10⁻³ when only color predictions are required.

Zoom In Zoom Out Reset image size

In this paper, we have learned a model that maps spectroscopic features (T_eff, , and ) to photometry, and at the same time describes the effect of dust on stellar photometry. Because this model is learned directly from the data, rather than using synthetic spectra and models of atmospheric, filter, and instrument transmission, it sidesteps potential problems in these spectral and transmission models. With a large number of input stars, we are able to average down the statistical errors in the photometry, and learn stellar magnitudes to better than the typical photometric uncertainties of individual stars in our training data. Our method is not restricted to any specific set of photometric passbands, and can be applied to develop accurate mappings from spectroscopic features to absolute magnitudes for whichever set of passbands one wishes to use for a given study. The only requirement is that the new photometric survey overlap sufficiently with spectroscopic surveys, so that it is possible to construct a training data set. Because our model has no concept of stellar spectra, filter transmission curves, atmospheric transmission, and other physical effects that underlie the mapping from spectroscopic features to photometry, it is unable to extrapolate to passbands that are not present in the training data. However, our model is not computationally expensive to train. Provided that training data is available, training a new model with additional passbands is a quick process.

Our trained model can be used to infer the spectroscopic features, distances, and extinctions of large numbers of stars with observed photometry—in the absence of spectral data. Using Bayes’ theorem, the posterior density of θ, μ, and E, given observed photometry , is given by

where is a normalizing constant. Because our model maps from spectroscopic features to absolute magnitudes and the reddening vector, R, it can be used to construct an accurate likelihood function, . The above inference additionally requires a prior on the distribution of stars and dust throughout the Milky Way, . There are many Galactic models that one could choose to construct such a prior, such as the Besançon model (Robin et al. 2012) or the SDSS tomography models (Ivezić et al. 2008; Jurić et al. 2008). More ambitiously, the Galactic model itself could be inferred hierarchically, using the distribution of a large number of observed stars and detailed knowledge of the selection function. However, the problem of inferring stellar parameters from photometry exceeds the scope of this paper, which is instead focused on developing a highly accurate mapping from spectroscopic features to stellar photometry. One promising area for the application of our data-driven stellar model is three-dimensional dust mapping, which often relies on hundreds of millions of stellar distances and reddenings inferred from photometric surveys (Andrae et al. 2018; Chen et al. 2019; Green et al. 2019; Lallement et al. 2019; Leike & Enßlin 2019).

In order to infer stellar parameters, we would like to construct models that cover as wide as possible a range of stellar parameters. With the increasing availability of spectroscopic data from ongoing surveys, such as LAMOST, APOGEE, and GALAH, and from upcoming surveys such as SDSS-V (Kollmeier et al. 2019), it will be possible to train models that span a wider range of spectroscopic features.

Our model also provides a clean and accurate way to convert stellar evolutionary models to photometric predictions. Stellar evolutionary models begin with intrinsic stellar properties, such as initial mass and chemical abundances, and predict the evolution of stellar interior structure, abundances, and atmospheric properties with age. In order to predict broadband photometry, synthetic atmospheric models are typically used, often with empirical corrections to better match observed photometry. Our models provide a simple alternative—and by construction accurate—method of mapping from the spectroscopic features θ output by stellar evolution models to stellar photometry.

Accurate stellar loci are also critical to developing color transformations between different photometric surveys. These color transformations are necessary for cross-calibration of photometric surveys (Scolnic et al. 2015; Finkbeiner et al. 2016). For example, Finkbeiner et al. (2016, “hypercalibration”) recalibrates SDSS photometry using PS1, and develops an empirical color transformation between the two surveys that is applicable to main-sequence stars in low-reddening regions. This procedure encounters difficulties in the SDSS u-band, however, due to the stronger metallicity dependence of stellar photometry in the near-ultraviolet. Accurate, metallicity-dependent color transformations derived using the methods presented in this work could therefore be of use in improving such cross-survey calibrations.

In this work, we have considered only the simplest possible extinction model, in which a single vector—with weak dependence on the spectroscopic features—describes the effect of dust on photometric magnitudes in the chosen passbands. However, the direction of this vector should depend on variation in the dust extinction spectrum itself, often parameterized by R_V (Fitzpatrick 1999). Here, we propose two possible avenues for future exploration of variation in the dust extinction spectrum. The first method is to model the residuals between the observed and predicted photometry, similarly to the analysis in S16. This is an after-the-fact analysis, which can be applied to the residuals generated by the machinery developed in this paper. The second method is to build a model of the variation in the extinction vector directly into our neural network. This could be achieved by modeling the extinction of star i as a sum of two or more components:

A separate dense neural network layer would represent each vector R_k, while the extinction of each star i would be described by K coefficients, E_ki, k = 1, …K. The vector R₁ would represent a mean extinction relation, while R₂ and higher-order vectors would describe variation in dust extinction properties. With appropriate regularizing conditions (e.g., orthogonality between the vectors ), such a model could be trained to obtain a more flexible and accurate description of dust extinction. Such a model could also be used as a basis to explore the spatial variation of dust properties throughout the Milky Way. We leave this extension of our model to future work.

In this work, we have presented a method for learning a mapping from spectroscopic features to absolute magnitudes. Our model additionally learns an extinction vector representing the effect of dust on stellar magnitudes in the chosen passbands. Our method can be easily adapted to an arbitrary set of photometric passbands or spectroscopic data sets, allowing inference of stellar parameters from a range of different photometric surveys.

We have demonstrated the effectiveness of our model on a data set consisting of spectroscopic features from LAMOST, APOGEE, and GALAH, parallaxes from Gaia, photometry from Gaia, PS1, 2MASS and WISE, and reddening estimates from SFD and Bayestar19. We obtain precision de-noised, zero-extinction stellar loci in the resulting three-dimensional magnitude space, along with a model of the reddening vector R and its (weak) dependence on spectroscopic features.

We provide our trained neural network model, our estimate of its region of validity, and a tutorial on how to use the model at doi:10.5281/zenodo.3902382.

This method is perhaps the simplest and most accurate way to convert stellar evolutionary models into photometric predictions. It also provides a simple method for constructing a highly accurate likelihood function of stellar photometry, given spectroscopic features, distance, and extinction. This is a key component of methods which infer stellar parameters, distances, and extinctions from photometric data. One immediate application of our method will be to three-dimensional dust mapping in the Milky Way, though we expect it to be of use in other fields, such as Galactic archeology and stellar population synthesis.

The computations in this paper were run on the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvard University. Jan Rybizki, Morgan Fouesneau, and René Andrae were funded by the DLR (German space agency) via grant 50 QG 1403. This work is partly based on data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Software: Astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), dustmaps (Green 2018), Keras (Chollet et al. 2015), Large Survey Database (Juric 2012), Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), SciPy (Jones et al. 2001), Tensorflow 2 (Abadi et al. 2016).

In specifying our model, we calculate a mean distance modulus from a noisy measured parallax: . As noted in the paper, we only attempt such a conversion when the parallax uncertainty is small. Here, we derive an estimate of the mean distance modulus in this regime.

Defining

and Taylor expanding the second logarithm in Equation (A2) in terms of x, we obtain

The nth moment of x is given by

where

In order to evaluate K_n, we assume that the likelihood of is Gaussian, and we Taylor expand the prior, , around . We additionally assume that , so that we can extend the lower limit of the integrand to −∞ with little change to the value of the integral. This turns K_n into a sum over Gaussian integrals, which can be exactly evaluated:

To second order in ,

Thus, the first two moments of x are given by

All higher-order moments of x depend on or higher-order terms. Plugging the moments of x into Equation (A4), we obtain

The lowest-order correction to the naive distance modulus thus depends on the slope of the prior on parallax at the location of the measured parallax. The sign of the correction depends on the choice of the parallax (or equivalently, of the distance) prior, and scales with . For the cutoff we make when calculating the distance modulus, , we can therefore expect the maximum correction on the order of several hundredths of a magnitude. As the sign of the correction depends on the exact choice of the distance prior, we do not attempt to correct for the effect in this work, and instead use the naive conversion from measured parallax to distance modulus. The errors incurred in this manner only affect absolute magnitudes, but not colors, which are independent of distance.

In this work, we assume that the effect of interstellar extinction can be described by a nearly universal vector, R, with only weak dependence on the spectroscopic features. As in Equation (2), the extinction of each star is proportional to this vector. This is an idealization that holds only when the wavelength-versus-dust-opacity relation does not change, and in the limit that the passbands we are working with are extremely narrow. However, the Gaia G passband covers a wide range of wavelengths, and therefore exhibits measurable nonlinear extinction (versus dust column density) and stellar-spectrum-dependent extinction. In the following, we show how these effects arise, and estimate the nonlinear extinction term.

The AB magnitude measured by a charge-coupled device (CCD) is given by

where is the transmission of the instrument (including the optics, filter, and CCD, and depending on convention, possibly the atmosphere) as a function of frequency, is the spectral flux density of the source, and the factors of correspond to the fact that CCDs count photons, rather than measuring incident energy. In the presence of interstellar dust with optical depth , the extinction is given by

where we are now using to denote the source spectrum in the absence of dust. For an infinitely narrow passband, , extinction is given by

In this limit, the vector R, describing the relative amount of extinction in each passband, can be calculated directly from the opacity of the dust at the central frequencies of the passbands, and does not depend on the spectroscopic features.

In order to explore the effect of passbandwidth on extinction, consider a Gaussian passband with width Δν:

The number of photons detected by the CCD is then given by

where is the opacity of the dust, and σ is the dust column density, such that optical depth is given by . Taylor expanding around ν = ν₀, we obtain

Extending the lower limit of the integral to −∞ allows us to evaluate it exactly, yielding

where k!! is the double factorial of k. As long as Δν ≪ ν₀, changing the lower limit of the integral impacts the result minimally. The extinction is given by

At second order in Δν, there are two separate effects to consider:

1. 1.  
   In the limit of low-optical depth, the ratio of extinction to optical depth depends on the source spectrum.
2. 2.  
   As optical depth increases, the nonlinear behavior of extinction with optical depth depends on the dust-opacity spectrum.

We can explore the second effect—the nonlinear dependence of extinction on dust column density—by splitting into a linear regime at low column density and nonlinear correction at larger optical depths:

Generically, ΔA ∝ A² for any dust-opacity law at second order in Δν. At this order in Δν, the quadratic coefficient is determined by the behavior of the dust-opacity relation, , with the shape of the source spectrum providing corrections at higher order in Δν. At second order in Δν, the effect of the source spectrum is to alter the slope of the extinction-versus-optical-depth relationship at low dust column density.

With a power-law dust opacity, , one obtains

For small dust optical depths, extinction increases proportionally to optical depth in all passbands. However, as the dust column increases, the widest passbands deviate more quickly from this linear relationship than the narrower passbands, so that the ratio of extinction in different passbands is no longer constant. The vector R is no longer sufficient to determine the relationship between extinction in the different passbands.

We can estimate the scale of this effect for Gaia G-band. Here, we use the Gaia G transmission curve calculated by Maíz Apellániz & Weiler (2018) and the Fitzpatrick (1999) extinction curve, improved by Indebetouw et al. (2005) in the near-infrared (Rodrigo et al. 2012; Rodrigo & Solano 2013). We set ν₀ and Δν to the transmission-weighted mean and standard deviation of frequency, respectively, and estimate γ from the shape of the extinction curve between ν₀ − Δν and ν₀ + Δν, obtaining ν₀ = 640 nm, Δν = 144 nm, and γ = 1.31. Plugging these values into Equation (B14), we estimate . This extinction nonlinearity therefore becomes comparable to our photometric uncertainty floor of 0.02 mag at an extinction of A_G ∼ 0.7 mag.

Numerically evaluating extinction using Equation (B2) yields similar results. We define ΔA to be the difference between the calculated extinction at a given optical depth and the extinction that would be obtained by linearly extrapolating extinction from the low-optical-depth limit. With a blackbody source spectrum at 6000 K, we obtain ΔA_G = −0.037 mag when A_G = 1 mag, closely matching our second-order perturbative result. The numerical result only varies by 20% for source temperatures between 4000 and 8000 K.