<figure><img src="image_1.png" /></figure>

THE OPTICAL-INFRARED EXTINCTION CURVE AND ITS VARIATION IN THE MILKY WAY

E. F. SCHLAFLY<sup>1,2,14</sup>, A. M. MEISNER<sup>2,3</sup>, A. M. STUTZ<sup>1</sup>, J. KAINULAINEN<sup>1</sup>, J. E. G. PEEK<sup>4</sup>, K. TCHERNYSHYOV<sup>5</sup>, H.-W. RIX<sup>1</sup>, D. P. FINKBEINER<sup>6,7</sup>, K. R. COVEY<sup>8</sup>, G. M. GREEN<sup>6</sup>, E. F. BELL<sup>9</sup>, W. S. BURGETT<sup>10</sup>, K. C. CHAMBERS<sup>11</sup>, P. W. DRAPER<sup>12</sup>, H. FLEWELLING<sup>11</sup>, K. W. HODAPP<sup>11</sup>, N. KAISER<sup>11</sup>, E. A. MAGNIER<sup>11</sup>, N. F. MARTIN<sup>1,13</sup>, N. METCALFE<sup>12</sup>, R. J. WAINSCOAT<sup>11</sup>, AND C. WATERS<sup>11</sup>

Received 2015 December 9; accepted 2016 March 4; published 2016 April 12

ABSTRACT

The dust extinction curve is a critical component of many observational programs and an important diagnostic of the physics of the interstellar medium. Here we present new measurements of the dust extinction curve and its variation toward tens of thousands of stars, a hundred-fold larger sample than in existing detailed studies. We use data from the APOGEE spectroscopic survey in combination with ten-band photometry from Pan-STARRS1, the Two Micron All-Sky Survey, and *Wide-field Infrared Survey Explorer*. We find that the extinction curve in the optical through infrared is well characterized by a one-parameter family of curves described by $R(V)$. The extinction curve is more uniform than suggested in past works, with $\sigma(R(V)) = 0.18$, and with less than one percent of sight lines having $R(V) > 4$. Our data and analysis have revealed two new aspects of Galactic extinction: first, we find significant, wide-area variations in $R(V)$ throughout the Galactic plane. These variations are on scales much larger than individual molecular clouds, indicating that $R(V)$ variations must trace much more than just grain growth in dense molecular environments. Indeed, we find no correlation between $R(V)$ and dust column density up to $E(B-V) \approx 2$. Second, we discover a strong relationship between $R(V)$ and the far-infrared dust emissivity.

Key words: dust, extinction – ISM: clouds – ISM: structure

1. INTRODUCTION

Dust is composed of small, solid grains of material. These grains are made of heavy elements formed from nuclear fusion in stars, blown into the interstellar medium (ISM) by stellar winds and explosions. The dust grains scatter and absorb light. Owing to the small size of the grains, dust preferentially scatters and absorbs blue light relative to red light in the optical through infrared. The resulting extinction as a function of wavelength is called the dust extinction curve (Draine 2003).

Effective parameterizations of the dust extinction curve in the ultraviolet (UV) and optical were developed by Fitzpatrick & Massa (1986, 1988). The work of Cardelli et al. (1989, CCM) showed that much of the variation could be described by a single parameter, $R(V) = A(V)/E(B-V)$, the total-to-selective extinction ratio, though especially in the UV this description is far from complete.

The shape of the extinction curve is a valuable diagnostic of the properties of the dust. Variation in $R(V)$ is sometimes attributed to variation in the size distribution of dust grains; dust with high $R(V)$ has a relatively gray, flat extinction curve in the optical, suggesting an abundance of large grains relative

to small. Variation in $R(V)$ may also be related to the formation of ice on dust grains or grain aggregation in dense environments ($E(B-V) > 1$) (Whittet et al. 1988; Ysard et al. 2013). Current observational evidence suggests a relationship between $R(V)$ and $E(B-V)$, though the strength of the correlation is not clear (Fitzpatrick & Massa 2007; Foster et al. 2013). Alternatively, variation in $R(V)$ may be driven by grain processing by UV photons or grain composition and chemistry (Jones et al. 2013; Mulas et al. 2013).

Many studies of the extinction curve have focused on relatively small samples of O and B stars or on particular molecular clouds. For example, the best atlases of Milky Way extinction curves are those of Valencic et al. (2004) and Fitzpatrick & Massa (2007), which have samples of a few hundred O and B stars. These stars are practical targets for studies of the extinction curve in the UV, owing to their intrinsic UV brightness. However, the small number of these stars and the atypical environments they inhabit make their use to characterize the variability of the extinction curve problematic. Additionally, the sparseness of appropriate bright O and B targets complicates morphological association of the observed extinction curve variations with known structures in the ISM, inhibiting efforts to uncover the underlying physical processes at work.

<footnote n="1"><sup>1</sup>Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="2"><sup>2</sup>Lawrence Berkeley National Laboratory, One Cyclotron Road, Berkeley, CA 94720, USA</footnote>

<footnote n="3"><sup>3</sup>Berkeley Center for Cosmological Physics, Berkeley, CA 94720, USA</footnote>

<footnote n="4"><sup>4</sup>Space Telescope Science Institute, 3700 San Martin Dr, Baltimore, MD 21218, USA</footnote>

<footnote n="5"><sup>5</sup>Department of Physics and Astronomy, The Johns Hopkins University, 3400 North Charles Street, Baltimore, MD 21218, USA</footnote>

<footnote n="6"><sup>6</sup>Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="7"><sup>7</sup>Department of Physics, Harvard University, 17 Oxford Street, Cambridge, MA 02138, USA</footnote>

<footnote n="8"><sup>8</sup>Department of Physics and Astronomy, Western Washington University, 516 High St., Bellingham, WA 98225, USA</footnote>

<footnote n="9"><sup>9</sup>Department of Astronomy, University of Michigan, 500 Church St., Ann Arbor, MI 48109, USA</footnote>

<footnote n="10"><sup>10</sup>GMTO Corporation, 251 S. Lake Ave., Suite 300, Pasadena, CA 91101, USA</footnote>

<footnote n="11"><sup>11</sup>Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA</footnote>

<footnote n="12"><sup>12</sup>Department of Physics, Durham University, South Road, Durham DH1 3LE, UK</footnote>

<footnote n="13"><sup>13</sup>Observatoire Astronomique de Strasbourg, CNRS, UMR 7550, 11 rue de l'Université, F-67000 Strasbourg, France</footnote>

<footnote n="14"><sup>14</sup>Hubble Fellow.</footnote>
Many studies of the extinction curve have focused on only a small range of wavelengths for practical reasons, for example, focusing either on only the optical (e.g., Schlafly et al. 2010) or only the near-infrared (NIR) (e.g., Wang & Jiang 2014) extinction curve. This means that individual parts of the extinction curve are known much better than the connections between these parts.

Perhaps unsurprisingly, then, current descriptions of the extinction curve are in substantial tension. For example, the CCM extinction curve contains no variation in the infrared, while the extinction curve of Fitzpatrick & Massa (2009) is more variable in the NIR than in the red-optical. Different authors have found both “universal” NIR extinction (e.g., Wang & Jiang 2014) and variable infrared extinction (e.g., Zasowski et al. 2009). Meanwhile optical studies essentially always find extinction curve variability, though often with the caveat that “most” extinction curves are compatible with an “average” Milky Way extinction curve (e.g., Krewski & Strobel 2012). That said, the detailed shape of the optical variation differs significantly between, for example, the extinction curves of CCM, Fitzpatrick & Massa (2007), and Maíz Apellániz et al. (2014).

Here we present a far more comprehensive, multiwavelength study of the Galactic dust extinction curve and its variation, combining data from spectroscopic and photometric surveys. The APOGEE survey has spectroscopically observed about 150,000 stars in the Galactic midplane, obtaining accurate temperatures, metallicities, and gravities (Majewski et al. 2015). The Pan-STARRS1 (PS1) survey has photometrically observed the entire sky north of declination -30°, providing optical photometry for essentially all of these stars. Infrared photometry from the Two Micron All-Sky Survey (2MASS; Skrutskie et al. 2006) and Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) complement the optical photometry, providing coverage over a factor of ten in wavelength. The APOGEE targets span Galactic longitudes from roughly 0° < *l* < 240°, with typical distances between 1 and 5 kpc. The broad coverage of the Galactic plane provides an excellent test for studying the extinction curve and its variation throughout the Milky Way.

We study extinction via a generalization of the “pair method,” where unextinguished “standard” stars are compared with extinguished stars of the same spectral type to assess the extinction to those stars. This technique was pioneered by Trumpler (1930), but large spectroscopic surveys have recently allowed it to be applied to hundreds of thousands of stars (e.g., Schlafly & Finkbeiner 2011; Yuan et al. 2013). It is impossible to apply the pair method directly to the APOGEE targets, because APOGEE has observed very few unextinguished stars; almost all of the sight lines observed are heavily extinguished. Moreover, the few unextinguished stars in the survey are typically low metallicity halo stars, while the extinguished stars are metal-rich disk stars, complicating the comparison of unextinguished standards with extinguished targets. We circumvent this difficulty by focusing on the *shape* of the extinction curve, leaving the extinction to any individual star relatively unconstrained: we compare stars of the same stellar types, but behind different amounts of extinction to determine the extinction curve.

An alternative approach to solving this issue using synthetic stellar spectra has been extensively explored in the literature (Fitzpatrick & Massa 2005, 2007, 2009; Schultheis et al. 2014,

2015; Clayton et al. 2015). In this approach, the problem of finding unreddened standard stars is avoided by using theoretical stellar spectra instead of observed spectra. We eschew this approach here in order to limit our exposure to any systematic errors in the theoretical spectra. In past work we have found such errors to be of the order of a few percent (Schlafly & Finkbeiner 2011), which is a few times larger than typical observational uncertainties.

This paper is divided into several sections. First, in Section 2, we discuss the observational data: the APOGEE spectroscopy, and PS1, 2MASS, and WISE photometry. In Section 3, we discuss our technique for studying the extinction curve and its variation with this data. In Section 4, we show the results of our analysis. In Section 5, we compare these results with the literature and discuss the consequences for our understanding of dust in the Milky Way. Finally, in Section 6, we conclude.

## 2. DATA

### 2.1. APOGEE

The APOGEE survey is a high-resolution ($R=22500$) spectroscopic, NIR H-band survey of the sky (Majewski et al. 2015). The APOGEE spectrograph (Wilson et al. 2010) is illuminated with 300 fibers from the Sloan Digital Sky Survey (SDSS) 2.5 m telescope (Gunn et al. 2006). The APOGEE survey is largely focused on obtaining abundances of distant giant stars in the Galactic disk, and so observes to a high signal-to-noise ratio of 100 per half-resolution element. The APOGEE spectroscopic pipeline (Nidever et al. 2015) and astrophysical parameter pipeline (García Pérez et al. 2015) measure temperatures, metallicities, and gravities for the stars, with typical uncertainties of <100 K, <0.05 dex, and <0.15 dex, respectively. The main APOGEE survey targets red clump and red giant stars based on their dereddened NIR color (Zasowski et al. 2013). The APOGEE spectroscopic parameters are determined from continuum-normalized spectra, and are therefore insensitive to the extinction of the source. This work uses data from SDSS-III Data Release 12 (Eisenstein et al. 2011; Alam et al. 2015; Holtzman et al. 2015), containing spectra of more than 150,000 stars.

### 2.2. WISE

The WISE is a NASA infrared space telescope that has surveyed the entire sky at 3.4 (W1), 4.6 (W2), 12 (W3), and 22 µm (W4) wavelengths (Wright et al. 2010). We use W1 and W2 photometry from the AllWISE data release (Cutri et al. 2013), which contains data from two complete sky coverage epochs and has identified more than 700 million objects. Approximately one third of all WISE detections of APOGEE targets have flags indicating problematic conditions (cc_flags), usually due to being in the halo of a nearby bright star or landing on a diffraction spike. We conservatively assign infinite uncertainties to the WISE photometry of these sources.

### 2.3. 2MASS

The 2MASS was a NIR survey of the entire sky in the *J* (1.25 µm), *H* (1.65 µm), and *K* (2.17 µm) bands, undertaken from 1997 to 2001 (Skrutskie et al. 2006). The 2MASS survey contains observations of more than 300 million objects, including all objects targeted by APOGEE.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** *E'* (roughly *E*(*B* − *V*)) to APOGEE targets. The best fit *E'* is shown as colored points. The background grayscale shows the Planck τ<sub>353</sub>-based extinction map (Planck Collaboration et al. 2014) for context, and ranges from 0 to 2.5 mag *E'*. Unsurprisingly, the most reddened APOGEE targets are at low latitudes in the inner Galaxy.</figcaption></figure>

## 2.4. Pan-STARRS1

We obtain optical photometry of APOGEE targets from the PS1 survey. The PS1 observations are made on a 1.8 m telescope on Haleakala (Hodapp et al. 2004). The telescope focal plane is outfitted with the 1.4 billion pixel GPC1 camera (Hodapp et al. 2004; Onaka et al. 2008; Tonry & Onaka 2009), which covers the 3° field of view of the telescope. Observations are performed in five broad passbands, covering about 400 nm to 1 μm (Stubbs et al. 2010). The effective wavelengths of the filters are roughly 480, 620, 750, 870, and 960 nm, for the $g_{\text{P1}}$, $r_{\text{P1}}$, $i_{\text{P1}}$, $z_{\text{P1}}$, and $y_{\text{P1}}$ filters, respectively. The data are automatically processed by the PS1 image processing pipeline (Magnier 2006), which analyzes images to deliver photometry, astrometry, and morphology (Magnier 2007; Magnier et al. 2008, 2013). The photometric calibration of the survey, both relative and absolute, is accurate to better than 1% (Schlafly et al. 2012; Tonry et al. 2012).

Many stars in APOGEE are too bright to have reliable photometric magnitudes in PS1. In our analysis, we mark the uncertainties of any measurements brighter than 14.0, 14.4, 14.4, 13.8, and 13.3 magnitudes in the $g_{\text{P1}}$, $r_{\text{P1}}$, $i_{\text{P1}}$, $z_{\text{P1}}$, and $y_{\text{P1}}$ bands as infinitely large to reflect the increased systematic uncertainties brightward of these limits.

## 2.5. Target Selection

This work uses only 37,000 stars of the full 150,000 stars in APOGEE DR12. The most significant restriction we adopt is to use only stars for which we estimate the PS1 $y_{\text{P1}}$ magnitude is fainter than 13th magnitude, estimated from the 2MASS colors and the Rayleigh–Jeans Color Excess estimated extinction (Zasowski et al. 2013). Stars not passing this cut will be saturated in many of the optical bands, limiting their use in characterizing the extinction curve. We also exclude any stars targeted as part of APOGEE “ancillary” programs, to include only stars from the main survey. Ancillary program stars may have selection criteria making them unsuitable for reddening studies. We further exclude any stars where accurate gravities could not be determined, or where the APOGEE flags indicate that stars had an abnormally high $\chi^2$ or rotation rate. Lastly, we exclude any stars whose PS1 or WISE position is more than 0.0′5 separated from its 2MASS position. These cuts reduce the more than 150,000 stars in APOGEE DR12 to under 40,000

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Histogram of $E'$ (roughly $E(B-V)$) to APOGEE targets. The typical $E'$ is roughly 0.65, though the distribution is broad and ranges from 0 to more than 5 mag.</figcaption></figure>

stars. Most of the stars are removed by the brightness cut (150,000 → 72,000). The requirement of a reliable gravity reduces the number further to 47,000, and the elimination of ancillary targets, stars with problematic flags, and astrometric separation reduces the number to 37,000, our full data set.

We show in Figure 1 the locations of the stars used in this work. Points are colored by their estimated reddenings $E'$ (roughly, $E(B-V)$), as determined in this work. The background image shows the Planck τ<sub>353</sub>-based extinction map (Planck Collaboration et al. 2014) for context. Figure 2 shows a histogram of the reddenings to these stars. The mean reddening is 0.65 mag with a standard deviation of 0.5, though the distribution extends to $E' > 5$ mag. Not all stars are detected in all photometric bands: for instance, when $E' > 2.5$ mag, most stars are not detected in $g_{\text{P1}}$.

## 3. METHOD

Our basic technique is to model the observed ten optical–infrared magnitudes of the APOGEE stars as a function of their temperature and metallicity, as well as their distance and reddening. The analysis has two parts. In the first part, we
assume the extinction curve is universal, completely non-varying throughout the Galaxy. We then determine the fixed shape of the extinction curve and the intrinsic colors of stars as a function of their temperatures and metallicities. In the second part, we use the previously determined intrinsic colors to determine the reddening to each star in each band. With these reddenings in hand, we are able to relax our initial assumption that the extinction curve has a single fixed shape, and study the extinction curve's variation throughout the Galaxy via a kind of principal component analysis.

This two-step approach allows us to initially fit a complicated function describing the intrinsic colors of stars as a function of temperature and metallicity, while adopting a simple description of the extinction. After we have intrinsic colors, we assume in the second step that these are known. This makes it possible to fit a more complicated model for the effect of reddening on the observed colors, where we allow the shape of the extinction curve to vary.

Our separation of the problem into two steps is not strictly optimal. To illustrate this point with an extreme example, consider a case where every star of a particular temperature and metallicity were behind an identical column of dust with a “peculiar” extinction curve. Then our initial fit would obtain incorrect intrinsic colors for stars of that temperature and metallicity. Our extinction curve variability fit would then find completely ordinary extinction for these stars, since the “peculiar” extinction would have been absorbed into the intrinsic colors and removed.

Fortunately, this situation is contrived: throughout the Galaxy we find a wide range of stellar types, and the dust column density varies on much smaller scales than the large scales over which the Galaxy’s stellar population varies. Moreover, by comparing our intrinsic colors with predictions from synthetic spectra, we verify that the intrinsic colors we obtain are reasonable.

### 3.1. Initial Fit

In our initial fit, we assume that the extinction curve is universal, described by a single, fixed reddening vector $R_0$ that is the same for all stars we consider. We then model the observed photometry $m$ as a function of the temperature and metallicity of the star, as well as the star’s distance and reddening. To be specific, we model the stars’ photometry as

$$m_i^m = f(T_i, [\text{Fe/H}]_i) + \mu_i + R_0 E'_i, \quad (1)$$

where $i$ indexes over stars, and each vector $m_i^m$ contains the $i$th star’s photometry in the 10 photometric bands. In this equation, $T_i$ and $[\text{Fe/H}]_i$ are the temperature and metallicity of the $i$th star, as determined by APOGEE. The parameters $E'_i$ and $\mu_i$ are the model distance modulus and dust column to the $i$th star. The function $f$ is an analytic function giving the intrinsic colors as a function of temperature and metallicity, with each element of $f$ corresponding to a particular band. We parameterize $f$ in each band as a fourth-order polynomial in $T$ and $[\text{Fe/H}]$, including cross terms, and additionally including terms proportional to $T^5$ and $T^6$. This means that 17 parameters are needed to describe $f$ in each band, for a total of 170 free parameters. Finally, $R_0$ is the fixed reddening vector giving the relative amount of extinction in each band, described by 10 parameters.

To completely determine the model, then, the following global parameters must be specified:

1. the intrinsic color function $f$, as determined by the 17 coefficients of the polynomial in each of the 10 bands, and

2. the reddening vector $R_0$, as determined by the 10 elements giving the relative extinctions in each band.

In addition to this, the following two parameters must be specified for each of the 37,000 stars:

1. the distance modulus $\mu_i$, and

2. the dust column $E'_i$.

The complete model then predicts the magnitudes of each star in each of the 10 bands according to Equation (1).

For a given choice of these model parameters, we evaluate $\chi$, the difference between the data and the model, scaled by the uncertainty

$$\chi_{i,j} = \frac{m_{i,j} - m_{i,j}^m}{\sigma_{i,j}}, \quad (2)$$

where $i$ indexes over the different stars, $j$ indexes over the different photometric bands, $m$ is the observed photometry, $\sigma$ is its uncertainty, and $m^m$ is given in Equation (1).

We determine the best fit model parameters by minimizing the sum of $\chi^2$. In practice, Equation (2) is vulnerable to outliers, and so we instead minimize

$$\sum_{i,j} \hat{\chi}_{i,j}^2 = \sum_{i,j} \left( \frac{\chi_{i,j}}{1 + \sqrt{|\chi_{i,j}|/5}} \right)^2, \quad (3)$$

smoothly switching from minimizing $\chi^2$ to minimizing $|\chi|$ at $\chi = 5$, philosophically similar to clipping at $5\sigma$.

The model requires a total of 74,180 parameters to be fit (ignoring, for a moment, 36 parameters which are subject to perfect degeneracies and are fixed; see Section 3.1.1). This is a substantial number, but they are constrained by roughly 300,000 photometric measurements, and so the model is well constrained.

We minimize Equation (3) to determine the best fit parameters using an alternating-least-squares algorithm. We first fix the global model parameters describing $f$ and $R_0$, and solve for the parameters $E'_i$ and $\mu_i$ for each star, one by one. We then fix these values of $E'_i$ and $\mu_i$, and solve for the global model parameters $f$ and $R_0$. With the improved global parameters, we then solve for $E'_i$ and $\mu_i$ again, iterating back and forth between solving for the per-star and global parameters until converged.

The fit has two important results. First, it gives us $f(T, [\text{Fe/H}])$, the intrinsic colors of the stars as a function of their APOGEE spectroscopic parameters. With these in hand, the observed reddenings are trivially computed as $m - f(T, [\text{Fe/H}])$, modulo a gray component of extinction and the distances to the sources. Second, the fit gives the mean extinction curve over the APOGEE footprint, as encoded in $R_0$.

### 3.1.1. Degeneracies

This model is subject to a number of perfect degeneracies which we eliminate by fixing certain parameters. It is
instructive to consider these degeneracies to see what signals we are sensitive to.

First, we are insensitive to any gray component of extinction ($R_0 \rightarrow R_0 + C$). Any gray component of extinction can be absorbed by appropriately adjusting $\mu$ for each star $i$ ($\mu_i \rightarrow \mu_i - CE_i'$). Likewise, we are insensitive to any change in the overall brightness of the stars ($f \rightarrow f + C$), since this can likewise be absorbed into $\mu$. To address these degeneracies, we fix one component of $R_0$ and fix $f$ to be 0 in the K band (18 parameters in total).

Second, we are insensitive to any change in the normalization of the extinction curve ($R_0 \rightarrow CR_0$). Any change in the normalization of the extinction curve can be compensated for by rescaling the reddenings of each star ($E_i' \rightarrow E_i'/C$). To remove this degeneracy, we fix a second component of $R_0$ (1 parameter).

Finally, we are insensitive to any change of $f$ along the reddening vector ($f(T, [Fe/H]) \rightarrow f + R_0g(T, [Fe/H])$). For such changes, we could modify the extinction $E_i' \rightarrow E_i' - g(T_i, [Fe/H]_i)$ for each star to cancel the effect. Accordingly, we can choose to set any color of choice in $f$ to theoretical expectations. We choose to set the $y_{P1}$ component of $f$ to a specific functional form picked to match expectations from synthetic spectra and the few low reddening stars in APOGEE; see Section 4.4 for further details. This fixes 17 free parameters. We note that if we had an adequate set of unreddened standard stars in APOGEE, we could avoid relying on theoretical expectations to fix these parameters.

These degeneracies mean that although $f$ is intended to represent stars' absolute magnitudes, it does not: we have no access to absolute magnitudes since we have no distances. Instead, it ultimately encodes only the intrinsic colors of stars as a function of their stellar parameters. Likewise, any additive offset to the reddening vector (that is, a gray component) is not measurable by our technique. This also means that our distance moduli $\mu$ for each star are combinations of the true distance moduli, errors in the gray component of our reddening vector, and errors in our absolute magnitudes. We do not use these "distance moduli" further in this work, and look forward to parallax measurements from Gaia, which will lift this degeneracy.

### 3.1.2. Limitations

Our model is correct in the limit that a star's intrinsic colors are polynomial functions of its temperature and metallicity, and that reddening in broad photometric bands can be described by a single, universal vector. The former assumption is expected to be valid up to the accuracy of the photometry in this work (~1%), though in principle a star's gravity, blended companions, rotation, and detailed abundances may have a small effect on the star's photometry.

The latter assumption—that reddening can be described by a single vector—is only true when the photometric bandpass is narrow and the extinction curve is universal. Neither condition applies. To mitigate the first problem, we could apply the results of Sale & Magorrian (2015). In this case, we could model the effect of reddening on magnitudes by

$$ \mathbf{m}^m = f(T, [Fe/H]) + \mu + R_0 E' + g(T, [Fe/H], \log g, E'), \quad (4) $$

where $g(T, [Fe/H], \log g, E')$ is a function describing the effect of reddening in the different bands for stars of different temperatures. The function $g$ can be fixed using existing extinction curves and synthetic stellar spectra as in Sale & Magorrian (2015). In this parameterization, the reddening vector $R_0$ would be a small perturbation to an existing extinction curve. We explored initially choosing $g$ to reproduce the extinction curve of Fitzpatrick & Massa (2009), assuming intrinsic stellar spectra given by the MARCS model grid (Gustafsson et al. 2008). However, when applying this model, we found no significant improvement in $\chi^2$. The primary effect of the more principled treatment was to scale the reddenings of highly reddened stars a few percent higher, but we are largely unconcerned here with the accuracy of the inferred monochromatic extinctions to individual stars, and so have instead applied the simpler treatment where the effect of reddening is linear and independent of the source spectrum.

We have chosen to describe $f$ as a fourth-order polynomial, plus terms proportional to $T^5$ and $T^6$. This places a limit on how well we can reproduce the intrinsic colors of stars. The choice of polynomial was driven by the desire to reproduce the intrinsic colors to 1% accuracy. We experimented with a number of different parameterizations of the intrinsic colors and examined the residuals for trends in temperature and metallicity to determine whether or not we had allowed $f$ sufficient freedom to describe the data. We note that the basic extinction curve results are largely insensitive to the order of the polynomial: the higher-order terms in the polynomial were motivated by reproducing the sharp curve of the intrinsic colors in the $g_{P1}$ band, but these cool stars compose only about 10% of the whole sample. The need for a higher order polynomial in $T$ than in $[Fe/H]$ is due to the greater dependence of the broadband photometry on temperature than on $[Fe/H]$, especially in the optical bands.

### 3.2. Extinction Curve Variation Fit

We are additionally interested in the variation of the extinction curve in the Galaxy, though we have assumed it to be universal in the previous step. Since the initial fit (Section 3.1) in concert with the APOGEE spectroscopic parameters gives us the intrinsic colors of each star, we can easily compute the reddenings of each star. Were the extinction curve universal, the observed reddenings would fall along a single line given by the reddening vector, with small dispersion due to photometric and spectroscopic uncertainties, and a small additional scatter owing to the different temperatures of the stars. Insofar as the extinction curve is in fact a single parameter family—characterized, for instance, by $R(V)$, as in Cardelli et al. (1989)—then this line will broaden into a 2D surface. Additional parameters will broaden the surface into higher-dimensional manifolds. In the limit that departures from a universal extinction curve are small, we can linearize the manifolds into linear subspaces. Accordingly, we can study the variation of the extinction curve by finding the low-dimensional subspaces that best explain the measured reddenings.

We find the best fit mean extinction curve and multi-parameter families of extinction curves by finding low-dimensional subspaces of the ten-dimensional space of observed reddenings that best explain the data. This procedure is essentially a weighted principal component analysis (PCA),
with separate weights ($\sigma^{-2}$) for each observation (Jolliffe 2002). We find these low-dimensional subspaces via the Heteroscedastic Matrix Factorization technique of Tsalmantza & Hogg (2012) (see also Gabriel & Zamir 1979; Roweis 1998; Tamuz et al. 2005). This technique, in contrast to classical PCA, appropriately accounts for the heteroscedastic uncertainty in the observations. In analogy with PCA, we call the vectors in these subspaces principal components, and order them according to the first subspace in which they appear.

The resulting principal components and the amount of variation in the data along each principal component describe the way in which the extinction curve varies and the significance of that variation.

We note that it is important to perform a heteroscedastic analysis rather than a classical PCA. Only roughly 30% of the stars we consider have photometry in all 10 bands, largely due to contaminated WISE photometry in the inner Galactic plane, but also due to saturation of unreddened stars in the PS1 bands. Moreover, the uncertainties in the PS1 bands are roughly half those in the infrared bands, leading to a significantly different weighting of the variabilities as compared with classical PCA.

To carry out the analysis quickly, we do repeated $\chi^2$ minimizations to take advantage of the bilinear nature of the problem (see Roweis 1998; Tsalmantza & Hogg 2012). However, $\chi^2$ minimization is vulnerable to outliers, so we again replace $\chi$ by a more robust version according to Equation (2). To keep the uncertainties diagonal, we perform the analysis in the 10 photometric bands. However, any gray component of extinction cannot be constrained by our technique. To address this, we force one vector in the low-dimensional subspace to be the “gray” reddening vector; the vectors in the reddening subspace of interest are constrained to be orthogonal to this vector.

We neglect the uncertainty in the model photometry stemming from uncertainty in the APOGEE temperatures and metallicities. For stars hotter than 4000 K, the uncertainties are close to aligned with the reddening vector and therefore contribute only to a small increase in noise in our estimates for the reddenings to individual stars. For colder stars, the temperature uncertainties can lead to significant dispersion in colors perpendicular to the reddening vector, which will be identified in this analysis as a reddening signal. To avoid this, we use only stars with $T > 4000$ K to determine the principal components.

## 4. RESULTS

The results of our analysis are:

1. the mean reddening vector,

2. the way it varies,

3. the reddenings of the APOGEE targets,

4. the intrinsic colors of APOGEE targets.

We perform the fit of Section 3.1 to obtain the mean reddening vector and the intrinsic colors of APOGEE targets. The model of Equation (1) proves extremely good at describing the colors of APOGEE targets in the optical through infrared. Figures 3–5 show the observed colors of APOGEE targets, their best fit model colors, and the distribution of residuals, respectively.

The model is an excellent fit to the data. As tabulated in Table 1, the scatter in the residuals in a given band is typically only marginally larger than the photometric uncertainties in that

band. Stars are observed up to $E' \approx 5$ in bands redward of $\dot{i}_{P1}$, but only up to about $E' \approx 2.5$ in the full 10 bands, as beyond this level of extinction the stars are either saturated in $y_{P1}$ or too faint in $g_{P1}$ (for these stars, $E(g - y) > 5.5$).

The description of reddening in terms of a single extinction curve is remarkably accurate in the optical and infrared. We compute extinctions in each band according to

$$ \mathbf{r} = \mathbf{m} - f(T, [\text{Fe/H}]), \qquad (5) $$

where $\mathbf{r}$ gives the observed extinctions, $\mathbf{m}$ is the observed magnitudes, and $f$ is the best fit function for the intrinsic colors of APOGEE targets as a function of their temperature $T$ and metallicity [Fe/H]. We note we have neglected $\mu$ in Equation (5): since this changes only the gray component of $\mathbf{r}$, we are insensitive to it in our analysis. Figure 6 shows the observed reddenings $\mathbf{r}$ in a variety of color combinations from the optical to infrared, with our mean reddening vector overplotted in red, and with points colored by their temperature. There is little signature of a correlation between residuals and temperature in Figure 6, suggesting that our fit to the intrinsic colors of APOGEE stars is accurate.

### 4.1. Mean Reddening Vector

The first principal component of our principal component fit (Section 3.2) is our best estimate of the mean reddening vector. It is effectively derived from the comparison of the photometry of stars of the same temperatures and metallicities, but different reddenings—i.e., via the “pair method,” with the caveat that neither star is likely to be unreddened. Table 2 shows our best fit values and their uncertainties. These are in good agreement with the literature; see Section 5.2.

We conservatively assess the uncertainty in the reddening vector $\mathbf{R}_0$ by splitting our target stars by temperature into ten equally sized subsamples and computing the root-mean-square dispersion in $\mathbf{R}_0$ over these ten subsamples. We divide the data by temperature because temperature has a more dramatic effect on the photometry than the metallicity or gravity, especially in the optical bands where the greatest variations in intrinsic color are present.

We treat the measurements in each filter here as independent. However, the measurements are covariant because of the free choice of mean and normalization for the vector, which is ultimately tied to the normalization and gray component degeneracies in our model (Section 3.1.1). It is difficult to completely determine the covariance given that we expect it to be dominated by systematics, and dividing the sample into much finer bins in temperature begins to introduce significant statistical uncertainty.

We note that the mean reddening vector we derive from our PCA (Section 3.2) and from our initial fit (Section 3.1) are in extremely good agreement, with typical differences of $<0.2\sigma$.

We also tabulate in Table 3 our measurements of the slope of the reddening vector (i.e., $E(a - b)/E(c - d)$ for different bands $a, b, c$, and $d$) in a variety of bands. These measurements have the advantage that they are independent of the normalization of the reddening vector and its gray component, though they also have non-trivial covariance because of shared photometric bands. We again determine uncertainties from the dispersion over subsamples of stars of different temperatures.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** The observed colors of APOGEE targets, colored by temperature, in the optical through infrared. In dense regions of the color–color diagrams, we have replaced individual points with colored bins giving the mean temperature of all points in that bin. The contours show the number density of points. The dashed box gives the region where the intrinsic colors lie, shown in more detail in Figure 19. The red arrow shows the reddening vector we measure. The observed colors are primarily determined by the stars’ reddening (most colors fall right along the reddening vector), though temperature also plays a role.</figcaption></figure>

## 4.2. Variation in Reddening

We assess the dimensionality of the space of reddenings by means of a version of PCA that appropriately considers the uncertainties in the measurements (Section 3.2). We show the first four principal components in Figure 7 (thick lines). Thin lines in Figure 7 show the four principal components

from different temperature subsamples of the data. Note that we force all principal components to be perpendicular to a “gray” reddening vector (see Section 3.2). Loadings for the first two principal components are tabulated in Table 2.

To assess the importance of these principal components, we compute the ratio of observed scatter in the component of
<figure><img src="image_5.png" /><figcaption>**Figure 4.** The best fit model colors of the APOGEE targets. The figure elements are the same as in Figure 3. These match the observed colors extremely well, modulo the reduced blurring due to noise. This is especially obvious in the reddest bands, where the observational signal is the smallest, since the *WISE* colors have much less sensitivity to temperature and reddening than the optical bands.</figcaption></figure>

reddening along each principal component to the expected scatter from photometric uncertainty alone. We term this ratio *W*, and tabulate it in Table 4.

In the limit that no signal is present in the data, we expect $W \approx 1$. Unsurprisingly, reddening, the first principal component, is detected at high significance ($W \approx 80$). A second principal component is detected at $W \approx 4$, which primarily acts to change the curvature of the extinction curve. We

identify this principal component with $R(V)$. The shape of the third principal component—largely, an offset between the 2MASS and *WISE* photometry—leads us to attribute it to small problems with the *WISE* photometry. The fourth and later principal components are essentially entirely consistent with noise. Therefore essentially all variation in the optical-infrared extinction curve detectable in our data is described by a single parameter.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** The residuals (data-model) in the optical through infrared, as a function of the best fit reddening $E'$ (roughly, $E(B-V)$). The solid lines show the median residual and the 16th and 84th percentiles of the distribution as a function of $E'$. The overall mean residual $\mu$ and rms scatter $\sigma$ are shown for each band in mmag, though we note that most of the sample has $0 < E' < 1$ mag. The dispersions are essentially the same as the photometric uncertainties (10–15 mmag in PS1 grizy<sub>P1</sub>, 22–24 mmag in 2MASS JHK and WISE W1 and W2), indicating that the model is an excellent fit to the data. At large $E'$, the residuals begin to increase significantly, due to a combination of photometric noise (in the optical bands) and variation in the extinction curve.</figcaption></figure>

<table><caption>**Table 1**
Typical Residuals and Photometric Uncertainties</caption><thead><tr><th>Filter</th><th>σ<sub>phot</sub>/mmag</th><th>σ<sub>resid</sub>/mmag</th></tr></thead><tbody><tr><td>g</td><td>12</td><td>11</td></tr><tr><td>r</td><td>12</td><td>10</td></tr><tr><td>i</td><td>12</td><td>17</td></tr><tr><td>z</td><td>12</td><td>16</td></tr><tr><td>y</td><td>11</td><td>13</td></tr><tr><td>J</td><td>24</td><td>23</td></tr><tr><td>H</td><td>26</td><td>23</td></tr><tr><td>K</td><td>22</td><td>25</td></tr><tr><td>W1</td><td>23</td><td>24</td></tr><tr><td>W2</td><td>22</td><td>26</td></tr></tbody></table>

**Note.** Photometric uncertainties versus root-mean-square dispersion of residuals between the data and our model. In most bands, the scatter in the residuals is dominated by the photometric uncertainties.

We use this result to parameterize the extinction to every star in terms of two principal components. We define the first principal component to be $\mathbf{R}_0$, and the second principal component to be $d\mathbf{R}/dx$, and tabulate these principal components in Table 2. We further tabulate uncertainties in $d\mathbf{R}/dx$ via the dispersion in our measurements over different temperature subsamples of the data. Then the extinctions $\mathbf{r}$ to the typical star can be expressed as

$$ \mathbf{r} = (\mathbf{R}_0 + x d\mathbf{R}/dx)CE', \quad (6) $$

where $x$ effectively determines the shape of the extinction curve toward a particular star, $E'$ is the star's extinction according to an average curve, and $C$ is a fixed overall scale factor. We choose $C$ so that $E' \approx E(B-V)$ for stars following a typical extinction curve; specifically, $C = (\mathbf{R}_0(g) - \mathbf{R}_0(r))^{-1} E(g-r)/E(B-V)$, where $E(g-r)/E(B-V) = 1.02$ is
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The measured reddenings for APOGEE targets in the optical through infrared, colored by the temperature of the stars. The figure elements are the same as in Figure 3. The reddenings are well described by a single fixed reddening vector, given by the red line. The residuals from the line are not correlated with temperature, indicating that our fit for the intrinsic colors of APOGEE stars is well determined.</figcaption></figure>

the expectation for a Fitzpatrick (1999) reddening law (taken from Schlafly & Finkbeiner 2011). Solving for $E'$ and $x$ for all targets with at least seven bands of photometry and $E' > 0.5$, we find that the root-mean-square deviation in $x$ is $\sigma(x) = 0.023$. In other words, roughly speaking, typical variations in the shape of the extinction curve lead to ~2% corrections in the reddening vector.

We note, however, that these results apply to typical stars in our sample: we do not consider the possibility that a small number of stars (say, 1%) may have significant variations in their reddenings that are not well described by a single parameter. In the context of the current work, these are hard to distinguish from cases where one of the surveys has provided spurious photometry.
<table><caption>**Table 2**
The Reddening Vector and its Variation</caption><thead><tr><th>Filter</th><th>λ</th><th>R<sub>0</sub></th><th>dR/dx</th></tr></thead><tbody><tr><td>g</td><td>5032</td><td>0.6537 ± 0.0014</td><td>-0.543 ± 0.012</td></tr><tr><td>r</td><td>6281</td><td>0.3906 ± 0.0012</td><td>0.034 ± 0.020</td></tr><tr><td>i</td><td>7572</td><td>0.2020 ± 0.0013</td><td>0.368 ± 0.017</td></tr><tr><td>z</td><td>8691</td><td>0.0787 ± 0.0013</td><td>0.423 ± 0.013</td></tr><tr><td>y</td><td>9636</td><td>-0.0048 ± 0.0017</td><td>0.382 ± 0.025</td></tr><tr><td>J</td><td>12377</td><td>-0.1421 ± 0.0028</td><td>0.141 ± 0.034</td></tr><tr><td>H</td><td>16382</td><td>-0.2366 ± 0.0012</td><td>-0.040 ± 0.025</td></tr><tr><td>K</td><td>21510</td><td>-0.2852 ± 0.0007</td><td>-0.135 ± 0.034</td></tr><tr><td>W1</td><td>32950</td><td>-0.3213 ± 0.0018</td><td>-0.269 ± 0.023</td></tr><tr><td>W2</td><td>44809</td><td>-0.3350 ± 0.0015</td><td>-0.363 ± 0.036</td></tr></tbody></table>

**Note.** The mean reddening vector $R_0$, and the wavelengths $\lambda$ for which we expect the monochromatic extinction to be nearest these values. The normalization and zero point of this vector is completely undetermined, and has been fixed by setting the mean of $R_0$ to 0 and the norm of $R_0$ to 1. We also tabulate $dR/dx$, which changes the shape of the extinction curve (Section 4.2). This vector is likewise mean 0 and norm 1, and we have fixed it to be perpendicular to $R_0$. In the language of principal component analysis, $R_0$ and $dR/dx$ are the (normalized) loadings of the first two principal components of our analysis.

<table><caption>**Table 3**
Reddening Vector Slopes and their Variation</caption><thead><tr><th>Filters</th><th>Slope</th><th>dSlope/dx</th></tr></thead><tbody><tr><td>g,r,i</td><td>1.395 ± 0.014</td><td>-0.59 ± 0.26</td></tr><tr><td>r,i,z</td><td>1.531 ± 0.013</td><td>-2.04 ± 0.36</td></tr><tr><td>i,z,y</td><td>1.477 ± 0.036</td><td>-1.36 ± 0.38</td></tr><tr><td>z,y,J</td><td>0.608 ± 0.010</td><td>-0.77 ± 0.29</td></tr><tr><td>y,J,H</td><td>1.454 ± 0.042</td><td>-0.25 ± 0.76</td></tr><tr><td>J,H,K</td><td>1.943 ± 0.020</td><td>-0.03 ± 0.48</td></tr><tr><td>H,K, W1</td><td>1.348 ± 0.042</td><td>-2.40 ± 2.17</td></tr><tr><td>K, W1, W2</td><td>2.627 ± 0.197</td><td>-8.18 ± 7.70</td></tr></tbody></table>

**Note.** The slope of the reddening vector in different bands, and the variation in that slope $dSlope/dx$. The filter combination a,b,c corresponds to the slope $E(a-b)/E(b-c)$. Uncertainties are from the root-mean-square dispersion over subsets of stars with different temperatures. These data are equivalent to those in Table 2, except for somewhat different uncertainties due to the covariance of the measurements.

Variation in the extinction curve is typically assessed via $R(V) = A(V)/E(B-V)$. We cannot however directly measure $A(V)$ in this work, because we are insensitive to any gray component of the extinction curve. However, a simple approximate proxy for $R(V)$ can be constructed from $(A(g) - A(W2) / A(g) - A(r))$ (see Section 5.3). Our mean reddening vector corresponds to a proxy $R(V)$ of 3.33, and changing $x$ changes $R(V)$. This is marginally different from the “standard” value of 3.1, but given the uncertainties on the standard value and the uncertainty in translating our proxy to $R(V)$, we do not see this as problematic.

To test the robustness of these results, we have compared our heteroscedastic PCA-like analysis of Section 3.2 with a traditional unweighted PCA analysis using the roughly 30% of stars for which photometry is available in all ten bands. We obtain qualitatively similar principal components for each of the first three principal components. Quantitatively, the second, $R(V)$-like principal component is different from our preferred principal component by up to $3\sigma$ in particular bands. We take this as remarkably good agreement given the significantly

<figure><img src="image_8.png" /></figure>

**Figure 7.** The first four principal components of the extinction curve (ignoring gray). The legend gives the amount of variation in the reddenings in the direction of each principal component, relative to expectations from photometric uncertainties alone. The top labels indicate the filter corresponding to each point. The first principal component is essentially the mean reddening vector. The second principal component is very similar to the effect of $R(V)$ in other formulations of the extinction curve; increasing $R(V)$ reduces the curvature of the extinction curve. The later principal components have $W \approx 1$ in this data: they are essentially not necessary to describe the observed spectrum of a single star, though formally over the whole APOGEE sample they are statistically significant. The thick lines show the best fit principal components from the full data set, while the thin lines show the principal components as determined from ten independent temperature subsamples.

<table><caption>**Table 4**
Principal Component Significance</caption><thead><tr><th>PC #</th><th>W</th><th>Note</th></tr></thead><tbody><tr><td>1</td><td>80.5</td><td>mean reddening vector (see Section 5.2)</td></tr><tr><td>2</td><td>3.5</td><td>R(V) (see Section 5.1)</td></tr><tr><td>3</td><td>1.2</td><td>mostly consistent with noise</td></tr><tr><td>4</td><td>1.0</td><td>mostly consistent with noise</td></tr></tbody></table>

**Note.** The significance of the principal components: the observed scatter in the reddenings along a particular component relative to the expected scatter from photometric uncertainty alone, $W$. We find that reddening (the first principal component) is detected, unsurprisingly, at extremely high significance in this data ($W \approx 80$). A second principal component, similar to $R(V)$ in other extinction curve prescriptions, has $W = 3.5$. The later principal components are detected at a significance of about $W = 1$: that is, the observed scatter in reddenings in these directions is almost completely consistent with photometric noise.

different weighting applied by our more correct analysis, where the optical photometry is weighted four times as heavily as the infrared data, due to its four times smaller variance.

### 4.3. Reddenings of APOGEE Stars

The model provides estimates of the reddening to each of the APOGEE targets, allowing the reddening to be mapped over the footprint, as shown in Figure 1. Observed reddenings range from $0 \le E' \le 5$ mag, though requiring a $g_{P1}$ band detection limits the range to roughly 2.5 mag, or an $r_{P1}$ band detection to roughly 3.5 mag.

The uncertainty of our reddening estimates to individual stars is constrained by the photometry to be less than 10 mmag.
The dominant source of noise, however, is in the APOGEE temperatures. Because the reddening vector and stellar locus are nearly covariant, uncertainty in temperature translates almost entirely into uncertainty in reddening (and hardly at all to uncertainty in $R(V)$ or any directions perpendicular to the reddening vector), except for stars with $T < 4000$ K. The temperature uncertainty translates into an $E'$ uncertainty of about 30 mmag.

### 4.4. Intrinsic Optical–Infrared Colors of Giants

Our analysis also gives the intrinsic optical–infrared colors of giant and subgiant stars as a function of their APOGEE parameters. Figure 8 shows color–color diagrams of the intrinsic colors our model predicts for each of the APOGEE stars. The colors look as expected: in $gri_{P1}$, the stellar locus features a prominent bend at $g_{P1} - r_{P1} \approx 1.2$; the other optical bands are essentially linear and relatively one-dimensional (little dependence of color on metallicity). In the infrared, metallicity begins to play a larger role relative to temperature than in the optical. The near collinearity of the reddening vector and the stellar locus in the optical largely disappears in the WISE bands, as recognized by Majewski et al. (2011), which motivated the APOGEE dereddening and target selection algorithm of Zasowski et al. (2013).

As discussed in Section 3.1.1, our intrinsic color measurements are subject to a perfect degeneracy. Since we do not assume we know where zero reddening lies, the entire set of colors can be shifted along the reddening vector without having any effect on the goodness of fit. Worse, the analysis essentially works by comparing the reddening of stars with similar temperatures and metallicities, and so in fact stars of different temperatures and metallicities can be shifted by different amounts along the reddening vector without affecting the analysis.

To address this shortcoming, we fix the intrinsic $y_{P1} - K$ color to be a function of our choice, intended to be the true $y_{P1} - K$ color of unreddened stars. We choose to force our intrinsic $y_{P1} - K$ colors to agree with the predictions of the MARCS model grid of Gustafsson et al. (2008), which were found to be good predictors of broadband colors by Edvardsson (2008) and Casagrande & VandenBerg (2014). It is challenging to fully assess the accuracy of these predictions with the APOGEE stars themselves. The most straightforward approach is to compare the $y_{P1} - K$ color of APOGEE targets in regions of low reddening with the predictions from model spectra. The results of such a comparison are shown in Figure 9.

The top panel of Figure 9 shows the predicted $y_{P1} - K$ color of APOGEE stars as a function of their temperature, with points colored by their metallicity. In the models, the $y_{P1} - K$ color has very little dependence on metallicity. The solid line shows our fit to these synthetic colors, at [Fe/H] = 0. The fit is extremely good; the root-mean-square difference between our fit and the actual prediction is 1.3 mmag. The observed colors, dereddened according to Schlegel et al. (1998, SFD), are shown in the bottom panel, using only stars for which $E(B-V)_{SFD} < 0.2$. We have again shown the [Fe/H] = 0 fit from the top panel, to illustrate the good agreement in general between the synthetic and observed magnitudes. Still, the fit is not perfect: there is a noticeable trend in $y_{P1} - K$ color with metallicity in the observed colors, which is absent or much reduced in the synthetic colors. Moreover, below 4000 K there are very few unreddened APOGEE targets, and none of these

are solar metallicity or above, making it impossible to test the accuracy of the fit in this region. Still, the fit is at least consistent with the limited data available there.

We can dramatically increase the number of stars available for this test by adopting a more permissive cut on $E(B-V)_{SFD}$, essentially using stars closer to the plane that are more likely to be distant solar metallicity giants. However, SFD is known to be problematic at low latitudes, especially in the inner Galaxy (e.g., Schlafly et al. 2014b). Extending Figure 9 to $E(B-V)_{SFD} < 0.5$ mag renders the observed metallicity trend invisible, presumably due to systematic overprediction of the reddening in the plane in SFD, possibly due to the fact that the APOGEE targets may not be behind the entire dust column.

Due to the degeneracies, our choice of $y_{P1} - K$ color does not affect the reddening vector we derive, or its variation. Accordingly, we simply fix $y_{P1} - K$ to the value we fit from the MARCS stellar models. Errors in this choice lead only to small variations in the $E'$ we infer to individual stars.

### 4.5. Quality of the Fit

We find in Section 4.2 that the extinction curve to our stars can be well parameterized in terms of two principal components. We can assess the quality of the fit by modeling the observed extinctions $m - f(T, [Fe/H])$ as a sum of a gray component, the mean reddening vector $R_0$, and $dR/dx$ (i.e., an $R(V)$-like component). We then compute the $\chi^2$ per degree of freedom for each star. We show the spatial distribution of $\chi^2$ per degree of freedom in Figure 10.

The mean $\chi^2$ per degree of freedom in Figure 10 is 0.77, slightly smaller than the expected value of 1, suggesting our uncertainties are slightly overestimated. As Figure 10 makes clear, however, the inner Galaxy and especially the Galactic bulge feature significantly larger $\chi^2$ per degree of freedom than the entire outer Galaxy, which is relatively featureless.

There are a few reasons for the large $\chi^2$ per degree of freedom in the inner Galaxy. The first is that we have included in the computation of $\chi^2$ only the contribution from the photometric uncertainties, and neglected the contribution of the spectroscopic uncertainties in *T* and [Fe/H]. For most stars, the spectroscopic uncertainties affect $\chi^2$ negligibly, because the reddening vector and stellar locus are well aligned, and any error in *T* can therefore be absorbed into the extinction in the model. For the coldest stars ($T \lesssim 4000$ K), however, the reddening vector becomes nearly perpendicular to the stellar locus, and colors become more sensitive to temperature. This can lead to uncertainties in *T* dominating the total uncertainties for these stars. In principle we could build the temperature uncertainty into the analysis, but we have found that doing so makes no difference for the reddening vector we derive, and so we have neglected this contribution. Because essentially all of the $T < 4000$ K stars reside in the inner Galaxy, this leads to elevated $\chi^2$ there.

However, even when excluding all cold stars we still find elevated $\chi^2$ per degree of freedom in the inner Galaxy. The primary driver seems to be residuals in the WISE bands, though additionally excluding these bands does not fully resolve the problem.

We have not been able to fully understand the cause of the large $\chi^2$ per degree of freedom in the inner Galaxy. One explanation could be that the extinction curve is more variable or differently variable in the inner Galaxy (e.g.,
<figure><img src="image_9.png" /><figcaption>**Figure 8.** The intrinsic colors of APOGEE stars, according to the best fit model of this work. The figure elements are as in Figure 3. The colors look much as expected: the $g_{P1} - r_{P1}$, $r_{P1} - i_{P1}$ diagram (top left), for instance, shows the expected sharp bend at $g_{P1} - r_{P1} = 1.2$ mag typical of PS1 filters. The colors in all filter combinations are determined primarily by temperature, though the distribution of points is also significantly broadened by colors' dependence on metallicity; iso-color lines show the variation in color due to varying metallicity at fixed temperature. The $r_{P1} - i_{P1}$, $i_{P1} - z_{P1}$ color-color diagram is found to be especially close to being a single-parameter family; meanwhile in the infrared, metallicity dramatically broadens the stellar locus.</figcaption></figure>

Nataf et al. 2013, 2015). However, applying our analysis only to stars with $|l| < 60^\circ$ and $|b| < 5^\circ$ produces similar principal components as when using our full data set. In particular, the third principal component continues to bear signs of a WISE-2MASS offset, which we find unlikely to have a physical origin.

We conclude that over the great majority of the sky, our two-component extinction curve model provides an excellent description of the data. In the inner Galaxy, this model clearly fails to account for the data completely. Nevertheless, even in this region, the two-component model describes the majority of the variation in the extinction curve.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The intrinsic $y_{p1} - K$ color of APOGEE stars, as synthesized from MARCS model spectra (top) and as observed for $E(B-V)_{SFD} < 0.2$ mag stars after dereddening according to SFD. The results of our fit to the MARCS predictions, for $[Fe/H] = 0$, is shown as the solid line in each panel. There is good overall agreement between the MARCS predictions and the observations, though the MARCS colors predict less variation in color with metallicity than is observed. Points are colored by metallicity, and a fit to the MARCS colors at $[Fe/H] = 0$ is shown in each panel as a solid line. The contours show the density of stars.</figcaption></figure>

# 5. DISCUSSION

In this section, we discuss the implications of our results and compare with similar measurements from the literature. We

1. discuss our one-parameter family of extinction curves,

2. compare our extinction curve with the literature,

3. link our observed variation with $R(V)$,

4. study the variation of the extinction curve through the Galaxy,

5. study correlations between the dust emission and extinction curves,

6. discuss measurements of the “gray” component of the extinction, and

7. compare our intrinsic colors with models.

## 5.1. *The Extinction Curve*

We find in Section 4.2 that the reddening to a star in the optical and infrared can be parameterized by two numbers:

essentially, the amount of reddening to the object, and a second parameter that slightly alters the direction of the reddening vector. This is the same conclusion that CCM made, and popularized the idea that the extinction curve can be parameterized by $R(V) = A(V)/E(B-V)$.

This work distinguishes itself from CCM in that that work focused largely on the UV, while we are concerned entirely with the optical-infrared extinction curve. Moreover, the UV extinction curve is in detail not a single-parameter family—the work of Fitzpatrick & Massa (1990) uses six parameters to describe the UV extinction curve, and the CCM curve was found to describe the majority of the variability among those parameters. However, we find that the bulk of the APOGEE data provide no support for using more than a single-parameter family to describe the extinction in the optical through infrared, with the possible exception of the inner Galaxy, where we can draw no firm conclusions (Section 4.5).

We find that the curvature of the extinction curve increases with decreasing $R(V)$ throughout the optical and infrared. This can be directly seen in Figure 11, where we show observed reddenings of stars, colored by $R'(V)$ as estimated from our PCA, for stars with photometry in at least 9 bands. In the optical bands, there is a clear association between the slope of the reddening vector and $R'(V)$. Meanwhile redward of J, the trend becomes harder to detect. However, we note that the work of Zasowski et al. (2009), which had access to stars of significantly larger reddening, found that the longest wavelengths had the most significant variations.

## 5.2. *The Extinction Curve Compared with Past Measurements*

In this subsection, we compare our mean extinction curve with past measurements of the optical-infrared extinction curve from broad-band photometry, and with standard parameterizations of the extinction curve. Our measurements rely on a large sample of highly reddened stars with precise, homogeneous stellar parameters and photometry, covering much of the Galactic plane. We therefore believe our measurements of the extinction curve supersede earlier works.

### 5.2.1. *Comparison with Photometric Measurements*

The slope of the reddening vector has been measured by a large number of authors in a variety of photometric bands in the optical through infrared. We show in Table 5 our $E(a-b)/E(c-d)$ in different combinations of bands, as compared with a selection of measurements from the literature.

The bands we study are not always a perfect match for the bands considered in other works. To translate measurements from one set of bands into another, we use an F99 extinction curve to compute predictions for the expected ratios $A(a)/A(b)$. Because the slope of the reddening vector also depends slightly on the intrinsic spectrum of a star and its reddening, we further additionally slightly correct the measurements of Schlafly & Finkbeiner (2011) and Yuan et al. (2013) to account for this difference. To make this adjustment, we assume in Table 5, we assume that these two sets of measurements studied lightly reddened ($E(B-V) \approx 0.1$)) blue main sequence stars (7000 K), while we study reddened ($E(B-V) \approx 0.65$)) red giants (4500 K). However, the choice of spectrum usually makes a difference of <0.01 in $E(a-b)/E(c-d)$.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** χ² per degree of freedom for APOGEE targets. Throughout most of the disk, we obtain a χ²/dof of 0.77, suggesting our uncertainties are very slightly overestimated. In the inner disk and especially toward the Galactic center, χ²/dof becomes significantly larger and more variable. Much of the contribution in the Galactic center comes from the WISE W1 and W2 bands, suggesting that crowding and blending may be problematic. This may also be a signal that the dust extinction curve is more variable in these directions, though we do not find conclusive signatures of that possibility.</figcaption></figure>

The agreement between our measurements and the literature is good. We agree with Schlafly & Finkbeiner (2011) in $E(g_{P1} - r_{P1})/E(r_{P1} - i_{P1})$ and $E(r_{P1} - i_{P1})/E(i_{P1} - z_{P1})$ to within about $1\sigma$. Our agreement with Yuan et al. (2013) is acceptable, though there is a more than $4\sigma$ disagreement in $E(r - i)/E(i - z)$ and in $E(J - H)/E(H - K)$. We also measure $E(K - W1)/E(W1 - W2)$ to be nearly twice as large as Yuan et al. (2013), though their measurement was uncertain due to the low reddenings available in their work. We are about $2\sigma$ discrepant from the work of Nataf et al. (2013), though we note that the spread in $E(V - I)/E(J - K)$ reported there is attributed to real variation in the extinction curve rather than uncertainty, making a $2\sigma$ offset very significant. We find in general poor agreement with the results of Davenport et al. (2014), possibly due to a systematic effect arising from adopting an intrinsic color relation best suited to dwarfs to find the reddenings of likely giants.

The work of Wang & Jiang (2014, W14) makes a measurement of $E(J - H)/E(H - K)$ from similar data to ours: APOGEE spectroscopy combined with 2MASS data. Our values, however, are in conflict: we find $E(J - H)/E(H - K) = 1.943$ as compared with their 1.780, a huge discrepancy given our estimated uncertainty of 0.019 and their quoted uncertainty of 0.008. This difference is especially disturbing given the largely identical data adopted by our work and that of W14. To determine the source of the discrepancy, we repeated our analysis, this time including only 2MASS data and excluding PS1 and WISE data. We could reproduce their results only if we also adopted their equations for the intrinsic colors of giant stars as a function of their temperature (Equations (1)–(3) in W14). This highlights the importance of the intrinsic colors to this analysis: because the reddest and most metal-rich APOGEE stars tend to lie at the greatest distances toward the Galactic center, and are accordingly the most reddened, a systematic trend in intrinsic color with temperature can masquerade as a different reddening vector. In our analysis, we allow the intrinsic colors to be fit simultaneously with the reddening vector, and obtain a fit with significantly better $\chi^2$ than we obtain with the W14 intrinsic color relation. We conclude that our value is much more likely to be correct, and note that it is in good agreement with the work of Zasowski et al. (2009).

### 5.2.2. Comparison with Existing Extinction Curve Parameterizations

We can additionally compare our results with parameterized extinction curves from the literature: we consider the curves of CCM, Fitzpatrick (1999, 2004, hereafter F99, FM04 respectively), Fitzpatrick & Massa (2009, FM09), and Maíz Apellániz et al. (2014, M14). The CCM and M14 extinction curves are not defined redward of 33333 Å; we extend these curves redward using the FM09 curve. Similarly, the FM09 extinction curve was developed using data redward of 6000 Å; we extend it blueward of 6000 Å using F99.

We note that each of these extinction curves is a one-parameter family. In the case of all but the FM09 extinction curve, the controlling parameter is referred to as $R(V)$; in the case of FM09, it is referred to as $\alpha$. FM09 has a second parameter, $R(V)$, that does not change the shape of the extinction curve redward of $V$; we ignore it here.

We compute predictions of the slope of the reddening vector we should observe from extinction curves by integrating the appropriate filter bandpasses over the MARCS synthetic spectrum of a 4500 K star with $\log g = 2.5$ and solar metallicity, typical of the APOGEE sample. Since the APOGEE sample has few unextinguished stars, the slope we observe is actually the slope of the reddening vector at $E(B - V) \approx 0.65$ rather than at $E(B - V) \approx 0$. To account for this, we compute $dm_b/dA$ for small variations of the extinction $dA$ about $E(B - V) = 0.65$, where $m_b$ is the observed magnitude in the bandpass $b$.

This procedure generates $dm_b/dA$ in each bandpass $b$, but in APOGEE we are insensitive to the normalization of the extinction and to any gray component of the extinction. So to compare with Table 2, we fit extinction curves to the measurements as

$$R_{0b,obs} = C \cdot dm_b/dA + D, \qquad (7)$$

in a least squares sense, using the uncertainties given in Table 2 and neglecting any covariance in the uncertainties.

Performing the fit to each family of extinction curves, the FM09 extinction curve is by far most consistent with our measurements ($\chi^2 = 25.5$). The F99 extinction curve is next best ($\chi^2 = 92.7$), followed by the FM04 extinction curve ($\chi^2 = 202.5$). The CCM and M14 extinction curves are strongly disfavored ($\chi^2 = 633.1$ and $\chi^2 = 861.0$). These $\chi^2$
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Reddenings of APOGEE stars in different color combinations, colored by our principal-component inferred $R'(V)$. The figure elements are as in Figure 3. The clearest changes in slope with $R'(V)$ are in the bluer bands ($g_{P1}$ through $z_{P1}$), while in the infrared the slope variation is relatively small.</figcaption></figure>

are all on 7 degrees of freedom, except for FM09, which is on 6 degrees of freedom, since we have excluded the $g_{P1}$ band, which is outside the region where the FM09 prescription was developed.

This simple description of which extinction curve works best hides many important details of the extinction curves. The M14 extinction curve, for instance, provides an extremely good fit at $R(V) = 3.7$ in the PS1 bands ($\chi^2 = 1.0$ on 2 degrees of

freedom). However, the IR extension of M14 follows CCM, which poorly matches our measurements, leading to high $\chi^2$ overall. We compare the behavior of the various extinction curves with that which we find from APOGEE in Figure 12, which shows the slope of the reddening vector in different photometric bands for a variety of extinction curves families.

The bottom panel of Figure 12 shows that we observe variation in the slope of the reddening vector across the optical
<table><caption>**Table 5**
Mean Reddening Vector Compared to Past Measurements</caption><thead><tr><th>Filters</th><th>This Work</th><th>Lit. Conv.</th><th>Orig. Filters</th><th>Orig. Value</th><th>References</th></tr></thead><tbody><tr><td>E(gP - rP)/E(rP - iP)</td><td>1.395 ± 0.013</td><td>1.449 ± 0.049</td><td>E(gS - rS)/E(rS - iS)</td><td>1.695 ± 0.057</td><td>Schlafly+2011</td></tr><tr><td>E(rP - iP)/E(iP - zP)</td><td>1.531 ± 0.012</td><td>1.578 ± 0.041</td><td>E(rS - iS)/E(iS - zS)</td><td>1.299 ± 0.034</td><td>Schlafly+2011</td></tr><tr><td>E(gP - rP)/E(rP - iP)</td><td>1.395 ± 0.013</td><td>1.410 ± 0.038</td><td>E(gS - rS)/E(rS - iS)</td><td>1.650 ± 0.044</td><td>Yuan+2013</td></tr><tr><td>E(rP - iP)/E(iP - zP)</td><td>1.531 ± 0.012</td><td>1.695 ± 0.041</td><td>E(rS - iS)/E(iS - zS)</td><td>1.395 ± 0.034</td><td>Yuan+2013</td></tr><tr><td>E(iP - zP)/E(zP - J2)</td><td>0.558 ± 0.008</td><td>0.558 ± 0.015</td><td>E(iS - zS)/E(zS - J2)</td><td>0.768 ± 0.020</td><td>Yuan+2013</td></tr><tr><td>E(zP - J2)/E(J2 - H2)</td><td>2.338 ± 0.076</td><td>2.434 ± 0.105</td><td>E(zS - J2)/E(J2 - H2)</td><td>2.154 ± 0.093</td><td>Yuan+2013</td></tr><tr><td>E(J2 - H2)/E(H2 - K2)</td><td>1.943 ± 0.019</td><td>1.627 ± 0.063</td><td>...</td><td>1.625 ± 0.063</td><td>Yuan+2013</td></tr><tr><td>E(H2 - K2)/E(K2 - W1)</td><td>1.348 ± 0.040</td><td>1.318 ± 0.093</td><td>...</td><td>1.333 ± 0.094</td><td>Yuan+2013</td></tr><tr><td>E(K2 - W1)/E(W1 - W2)</td><td>2.627 ± 0.187</td><td>4.653 ± 1.440</td><td>...</td><td>4.615 ± 1.428</td><td>Yuan+2013</td></tr><tr><td>E(gP - iP)/E(J2 - K2)</td><td>3.157 ± 0.066</td><td>3.676 ± 0.214</td><td>E(VL - IL)/E(J2 - K2)</td><td>2.913 ± 0.170</td><td>Nataf+2013</td></tr><tr><td>E(gP - rP)/E(rP - iP)</td><td>1.395 ± 0.013</td><td>1.443</td><td>E(gS - rS)/E(rS - iS)</td><td>1.625</td><td>Davenport+2014</td></tr><tr><td>E(rP - iP)/E(iP - zP)</td><td>1.531 ± 0.012</td><td>1.283</td><td>E(rS - iS)/E(iS - zS)</td><td>1.043</td><td>Davenport+2014</td></tr><tr><td>E(iP - zP)/E(zP - J2)</td><td>0.558 ± 0.008</td><td>0.733</td><td>E(iS - zS)/E(zS - J2)</td><td>1.000</td><td>Davenport+2014</td></tr><tr><td>E(zP - J2)/E(J2 - H2)</td><td>2.338 ± 0.076</td><td>2.852</td><td>E(zS - J2)/E(J2 - H2)</td><td>2.556</td><td>Davenport+2014</td></tr><tr><td>E(J2 - H2)/E(H2 - K2)</td><td>1.943 ± 0.019</td><td>1.500</td><td>...</td><td>1.500</td><td>Davenport+2014</td></tr><tr><td>E(H2 - K2)/E(K2 - W1)</td><td>1.348 ± 0.040</td><td>1.000</td><td>...</td><td>1.000</td><td>Davenport+2014</td></tr><tr><td>E(K2 - W1)/E(W1 - W2)</td><td>2.627 ± 0.187</td><td>1.500</td><td>...</td><td>1.500</td><td>Davenport+2014</td></tr><tr><td>E(J2 - H2)/E(H2 - K2)</td><td>1.943 ± 0.019</td><td>2.000 ± 0.050</td><td>...</td><td>2.000 ± 0.050</td><td>Zasowski+2009</td></tr><tr><td>E(J2 - H2)/E(H2 - K2)</td><td>1.943 ± 0.019</td><td>1.780 ± 0.008</td><td>...</td><td>1.780 ± 0.008</td><td>Wang+2014</td></tr><tr><td>E(J2 - H2)/E(H2 - K2)</td><td>1.943 ± 0.019</td><td>1.778 ± 0.154</td><td>...</td><td>1.778 ± 0.154</td><td>Indebetouw+2005</td></tr></tbody></table>

**Note.** The extinction curve of this work compared with the literature. The first column gives the reddening vector slope of interest, and the second column gives our measurement of it. The third column gives measurements and uncertainties from the literature, converted to be in the same filters as our measurements, when applicable. The fourth column gives the original filters the literature measurement was made in, if different from ours, and the fifth column gives the original values of the measurement from the literature. Finally, the sixth column gives the literature reference. Bandpasses from the Pan-STARRS1, SDSS, 2MASS, and Landolt systems are denoted with the characters P, S, 2, and L, respectively. The measurements of Schlafly & Finkbeiner (2011) and Yuan et al. (2013) have been further adjusted to account for the fact that these works studied intrinsically hotter, less reddened stars than the stars targeted in APOGEE; see text.

and infrared, with significantly reduced variation in the NIR relative to the optical. Associating this variation with $R(V)$ according to Equation (8), we find that the slope of the reddening vector increases with decreasing $R(V)$ throughout the optical and infrared. The general trend of the variation is qualitatively similar to that of all curves except the FM09 and F99 curves: the optical bands show more variability than the infrared bands. However, in detail none of the extinction curves come especially close to our observations.

Restricting to bands blueward of J, the M14 family provides both a good description of the mean extinction curve and its variation with $R(V)$. It however qualitatively disagrees with our measurement for the amount of variation in $E(g - r)/E(r - i)$ relative to $E(r - i)/E(i - z)$: we find more significant variation in the latter than the former.

We summarize these results in Figure 13. The upper panel compares the predictions for the slope of the reddening vector in different filter combinations for various extinction curves (colors) as compared with our measurements (thick gray line). The bottom panel shows the derivative of the slope with respect to $R(V)$ (or $\alpha$ in the case of FM09, with an arbitrary rescaling) for different filter combinations. A sense for the uncertainty in our measurements is given by the thin gray lines, which show our measurements on 10 different temperature subsamples of the APOGEE data. The top panel shows that several extinction curves provide acceptable descriptions of the reddening in the optical (M14, F99, FM04, FM09), but only FM09 also provides an acceptable description in the infrared.

We note that Figure 12 highlights the disagreement between extinction curves as much as possible. The slopes of the reddening vectors in Figure 12 are related to second derivatives of the extinction curve. An alternative would be to plot

$E(\lambda - r)/E(g - r)$. In such plots, all of the various extinction curves have similar behavior, since all curves predict that $E(g - \text{infrared})/E(g - r)$ increase with $R(V)$, and all curves look similar in the mean. So the mixed agreement of Figure 12 is possible despite the fact that ultimately all of these curves share broad similarities in their description of extinction.

For completeness, we also consider the O'Donnell (1994) extinction curve. This curve is identical to the CCM extinction curve in the infrared, but refines CCM in the optical. We find that the curve fits our data about as poorly as CCM and M14 do ($\chi^2 = 611$ on 7 degrees of freedom), due to the mismatch in the infrared. Restricting to the PS1 bands improves the fit only somewhat, obtaining $\chi^2 = 143$ on 2 degrees of freedom.

### 5.3. Linking Our Extinction Curve with $R(V)$

It is useful to link our description of the extinction curve with typical descriptions in terms of $R(V) = A(V)/E(B - V)$. We cannot directly measure $R(V)$ because it depends on the gray component of the extinction. However, since $B - V$ is similar to $g_{P1} - r_{P1}$ in that they are both optical colors covering similar wavelengths, and since $A(V)$ is similar to $E(g_{P1} - W2)$, in that $A(W2) \ll A(g_{P1})$, we are motivated to look for a linear relationship between $R(V)$ and $E(g_{P1} - W2)/E(g_{P1} - r_{P1})$. Figure 14 shows our calculations for this quantity as a function of $R(V)$ for a variety of extinction curves at $A(V) = 2$. In all cases, the relationship is not far from linear, though the F99 and FM04 predictions for $E(g_{P1} - W2)/E(g_{P1} - r_{P1})$ differ from the CCM predictions by as much as 0.5. All extinction curves predict similar slopes, so we fit a line to the F99 predictions to
<figure><img src="image_13.png" /><figcaption>Figure 12. Different parameterizations of the extinction curve compared with our measurements from APOGEE. We consider the CCM, M14, F99, FM04, and FM09 extinction curves. The y-axis shows the slope of the reddening vector in a particular combination of filters, as indicated on the x-axis. The filter combination g, r, i corresponds to the slope $E(g - r)/E(r - i)$, and analogously for the other filter combinations. Different colors correspond to different values of $R(V)$ or $\alpha$. The blue line with crosses shows the APOGEE measurements. The bottom panel shows our APOGEE-determined extinction curve. The FM09 extinction curve best matches our mean extinction curve, though its prediction for the variation in the extinction curve is significantly different from what we measure.</figcaption></figure>

obtain

$$R'(V) = 1.2E(g_{P1} - W2)/E(g_{P1} - r_{P1}) - 1.18, \quad (8)$$

and use this as a proxy for $R(V)$. In terms of our description of the reddening vector in Equation (6), this roughly corresponds to $R(V) = 3.3 + 9.1x$.

<figure><img src="image_14.png" /><figcaption>Figure 13. Comparison of various extinction curves and our measurements from APOGEE. The top panel shows predictions for the slope of the extinction curve in different filter combinations (see Figure 12 caption for details), and the bottom panel shows the derivative of the slope with respect to $R(V)$ (or $\alpha$ in the case of FM09, with an arbitrary rescaling). The black line with error bars shows our APOGEE measurements, while the thin gray lines show our measurements for subsamples of the data with different temperatures. The FM09 extinction curve provides the best match, but quantitatively no extinction curve agrees especially well with our measurements of the variation of the shape of the extinction curve.</figcaption></figure>

### 5.4. R(V) Variation in the Galaxy

We have identified a single $R(V)$-like parameter which describes the shape of the extinction curve in the optical and UV. In this subsection, we make measurements of this parameter for every star in our sample, and use them to study the extent to which $R(V)$ varies in the Galaxy.

We use the $R(V)$ proxy $R'(V)$ of Equation (8) to determine the typical extent to which $R(V)$ varies in the Galactic plane. Figure 15 shows the distribution of $R'(V)$ for our sample, for stars with $g_{P1}$ photometric uncertainty of less than 0.1 mag, $E' > 0.5$, and $r_{P1}$ and W2 photometry. The distribution is well described by a Gaussian with a mean of 3.32 and a standard deviation of 0.18. This is a somewhat tighter distribution of $R(V)$ than found by Fitzpatrick & Massa (2007) ($\sigma = 0.27$), and the same size as inferred by Schlafly et al. (2010). For $E' > 0.5$, the typical uncertainty in $R'(V)$ is less than 0.1; the observed scatter is dominated by the intrinsic width of the $R'(V)$ distribution. Moreover, the tail to large $R'(V)$ is much less pronounced in Figure 15 than in the work of Fitzpatrick & Massa (2007), possibly owing to the very different populations of stars probed in the two works (O stars versus background
<figure><img src="image_15.png" /><figcaption>**Figure 14.** The ratio $E(g_{P1} - W2)/E(g_{P1} - r_{P1})$ for a variety of extinction curves at $A(V) = 2$, as a function of $R(V)$. The ratio is very close to linear in $R(V)$, motivating us to adopt the proxy $R(V) \approx 1.2E(g_{P1} - W2)/E(g_{P1} - r_{P1}) - 1.18$ as a simple proxy for $R(V)$, based on a fit to the Fitzpatrick (1999) extinction curve.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Distribution of $R'(V) \approx R(V)$ according to the equation $R'(V) = 1.2E(g_{P1} - W2)/E(g_{P1} - r_{P1}) - 1.18$ (Equation (8)). The width of the distribution of $R'(V)$ among these stars is remarkably small (0.18). Moreover, notably absent is a significant tail in this distribution toward large $R'(V)$: only 0.8% of the sample has $R'(V) > 4$, in comparison with 9.5% in the sample of Fitzpatrick & Massa (2007).</figcaption></figure>

giants). We find, for instance, that only 0.8% of our stars have $R'(V) > 4$, in comparison with 9.5% in the work of Fitzpatrick & Massa (2007). Moreover, many of these stars are simply poor fits, presumably owing to spurious photometry from one of the surveys; the true fraction of high $R(V)$ sight lines is presumably still smaller. Our results should much better describe the variation in $R(V)$ expected along a typical line of sight.

Dust properties are expected to change in different environments: a dust grain in a diffuse, atomic cloud is subject to a very different radiation field and collision rate than a dust grain in a dense molecular cloud. For example, Whittet et al. (1988) found signatures of water-ice in infrared spectra of stars in Taurus, along sight lines with $E(B-V) \gtrsim 1$, but sight lines with lower $E(B-V)$ were free of ice. The work of Ysard et al.

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Distribution of $R'(V)$ vs. $E'$, which are roughly equivalent to $R(V)$ and $E(B-V)$. The grayscale shows the density of points, and the solid lines show the 16th, 50th, and 84th percentiles of the distribution at each $E'$. Below $E' = 0.5$ the uncertainty begins to increase rapidly, but outside this region $R'(V)$ is remarkably independent of $E'$, described by a simple Gaussian with a standard deviation of 0.2. There is a slight tendency at the highest $E'$ where we have $g_{P1}$ measurements ($E(B-V) = 2.5$) for $R'(V)$ to be slightly higher than average (by less than 0.1); this may however be due to a systematic bias, in that there is less g band extinction at higher $R'(V)$ and fixed $E'$.</figcaption></figure>

(2013) finds signs of grain agglomeration at similar dust column densities. The extinction curve might then be predicted to change from $E' < 1$ to $E' > 1$ mag. We consider this possibility in Figure 16.

Figure 16 shows that we find that the distribution of $R'(V)$ is remarkably independent of $E'$. When $E' < 0.5$ mag, the uncertainty in $R'(V)$ increases significantly, though to as low as $E' = 0.2$ mag we find no significant change in $R'(V)$. At the highest $E'$, there is a slight tendency for $R'(V)$ to increase, though the amount is small (0.1) and high $R'(V)$ stars are easier to observe in $g_{P1}$, since $g_{P1}$ band extinction decreases with $R'(V)$ at fixed $E'$. We conclude that there is no trend in $R'(V)$ with $E'$ for $E' < 2$ mag. There is some suggestion of a trend for $E' > 2$ mag, but we are not confident of its significance because our $g_{P1}$ band photometry is insufficiently deep.

We interpret Figure 16 as indicating that there is little change in dust extinction curve properties at $E' = 1$ mag, despite the formation of ice mantles found by Whittet et al. (1988). However, we note that $E'$ is only a rough proxy for dust volume density, as high $E'$ can either indicate individual dense clouds or a number of diffuse clouds along the line of sight. Despite this, the majority of our sight lines are in the outer Galaxy where often much of the total dust column is found in a single cloud, so we expect $E'$ to be an acceptable proxy for dust density there. This ambiguity will be resolved by the APOGEE Reddening Survey, which specifically targets giants in the background of dense regions of local molecular clouds.

We can also map the variation in the dust extinction curve over the sky. Figure 17 shows the spatial distribution of $R'(V)$ for the APOGEE targets with $E' > 0.3$ mag. Large, coherent trends in $R'(V)$ are readily detected. We detect regions with $R'(V)$ as low as 2.9 and as high as 3.9. The most obviously detected cloud with atypical $R'(V)$ in the APOGEE footprint is the Rosette Nebula, at $(l, b) = (206^\circ, -2^\circ)$, which is found to have an $R'(V)$ of about 4. The extinction curve in the Rosette Nebula was formerly studied by Fernandes et al. (2012); they
<figure><img src="image_18.png" /><figcaption>**Figure 18.** *R'(V)* vs. Planck Collaboration et al. (2014) *β* and Meisner & Finkbeiner (2015) *f*₁, averaged on 1° scales. We find a strong correlation between the two quantities, suggesting that variations in the optical-infrared extinction curve and FIR SED have a related origin in dust physics.</figcaption></figure>

the two parameterizations are tracking the largely same variability in the dust SED.

The ideal comparison between dust emission spectra and extinction curves would be performed in high latitude molecular clouds where the stars used for tracing *R(V)* are behind the entire dust column and where the dust emission is dominated by a single cloud. The majority of the APOGEE targets, however, are at low Galactic latitudes where both of these conditions are violated. At these latitudes, even relatively simple measurements like FIR optical depth are problematic (e.g., Schlafly et al. 2014b). Ideal sight lines for comparison would also contain no significant CO, which complicates the SED modeling, and would be at high ecliptic latitudes, where zodiacal light is a small contributor to the total FIR emission. However, in Figure 18 we simply include all sight lines. We therefore expect that the underlying relationship between *R'(V)* and *β* or *f*₁ may be significantly stronger, though we currently lack adequate coverage of high latitude clouds in APOGEE to confirm this hypothesis. The correlation we observe persists regardless of the cuts on CO emission or ecliptic latitude we impose.

## 5.6. The “Gray” Component of the Extinction

Our analysis is only sensitive to the colors of stars. Any component of the extinction curve that uniformly extinguishes light across the optical and infrared is undetectable in this analysis. This insensitivity stems from our ignorance of the distance to any of the sources we measure. This is a serious limitation, because many practical applications of the extinction curve require knowledge of, for instance, *A(J)/A(H)*, which we are unable to measure.

Many measurements from the literature (e.g., CCM) determine *A(V)* from reddenings using a fixed extinction curve in the infrared. This is essentially the solution we adopt in the Appendix, but we are hesitant to employ this procedure, since we seek to measure the variation of the extinction curve.

The simplest solution to this problem would be to adopt measurements of the gray component of the extinction curve from studies of globular clusters, the Galactic bulge (e.g., Stutz et al. 1999; Nataf et al. 2013, 2015), or external galaxies, like the SMC, LMC, or Andromeda (e.g., Maíz Apellániz et al. 2014; De Marchi et al. 2015), where the distances to all of the stars are known. Because the variations in *R'(V)* we measure are presumably linked to variations in the gray component of the extinction curve, we need measurements of the gray component over a wide range of *R'(V)*, and preferably not limited to a single star-forming region. This may be possible in Andromeda with the PHAT survey (Dalcanton et al. 2012), but mapping the gray component of extinction throughout the Milky Way will have to wait for parallaxes from the Gaia mission (Perryman et al. 2001).

## 5.7. The Intrinsic Colors of Giants Compared with Synthetic Models

Our model produces estimates of the intrinsic colors of giant stars as a function of their temperature and metallicity. Because faint, metal-rich giants are almost always located far away and in the disk, they are typically significantly reddened. Our colors may therefore be some of the best empirical estimates of the intrinsic colors of these stars in the combined PS1, 2MASS, and WISE bands, useful for informing models. We note however two limitations: first, we are ultimately tied to the $y_{P1} - K$ color from the MARCS synthetic stellar grid, and there are hints of problems of around a few hundredths with those models (Section 4.4). Second, we are projecting observed reddened colors back to intrinsic colors across a typical $A(V) \approx 2$ mag linearly along the reddening vector; we should be considering the full non-linear effect of reddening on magnitudes here, which may make a difference of a couple hundredths in the optical.

Figure 19 shows the intrinsic colors we derived, as compared with observed, SFD-derreddened colors and with synthetic colors from the MARCS grid, for stars with $E(B-V)_{SFD} < 0.2$. Observed colors are shown by circles, colored by their temperature. The model colors of Section 3.1 are shown by the solid black lines; the three lines correspond to $[Fe/H] = (-0.75, -0.25, 0.25)$ and the color along the line corresponds to temperature. The model colors match the observed colors well (up to photometric noise in the observed colors, most obvious at long wavelengths). Synthetic MARCS colors are shown with dashed lines. Generally, the synthetic and model colors are in close agreement, though differences of up to several hundredths are present. Offsets of similar size
<figure><img src="image_19.png" /><figcaption>**Figure 19.** Observed, SFD-dereddened colors of APOGEE stars, compared with the model intrinsic colors we determine and synthetic colors from the MARCS spectral grid. Points are colored by their corresponding temperature, and in dense regions we have replaced the points with a colored bin giving the average temperature of all points in that bin. The contours show the density of points. The three solid lines show the model colors for metallicities of −0.75, −0.25, and 0.25, while the three dashed lines show the synthetic MARCS colors for the same metallicities. The color along each line shows the temperature of the model. In general, there is good agreement between the observed dereddened colors (points), the model intrinsic colors (solid lines), and the synthetic colors (dashed lines). The biggest differences appear for the coldest stars, where the synthetic models are expected to have difficulties due to the formation of molecules.</figcaption></figure>

were also found by Schlafly & Finkbeiner (2011) when comparing observed colors of stars with MARCS models, and may be partially due to errors in the APOGEE temperature and metallicity scales.

The largest differences between the model colors and the synthetic colors occur for the coolest stars, $T < 4000$ K. For

these stars, the onset of molecule formation makes their synthetic colors especially uncertain, so color differences are expected.

In general, our model colors provide a better match to the observed, SFD-dereddened colors than the synthetic colors. However, very few of the coolest $T < 4000$ K stars are present
in the $E(B-V)_{SFD} < 0.2$ sample where we believe the SFD-dereddening to be adequate, making it hard to assess which of the two is actually more accurate in this region. That said, we expect our technique to perform nearly as well for cold stars as warm stars, subject only to the accuracy of our fit in Figure 9, which itself depends on an accurate MARCS $y_{P1}-K$ synthetic color.

## 6. CONCLUSION

We present sensitive measurements of optical-infrared reddenings to 37,000 stars in the Galactic disk, enabled by APOGEE spectroscopy and PS1, 2MASS, and WISE photometry. The typical star has a reddening of 0.65 mag $E(B-V)$, and the stars probe the dust over much of the disk within a few kiloparsecs, making for a uniquely powerful set of reddening measurements.

We use these reddening measurements to determine the shape of the extinction curve in the Milky Way and its variation. We draw the following conclusions:

1. We make new measurements of the mean extinction curve in the optical through infrared. Agreement is good with past measurements and the extinction curve of FM09, but other extinction curves (CCM, F99, FM04, M14) provide poor matches to the full optical-infrared extinction curve.

2. We find that the shape of the extinction curve in the optical and infrared can be well characterized by a single parameter, for instance, $R(V)$. The curvature of the extinction curve increases with decreasing $R(V)$ throughout the optical and infrared, with smaller variation in the NIR than in the optical.

3. The shape of the extinction curve is surprisingly uniform, with $\sigma(R(V)) \approx 0.18$, and fewer than 1% of sightlines having $R(V) > 4$.

4. The variation in $R(V)$ that does exist is uncorrelated with column density for $E(B-V) < 2$.

5. The variations in $R(V)$ we observe are spatially coherent on large scales (>30°), suggesting that most of the observed variation in $R(V)$ is driven by processes that act on large scales. In particular, the lack of correlation with column density suggests that the variation is tracing much more than grain growth in dense molecular clouds.

6. Finally, we discover a previously unknown, strong correlation between the thermal dust SED and the shape of the extinction curve.

Our work leaves at least two important questions unanswered. First, we are unable to measure any gray component of the extinction, or the variation in this gray component with the variation we observe in the extinction curve. Addressing this question is crucial to providing measurements of the full extinction curve $A(\lambda)$, rather than the reddening curve $E(\lambda - \lambda_0)$. Up to now, extinctions have been much more challenging to measure than reddenings, but the upcoming release of data from Gaia will resolve this long-standing problem.

Second, we have discovered large, coherent variations in the shape of the dust extinction curve, and a relative absence of expected small-scale variations in dense regions. These variations are strongly negatively correlated with maps of the FIR spectral index of the dust SED. We are aware of no adequate theoretical framework for understanding the large-scale variations.

In future work, we plan to better characterize ISM structures leading to these signals, in the hopes of providing observational clues to the source of the variations. Recent improvements in 3D dust mapping (Schlafly et al. 2014a; Green et al. 2015) will facilitate this effort.

The forthcoming APOGEE Reddening Survey, part of APOGEE-II, will also help to address this issue. The survey targets bright red giants in the background of the densest parts of several nearby molecular clouds: Orion, Perseus, Taurus, and Monoceros-R2. APOGEE temperatures, metallicities, and gravities for these stars combined with PS1, 2MASS, and WISE photometry will allow us to study the shape of the extinction curve in the densest parts of local molecular clouds, indicating to what extent the extinction curve varies in dense regions.

ES acknowledges support for this work provided by NASA through Hubble Fellowship grant HST-HF2-51367.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555, and funding from Sonderforschungsbereich SFB 881 "The Milky Way System" (subproject A3) of the German Research Foundation (DFG). DF acknowledges support of NASA grant NNX10AD69G. GMG and DPF are partially supported by NSF grant AST-1312891. The work of JK was supported by the Deutsche Forschungsgemeinschaft priority program 1573 ("Physics of the ISM"). NFM gratefully acknowledges the CNRS for support through PICS project PICS06183.

Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web Site is http://www.sdss.org/.

The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max Planck Institute for Astronomy (MPIA), the Max Planck Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.

This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.
The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen's University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE).

APPENDIX

We have made new measurements of the extinction curve in the optical through infrared with an unprecedentedly sensitive sample of reddening targets, and find substantial disagreement with existing extinction curves. Accordingly, we wish to provide a new extinction curve consistent with our data.

Unfortunately, our measurements are lacking in three important ways. First, we measure the extinction curve only in broad photometric bands, and are largely insensitive to the shape of the curve within those bands. Second, we are insensitive to any gray component of the extinction curve—but important quantities like $A(\lambda_1)/A(\lambda_2)$ require knowledge of the gray component. Third, we are only able to measure the extinction curve from the $g_{P1}$ to W2 bands.

To address the first of these problems, we simply interpolate between the broad photometric bands we measure with a cubic spline, acknowledging that we are ignorant of the detailed shape. To address the second, we use the measurement $A(H)/A(K) = 1.55$ from Indebetouw et al. (2005) to fix the gray component of the extinction curve. This latter procedure is problematic if there is significant variation in the infrared extinction curve, but at least that is the least variable part of the extinction curve that we observe. Despite this limitation, the extinction curves we derive by this technique look reasonable, and this procedure provides a simple way to fix the gray component observationally. To address the third problem, we can only caution the reader that the extinction curve is unreliable outside 5000–45000 Å, and that at the edges of this range the slope is relatively uncertain. At the edges of the spline we set the third derivative of the spline to zero, an arbitrary choice.

We construct the extinction curve by first determining the monochromatic wavelengths to which the broad band measurements of Table 2 apply. In analogy with the definition of isophotal wavelengths, we define isoextinction wavelengths $\lambda_b^e$

$$A(\lambda_b^e)/A_0 = \frac{d}{dA_0} \left( -2.5 \log \int d\lambda F(\lambda) 10^{-A(\lambda)/2.5} T_b(\lambda) \right), \quad (9)$$

where $A(\lambda)$ is the extinction at the wavelength $\lambda$, $A_0$ is the extinction at some reference wavelength, $F(\lambda)$ is the flux from a star, and $T_b(\lambda)$ is the total system throughput in the band b.

<figure><img src="image_20.png" /><figcaption>Figure 20. Extinction curves for $x = 0.04, 0.0$ and $-0.04$, according to our APOGEE measurements. Curves are normalized at 5420 Å, and have gray components fixed according to $A(H)/A(K) = 1.55$ from Indebetouw et al. (2005).</figcaption></figure>

Roughly, the isoextinction wavelengths are the wavelengths at which the monochromatic extinction equals the rate of change of extinction, in magnitudes, in the broad photometric band b.

Given the isoextinction wavelengths, a smooth, monochromatic extinction curve reproducing our measurements is given by a cubic spline connecting our measured $dA_b/dA_0$, as given in Table 2. The full procedure is then:

1. Produce an extinction vector as $A = R_0 + xdR/dx$, with $R_0$ and $dR/dx$ taken from Table 2.

2. Fix the gray component, by sending $A \to A + C$, so $A(H)/A(K) = 1.55$ (Indebetouw et al. 2005).

3. The extinction curve $A(\lambda)$ is given by a cubic spline passing through $A$ at the wavelengths $\lambda_b^e$ (Table 2), additionally imposing that the third derivative of the spline is zero at the boundaries.

Code implementing the above procedure is available at our web site$^{15}$. Figure 20 shows three example extinction curves, for $x = 0.04, 0.0$, and $-0.04$, corresponding roughly to $R(V) = 3.6, 3.3$, and 3.0.

This procedure is somewhat circular. The isoextinction wavelengths depend on the shape of the extinction curve $A(\lambda)$, but we use these wavelengths to determine the extinction curve. We resolve this circularity by initially setting $A(\lambda)$ to be the extinction curve of F09, solving for the isoextinction wavelengths, and using those wavelengths to construct our own curve. This process is then iterated with the new extinction curve until the isoextinction wavelengths have converged.

In detail, the isoextinction wavelengths depend on the amount of extinction and the spectrum of the source. We find the isoextinction wavelengths corresponding to a 4500 K star with [Fe/H] = 0 and $\log g = 2.5$ at a reddening $g - r = 0.65$, roughly typical of our sample. The isoextinction wavelengths also depend on the shape of the extinction curve. However, for changes of $R(V)$ of 0.3, the wavelengths change by less than one part in a thousand, so we neglect this variation, and always use the isoextinction wavelengths corresponding to our mean extinction curve.

<footnote n="15">15 http://faun.rc.fas.harvard.edu/eschlafly/apored/extcurve_s16.py</footnote>
One can compute $R(V) = A(V)/E(B-V)$ directly from our extinction curve. For $x=0$, the curve gives $R(V) = 3.6$, significantly larger than traditionally associated with the diffuse ISM. However, the central wavelength of the Landolt B band is roughly 500 nm blueward of the central wavelength of the PS1 g band, our bluest band, so computing $R(V)$ via this method requires extrapolating our extinction curve. If instead the extinction curve of F99 and this work are spliced together, we obtain an $R(V)$ lower by a few tenths; direct computation of $R(V)$ is sensitive to how the extinction curve is extrapolated.

## REFERENCES

Alam, S., Albareti, F. D., Allende Prieto, C., et al. 2015, arXiv:1501.00963

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Casagrande, L., & VandenBerg, D. A. 2014, MNRAS, 444, 392

Clayton, G. C., Gordon, K. D., Bianchi, L. C., et al. 2015, ApJ, 815, 14

Cutri, R. M., Wright, E. L., Conrow, T., et al. 2013, Explanatory Supplement to the AllWISE Data Release Products, Tech. Rep.

Dalcanton, J. J., Williams, B. F., Lang, D., et al. 2012, ApJS, 200, 18

Davenport, J. R. A., Ivezić, Ž., Becker, A. C., et al. 2014, MNRAS, 440, 3430

De Marchi, G., Panagia, N., Sabbi, E., et al. 2015, arXiv:1510.08436

Draine, B. T. 2003, ARA&A, 41, 241

Edvardsson, B. 2008, PhST, 133, 014011

Eisenstein, D. J., Weinberg, D. H., Agol, E., et al. 2011, AJ, 142, 72

Fernandes, B., Gregorio-Hetem, J., & Hetem, A. 2012, A&A, 541, A95

Finkbeiner, D. P., Davis, M., & Schlegel, D. J. 1999, ApJ, 524, 867

Fitzpatrick, E. L. 1999, PASP, 111, 63

Fitzpatrick, E. L. 2004, in ASP Conf. Ser. 309, Astrophysics of Dust, ed. A. N. Witt, G. C. Clayton, & B. T. Draine (San Francisco, CA: ASP), 33

Fitzpatrick, E. L., & Massa, D. 1986, ApJ, 307, 286

Fitzpatrick, E. L., & Massa, D. 1988, ApJ, 328, 734

Fitzpatrick, E. L., & Massa, D. 1990, ApJS, 72, 163

Fitzpatrick, E. L., & Massa, D. 2005, AJ, 130, 1127

Fitzpatrick, E. L., & Massa, D. 2007, ApJ, 663, 320

Fitzpatrick, E. L., & Massa, D. 2009, ApJ, 699, 1209

Foster, J. B., Mandel, K. S., Pineda, J. E., et al. 2013, MNRAS, 428, 1606

Gabriel, K. R., & Zamir, S. 1979, Technometrics, 21, 489

García Pérez, A. E., Allende Prieto, C., Holtzman, J. A., et al. 2015, arXiv:1510.07635

Green, G. M., Schlafly, E. F., Finkbeiner, D. P., et al. 2015, ApJ, 810, 25

Gunn, J. E., Siegmund, W. A., Mannery, E. J., et al. 2006, AJ, 131, 2332

Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951

Hodapp, K. W., Kaiser, N., Aussel, H., et al. 2004, AN, 325, 636

Holtzman, J. A., Shetrone, M., Johnson, J. A., et al. 2015, arXiv:1501.04110

Indebetouw, R., Mathis, J. S., Babler, B. L., et al. 2005, ApJ, 619, 931

Jolliffe, I. T. 2002, Principal Component Analysis (2nd ed., Berlin: Springer)

Jones, A. P., Fanciullo, L., Köhler, M., et al. 2013, A&A, 558, A62

Krelowski, J., & Strobel, A. 2012, AN, 333, 60

Magnier, E. 2006, in The Advanced Maui Optical and Space Surveillance Technologies Conf., ed. S. Ryan (Kihei, HI: The Maui Economic Developer Board), E50

Magnier, E. 2007, in ASP Conf. Ser. 364, The Future of Photometric, Spectrophotometric and Polarimetric Standardization, ed. C. Sterken (San Francisco, CA: ASP), 153

Magnier, E. A., Liu, M., Monet, D. G., & Chambers, K. C. 2008, in IAU Symp. 248, A Giant Step: from Milli- to Micro-arcsecond Astrometry, ed. W. J. Jin, I. Platais, & M. A. C. Perryman (Cambridge: Cambridge Univ. Press), 553

Magnier, E. A., Schlafly, E., Finkbeiner, D., et al. 2013, ApJS, 205, 20

Mańz Apellániz, J., Evans, C. J., Barbá, R. H., et al. 2014, A&A, 564, A63

Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2015, arXiv:1509.05420

Majewski, S. R., Zasowski, G., & Nidever, D. L. 2011, ApJ, 739, 25

Meisner, A. M., & Finkbeiner, D. P. 2015, ApJ, 798, 88

Mulas, G., Zonca, A., Casu, S., & Cecchi-Pestellini, C. 2013, ApJS, 207, 7

Nataf, D. M., Gonzalez, O. A., Casagrande, L., et al. 2015, arXiv:1510.01321

Nataf, D. M., Gould, A., Fouqué, P., et al. 2013, ApJ, 769, 88

Nidever, D. L., Holtzman, J. A., Allende Prieto, C., et al. 2015, arXiv:1501.03742

O'Donnell, J. E. 1994, ApJ, 422, 158

Onaka, P., Tonry, J. L., Isani, S., et al. 2008, Proc. SPIE, 7014, 70140D

Perryman, M. A. C., de Boer, K. S., Gilmore, G., et al. 2001, A&A, 369, 339

Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2014, A&A, 571, A11

Roweis, S. 1998, in Proc. 1997 Conf. on Advances in Neural Information Processing Systems 10, NIPS '97, ed. M. I. Jordan, M. J. Kearns, & S. A. Solla (Cambridge, MA: MIT Press), 626

Sale, S. E., & Magorrian, J. 2015, MNRAS, 448, 1738

Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schlafly, E. F., Finkbeiner, D. P., Jurić, M., et al. 2012, ApJ, 756, 158

Schlafly, E. F., Finkbeiner, D. P., Schlegel, D. J., et al. 2010, ApJ, 725, 1175

Schlafly, E. F., Green, G., Finkbeiner, D. P., et al. 2014a, ApJ, 786, 29

Schlafly, E. F., Green, G., Finkbeiner, D. P., et al. 2014b, ApJ, 789, 15

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

Schultheis, M., Kordopatis, G., Recio-Blanco, A., et al. 2015, A&A, 577, A77

Schultheis, M., Zasowski, G., Allende Prieto, C., et al. 2014, AJ, 148, 24

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Stubbs, C. W., Doherty, P., Cramer, C., et al. 2010, ApJS, 191, 376

Stutz, A., Popowski, P., & Gould, A. 1999, ApJ, 521, 206

Tamuz, O., Mazeh, T., & Zucker, S. 2005, MNRAS, 356, 1466

Tonry, J., & Onaka, P. 2009, in Advanced Maui Optical and Space Surveillance Technologies Conf., ed. S. Ryan (Kihei, HI: The Maui Economic Developer Board), E40

Tonry, J. L., Stubbs, C. W., Lykke, K. R., et al. 2012, ApJ, 750, 99

Trumpler, R. J. 1930, PASP, 42, 214

Tsalmantza, P., & Hogg, D. W. 2012, ApJ, 753, 122

Valencic, L. A., Clayton, G. C., & Gordon, K. D. 2004, ApJ, 616, 912

Wang, S., & Jiang, B. W. 2014, ApJL, 788, L12

Whittet, D. C. B., Bode, M. F., Longmore, A. J., et al. 1988, MNRAS, 233, 321

Wilson, J. C., Hearty, F., Skrutskie, M. F., et al. 2010, Proc. SPIE, 7735, 1

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868

Ysard, N., Abergel, A., Ristorcelli, I., et al. 2013, A&A, 559, A133

Yuan, H. B., Liu, X. W., & Xiang, M. S. 2013, MNRAS, 430, 2188

Zasowski, G., Johnson, J. A., Frinchaboy, P. M., et al. 2013, AJ, 146, 81

Zasowski, G., Majewski, S. R., Indebetouw, R., et al. 2009, ApJ, 707, 510