The Near-infrared Extinction Law at High and Low Galactic Latitudes

The Milky Way dust extinction curve in the near-infrared (NIR) follows a power-law form, but the value of the slope, β NIR, is debated. Systematic variations in the slope of the Milky Way UV extinction curve are known to be correlated with variations in the optical slope (through R V ), but whether such a dependence extends to the NIR is unclear. Finally, because of low dust column densities, the NIR extinction law is poorly understood at high Galactic latitudes where most extragalactic work takes place. In this paper, we construct extinction curves from 56,649 stars with Sloan Digital Sky Survey (SDSS) and Two Micron All Sky Survey photometry, based on stellar parameters from SDSS spectra. We use dust maps to identify dust-free stars, from which we calibrate the relation between stellar parameters and intrinsic colors. Furthermore, to probe the low-dust regime at high latitudes, we use aggregate curves based on many stars. We find no systematic variation of β NIR across low-to-moderate dust columns (0.02 < E(B − V) ≲ 1), and report average β NIR = 1.85 ± 0.01, in agreement with the law in the 2019 Fitzpatrick et al. study, but steeper than the Cardelli et al. and 1999 Fitzpatrick laws. Star-to-star scatter in β NIR is relatively small (σ(β NIR) = 0.13). We also find no intrinsic correlation between β NIR and R V (there is an apparent correlation that is the result of the correlated uncertainties in the two values). These results hold for typical sightlines; we do not probe very dusty regions near the Galactic Center, nor rare sightlines with R V > 4. Finally, we find R H = 0.345 ± 0.007 and comment on its bearing on Cepheid calibrations and the determination of H 0.


Introduction
Determinations of the dust extinction law in our galaxy, and its variations along different lines of sight, are relevant for many areas of astronomy (see the recent reviews of Salim &Narayanan 2020 andHensley &Draine 2021).The result from Cardelli et al. (1989, hereafter CCM) that UV and optical extinction curves in the Milky Way belong to a singleparameter family was a watershed development toward systematizing observed variations.That single parameter, R V , is essentially the slope of the extinction curve in the optical region (1/R V = A B /A V − 1).However, in CCM, this R V dependence did not extend to the near-infrared (NIR) part of the extinction curve, where they fit a power law of the form A λ ∝ λ −βNIR to the data from Rieke & Lebofsky (1985), yielding β NIR = 1.61.CCM find this slope is applicable regardless of R V .Indeed, the existence of such a "universal" NIR extinction law had been found in several prior and contemporary studies (e.g., Jones & Hyland 1980;Koornneef 1983;Smith 1987).
Though the power-law nature of the NIR extinction curve has been concretely established by numerous studies, the value of the slope varies from study to study.Early studies prior to and after Rieke & Lebofsky (1985) found values in the range between about 1.6 to 1.85 (see the reviews of Draine 1989, Mathis 1990, and Draine 2003).Martin & Whittet (1990) found one of the highest values of the time, β NIR = 1.84; this value was subsequently used in the parameterization developed by Fitzpatrick & Massa (2007).More recent studies of the NIR portion of the curve have found higher values for β NIR (up to 2.64; Gosling et al. 2009).While there is clearly a wide spread in β NIR determined between studies, the origin of the range is unclear.One possible explanation is variance with the environment or amount of dust along sightlines, but this has not been concretely shown.While different studies have found differing slopes, few studies claimed that there is a systematic change in this slope.The general assumption still is that it is universal, and evidence for dependence on R V has been limited to very small-scale studies (Fitzpatrick & Massa 2009;Decleir et al. 2022).Many of the recent studies finding high β NIR have focused on the extreme-extinction region in the Galactic Center (e.g., Nishiyama et al. 2009;Gosling et al. 2009;Schödel et al. 2010;Hosek et al. 2018;Nogueras-Lara et al. 2021;Sanders et al. 2022), which, if the curve is dependent on the environment, could explain higher values at least in part.
What essentially all previous studies of NIR extinction laws have had in common is their focus in the plane of the Galaxy, i.e., low Galactic latitude regions of relatively high dust content.It is inherently more difficult to study dust in regions where there is very little, especially in the NIR.The effects of the dust on incoming light are subtle and the photometry must therefore be very precise.An additional problem for the determination of dust extinction curves at high latitudes arises for the empirical pair method (e.g., Fitzpatrick & Massa 1990) because the comparison stars rarely have zero dust along their sightlines and may have a comparable amount of dust to target high-latitude sightlines.This is not an issue where models are used to constrain extinction (see, e.g., Fitzpatrick & Massa 2005), but the requirement of high accuracy remains.
Although it is generally assumed that the nature of the dust (and therefore the extinction curve) is the same at high latitudes as it is at low latitudes, this has not been confirmed in the NIR.Since extragalactic astronomy and cosmology is largely con-ducted at high latitudes, it is important to constrain extinction in that low-dust regime, especially in the era of JWST (Gardner et al. 2006) and LSST (Ivezić et al. 2019), in addition to the upcoming Euclid (Racca et al. 2016) and Nancy Grace Roman (Akeson et al. 2019) missions, which will require extinction corrections of high accuracy.
Among the most prominent open questions in astronomy and cosmology is the accurate value of the Hubble constant (H 0 ).The most recent result based on Cepheids and Type Ia supernovae (Riess et al. 2022) is in 5σ tension with the result derived from cosmic microwave background anisotropies (Planck Collaboration et al. 2020).The local determination of H 0 uses Cepheid variables, which follow a tight relation between period and luminosity (Leavitt & Pickering 1912), allowing for very precise distance measurements.The observed magnitudes of the Cepheids must be corrected for dust extinction along the sightlines.There have been extensive studies into Cepheid calibrations and SN host galaxy dust, with no small amount of attention being paid to the dust extinction (e.g., Riess et al. 2016Riess et al. , 2021;;Mandel et al. 2017;Follin & Knox 2018;Brout & Scolnic 2021;Meldorf et al. 2023).Precise NIR extinction values are crucial for Cepheid calibrations specifically (the first rung in the distance ladder), and knowledge of the appropriate NIR extinction law is required.
In recent years, there has been a proliferation of extinction studies which use large samples of stars taken from surveys (e.g., Peek & Graves 2010;Schlafly et al. 2010;Schlafly & Finkbeiner 2011;Berry et al. 2012).Given such a sample, one can apply several different methods to derive intrinsic colors for stars, which can then be subtracted from observed colors to obtain reddenings.Some of these methods employ template spectra or atmospheric models to generate simulated comparison stars (e.g., Jones et al. 2011;Schlafly & Finkbeiner 2011), while others use a set of stars with very low extinction to construct a relation between their photometry and stellar parameters, usually T eff , log g, and [Fe/H] (e.g., Yuan et al. 2013Yuan et al. , 2015;;Li et al. 2018;Sun et al. 2018;Wang & Chen 2019).We follow a procedure of the latter kind to determine intrinsic colors for our sample of stars, with the goal of investigating NIR extinction curves at both high and low latitudes.
In this study, we attempt to clarify the potential dependence of β NIR on the characteristics of individual sightlines.Namely, we look for a correlation between β NIR and the dust column density along a sightline, and separately for a correlation between β NIR and R V .We construct extinction curves across a wide range of dust columns, made possible specifically at high latitudes by a large sample size.In the low-dust (high-latitude) regime, NIR extinction is most poorly constrained, simply due to the very small effect of dust on the incoming light.Because of this, we use aggregates of large numbers of stars from optical and NIR surveys (SDSS, 2MASS) to increase our sensitivity to the effects of the dust.
Section 2 describes our optical and NIR data sources.Section 3 details the methods we use for deriving intrinsic colors and subsequently constructing extinction curves.In Section 4, we present the results of our study, relating β NIR to the dust column and also to R V .We additionally comment on the bearing of our results on the Cepheid distance ladder calibration in the pursuit of a precise determination of the Hubble constant, H 0 .In Section 5, we discuss our results in the context of prior results from the literature.

Data and Sample
We incorporate both optical and near-IR photometry of Galactic stars into our analysis, along with stellar parameters derived from optical spectra.

Data Sources
We draw our optical photometry (ugriz) and spectroscopy from the Sloan Digital Sky Survey (SDSS; York et al. 2000) Data Release 16 (Ahumada et al. 2020), in particular its Legacy (henceforth referred to as SDSS) and SEGUE components, which collectively cover 14,555 square degrees.We use PSF magnitudes, which are most appropriate for stars.
The spectra from which our stellar parameters are derived come from any of three separate (but closely related) spectroscopic campaigns: the SDSS Legacy Survey (reported in the database as SDSS), SEGUE-1 (Yanny et al. 2009), and SEGUE-2 (Eisenstein et al. 2011).These surveys all employed the same twin spectrographs with wavelength coverage from 3900 Å to 9000 Å at a moderate resolution of R ∼ 1800.We take the stellar parameters T eff , log g, and [Fe/H] as reported in the DR16 sppParams table, derived using the SEGUE Stellar Parameter Pipeline (SSPP; Lee et al. 2008).Lee et al. (2008) derive these three parameters using multiple methods, with the weighted averages reported as adop.The colorindependent, spectroscopically-determined weighted averages are reported as spec.We choose the spec parameters instead of the recommended adop parameters because the former are independent from photometry.We also take the spectroscopic subClass of the stars.
For JHK near-infrared photometry, we use the Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006) Point Source Catalog (PSC), which covers the SDSS footprint in its entirety.We use the "default" magnitudes as reported in the PSC.Note that we will use K to refer to the 2MASS K s band henceforth.We also use UKIDSS near-infrared photometry (Lawrence et al. 2007), which covers only a portion of the SDSS footprint (but is deeper than 2MASS) to perform checks on 2MASS.

Sample Selection
We queried SDSS within the DR16 context in CasJobs to obtain the base optical sample from SDSS.We require primary objects classified as stars (class='STAR') and with a clean photometry flag (clean=1).We also require that valid T eff , log g, and [Fe/H] values are present, and that the spectroscopic flag zWarning is set to 0 or 16, as recommended by SDSS documentation.
From this original optical sample of ∼342,000 stars with spectra, we make several further quality cuts on the photometric data.We make a cut on the quality of field (score< 0.7), which removes a small clump of poor-quality fields.Finally, we require that magnitudes in all SDSS bands are present for all stars, and that their photometric uncertainties are less than 0.1 mag.After these cuts, our base optical sample contains ∼314,000 stars.
We make use of the TwoMass table within CasJobs to get JHK photometry for our sources.Since 2MASS is less deep than SDSS, the stars for which we recover JHK photometry number about 60% of the optical sample.Comparing 2MASS and UKIDSS, we find that 2MASS tends to overestimate the brightness of an object at faint magnitudes in K band.To eliminate this issue, we remove all stars with K > 14.8 from all analyses that pertain to the near-IR.In a large-survey-based endeavor such as the current study, it may not be possible to eliminate all photometric anomalies using only catalog flags and simple cuts.We therefore inspect color-color plots in order to eliminate any remaining outliers.In the optical, we find i to be the most potentially problematic band, so we eliminate outliers beyond 6 robust σ away from the running median in a plot of r − i vs. g − z.We also find it necessary to remove any stars with (u − g) < 0.7.We make an additional cut determined using a similar method to the r − i cut above, but this time using a plot of r − J vs. g − z.We have verified that our outlier cuts are not removing any real measurements arising from atypical R V values.
For the calibration sample (see Section 3), we use only stars with E(B − V ) SFD < 0.0125.In addition, we do not require that the stars have valid measurements in JHK (as we do for the program sample); rather, they only need valid measurements in r and the band in question for calibration.This allows us to maximize our calibration sample to allow for better precision in intrinsic color estimates.All other cuts still apply for this sample.
The final (optical+NIR) sample from which we draw reddened (program) stars for this study contains 56,649.Figure 1 shows the distribution of r magnitudes for the stars in our final program sample, along with their median.The distribution is perhaps more symmetrical than expected for a general selection of Milky Way stars; this is likely due to the matching with 2MASS, which removes many of the fainter stars.Since we are not concerned with completeness in this study, these effects are not important.
Figure 2 characterizes our sample with respect to the absolute value of Galactic latitude and the log of reddening E(u − z).
We consider E(u − z) a proxy for the dust column (see Section 3.2).We cover a wide range of absolute latitudes, but note that very few stars lie below |b| ≈ 7 • .There is a shift in the latitude distribution around E(u − z) = 0.4 (or equivalently E(B − V ) = 0.1), which we will use to distinguish between dust columns typical of low and high latitude.The banding in the figure is a result of the SEGUE discrete data-taking pattern.
The numbers of stars in each spectral subclass from our final program sample are shown in is also shown.Most of our program stars come from SEGUE, and fall generally within the F type.
3. Methods In this study, we endeavor to obtain parameterized nearinfrared extinction curves across a wide range of dust column densities.To do this, we need (1) a way to estimate intrinsic colors, (2) color excesses (reddenings) calculated using (1), and (3) a model to which the color excesses can be fitted to obtain a parametric absolute extinction curve (A(λ)).These steps are presented in the sections to follow.

Intrinsic Color Estimates
In order to construct an extinction curve for the dust along a given sightline, estimates of the intrinsic (unreddened) colors of the object are required.We make use of stellar parameters determined from SDSS spectroscopy to estimate the intrinsic colors of stars in our sample.In general terms, the spectral type of a star is related to its color, but there is still a significant spread of colors even within one subclass.We show an example survey+subclass combination (SEGUE2 F9) in Figure 3.There is a general trend such that the observed colors become redder with increasing E(B − V ) SFD , the reddening as determined by Schlegel et al. (1998, SFD).However, there is quite a large range (∼0.25) in observed g − r color, even near E(B − V ) SFD = 0, i.e., where the SFD dust map tells us there is no dust.Thus, as expected, a pair-matching technique for constructing extinction curves, based on spectral type alone, is rather imprecise for later spectral types.Fortunately, with stellar parameters (log g, [Fe/H], and T eff ) in addition to the spectral class, we can more closely relate a star's intrinsic color to its physical properties than if spectral type alone were used.
In order to calibrate the relation between the intrinsic color and stellar parameters (including spectral type), we select the stars from our sample with very low reddening, E(B−V ) SFD < 0.0125.Note that the E(B − V ) SFD value is the upper limit because the reddening map includes measures the full column of dust in the direction of the star, and the star may not have all of this dust in front of it.The stars are then grouped by spectral subclass and survey (e.g., a group may be SEGUE2 F9).In principle, grouping by survey should not be necessary, but we maintain it to remove any potential systematics.
Within each survey+subclass grouping, we carry out calibrations for λ−r colors where λ refers to ugizJHK.We use r as the basis for colors because the photometry in that band is the most precise of all bands.Note that although the focus of this paper is on the NIR extinction curve, we still need reddenings in optical bands to determine E(u − z) (our proxy for dust column density) and R V .We require there to be at least 10 calibrator stars (E(B − V ) SFD < 0.0125) in a group for the regression to take place (see Section 2.2 for details on cuts).The intrinsic colors are fit using Several other studies have used stellar parameters of large number of stars to infer intrinsic colors, and we highlight some of the differences in our approach.Wang & Chen (2019) used log T eff and metallicity [M/H] (along with their second-order terms), as well as log g.They, however, did not split by spectral class (their sample contained only clump giants).For our sample where we split into spectral classes, a second-order term in T eff is not required.Furthermore, Wang & Chen (2019) consider as dust free the bluest 5% of stars at a given T eff , whereas we use SFD dust maps.Note also in our formulation the inclusion of an E(B−V ) SFD term, the purpose of which is to account for the fact that, strictly speaking, these calibrator stars have some (very small) reddening.Schlafly & Finkbeiner (2011) also use SSPP parameters as the basis for their intrinsic color prediction, but the correspondence between the parameters and colors comes from spectrophotometric models of stars, whereas in our case it is based on empirically selected dust-free stars.
We use the Huber Regressor as implemented in the Python package scikit-learn (Pedregosa et al. 2011) to determine the coefficients in Equation 1.This regression model employs the Huber loss function (Huber 1964), which allows it to be less sensitive to outliers than more standard loss functions, while not ignoring them altogether.The weights are set to 1/σ 2 , where σ is the photometric error on the color.We set the parameter ϵ = 3.5, which is similar to 3.5σ outlier rejection.
The regression produces one set of 7 coefficients for each color within each survey+subclass grouping.We calculate the standard deviation of the residuals (σ res ) in each case as an assessment of the constraining power of that group or stars (which is some combination of fit quality and photometry errors).Figure 5 gives an idea of how σ res varies across the range of colors and groups.Each line corresponds to a survey+subclass group.Residuals are the smallest in the optical bands, and in each color the majority of groups have similar residuals.In certain cases, the residuals are large enough that the corresponding group of stars warranted exclusion from the analysis.While there is no absolute sense in which σ res can be used to exclude a group, we used a relative assessment within each color.For a given color λ − r, we plot a cumulative distribution function (CDF) of the number of program stars (the ones to which the calibration will be applied) ranked by σ res .This allows us to eliminate groups which would provide only a modest increase in final program sample size, but at the cost of relatively large σ res .We make note of a cutoff σ res value for each λ − r, and exclude that sur-vey+subclass group in that color from further analysis.The gray-shaded region in Figure 5 indicates where we remove the grouping for a given color based on the CDF method described above.
In Figure 6  modified the model by replacing T eff with log T eff , and added a second-order term in that parameter.These changes (especially the addition of the second-order term) improved the universal calibration beyond its performance using Equation 1 as-is.
However, there is still some curvature to the locus, with both A0 and K5 stars having appreciably poor predictions.There is also a larger spread in the middle of the locus than in the right panel, which shows our adopted type-specific calibration.Aside from a very small number of outliers, our method appears to work remarkably well across a range of observed g − r.As mentioned, we also split by "survey" (one of three spectroscopic campaigns under the umbrella of SDSS: the Legacy component, which we simply refer to as SDSS, along with SEGUE1 and SEGUE2).We find that given a fixed subclass, there is still subtle variation in the calibration among these; however, the effect is much smaller than the one from the subclasses themselves illustrated in Figure 6.

Color Excesses (Reddenings)
For each program star, i.e., those with E(B − V ) SFD > 0.0125 in survey+subclass groupings which were not excluded by the CDF cutoff (see Section 3.1) we use its stellar parameters with Equation 1 to determine intrinsic color (λ − r).We set E(B − V ) = 0, since we want the predicted color of the star if it had exactly zero dust along it sightline.We then subtract these predicted intrinsic λ − r from the observed colors to produce a set of color excesses: We also calculate E(u − z) as which we consider to be a proxy for the dust column.If extinction curves were identical across the sky (i.e., R V = const.),any color excess would be equally well correlated with the dust column, and one would pick whichever reddening is best measured, or is closest to the traditional E(B − V ).However, the extinction curve is not universal; Butler & Salim (2021) showed that A u correlates with N (H), and therefore the dust column, better than any other optical/NIR extinction.This, combined with the fact that z is the optical band least affected by dust, leads us to choose E(u − z) as the best available proxy for the dust column.

NIR Extinction Curves
There are two principal ways in which we construct curves: (1) for individual stars that have higher reddening (low Galactic latitude regime), and (2) for aggregated stars.In the lowdust/high-latitude regime, many individual stars have very uncertain color excesses, making curve derivation meaningless.For that reason, in method (2), we instead use the power of our large sample to aggregate the stars by averaging their color excesses in bins of dust along their sightline (namely, its proxy E(u − z)).
To determine the parameterized NIR extinction curve, we choose a power law, which is known to provide a good description of the extinction curve between 0.75 and ∼5 µm (e.g., Martin & Whittet 1990;Fitzpatrick & Massa 2009;Wang & Chen 2019).In our dataset this range includes bands i through K.We switch our basis (or "anchor") from r to i: , because i, unlike r, belongs to this power-law region of the extinction curve.The power law Figure 6.Observed SDSS g − r color vs. g − r color predicted from stellar parameters, using Equation 1 and the method described in Section 3.1.Shown are dust-free stars from which calibrations are produced, i.e., for which the observed color is essentially the intrinsic color.Left: The case where calibrations are based solely on stellar parameters (not used in this work).We used a log T eff instead of T eff here, along with an added second-order term in T eff .See text near the end of Section 3.1.Right: The calibrations used in this work, derived from stellar parameters and done individually for each survey+subclass group.formulation of the reddening curve is where λ i = 7510 Å.The two free parameters are the slope of the power law (β NIR ) and the extinction in i band (A i ), the latter of which allows for any reddening E(λ − i) to be transformed into an absolute extinction A λ .Technically, the extinctions in Equation 2 should be monochromatic, whereas we use bandpass extinctions.However, in the optical and NIR, extinctions within bandpasses are close to monochromatic extinctions at their central wavelengths.
To generate curves from aggregated stars, we create 22 bins in E(u − z), the edges of which are given in the titles of the panels in Figure 7.The bins become progressively wider with increasing E(u − z) due to the majority of stars being at low dust.For 18% of the stars in our sample, the E(u − z) value we determine is negative, which can be attributed to the inherent uncertainties in the determination of intrinsic colors (and therefore color excesses) of stars which have very small amounts of dust along their sightlines.Such stars are not included in the analysis.
In each E(u − z) bin, we fit Equation 2 to the median of the E(λ − i) values of the individual stars.We restrict our sample in a given bin to only those stars which have a valid magnitude and error measurement in all four bands zJHK (in addition to passing the cuts in Section 2.2).
The fits are carried out using the unweighted non-linear least squares method, as implemented in the curve_fit function within the scipy (Virtanen et al. 2020) package in Python.We return A i and β values along with their errors.We choose to perform the fits without weights because the proportionally much smaller error on E(λ − z) disproportionately affects the fits (as compared to its NIR counterparts).We determine the effective errors of E(λ−i) asking that they agree with the power-law model.To do so, we calculate a χ 2 red value for each fit, and rescale the nominal errors of the median E(λ−i) such that they produce χ 2 red ≡ 1.The nominal error of the median where n is the total number of stars in the bin.
To derive parameterized extinction curves for individual stars, we use much the same method as described above.We also derive R V values for both aggregated and individual stars.To summarize, we derive NIR extinction curves parameterized by the power-law slope β NIR (with normalization A i ), as well as R V , for stars aggregated in 22 reddening bins spanning 0 < E(u−z) < 4, and for 2700 individual stars with E(u−z) > 1.We use the Fitzpatrick (1999) extinction law at R V = 3.1 to convert between reddenings in SDSS colors to reddenings in Johnson colors, finding E(B − V ) = 0.86E(g − r) and This conversion is somewhat sensitive to the choice of extinction law.We use the Fitzpatrick (1999) extinction law rather than e.g.CCM because the former has been shown to correctly predict E(g − r) (Schlafly et al. 2010), whereas the CCM prediction is about 20% too small.For diffuse dust, extinction curves should not depend on the dust column alone, unless the dust composition and grain size distribution are different (which in principle they can be between high and low latitudes).However, this expectation has not been tested in the NIR at high and low latitudes simultaneously.The panels in Figure 7 show aggregate NIR reddening as a function of inverse wavelength, where each panel is a different bin of the dust column, E(u − z).The dust column increases to the right and downward within the set of panels.In each panel we show the best fitting power-law curves and their  Ai = 0.24 βNIR = 1.84   parameters.Also shown for reference are standard Fitzpatrick (1999) extinction curves with R V = 3.1, scaled to match the mean E(u − z) of the bin.The Fitzpatrick (1999) curve is not defined as a strict power law in any wavelength regime.In some bins it matches the observations better than in others, without any particular pattern.Note that the range of E(λ−i) increases with each succeeding bin.We can see that the power law represents a good description of reddening regardless of dust column density.
In this study we are particularly interested in how β NIR may vary with the dust column, which is shown in Figure 8.The purple line represents the β NIR values determined by the fits to aggregated (binned) stars, as described above.We also include a secondary x-axis with E(B − V ), for convenience and ease of comparison with other studies.This E(B − V ) is related to E(u − z) by the simple scaling E(B − V ) = 0.26E(u − z), determined from all stars in our sample.Note that the reddening is presented on a log scale.A large portion of the horizontal span is occupied by rather minuscule reddening.
It is clear from Figure 8 that there is no strong correlation between β NIR and the amount of dust along a sightline.As a further test, we show, in the same bins, the median power-law slopes of individual stars with E(u − z) > 1 (see Section 3.3).Medians of β NIR based on extinction curves for individual stars (orange line in Figure 8) agree with those from the aggregated stars.
Since we find no correlation between β NIR and the dust column, we find it appropriate to report a single value for our study: β NIR = 1.85 ± 0.01.This value is determined by taking a weighted average of the β NIR values from the 22 bins.We report it along with separate values for high-and low-latitude regimes (corresponding to low and high values of E(u − z), respectively) in Table 1, which are again consistent with a single overall value.4.2.NIR Extinction Curve vs. R V UV and optical extinction curves are known to vary as a function of R V , but whether that extends to the NIR regime has not been firmly established.The relation we find between β NIR and R V is plotted in the leftmost panel of Figure 9.Note that in this analysis we use individual (high reddening) extinction curves, since little variation in R V (and β NIR ) is seen for aggregated stars.There is a strong apparent correlation between β NIR and R V .At the average Milky Way value of R V = 3.1, the points reassuringly passes quite near the commonly used β NIR = 1.84 (Martin & Whittet 1990;Fitzpatrick & Massa 2007), but the ranges in both β NIR and R V are quite large.Although the correlation between β NIR and R V appears quite strong based on formal values, we will show that it is not a true reflection of the intrinsic relationship between the two parameters.
The β NIR and R V values we determine for individual sightlines have associated uncertainties which stem from photometric and spectroscopic observations as well as our intrinsic color estimation process (Section 1).Larson & Whittet (2005) find that, for low-extinction sightlines, relatively small photometric errors can produce very large errors in R V .Furthermore, although β NIR and R V are independent quantities, the former being the slope of the extinction curve in the NIR whereas the latter being related to the slope in the optical (1/R V = A B /A V − 1), they both depend on the determination of the absolute extinction, which via Equation 2 is correlated with β NIR .Indeed, Fitzpatrick & Massa (2009)  which are not independent, and follow the general shape of the correlation we see in Figure 9.To test the extent to which the apparent correlation is driven by correlated errors, we perform an analysis using mock observations.The idea of the mock analysis is the following.We will assume that the actual extinction curves have identical β NIR irrespective of R V .We will then add realistic noise to the reddenings that these fixed-beta model curves predict and see what values of β NIR (and R V ) we recover.It is therefore important to distinguish between true (or input) β NIR and R V and the ones we derive ("observe") due to the uncertainties in reddening.
The mock analysis uses the same sample size as we have for individual stars.We generate curves by assuming the Fitzpatrick (1999) law, which essentially has an invariant NIR extinction curve, and therefore no dependence on R V .The input (intrinsic) R V values for the mock curves are drawn from a Gaussian distribution with mean R V = 3.1.The width of the distribution of intrinsic R V values is much smaller than the observed range and will be discussed shortly.
For each real star, we thus have its mock E(λ − i) values for λ = (z, J, H, K).We perturb these color excesses by an amount derived from the uncertainties of the real stars.We then employ the same method as was used in Section 3.3 to fit the power-law extinction model (Equation 2) and obtain "observed" β NIR and A i values.The output R V values of the mock stars are derived in the same way as in the observed case.A  slight "zero-point" adjustment is made to the derived mock β NIR in order to take into account that Fitzpatrick (1999) assumes a shallower beta than our overall value and is not a perfect powerlaw over the izJHK range.
As can be seen in the middle panel of Figure 9, the shape of the mock distribution bears striking resemblance to the shape of the observed distribution (left panel).This similarity strongly indicates that the apparent correlation we find between β NIR and R V for observed stars is an artifact of correlated uncertainties, and not the nature of the intrinsic relationship between the parameters.In order for the mock to most closely match the observed distribution, we have allowed the input R V to have a range of values drawn from a Gaussian distribution with a standard deviation of 0.28.Note that this intrinsic spread in R V is much smaller than the spread of mock "observed" R V values (the ones that result from taking observational uncertainties into account).
In reality, both the intrinsic R V and intrinsic β NIR can vary.Armed with the information that the intrinsic span in R V is small, we refine the value of intrinsic σ(R V ) and determine σ(β NIR ) in the following way.We add various amounts of scatter in both R V and β NIR to the mock distribution after it has been generated using a fixed R V = 3.1.When the scatter in the output R V values in the upper tail of the mock distribution (where the horizontal scatter is dominated by R V ) equals the scatter in the same part of the observed distribution, and the same is true for the scatter in β NIR in the lower tail (where the vertical scatter is dominated by β NIR ), we take those scatters to be the intrinsic.In this way, we get σ(R V ) = 0.24 and σ(β NIR ) = 0.13.We schematically plot these derived intrinsic distributions using 2D Gaussian contours in the right panel of Figure 9.We tested a scenario where the intrinsic β NIR and R V were linearly correlated, but found that values with correlation coefficients other than 0 were taking us away from the observed distribution.We therefore conclude that intrinsically β NIR and R V not only have a much smaller range than what the nominal values suggested, but also seem to be uncorrelated.
We also tested the correlation found in Fitzpatrick & Massa (2009) as opposed to a simple linear correlation.The run which matched the observed distribution the best had a similar spread in R V , but deviated ∼10% more in β NIR spread than the best match with no intrinsic correlation assumed.Our sightlines, with a relatively small range of R V values, thus disfavor a Fitzpatrick & Massa (2009) (or any other) correlation.

R V vs. the Dust Column
In addition to testing for a correlation between R V and β NIR , we also investigated whether R V might be dependent on the dust column density along sightlines, represented here by E(u − z).The resultant plot, Figure 10, is analogous to Figure 8, simply replacing β NIR on the y axis with R V .We again use the aggregated stars with the same bins as in Figure 7.The R V values are determined using the same conversion as described at the end of Section 3, this time employing the median E(g − r) in each bin (alongside median E(r − i) and fitted A i ), resulting in a median-based R V .
One may in principle not expect R V to vary with the dust column, since R V is generally thought to track the size distribution of grains.However, in denser regions, the size distribution may change (toward larger grains) due to higher probability of accretion and coagulation (e.g., Draine 2003).It is important to note that density, when applied to the discussion of size distribution, is physical as opposed to a column density, 0.01 0.1 1 since physical density is what drives grain interaction.Column density may happen to track physical density, but for a given column density, one can imagine several different dust distribution scenarios along a line of sight (e.g., clouds interspersed with relative voids versus a constant diffuse distribution).
Since we do not probe the densest regions in this work (at least in terms of column density), we may not expect to see any rise in R V with the dust column.At the lower dust columns associated with high Galactic latitudes, the expectation is less clear.Schlafly et al. (2016) produced a relationship between their R ′ V and E ′ (roughly equivalent to R V and E(B − V )) in their Figure 16.We can compare the trends in that figure directly with those in our Figure 10.At the highest E(B − V ), they find a rather sharp rise in R V (as may be expected per the above discussion), and we do as well in our highest-E(u − z) bin.However, the maximum E(B − V ) we probe is significantly lower than theirs (∼1 vs. ∼2.5).There is evidence for a slight rise around E(B − V ) = 1 in their results, but not of nearly the magnitude we find.Since our densest bin has only 18 stars, we take the sharp rise there with a grain of salt (although there may be some expectation for a change around E(B − V ) = 1, as Whittet et al. 1988 find the evidence that ice mantles form in almost all cases above that point).Schlafly et al. (2016) note a large increase in uncertainty below E(B − V ) = 0.5, but seem to show a rise below roughly that point; we find increased R V below E(B − V ) ∼ 0.1, albeit with large uncertainties at the very lowest dust columns.
Overall, given the low number of stars in the highest-dust bin and the relative unreliability of the results in the lowest 2-3 bins, we seem to find rough agreement with the established MW standard R V = 3.1 across a range of dust columns.There 0.01 0.1 1 vs. a proxy for the dust column (E(u−z)) using the same bins and fits as Figure 8.The dotted line indicates the adopted value of R H = 0.386 used by the SH0ES team in Cepheid dereddening en route to the derivation of H 0 (e.g., Riess et al. 2016Riess et al. , 2022)).The upper secondary is some evidence for a rise toward high latitudes, but it is fairly weak.

The NIR Extinction Curve and Cepheid Calibration
Recent Cepheid period-luminosity calibrations have been performed in the near-IR, in particular the H band. Whether the calibration is performed using Wesenheit magnitudes (Madore 1982): or unreddened magnitudes: one needs to know R H , an extinction curve-dependent parameter, precisely.Since modern measurements of the Hubble constant are based on Hubble Space Telescope photometry, H, V , and I here refer to the Hubble filters F160W, F555W, and F814W, respectively.The parameter R H (elsewhere called simply R or sometimes R E ) is of similar makeup to the more familiar R V : The most recent studies from SH0ES (Riess et al. 2016(Riess et al. , 2022) ) use for their baseline fit the value R H = 0.386, derived from the Fitzpatrick (1999) reddening law assuming a somewhat higher than standard R V = 3.3.Given that law, R H only changes by ∼0.005 between R V = 3.3 and R V = 3.1, so the choice of R V is not of great importance.The use of a constant R H for any Cepheid is equivalent to assuming a universal power-law slope in the NIR.Therefore, based on the analyses vs. a representation of the optical slope of the extinction curve (R V ) for three cases.In each case, standard literature values of R V = 3.3 and R H = 0.386 (Riess et al. 2016(Riess et al. , 2022) ) are included as brown dashed lines.These panels are directly analogous to the three β NIR panels, but show R H instead. Left: Observed individual stars with E(u − z) > 1. Center: Analogous to the left panel, but using a comparably sized sample of mock stars generated from the Fitzpatrick (1999) family of curves, which are consistent with a fixed β NIR (and roughly fixed R H ). See Section 4.2 for details.Right: Values calculated from the Fitzpatrick et al. (2019, FM19, green), Fitzpatrick (1999, F99, orange), CCM (purple), and Fitzpatrick & Massa (2007, FM07, cyan, R V = 3.1) literature extinction laws.The gray shaded ellipses represent 1σ, 2σ, and 3σ contours of a two-dimensional Gaussian with σ(R V ) = 0.24 and σ(R H ) = 0.10.These are intrinsic values derived from mock comparisons to observed data (see Section 4.4).
in Sections 4.1 and 4.2, we would expect R H to be constant as well, but we still need to verify that and determine its value.With a straightforward reparameterization of our results, we can investigate how R H may vary with sightline parameters.To determine R H from our framework, we again use Fitzpatrick (1999) to construct transformations between reddenings, here converting between SDSS bands and HST bands.We find We obtain A H to be used in Equation 5 by adjusting the A H as defined in 2MASS photometry as where we use β NIR = 1.85 as determined in section 4.1.
In Figure 11, we plot the relation between R H and E(u − z).With the exception of some large spikes at very low dust where we do not constrain it well, R H is essentially invariant with the dust column.We find the mean R H = 0.345 ± 0.007 across the last 10 (low-latitude) bins, whose E(V − I) range corresponds to that of the MW Cepheid sample presented in Riess et al. (2021Riess et al. ( , 2022)).
We can also use the same mock analysis framework from Section 4.2 to investigate potential variation of R H with R V .The left panel of Figure 12 indicates that, for individual sightlines (E(B − V ) > 0.25), there is quite a steep correlation between R H and R V .However, we must again test whether the correlation is intrinsic by comparing to a mock dataset.The set of mock stars we use here is exactly the same as is shown in the middle panel of Figure 9.This time, we simply calculate R H using the method described above and plot those values vs. R V instead.Indeed, we see a similar correlation in steepness and direction as in the observed case.The observed trend of R H with R V again appears to be a result of uncertainties in measurement.
To establish an intrinsic scatter in R H , we propagate the intrinsic scatter in β NIR (Section 4.2) and obtain σ(R H ) = 0.10, depicted in the vertical direction of the 2D Gaussian shown in the right panel of Figure 12.

Recipe for Use
We recommend the following simple recipe by which to make use of our results in correcting NIR observations: 1. Produce an extinction curve (A λ ) based on some reddening (e.g., E(B − V ) from SFD or Schlafly & Finkbeiner 2011) using a standard curve of choice.This curve can be R V -dependent, as in CCM, Fitzpatrick (1999), andFitzpatrick et al. (2019), or the fixed R V = 3.1 curve of Fitzpatrick & Massa (2007).The Python packages extinction (Barbary 2016) and dust_extinction (Gordon 2022) may be useful for this, or CCM_UNRED and F99_UNRED in IDL.Note that CCM_UNRED may actually use O'Donnell (1994) parameters as default.
2. Beyond λ = 7500 Å, replace the standard curve A λ with where the relevant values of β NIR are given in Table 1.

Standard Extinction Laws
There is a large number of studies, especially in recent years, that have investigated the NIR extinction law.Nonetheless, the NIR extinction curves that are most commonly used still come from the general UV through NIR extinction curves, especially the ones which have gained traction over the years and are referred to as the "standard" Milky Way curves.We start by comparing our average value for the power-law slope of the NIR extinction curve, β NIR = 1.85 ± 0.01, to other NIR power-law slopes associated with these standard curves.
The earliest standard curve still widely used is from CCM, which assumes a shallower curve than the one we find: β NIR = 1.61, originally found by Rieke & Lebofsky (1985).Rieke & Lebofsky (1985) investigated the extinction law in bands J through N (1 to 10 µm) toward a handful of M supergiants in the Galactic Center (GC), as well as toward o Sco (an A5 II star at b = 18 • ) and Cyg OB2 #12 (a B supergiant within an OB association in the Galactic plane).They found no variance in the NIR law between these sightlines.So even though CCM curves are R V -dependent in the UV and optical, they are universal beyond 0.9 µm.
The parameterization of another widely used law, from Fitzpatrick (1999), does not implement a strict power law in the NIR regime; rather, it makes use of fixed cubic spline anchor points at 1.2 and 2.7 µm, which we find to also be close to β NIR = 1.6 and thus shallower than our typical slope.The sightlines on which the law is based are described in Fitzpatrick & Massa (1990).Their program sample consists of 45 OB stars within ∼20 • of the Galactic plane, measured in the NIR regime in bands I through M (out to ∼5 µm).Again, although this law has R V dependence in the UV and optical, it is essentially fixed in the NIR.Fitzpatrick & Massa (2007) use 328 near-main-sequence OB stars to revisit the UV through NIR extinction law, but using the synthetic stellar model method (Fitzpatrick & Massa 2005) rather than the empirical pair method used in previous studies.For >1 µm, however, they just assume β NIR = 1.84 from Martin & Whittet (1990), a value in very close agreement with our result.Martin & Whittet (1990) investigate both diffuse and highly reddened sightlines.The diffuse data come from previous compilations (Savage & Mathis 1979;Whittet 1988) of many sightlines generally in the Galactic plane, and the highly reddened sightlines are ρ Oph and the Tr 14/16 OB stars.The Martin & Whittet (1990) study finds their data to be consistent with a single-β NIR universal curve between I and at least 5 µm for both the diffuse and highly reddened regimes.
In the most recent update to their work, Fitzpatrick et al. (2019) produce a new UV through NIR parameterization based on a sample of 72 OB stars.They have spectrophotometric data up to 1 µm, beyond which 2MASS photometry is used.The extinction curves are derived using synthetic stellar models, now in the NIR as well.In that method, the observed SEDs of reddened near-main-sequence B (and late O) stars are modeled in order to simultaneously derive estimates of the physical properties of the star along with the shape of its extinction curve.Fitzpatrick et al. (2019) derive and present their extinction curves as a series of finely spaced points in the region where they have spectra, out to 1 µm.Beyond that, they present one point for each 2MASS band (JHK), and a point representing λ = ∞.One can then interpolate to gain a continuous extinction curve.In the NIR, we find that their curve can be approximated by a power-law with β NIR = 1.86 at R V = 3.1, again in excellent agreement with our result.
To conclude, two of the earlier standard laws (CCM; Fitzpatrick 1999) featured shallower NIR laws (β NIR ≈ 1.6), whereas the two more recent laws feature somewhat steeper laws (β NIR ≈ 1.9).Our results agree with the latter two.

NIR studies of moderately dusty sightlines
In the intervening period between the Rieke & Lebofsky (1985) study and the review on dust grains by Draine (2003), several different values for β NIR were found, all larger than 1.61, with an upper bound around 1.8 (e.g., Martin & Whittet 1990).Since the turn of the century, derived values of β NIR have generally crept even higher.However, relatively few NIR studies have focused their attention toward diffuse regions of the Galaxy; most investigate the Galactic Center or dense clouds (such as the Coalsack).Those studies tend to find large values of β NIR (see Section 5.1.3).Even so, some recent studies find diffuse sightlines to have relatively steep NIR slopes: 2.14 (Stead & Hoare 2009), 1.95 (Wang & Jiang 2014), and 2.07 (Wang & Chen 2019).
The Wang & Chen (2019) study is perhaps the most similar in methods to ours in the literature to date.They select 61,111 red clump stars along diffuse Galactic sightlines, collecting photometry from a plethora of surveys from the optical to the MIR (including SDSS and 2MASS).To derive intrinsic colors, they employ stellar parameters from APOGEE (Eisenstein et al. 2011) in a similar fashion to our method in Section 3.They include a second-order term in T eff , but not log g.Their β NIR = 2.07 is higher than ours with a small uncertainty (0.03), making it in disagreement with us.Indebetouw et al. (2005) presents a bit of an outlier as far as recent diffuse β NIR : 1.65.However, they do not explicitly compute this value themselves; rather, it is quoted in other papers that way as derived from color excess ratios.Stead & Hoare (2009) present the fact that determining β NIR in this way depends sensitively on the wavelengths chosen for the filters.Using isophotal wavelengths, they find β NIR = 1.81 for Indebetouw et al. ( 2005); using effective wavelengths instead, they find β NIR = 2.05.It is possible that this effect impacts β NIR values reported for other studies which calculate color excess ratios based entirely within JHK.Stead & Hoare (2009) suggest that the use of isophotal wavelengths could be the reason for shallower NIR slopes in earlier studies.Fritz et al. (2011), however, find that the use of effective wavelengths actually slightly overestimates β NIR .In our study, since we do not use color excess ratios to determine β NIR , we are much less sensitive to the choice of wavelength for JHK.Notably, Wang & Chen (2019) also derive β NIR using color excess ratios, so they may be susceptible to similar issues.
While not large in sample size, the recent study of Decleir et al. (2022) uses very detailed spectrophotometric data that allows precise determination of the NIR extinction law.The unweighted average of their 13 stars yields β NIR = 1.77, supporting the canonical value.

NIR Slope vs. the Dust Column
All of the "standard" laws discussed in Section 5.1.1 are derived from sightlines in the plane of the Galaxy.In that lowlatitude regime, we find β NIR = 1.84.Since in this study we are particularly interested in the less explored high-latitude regime, we determine a separate value for it: β NIR = 1.91.The results between the two regimes are in better than 2σ agreement, so we find it appropriate to claim there are no differences in the character of the NIR power-law extinction between the low-and high-latitude regimes.In other words, β NIR values from studies based solely on sightlines confined to the Galactic plane should also apply at high latitude, and our agreement with canonical value supports that.
Many of the most recent studies of NIR extinction focus on sightlines toward the Galactic Center.These studies find preferentially larger values of β NIR than those which are focused away from the GC.This could be indicative of different NIR extinction behavior associated with "extreme" dust.In that region, E(B − V ) ∼ 10, far beyond what we probe in this study (even for our most reddened sightlines).Vexingly, though, there has not been uniform agreement across studies as to the value of β NIR even toward the GC alone.For example, Rieke & Lebofsky (1985) pull most of their sightlines from the GC, but they find β NIR = 1.61, one of the smallest values in regular use.Modern studies are more consistent, with β NIR values ranging from ∼2 to >2.5 in that direction, depending on the study.Although the methods for constructing extinction curves are relatively well-understood, it would appear that something still lies hidden in the process which leads to a wide dispersion of results.
There are a few studies which investigate both the GC and other Galactic plane regions together, thus standardizing technique and allowing for an internally consistent comparison between NIR extinction at different Galactic longitudes.One of these, Zasowski et al. (2009), uses 2MASS NIR and Spitzer MIR photometry to investigate the IR extinction law across 150 • of longitude in the Galactic plane.They report their results in terms of color excess ratios, which can be converted to β NIR by assuming the underlying dependence is a power law.
Thus their E(H − 5.While there is evidence that the slope of the NIR extinction curve takes a canonical value everywhere in the sky except near the center of the Galaxy, some studies find steeper values elsewhere in the plane, not allowing us to reach a final conclusion as to where the change between steep and less steep laws occurs.

NIR Slope vs. RV
In Section 4.2, we show that the apparent correlation we find between β NIR and R V can be explained by covariant uncertainties in those parameters, and that the intrinsic variation in β NIR has a standard deviation of only about 0.1.In other words, we confirm a universal NIR law in the sense that there are no R Vdependent trends as there are in the UV and optical regimes.
The strongest evidence for an R V dependence so far was found by Fitzpatrick & Massa (2009).They find that, for a selection of OB stars, β NIR does vary with R V .The orange points in the right panel of Figure 9 are taken from that study, with associated fit to the trend from Salim & Narayanan (2020, their Fig.4D).The result from Fitzpatrick & Massa (2009) is in stark contrast to those presented/assumed in previous parameterizations of the optical+NIR extinction curve (e.g., CCM, Fitzpatrick 1999;Fitzpatrick & Massa 2007), which provide extinction laws dependent on the single parameter R V , except in the NIR, where the curve becomes invariant.Though only 14 sightlines are studied in Fitzpatrick & Massa (2009), they show a relatively robust correlation between β NIR and R V .They suggest that it would be "surprising" for extinction to follow a constant power-law everywhere across a large wavelength range.Fitzpatrick & Massa (2009) attribute the ability to see the correlation to a large (intrinsic) range of R V values in their sample.Fitzpatrick et al. (2019) provide an updated R V -dependent Milky Way extinction curve in which β NIR is no longer assumed to be mostly constant as was the case in CCM and Fitzpatrick (1999).The resulting relationship between β NIR and R V (shown in the right panel of Figure 9) is nonetheless quite a bit shallower than the one from Fitzpatrick & Massa (2009), only spanning a range of ∼0.5 in β NIR on 2 < R V < 6.Such a weak dependence may not be detectable in our analysis, where we find the R V range is quite small.This relationship and its origin are not specifically discussed in the study, but are inherent in the presented R V -dependent curves.
In a recent study, Decleir et al. (2022) use NIR spectra to construct extinction curves along 15 sightlines via the pair method.Their sample is similar to that of Gordon et al. (2021), including 5 of the same stars.The sightlines are quite dust-heavy compared to ours, with two dense enough to show an ice extinction feature.Decleir et al. (2022) use IRTF/SpeX spectra covering the NIR and find a moderate negative correlation between R V and β NIR (ρ = −0.68).In the right panel of Figure 9, we show the diffuse sightlines from Decleir et al. (2022) in purple (excluding two with large R V uncertainties).They follow a similar relation to the Fitzpatrick & Massa (2009, orange).
Analysis of our sample reveals that the R V range is quite small, so any dependence of β NIR on R V would be relatively subtle.Nonetheless we test, again using the mock analysis, different degrees of linear correlation (ρ) between β NIR and R V .We find that ρ departing from zero yields a poorer match between the mock and actual measurements.From this we conclude that, for the general, low-latitude sightlines we have in our sample of individual curves, there is no significant correlation between β NIR and R V .Our findings contrast the studies discussed above that use "representative" samples, which select sightlines with a wide range of values in parameter space (especially R V ).Within these samples, it would appear that the inclusion of sightlines with R V significantly beyond 4 may lend itself to the finding of a correlation between R V and β NIR .
There have been a few studies that explore the relationship between the NIR extinction law and R V using statistical samples.One such study, Schlafly et al. (2016), derived the extinction law in the Galactic plane using stars with APOGEE spectra and optical (Pan-STARRS1) and NIR (2MASS and WISE) photometry.They find no correlation between E(y −J)/E(J −H) or E(J − H)/E(H − K) and R V (their Figure 12, bottom panel).Like us, they find R V to have a small range (σ(R V ) = 0.18), so sightlines that may potentially show a non-universal NIR law are uncommon.A similar conclusion of a lack of dependence between the NIR extinction curve slope and R V can be gleaned from the Nataf et al. (2016) study of extinction near the GC.Their Figure 8 shows that E(J − K s )/E(I − J) does not depend on R I = A I /E(V − I), which should correlate with R V .
The landscape of results pertaining to a dependence of NIR slope on R V is confusing given that some studies do and others do not find a dependence.Of the two that do (Fitzpatrick & Massa 2009;Decleir et al. 2022), the key seems to be having sightlines with R V > 4, which seem to be rare in general samples of stars like ours.One the other hand, according to Martin & Whittet (1990), the ρ Oph molecular cloud has R V = 4.3, but its NIR slope (1.85 ± 0.09) agrees with the β NIR found in the same study for diffuse R V = 3.1 sightlines (1.84 ± 0.03).It should be pointed out, however, that ρ Oph stars have 1.1 < E(B − V ) < 2.9, where as the two stars with R V > 4 in Decleir et al. (2022) have E(B − V ) ∼ 0.5.Similarly, the four stars with R V > 4 in Fitzpatrick & Massa (2009) also have lower reddening (0.5 < E(B − V ) < 0.9) than ρ Oph.Could it be that in the case of ρ Oph, high R V coupled with high dust column somehow conspire to yield a "normal" β NIR ?
Both Fitzpatrick & Massa (2009) and Decleir et al. (2022) are based on very small samples and Fitzpatrick & Massa (2009) may include some problematic stars (Maíz Apellániz & Barbá 2018).The hope is that the sample of stars for which β NIR and R V can be determined well enough to not be affected by covariances will be significantly expanded in the future so that the question can be definitively answered.

σ(R V )
Much effort has been put into the determination of R V in the literature since CCM codified it as the single value upon which the extinction curve in the UV and optical relies.The value determined in that study, R V = 3.1, is the one still generally applied today.Many subsequent studies (e.g., Fitzpatrick 1999;Fitzpatrick & Massa 2007;Schlafly & Finkbeiner 2011) have found general agreement with this value, either by direct measurement of R V or by comparing their extinction values to those predicted by curves of different R V .As noted above, we adopt R V = 3.1 for the generation of the mock stars in our analysis.We find no evidence to contradict that value from our data.
Of greater interest to our study is the intrinsic spread in R V , σ(R V ).While R V = 3.1 is characteristic of the Milky Way average extinction curve, it is important to know the chance that any given sightline may vary from that value significantly.We derive the spread using a mock analysis in two different ways (see Section 4.2), yielding the similar results σ = 0.28 and 0.24.Fitzpatrick & Massa (2007) find σ = 0.27 from their sample of 328 extinction curves along sightlines to OB stars.Schlafly et al. (2016) find a narrower spread with σ = 0.18, using 37,000 stars from the APOGEE survey largely confined to the Galactic plane (roughly the same width as found earlier by Schlafly et al. 2010).The characteristics of the distributions found by Fitzpatrick & Massa (2007) and Schlafly et al. (2016) are somewhat different.Fitzpatrick & Massa (2007) found quite a strong high-R V tail to their distribution, such that the Gaussian which yielded their spread was only fit to a core region around R V = 3.0; the tail dragged the mean of the entire sample up to R V = 3.22.Schlafly et al. (2016) find an almost symmetric Gaussian distribution about R V = 3.32, with a very small tail; they posit that the difference could be due to the different populations of stars probed (OB stars vs. background giants).They further note that many high-R V stars in Fitzpatrick & Massa (2007) are relatively poor fits, and claim that there may be an even smaller number of high-R V sightlines in reality (i.e., essentially no tail).We assumed a Gaussian distribution for intrinsic R V in our mock analysis partially due to this finding.In the future, it may be useful to carry out a more detailed analysis which allows for tailed and/or non-Gaussian distributions of intrinsic R V .

Cepheid Calibration
In Section 4.4, we find a new value for the NIR extinction parameter used to correct Cepheid variable magnitudes: R H = 0.345±0.007.Our value is somewhat lower than the Fitzpatrick (1999) R H = 0.386 used in Riess et al. (2016Riess et al. ( , 2021Riess et al. ( , 2022)).In Riess et al. (2022) an additional analysis was performed to determine, rather than assume, the value of R H by searching for the value that improves the overall fit the best, yielding R H = 0.363 ± 0.038 for the MW Cepheid sample.This value comes somewhat closer to what we find (but is still consistent with the Fitzpatrick (1999) value as well).
In the right panel of Figure 12, we show the relation between R H and R V which comes directly from the literature extinction laws CCM (purple), Fitzpatrick (1999, orange), andFitzpatrick et al. (2019, green).For Fitzpatrick (1999), there is almost no dependence of R H on R V .We also show a cyan point at R H = 0.336 for Fitzpatrick & Massa (2007), in which R V = 3.1.At R V = 3.1, Fitzpatrick et al. (2019) yields R H = 0.411.Our derived intrinsic values (gray contours) are in general agreement with these standard curves in the relevant R V range.Riess et al. (2022) find the relation ∆(5 log H 0 ) = 0.08∆R H , and therefore ∆H 0 /H 0 = 0.04∆R H . Adjusting from their baseline value of H 0 = 73.04km s −1 Mpc −1 , our R H = 0.345 would yield H 0 = 73.16km s −1 Mpc −1 , a change of ∼0.1σ.
We note that the recent study Mörtsell et al. (2022) finds a very strong change in the fitted value of H 0 with increasingly stringent cuts on the reddening E(V − I) of the Cepheids included in their sample.Given that we find no correlation between R H and E(u−z) (or by extension E(V −I)), our results would not support such a change in the final H 0 value.

Conclusions
In this work, we used ∼50,000 stars with SDSS spectra and SDSS+2MASS photometry to investigate the nature of the Milky Way extinction curve in the near-infrared (out to K band).We examine whether the curve is universal or varies with the amount of dust or with R V .
1. We report an overall mean power-law slope of the nearinfrared extinction curve β NIR = 1.85 ± 0.01.We recommend this β NIR for correcting the magnitudes of extragalactic objects, MW Cepheids, and other objects with low to moderate reddening (0.02 < E(B − V ) ≲ 1).
2. We find no correlation between β NIR and the dust column over a wide range of dust column densities, from very low ones at high latitudes (E(B − V ) ∼ 0.02) to moderately high ones at lower latitudes (E(B − V ) ∼ 0.8).However, we do not explore highly reddened regions such as the Galactic Center, where some studies find steeper curves (β NIR ≳ 2).
3. Our β NIR confirms the NIR law used in Fitzpatrick & Massa (2007) and Fitzpatrick et al. (2019).It is however steeper than in either Cardelli et al. (1989) or Fitzpatrick (1999) (β NIR ∼ 1.6), so we do not recommend those two for the NIR.We do not determine which standard law may be favored at shorter wavelengths.We propose replacing the λ > 7500 Å portion of whichever law one uses with 4. We find no strong evidence for a correlation between R V and the dust column, although there is weak evidence for a rise beyond R V = 3.1 below E(B − V ) ≈ 0.1.

5.
There is an apparent correlation between β NIR and R V (i.e., between NIR and optical slopes, since R −1 V = A B /A V − 1), which we entirely explain (using a mock analysis) as arising from covariant uncertainties in the two parameters.We conclude that an intrinsic correlation between β NIR and R V must be small or nonexistent in our sightlines.Some recent studies that do find a correlation include sightlines with R V > 4, which seem very rare in our general sample.
6.While the apparent star-to-star scatter in β NIR is quite large, we find (again using the mock analysis) that this is the result of observational uncertainties, and that the intrinsic star-to-star scatter in β NIR is relatively small (σ(β NIR ) = 0.13).
7. With a similar mock analysis, we derive a relatively small intrinsic spread in R V of σ(R V ) = 0.24, confirming the results of Schlafly et al. (2016) and other large-sample studies.
8. We derive a new value for the reddening parameter used in Cepheid calibrations in the measurement of H 0 : R H = 0.345 ± 0.007.Like β NIR , R H does not appear to intrinsically depend on the amount of dust or on R V .
While we have attempted to not only present our results but also arrive at a consistent picture, the landscape of literature results regarding NIR extinction is somewhat difficult to navigate.One thing that may help in the future is to consider using mock analysis and forward modeling to test whether any correlation between extinction parameters is robust.

Figure 1 .
Figure 1.Distribution of observed SDSS r magnitudes for the 56,649 stars used in our final program sample for the main analysis.The dashed red vertical line is the median, r = 15.8.

Figure 2 .
Figure 2. A 2D histogram of the absolute value of Galactic latitude vs. E(u − z), a proxy for dust.All of the 56,649 stars used in our final program sample for the main analysis are included, with the exception of a small number which have very low E(u − z).The dashed vertical line indicates the approximate separation between high-latitude (low-dust; left) and low-latitude (high-dust; right) regimes.We also provide a secondary x-axis to show the more familiar E(B − V ) value, determined via a simple ratio from E(u − z) (see Section 3.3).

Figure 3 .
Figure 3. 2D histogram of observed SDSS g − r color vs. E(B − V ) from Schlegel et al. (1998) for SEGUE2 F9 stars in our original sample.The stars to the left of the red line indicating E(B − V ) < 0.0125 (essentially unreddened stars) were used to determine the intrinsic color calibration by fitting to Equation 1.

Figure 4 .
Figure 4. Distribution of SDSS spectral subclasses (SUBCLASS) for the 56,649 stars used in our final program sample for the main analysis.The stacked colored bars indicate the relative numbers of stars in each spectral subclass which were measured by the SDSS (Legacy), SEGUE1, and SEGUE2 spectroscopic campaigns.

Figure 5 .
Figure 5.Standard deviation of the residuals vs. color for calibrator stars (those with E(B − V ) SFD < 0.0125).For each color, all survey/spectral subclass combinations in our original sample with at least 10 stars at E(u − z) are shown.The residuals come from regression to Equation 1.

Figure 7 .
Figure 7. E(λ − i) vs. λ −1 in 22 bins of E(u − z).The magenta dashed curves represent best fits to the two-parameter model described by Equation 2. The gray curves represent a standard Fitzpatrick (1999) extinction law with R V = 3.1 based on the mean E(u − z) of the bins.Error bars on the points are rescaled such that χ 2 red ≡ 1.The central wavelengths for each of the four filters incorporated in the fits are labeled along the x-axis.Fitted β NIR and A i values are shown near the top of each panel, above the E(λ − i) = 0 dotted line.

Figure 8 .
Figure 8. Power-law slope of the NIR extinction curve (β NIR ) vs. a proxy for the dust column (E(u − z)) for the fits described in Section 3.3.The purple line and associated shaded error region shows the values determined from the 22-bin aggregated fits as shown in Figure 7.The orange line shows individual stars, which were then binned after fitting for β NIR and A i .Two horizontal dotted lines are included, showing β NIR values from impactful studies: 1.84, from Fitzpatrick & Massa (2007) (originally from Martin & Whittet 1990), and 1.61, from Cardelli et al. (1989) (originally from Rieke & Lebofsky 1985).Note that Fitzpatrick (1999) find β NIR ≈ 1.62.The upper secondary x-axis is E(B − V ) = 0.26E(u − z).The arrows indicate an approximate separation between the high-and low-latitude regimes in our study, to either side of E(B− V ) ∼ 0.1.See Figure 2.

Figure 9 .
Figure9.Power-law slope of the NIR extinction curve (β NIR ) vs. a representation of the optical slope of the extinction curve (R V ) for three cases.In each case, standard literature values of R V = 3.1 (e.g.,Cardelli et al. 1989) and β NIR = 1.84(Martin & Whittet 1990; Fitzpatrick & Massa 2007)  are included as brown dashed lines.Left: Observed individual stars with E(u − z) > 1. Center: Analogous to the left panel, but using a comparably sized sample of mock stars generated from the Fitzpatrick (1999) family of curves, which are consistent with a fixed β NIR .See Section 4.2 for details.Right: Selected literature values for sightlines with measured β NIR and R V .The orange points and associated fit come fromFitzpatrick & Massa (2009).The purple points come fromDecleir et al. (2022).The gray shaded ellipses represent 1σ, 2σ, and 3σ contours of a two-dimensional Gaussian with σ(R V ) = 0.24 and σ(β NIR ) = 0.13.These are intrinsic values derived from mock comparisons to observed data in the tails of the distribution (see Section 4.2).

Figure 10
Figure10.R V ≡ A V /E(B − V ) vs. a proxy for the dust column (E(u − z)) using the same bins and fits as Figure8.The dotted line indicates the accepted MW value R V = 3.1 (e.g.,Cardelli et al. 1989).The upper secondary xaxis is E(B − V ) = 0.26E(u − z).The arrows indicate an approximate separation between the high-and low-latitude regimes in our study, to either side of E(B − V ) ∼ 0.1.See Figure2.
8)/E(H − K) values yield β NIR = 2.2 for |l| < 25 • and β NIR = 1.8 for |l| ∼ 70 • .The latter result stands in excellent agreement with our results, whereas the former result agrees equally well with the most recent estimate for the GC inner 3 × 3 deg 2 extinction law from Sanders et al. (2022, also β NIR = 2.2).Although the dependence of color excess ratio on longitude in Zasowski et al. (2009) unarguably leads to the conclusion of a different law near and away from the central region of the Galaxy, the more recent study of Maíz Apellániz et al. (2020) find the same β NIR = 2.26 near the Galactic Center and other regions of the plane.

Table 1 .
Values for the slope of the NIR extinction curve.