SpeX near-infrared spectroscopic extinction curves in the Milky Way

Interstellar dust extinction curves provide valuable information about dust properties, including the composition and size of the dust grains, and are essential to correct observations for the effects of interstellar dust. In this work, we measure a representative sample of near-infrared (NIR; 0.8-5.5 $\mu$m) spectroscopic extinction curves for the first time, enabling us to investigate the extinction at wavelengths where it is usually only measured in broad photometric bands. We use IRTF/SpeX spectra of a sample of reddened and comparison stars to measure 15 extinction curves with the pair method. Our sample spans A(V) values from 0.78 to 5.65 and R(V) values from 2.43 to 5.33. We confirm that the NIR extinction curves are well fit by a power law, with indices and amplitudes differing from sight line to sight line. Our average diffuse NIR extinction curve can be represented by a single power law with index $\alpha = 1.7$, but because of the sight line-to-sight line variations, the shape of any average curve will depend on the parental sample. We find that most of the variation in our sample can be linked to the ratio of total-to-selective extinction R(V), a rough measurement of the average dust grain size. Two sight lines in our sample clearly show the ice extinction feature at 3 $\mu$m, which can be fitted by a modified Drude profile. We find tentative ice detections with slightly over 3$\sigma$ significance in two other sight lines. In our average diffuse extinction curve, we measure a 3$\sigma$ upper limit of A(ice)/A(V) = 0.0021 for this ice feature.


INTRODUCTION
Interstellar dust plays a significant role in several physical and chemical processes in the interstellar medium (ISM).The dust grains regulate the temperature of the gas, catalyze the formation of molecular hydrogen, and provide a reservoir of heavy elements.Furthermore, dust absorbs and scatters starlight at ultravi-Corresponding author: Marjorie Decleir mdecleir@stsci.eduolet (UV), optical, near-and mid-infrared (NIR-MIR) wavelengths, and re-emits this energy at infrared (IR) wavelengths.This reprocessing of electromagnetic radiation highly affects our ability to study the universe at these wavelengths, and alters the observed spectral energy distribution (SED) of celestial objects.It is therefore of utmost importance to understand the interplay between dust and starlight.
The extinction of starlight by dust, i.e. the combined effect of absorption and scattering, is wavelength dependent and can be described by an extinction curve.Extinction curves are not only required to correct the observations of a range of astrophysical objects for the effects of dust, but also provide insights into the properties of the interstellar dust: the continuum shape contains details about the size distribution of the grains, while extinction features reveal the composition of the dust grains.Over the past decades, several techniques have been developed to measure extinction curves, which can broadly be categorized in two types: those using individual sightlines, and those using an ensemble of stars.Both types are complimentary and come with their own advantages and disadvantages.
A commonly used method of the first type is the pair method (Stecher 1965;Massa et al. 1983): the observed SED of a reddened star is compared to (i.e.divided by) the observed SED of an unreddened comparison star with similar stellar properties.The difference in their SEDs is then attributed to the dust in the lineof-sight towards the reddened star.This method has been used extensively to measure extinction curves in samples of diffuse and dense sightlines in the Milky Way (e.g., Rieke & Lebofsky 1985;Fitzpatrick & Massa 1986, 1988;Cardelli et al. 1989;Martin & Whittet 1990;Clayton et al. 2003;Valencic et al. 2004;Gordon et al. 2009Gordon et al. , 2021)), as well as in the Magellanic Clouds (e.g., Gordon et al. 2003).The main benefit of using comparison stars is that it does not require any modeling or absolute flux calibrations.One caveat is that the extinction measurements are relative (to the comparison star), and often only a limited number of comparison stars are available which makes it more difficult to exactly match the stellar properties of the reddened star.It is also possible to use stellar atmosphere models instead of observed comparison stars with which to compare the SED (or the colors) of the reddened star.This has been done in the Milky Way (by e.g., Fitzpatrick & Massa 2005, 2007, 2009;Fitzpatrick et al. 2019;Massa et al. 2020), in M31 (Clayton et al. 2015), and in the Small Magellanic Cloud (Maíz Apellániz & Rubio 2012).This method makes it easier to closely match the properties of the reddened star, but its accuracy relies on a precise absolute flux calibration and accurate stellar atmosphere models.
The second technique to measure extinction uses an ensemble of stars, rather than individual sightlines.Nishiyama et al. (2006) and Alonso-García et al. (2017), for example, used the positions of red clump (RC) stars in color-magnitude diagrams (CMD) as a tracer of the extinction and reddening toward the Galactic center (GC).This technique is referred to as the RC method, and a variant was used by Nishiyama et al. (2009) to determine the extinction curve toward the GC using magnitudes of bulge RC stars and colors of red giant branch (RGB) stars.Nogueras-Lara et al. (2018, 2019) explored different methods based on RC stars to measure the extinction toward the GC, including using a combination of stellar atmosphere models and an extinction grid, using a fixed extinction, using a color-color diagram (CCD), and using a CMD.Indebetouw et al. (2005), Stead & Hoare (2009) and Maíz Apellániz et al. (2020) also used the position of RC stars in a CCD to measure the extinction in the Galactic plane.Finally, Fritz et al. (2011) derived the extinction curve toward the GC using nebular hydrogen emission lines.
The ensemble methods listed above are very valuable to measure extinction toward the GC in the NIR, especially given the large number of stars and strong extinction effects.However, it is not possible to measure UV and optical extinction, because the GC cannot be detected at those wavelengths due to its very high extinction.To obtain a detailed understanding of the interstellar dust properties (e.g.size, composition), it is critical to combine multi-wavelength extinction measurements from UV to MIR.Furthermore, dust depletions, which are measured from UV absorption line spectra of bright nearby stars, provide unique insights into the dust properties.Since UV-optical extinction measurements and depletions can only be measured in the local ISM, we focus on the local ISM (i.e.within 3 kpc) in this paper.
All of the above-mentioned IR studies (except for Fritz et al. (2011)) are based on a limited number of broadband photometric data points (usually a subset of the IJHKLM bands), in contrast to most of the UV studies listed above and the recent optical measurements of Fitzpatrick et al. (2019) and Massa et al. (2020) which use spectra.As, for example, suggested by Maíz Apellániz et al. (2020), there is a limit to what can be studied about the extinction curve with photometry.Recently, Gordon et al. (2021) characterized the MIR extinction curve based on Spitzer photometry (3.6-24 µm) and spectra (5-37 µm).They showed that the average diffuse Milky Way extinction curve at these wavelengths can be represented by a power law with index α = 1.48 (and two modified Drude profiles for the silicate features).However, they also found large variations in the shape of the extinction curve between different sightlines, proving that the IR extinction is not uniform within the Galaxy.
With the work presented in this paper, we fill the gap between the spectral UV/optical extinction studies and the recent MIR results from Gordon et al. (2021), by measuring extinction curves at wavelengths between 0.8 and 5.5 µm, for the first time using NIR spectra (instead of broadband data).We utilize SpeX spectra for a sample of 25 reddened stars and 15 comparison stars, and measure the extinction towards the reddened stars with the pair method.With this data set, we are not only able to characterize the shape of an average diffuse Milky Way NIR extinction curve at spectroscopic resolution, but also to study the variations between the different sightlines in our sample.
Section 2 describes the sample, and the processing of the SpeX spectra.In Section 3 we explain how we measured extinction curves from these data.The fitting of the extinction curves is outlined in Section 4. Section 5 presents and discusses the results of the fitting, the correlation between the fitting parameters, the average diffuse Milky Way extinction curve with a comparison to other studies and dust grain models, the correlation of the sightline variations with R(V ), and the observed extinction features.Finally, Section 6 summarizes this work.

Sample
Our sample consists of 15 comparison and 25 reddened Milky Way OB stars.They are listed in Table 1 with their spectral type, B and V-band photometry, and distance, all obtained from the literature (see references in the table).OB stars are particularly suited for determining extinction curves, because their spectra exhibit fewer stellar lines compared to later type stars.Furthermore, these stars are luminous at UV-MIR wavelengths, which are the wavelengths of interest to study dust extinction.The comparison stars were selected to have very little dust along their line-of-sight, as can be seen from their estimated reddening, E(B −V ), listed in Table 1.To obtain these reddening values, we calculated the observed (B − V )-color for every comparison star from its B and V-band photometry (Table 1), and subtracted the intrinsic (B − V )-color from Table 1 in Fitzgerald (1970) for the corresponding spectral type and luminosity class (also given in our Table 1).Note that these values are not used in our analysis, but just given as a reference. 1he sample of reddened sightlines was chosen so that it represents a large range in R(V ) values (2.4-5.3), which probe the dust grain size along the line-of-sight.The E(B − V ) values range from 0.2 to 1.6, and A(V ) values from 0.8 to 5.6.2

SpeX NIR spectra
The NIR spectra have been obtained during an observational campaign at the 3.2 m NASA Infrared Telescope Facility (IRTF) on Mauna Kea.Observations were obtained with SpeX, a medium-resolution 0.8 5.5 µm spectrograph, in two modes: the Short wavelengths Cross-Dispersed mode (SXD, ∼0.8 2.4 µm), and the Long wavelengths Cross-Dispersed mode (LXD, ∼1.9 5.5 µm) (Rayner et al. 2003).The observations have been performed over several nights in 2005, 2006 and 2007.In addition, we used 7 spectra from the IRTF spectral library, observed over multiple years, starting in 2000 and ending in 2003.
Table 1.Sample stars with their spectral type, B and V-band photometry, distance, and corresponding references.For the comparison stars, also an estimate of their reddening, E(B − V ), is given. [mag] [mag] ref.
[pc] ref.Note-Stars with an * were used to measure extinction curves, as explained in Section 3.
We reduced and processed the spectra with Spextool (SPectral EXtraction TOOL) (v.4.1), an IDL-based data reduction package written by Cushing et al. (2004).The data reduction started with extracting the spectra from the raw data, and the flat fielding and wavelength calibration of these spectra, using the xspextool GUI.Subsequently, multiple spectra of the same star were combined using the xcombspec GUI.Next, the stellar spectrum was corrected for the telluric absorption and the instrument throughput with the xtellcor program (Vacca et al. 2003).For this step, observations of a standard A0V star, and a high resolution model of Vega were used.After this, the different orders of the corrected spectrum were scaled and merged using the task xmergeorders, which results in one SXD and one LXD spectrum per star.It has to be noted that some wavelength regions of the spectra still suffered significantly from telluric absorption effects, despite the corrections applied with Spextool.These regions are shaded in red in Fig. 1, and were masked and excluded from all further analyses (such as measuring and fitting the extinction curves).
All spectra were then placed on the same wavelength grid (between 0.8 and 2.45 µm for SXD, and between 1.9 and 5.5 µm for LXD) with a resolution of 2000, using the measure extinction python package (Gordon & Decleir 2021).Finally, with the same package, the spectra were calibrated based on photometric data points.First, the SXD spectrum was scaled to match J, H, K S photometry from 2MASS (Skrutskie et al. 2006), obtained from the IRSA 2MASS All-Sky Point Source Catalog (Skrutskie et al. 2019).Subsequently, the LXD part of the spectrum was scaled to align with the SXD spectrum, resulting in one smooth NIR spectrum. 3The final spectra of all stars are electronically available (Decleir 2022a), and can be found in Figs.2-4.They are plotted multiplied by λ4 to remove the strongly decreasing Rayleigh-Jeans tail in this wavelength range.When plotted in this way, the spectrum of an unreddened (comparison) star flattens towards higher wavelengths.
The flux uncertainties that Spextool yields only include the photon noise and read noise.However, there are several other sources of uncertainty that need to be taken into account.For example, scaling and merging the different orders in the spectrum introduces some uncertainty, especially in those wavelength regions suffer-ing from telluric absorption where it is difficult to measure the scaling factor based on the overlap of two adjacent orders.In addition, scaling the SXD spectrum based on the 2MASS photometry, and manually scaling the LXD spectrum to match the SXD spectrum, induces more uncertainty.Since it is very hard to quantify these uncertainties, we quadratically added a 1% uncertainty to the photon and read noise to account for these uncertainties.
Finally, we excluded data points with a signal-to-noise ratio (SNR) below 10 from all further analyses, limiting the wavelength range to ∼5.2 µm for most stars.The median SNRs of the stellar spectra used in this work (indicated with an * in Table 1) are in the range 107-141 below 2.5 µm, 38-137 between 2.8 and 4 µm, and 18-72 beyond 4.5 µm, at a resolution of about 2000.

MEASURING NIR EXTINCTION CURVES
As mentioned in the introduction, we measured NIR extinction curves using the pair method: The spectrum of a reddened star was divided by the spectrum of an unreddened comparison star with similar stellar properties.We want to point out that it is not possible to use stellar atmosphere models to measure NIR extinction curves with the pair method.The NLTE (non-local thermodynamic equilibrium) TLUSTY models (Lanz & Hubeny 2003), which are needed for OB stars, are not complete beyond 0.8 µm.Many stellar lines (including the upper Paschen lines) are missing (Massa et al. 2020), and also the continuum level is not reliable (priv.comm.I. Hubeny).Usually, reddened and comparison stars are matched based on their spectral type, and if possible also on their luminosity class.The spectral types and luminosity classes listed in Table 1 were taken from SIMBAD (Wenger et al. 2000).It has to be noted, however, that for several stars, multiple spectral types are listed on SIMBAD (obtained from different references).Furthermore, spectral types might be different when derived using different methods or from different wavelength regions (see e.g., Hanson et al. 1996;Smith Neubig & Bruhweiler 1997).For these reasons, we did not use the literature spectral types in Table 1 to match comparison and reddened stars, but merely list them as a reference. 4Instead, we validated empirically which comparison star is the best match.We tried all comparison stars for every reddened star in our sample and retained the one that results in the smoothest extinction curve, i.e. that cancels out the stellar (hydro-  2004) with the atmospheric transmission tool ATRAN (Lord 1992).The red shading indicates wavelength regions where the atmospheric transmission is very low, and the spectra are significantly affected by the telluric absorption.These regions were masked in the spectra.
gen) lines and jumps as well as possible.Tables 2 and 3 list the comparison star used to measure the extinction curve for every reddened star.Given the limited number of available comparison stars (only 15), the spectral match is not always perfect, and remaining stellar lines and jumps might be visible in some extinction curves (see Fig. 5).
The extinction curves were calculated with the measure extinction package (Gordon & Decleir 2021) using the following procedure.The absolute dust extinction A(λ) at a wavelength λ is given by: where F is the spectral flux density, d is the distance to the star, "red" refers to the reddened star, and "comp" refers to the comparison star.Extinction is usually measured relative to a reference wavelength measurement in order to avoid the need of an accurate distance to the star.Most commonly, the V-band measurement is used as the reference.Thus, the relative dust extinction E(λ−V ) (also called differential extinction, dust reddening or color excess) at wavelength λ can be calculated as: The V-band magnitudes of all stars are given in Table 1.
In order to compare between extinction curves from different sightlines, they must be normalized to the total level of extinction in that sightline, e.g.represented by A(V ), the absolute extinction in the V-band.The differential measurement E(λ − V ) can be converted into an absolute normalized extinction measurement A(λ)/A(V ) from Eq. 2: This conversion requires knowledge of A(V ), which we measured by fitting a power law to the NIR extinction curve, as discussed in the next section.Once A(V ) is known, one can also compute the total-to-selective extinction R(V ): with the color excess E(B − V ) calculated as ) with m(B) the apparent B-band magnitude and m(V ) the apparent V-band magnitude listed in Table 1.The obtained values for E(B − V ) are given in Tables 2 and  1 2 3 4 5 [ m] 3. 5 R(V ) probes the dust grain size along the lineof-sight, with larger values corresponding to sightlines dominated by larger dust grains (e.g., Fitzpatrick 1999).On average, R(V ) = 3.1 in the Milky Way diffuse ISM (Cardelli et al. 1989).
The normalized extinction curves are shown in Fig. 5, and are electronically available (Decleir 2022a).We were able to measure extinction curves for 15 reddened stars of our sample (with spectra shown in Fig. 3).The other 10 reddened stars were not suitable to measure extinction curves, due to peculiarities in their spectrum (plotted in Fig. 4).Four of these stars (HD014422, HD034921, HD206773 and HD052721) show strong emission lines in their spectrum (as also indicated in their spectral type with the letter "e").HD166734 and HD294264 have a strongly rising spectrum towards longer wavelengths.Emission lines and/or very steeply rising spectra are most likely signatures of a stellar wind or circumstellar disk.HD034921, HD206773 and HD166734 have also been classified as stars with clear wind signatures by Gordon et al. (2021).As an additional check, we used a similar approach as in Fig. 4 of Gordon et al. (2021) to separate stars with wind signatures, based on their IR colors.Fig. 6 shows K S −WISE 4 versus J − K S for the reddened and comparison stars in our sample that have WISE photometry (Wright et al. 2010), obtained from the IRSA AllWISE Source Catalog (Wright et al. 2019).The J − K S color is a measure of the reddening and high K S −WISE 4 values ( 1) indicate stellar wind signatures.The 4 remaining stars (HD037020, HD037023, HD037022 and HD014250) have either an "n" (indicating nebulous absorption) or a "p" (indicating another unspecified peculiarity) in their spectral type.All these peculiarities are unique to the specific star and complicate the extinction measurement using the pair method, because a comparison star with the same characteristics (e.g.same wind emission) would be needed.We thus excluded these sightlines from all further analyses.We indicated with an * in Table 1 which stars we used to measure extinction curves.These are plotted as green diamonds in Fig. 6, and do not show signatures of stellar winds.
The median SNRs of the normalized extinction curves are in the range 2.7-27 below 2.5 µm, 0.4-5.8 between 5 Several previous extinction studies pointed out that using monochromatic quantities to characterize an extinction curve is more appropriate than using the band-integrated equivalents we use here.However, as explained in App.A of Maíz Apellániz & Barbá (2018), and as visible in Fig. 3 of Maíz Apellániz ( 2013) and Table 4 of Fitzpatrick et al. (2019), these effects are very small for low-extinction OB stars, which we use in our sample.
2.8 and 4 µm, and 0.3-1.5 beyond 4.5 µm, at a resolution of about 2000.We would like to note here that some of the 2MASS photometric data points in Fig. 5 seem to deviate slightly from the spectral extinction curves.This is likely due to saturation issues in the J-band for some of the (bright) stars in our sample, resulting in a lower quality (as flagged in the 2MASS catalog) and hence in larger error bars on the magnitudes, and consequently on the extinctions in those bands.Given that all three bands (J, H and K S ) were used to scale the (SXD) spectrum (as discussed in Section 2.2), a small offset in one of the bands could slightly shift the spectrum up or down, causing the other bands to look slightly deviant from the spectrum as well.However, because we are using all three bands for the scaling, an issue with one band does not significantly impact the extinction measurements and fitting results.Furthermore, all JHK S photometric data points (for all stars) are within 3σ from the spectrum.Finally, when available, we tested with other JHK photometry, but this did not significantly affect the final fitting results.Hence, we decided to stick with the 2MASS photometry for all sightlines to be as uniform as possible.

Continuum extinction
In previous works, the continuum NIR extinction was usually described with a power law (see e.g., Draine 1989;Cardelli et al. 1989;Martin & Whittet 1990), as it provides a convenient analytical representation of the extinction curve.We followed their example and fitted the continuum extinction between 0.8 and 5.5 µm with a power law.From Eq. 2, the measured differential extinction can be written as: where where S is the amplitude (by definition equal to the normalized extinction at 1 µm) and α is the index of the power law.At infinite wavelengths, Eq. 7 reduces to E(∞ − V ) = −A(V ).In other words, fitting the observed extinction curve with the function in Eqs.7 and 8, gives us a direct measurement of A(V ).We used a combination of the PowerLaw1D Astropy (Astropy Collaboration et al. 2013) model (for Eq. 8), and the AxAvToExv conversion function from the dust extinction python package (Gordon et al. 2022) to implement the conversion in Eq. 7.
The fitting was done in two steps.First, the Levenberg-Marquardt algorithm was used to obtain prelimi-2 3 4 5 [ m] nary fit results for the three parameters S, α and A(V ), using the LevMarLSQFitter from Astropy.These results were then used as initial guesses for the Markov chain Monte Carlo (MCMC) fitting based on the Emcee python tool (Foreman-Mackey et al. 2013).We used 6 walkers, each with 9000 steps after a burn in of 1000 steps to sample the parameter space.We experimented with more walkers, more steps and larger burn fractions, but the results were the same.The inverse of the squared uncertainties on the extinction were used as weights in the fitting, so that wavelength regions with larger uncertainties contribute less to the fit.The fitted curves are shown as red dashed lines in Fig. 5.The uncertainties on the fitted parameters given by the MCMC fitting only include random noise, computed as the difference between the 84 th and 50 th percentile (upper uncertainty), or between the 50 th and 16 th percentile (lower uncertainty) of the posterior distribution function.However, there is also a systematic uncertainty related to the fact that the extinction curves are with and without strong stellar winds.The red circles represent reddened stars that were not used to measure extinction curves (see Fig. 4), because they have winds.The purple square corresponds to HD014250, which was not used because it has other peculiarities in its spectrum (as explained in Section 3).The green diamonds represent the reddened stars that were used to measure extinction curves.They do not show signatures of stellar winds.The green diamond above the line corresponds to HD029647, which is not windy (as confirmed by Gordon et al. (2021)), but its WISE photometry possibly suffers from extended source and/or scattered moonlight contamination.
measured relative to the V-band extinction.To assess the effect of this uncertainty, we first combined the uncertainties on the V-band magnitude of the reddened and the comparison star by adding them quadratically and taking the square root.Subsequently, we ran the MCMC fitting three times: using the measured extinction, using the extinction subtracted by the combined V-band uncertainty, and using the extinction summed with the combined V-band uncertainty.For all parameters, we then computed the systematic uncertainty as half of the difference between the maximum and minimum of the three fitted values.Finally, we added this systematic uncertainty in quadrature to the random uncertainty, and took the square root to obtain the total uncertainty on each fitted value.
The code that was used to perform the fitting, compute the uncertainties, and analyze and plot the results, is available as part of the spex nir extinction package on GitHub (Decleir 2022b).

Extinction features
Strong extinction features can influence the shape of the fitted power law, because the procedure described above will try to fit the features and the continuum together with one power law.Therefore, for sightlines with strong features, it is important to explicitly consider these features in the fitting by, e.g., adding a Drude profile to the power law, and fit the data with the combined function, as such constraining the free parameters of the power law and the Drude function simultaneously.
From a visual inspection of the extinction curves, only two sightlines in our sample seem to have an obvious extinction feature around 3 µm: HD283809 and HD029647, as can be seen in Fig. 7.This feature is caused by water ice along the line-of-sight, as will be discussed in more detail in Section 5.4.Since it is generally accepted that the detection of ice indicates the presence of a dense molecular component along the line-of-sight (e.g., Whittet et al. 1997;Boogert et al. 2015), we consider these two sightlines as "dense".Note that we did not use the total V-band extinction A(V ) to divide our sample into diffuse and dense sightlines, as it is known that a high A(V ) can also be measured in sightlines without any dense material (e.g., in the diffuse sightline towards Cyg OB2 no. 12 which has A(V ) ∼ 10, Whittet (2015)).The extinction curve of HD283809 also shows a weaker bump around 3.4 µm, which is likely caused by hydrocarbons, as will be discussed in Section 5.4.This feature is, however, not obvious in the extinction curve towards HD029647.There also seems to be a small systematic residual around 2.3 µm for HD283809, which is not visible for HD029647.However, it is possible that this residual is due to the transition between the SXD and LXD spectrum around 2.4 µm.
One approach to study these extinction features would be to compare them with laboratory measurements, as was e.g.done by Thi et al. (2006) and Boogert et al. (2011).However, since one of the goals of this work is to obtain a functional form for the NIR extinction curve, we opted to directly fit functional profiles to our data, without making any assumptions on the material causing the features.
It has to be noted here that it is also possible to constrain an extinction feature directly from the spectrum, as has been done in the past (e.g., Pendleton & Allamandola 2002;Thi et al. 2006).However, this requires knowledge of the underlying continuum flux, which is often approximated by a straight line between two points left and right from the feature.Similarly, fitting extinction features from a continuum-subtracted extinction curve requires assumptions on the underlying continuum extinction and the extent of the features, which  3.
introduces uncertainties on the measurement of the extinction features (e.g., Boogert et al. 2011).Therefore, we have opted to fit the extinction continuum and features simultaneously in the extinction curves.
The strong telluric absorption (see Fig. 1) significantly increases the noise between 2.5 and 3.5 µm, which makes it hard to obtain any meaningful constraints on the weaker feature around 3.4 µm.Furthermore, we found that the impact of this weak bump on the continuum fit is negligible.We thus only considered the strong ice feature in the fitting.We fitted the extinction curve of both dense sightlines with a combination of a power law and a modified Drude profile D(λ): with with B the strength, λ 0 the central wavelength, and γ the width of the profile.Because of the asymmetric appearance of the feature, we modified the standard Drude function to allow for extra asymmetry. 6The modified profile is obtained by replacing the standard width in Eq. 10 by a width γ(λ) that depends on the wavelength: with γ 0 the standard width, and a the asymmetry parameter.This modification is based on the work of Stancik & Brauns (2008), who introduced asymmetric Gaussian and Lorentzian profiles with widths varying across the profiles to model infrared absorption profiles.It was also successfully used by Gordon et al. (2021) to fit the 10 and 20 µm silicate features in MIR extinction curves.
The fitting results (i.e. the 50 th percentiles of the posterior distribution functions) and their uncertainties (including both random and systematic uncertainties, as explained in Section 4.1) for the two dense sightlines are given in Table 3.To verify whether an ice feature is present in the extinction curves of the other sightlines, we also fitted them with a combination of a power law and a modified Drude profile as described above, but with the central wavelength, width and asymmetry fixed to the average values of the two dense sightlines.For only 4 sightlines, the fitted strength of the Drude profile is larger than 3 times its uncertainty (i.e. a > 3σ detection): HD038087, HD156247, HD183143 and HD229238.However, an investigation of the extinction curves and residuals of HD038087 and HD156247, shows that these sightlines significantly suffer from telluric absorption around 3 µm.In sightlines HD183143 and HD229238, there is a very tentative detection of a weak ice feature with strengths B = 0.0026 ± 0.0008 (3.3σ) and B = 0.0042 ± 0.0012 (3.5σ) respectively (as compared to a > 50σ detection for the two dense sightlines).These sightlines have the highest A(V ) values of the diffuse sample (A(V ) = 3.86 and A(V ) = 2.99, respectively).Given the non-detection (or very tentative detection) of ice in all sightlines other than HD283809 and HD029647, these can be considered "diffuse" for the purpose of this paper.The fitting results and their uncertainties for the 13 diffuse sightlines (without fitting any features) are given in Table 2.

Extinction curve parameter trends
The fitting results in Tables 2 and 3 show a wide range in A(V ) (0.78-5.65) and R(V ) values (2.43-5.33).As a cross-check, in Fig. 8 we compared our obtained A(V ) values with values reported in the literature for the same sightlines.A(V ) values were collected from Cardelli et al. (1989), Valencic et al. (2004), Gordon et al. (2009) and Gordon et al. (2021), which all used the pair method to measure extinction curves.Cardelli et al. (1989) fitted the Rieke & Lebofsky (1985) extinction curve to the R, I, J, H, K, and L color excesses E(λ−V ) with A(V ) as the only free parameter.Valencic et al. (2004) estimated R(V ) from the J, H, and K color excesses as described in Fitzpatrick (1999), and calculated A(V ) using R(V ) and E(B − V ). Gordon et al. (2009) extrapolated the J, H, and K E(λ − V ) curves to infinite wavelength using the Rieke et al. (1989) IR extinction curve to derive A(V ) values.Finally, Gordon et al. (2021) derived A(V ) by fitting a combination of a power law and two modified Drude functions (for the silicate features) to their measured NIR/MIR E(λ − V ) extinction curves.The agreement between our A(V ) values and the literature values is very good for all sightlines, even though they were derived using different methods and different data.
It is interesting to see if there are any correlations between the different parameters describing the extinction curve: the amplitude S and the index α of the power law, the V-band extinction A(V ) and the total-to-selective extinction R(V ).Fig. 9 shows scatter plots of these parameters.We found a clear anti-correlation (with Spearman's rank correlation coefficient ρ = −0.91) between the amplitude and the index of the power law.This is straightforward to explain considering the fact that the amplitude is equal to the normalized extinction at 1 µm, and that the extinction curves are normalized to A(V ).A steeper extinction curve (i.e.larger α) implies a faster decrease in extinction at wavelengths beyond the V-band compared to a flatter extinction curve, resulting in a lower normalized extinction at 1 µm (i.e.smaller S).R(V ) correlates with the amplitude of the power law (ρ = 0.84), and anti-correlates with the power law index (ρ = −0.68).This shows that R(V ) is linked to the slope of the NIR extinction curve, with higher R(V ) values corresponding to flatter curves (i.e.smaller α and larger S), as is also the case in the UV and op-  2 and 3) to values reported in the literature for the same sightlines (taken from Cardelli et al. 1989, Valencic et al. 2004, Gordon et al. 2009and Gordon et al. 2021.).There is a good agreement for all sightlines.
tical (e.g., Cardelli et al. 1989;Fitzpatrick 1999).The R(V )-dependence of the shape of the extinction curve is investigated and discussed in more detail in Section 5.3.

Measurement and fitting
We averaged the extinction curves of the 13 diffuse sightlines in our sample.The average was only computed in those wavelength regions where at least 5 curves have data, as such excluding the wavelengths with limited and noisy extinction data (e.g.beyond 5 µm).The average NIR curve is plotted in the right panel of Fig. 10, together with some average Milky Way extinction curves from the literature (Rieke & Lebofsky 1985;Martin & Whittet 1990;Indebetouw et al. 2005;Gordon et al. 2021).Our average curve is very close to that from Martin & Whittet (1990), while the other literature curves are slightly flatter.Note that all these average curves are based on a different sample of sightlines and some differences are to be expected.We also compared with average curves from Rieke et al. (1989) and Chiar & Tielens (2006), which are not shown in the plot because they are very similar to the other literature curves.A more detailed comparison to previous studies is given in Section 5.2.2.Given that our average curve is close to most diffuse literature curves, we can Note-The measured extinction values in this table are relative to the comparison star.The reported uncertainties include both random and systematic uncertainties, as explained in Section 4.1.
Table 3. MCMC fitting results for the 2 dense extinction curves.conclude that our sample of 13 diffuse sightlines provides a reasonable representation of the average diffuse Milky Way extinction.As an additional check, we also measured the average UV extinction curve for our sample of sightlines using International Ultraviolet Explorer (IUE) spectra (taken from Valencic et al. 2004).The average UV curve for our sample is plotted in the left panel of Fig. 10, and is very similar in shape to the Cardelli et al. (1989) and Fitzpatrick et al. (2019) average R(V ) = 3.1 diffuse UV curves, hence confirming that our sample of 13 diffuse sightlines is representative of the average diffuse Milky Way extinction.Furthermore, we find R(V ) = 3.12 ± 0.05 for our average curve, which is consistent with the average diffuse Milky Way value (Cardelli et al. 1989).The measured average NIR extinction curve is electronically available (Decleir 2022a), and a binned version of the curve and its uncertainty is given in Table 4 (columns 2 and 3) for wavelengths between 0.  2 and 3).The magenta squares indicate the dense sightlines, the black circles indicate the diffuse sightlines, and the red star shows the average diffuse extinction curve (see Section 5.2.1).The Spearman's rank correlation coefficient is shown in every plot.
Table 4. Average diffuse Milky Way extinction curve and parameters of the linear relationship between extinction A(λ)/A(V ) and 1/R(V ).This table is also available as part of the D22 MWAvg model in the dust extinction python package (Gordon et al. 2022).
band average extinction R(V )-dependent ext.
or measured fit       We also fitted the average curve with a power law and obtained:7 The fitted parameters are given in Table 2, and indicated with a red star in Fig. 9.The fitted model is shown in red in Fig. 10, and listed at specific wavelengths and in photometric bands in Table 4 (column 4).The fitted model is also available as part of the D22 MWAvg model in the dust extinction python package (Gordon et al. 2022).
The residuals of the fitting (i.e. the data subtracted by the fit) are plotted in Fig. 11, and are small (< 0.02) at most wavelengths.The horizontal magenta lines represent the sigma-clipped standard deviation of the residuals in different wavelength regions: 0.0050 (0.8-1.4 µm), 0.0032 (1.4-1.8 µm), 0.0030 (1.9-2.5 µm), 0.0047 (2.8-4 µm) and 0.016 (4.5-5.0 µm).These values can be considered as the upper limits on the strengths of any features in the average extinction curve.At several specific wavelengths hydrogen lines are visible in the residuals, due to small spectral mismatches between the reddened stars and their comparison star.These are indicated with vertical blue lines.We also indicated peaks in the residuals caused by the Paschen (Pa) and Brackett (Br) jump mismatches.At the top, the same atmospheric transmission model is plotted as in the residuals can most likely be attributed to the telluric absorption.
To verify whether the 3 µm water ice feature is present in the average diffuse extinction curve, we also fitted it with a combination of a power law and a modified Drude profile with fixed central wavelength, width and asymmetry.The resulting strength of the Drude profile is B = 0.0019 ± 0.0007.This is below the 3σ detection threshold of A(ice)/A(V ) = 0.0021.However, the strength of the feature is very close to what Potapov et al. (2021) have found for the diffuse sightline Cyg OB2 no. 12, which adds weight to the possible detection of ice in our diffuse average extinction curve.

Comparison to other studies and dust models
Several studies in the past have used a variety of methods to measure NIR extinction curves, and have targeted different Galactic environments.Furthermore, they all had a different way of presenting the extinction measurements, e.g. by giving color excess ratios, reporting absolute or relative extinctions, adopting a different normalization, using various photometric systems, etc.This complicates a comparison between different studies.One thing most of these studies seem to have in common, is that they approximated the NIR extinction curve by a power law, i.e.A(λ) ∼ λ −α .Therefore, in this section, we use the power law index α to compare our result with previous measurements.Table 5 compares the fitted power law index for our average diffuse NIR extinction curve to indices found in the literature.
Our value for α is in reasonable agreement with that from Draine (1989), and between the results of Cardelli et al. (1989) and Martin & Whittet (1990).Later measurements toward the GC have a larger power law index compared to our result.It has to be noted that those measurements probe highly extinguished (A(V ) 8) regions in the inner Galactic disk and bulge, at about 8 kpc distance.On the contrary, we measure extinction in the local (closer than 3 kpc, see distances in Table 1), low-extinction (A(V ) < 3.9), diffuse ISM.The very long sightlines toward the distant GC have a higher probability of encountering molecular clouds along the lineof-sight compared to our local measurements.Those sightlines thus likely contain a mix of diffuse and dense material, which can have an effect on the shape of the extinction curve.The presence of dense material towards the GC is confirmed by the detection of ice features (see e.g., Fritz et al. 2011), which are likely only (clearly) visible in sightlines through molecular clouds and not in purely diffuse extinction regions (Whittet et al. 1997, and see our discussion in Section 4.2).On the other hand, the detection of the strong aliphatic hydrocarbon feature at 3.4 µm shows that those sightlines also contain diffuse dust (Fritz et al. 2011).We are confident that our average curve probes purely diffuse extinction, given that no strong ice feature has been detected and that our average UV extinction curve is close to diffuse curves from the literature (see Fig. 10).Fritz et al. (2011) also suggested that it is likely that there is a transition of a mostly flatter NIR extinction curve in the solar neighborhood to a steeper one in most parts of the Galactic disk.
We want to emphasize that we find a quite large range of indices within our sample of diffuse sightlines (between 1.36 and 2.20, see Table 2), showing that there are real sightline-to-sightline variations in the shape of the NIR extinction curve within the local diffuse ISM.Nishiyama et al. (2006Nishiyama et al. ( , 2009)), Fritz et al. (2011), andAlonso-García et al. (2017) also found that the NIR extinction curve changes from one sightline to another.As a consequence, the shape of any average extinction curve will depend on the sightlines that were used to measure that average, and it is not surprising that different studies, using different samples, find different average curves.As can be seen from Fig. 9, and as briefly discussed in Section 5.1, the variation in α seems to be anti-correlated with R(V ).A more detailed discussion on the R(V )-dependence of the NIR extinction curve is given in Section 5.3.
As mentioned in the introduction, some studies toward the GC (and Galactic plane) reported a dependence of the power law index on the wavelength region.For example, Indebetouw et al. (2005) and Nishiyama et al. (2009) saw a flattening of the extinction curve at wavelengths beyond ∼ 3 µm.Nogueras-Lara et al. ( 2019) measured a smaller index α in HK S than in JH (see Table 5).Fritz et al. (2011) obtained a flatter extinction curve beyond ∼ 3.7 µm.They also found a smaller α between 2.2 and 2.8 µm than between 1.2 and 2.2 µm (see Table 5), but stated that this change in slope is not significant.We demonstrate that our local average diffuse extinction curve between 0.8 and 4 µm can be represented by a single power law, since there are no long-range trends visible in the residuals in Fig. 11.We argue that it is possible that the mix of diffuse and dense material in sightlines toward the GC and in the Galactic plane imposes the need for multiple power law indices across the NIR.
In Fig. 12, we compare our average diffuse NIR extinction curve to the Draine (2003a,b), Zubko et al. (2004), Compiègne et al. (2011) and Jones et al. (2013) diffuse ISM dust grain models.The Zubko et al. (2004) model seems to agree with our average curve below 1 µm, while the other models correspond better beyond 3 µm.

The R(V)-dependent NIR extinction curve
Although the average extinction curve described in the previous section can be very useful to get an idea of the overall shape of the diffuse NIR extinction curve, we found clear variations around this average for different sightlines, as can be seen from the fitting results in Table 2 and in Fig. 9.Because of these strong sightline-to-  Draine (1989) 1.75 multiple a 0.9-5 µm (I, J, H, K, L, M) Galactic a Cardelli et al. (1989) 1.61 pair method 0.9-3.3µm (I, J, H, K, L) local ISM Martin & Whittet (1990) 1.84 ± 0.03 pair method 0.9-4.8µm (I, J, H, K, L, M) local diffuse ISM Nishiyama et al. (2006) 1.99 ± 0.02 RC CMD 1.2-2.2µm (J, H, K S ) Galactic center Nishiyama et al. (2009) 2.0 RC+RGB CMD 1.2-2.2µm (J, H, K S ) Galactic center Stead & Hoare (2009) 2.14 sightline variations, the shape of any average extinction curve will depend on the sample of sightlines used to measure that average, so one should be cautious when using an average extinction curve.Interestingly, several previous studies have shown that these variations are largely dependent upon a single parameter, often chosen to be R(V ).Cardelli et al. (1989) presented a linear relationship between the extinction curve and 1/R(V ), based on photometric extinction measurements between 0.1 and 3.4 µm.For the NIR region, they found that the slopes and intercepts of this linear relationship can be described by a power law.More recently, Fitzpatrick et al. (2019) presented a linear relationship between the extinction curve and R(V ) based on a combination of spectral UV and optical extinction measurements and 2MASS JHK S measurements.
In this section, we investigate our sample of observed NIR extinction curves as a function of 1/R(V ).8 Fig. 13 shows A(λ)/A(V ) vs. 1/R(V )−1/3.1 at 5 different wavelengths (chosen in regions where the atmospheric absorption is not significantly affecting the extinction measurement).We indeed see a linear trend at all wavelengths.It has to be noted, however, that at the longest wavelength (4.7 µm, bottom plot), the range in R(V ) values is limited due to the limited number of sightlines with data at wavelengths beyond 4 µm.
To quantify this correlation, using the LinearLSQFitter from Astropy, we fitted a linear function (Linear1D Astropy model) at every wavelength in the SpeX range where at least 5 sightlines have data (∼2800 data points between 0.8 and 5.1 µm): By choosing 1/R(V ) − 1/3.1 as the abscissa, the intercepts a(λ) correspond to the average R(V ) = 3.1 Milky Way extinction curve.The slopes b(λ) illustrate the R(V )-dependence.The fits are shown as green lines in Fig. 13.We also added the R(V )-dependent relationships from Cardelli et al. (1989) and Fitzpatrick et al. (2019) (where available) for comparison.At λ = 0.9 µm these have a similar slope to our fitted line, while at longer wavelengths our fit is steeper than the literature lines, which means that we find a stronger dependence on R(V ).We also calculated the standard deviation of the individual sightlines from the fitted relationship, which reflects both measurement uncertainties and real deviations from the R(V )-dependent relationship.Fig. 14 shows the obtained intercepts, slopes and standard deviations as a function of wavelength.It is clear that the R(V )-dependence of the extinction slowly decreases towards longer wavelengths, but it does not entirely disappear within our wavelength range.We want to note here that the fact that A(λ)/A(V ) and 1/R(V ) = E(B − V )/A(V ) have a common factor, namely A(V ), could possibly lead to an artificial correlation between both quantities and between their uncertainties.However, we found that the correlation coefficients between both quantities, due to their common factor, are small, ranging from -0.02 to 0.12, with a median value over all sightlines and all NIR wavelengths of 0.02.More importantly, the correlation coefficients are positive for most wavelengths and most sightlines.In contrast, the observed relationship between both quantities is a negative linear trend (i.e. an anti-correlation) as b(λ) < 0 for all NIR wavelengths (see Figs. 13 and 14), and could thus not be caused by the (small) positive artificial correlation due to their common factor A(V ).9 In order to provide a smooth representation of the R(V )-dependence for use in future extinction studies, using the LevMarLSQFitter from Astropy, we fitted the intercepts with a power law (PowerLaw1D Astropy model) and found: As mentioned above, the intercepts correspond to the average R(V ) = 3.1 Milky Way extinction curve.Hence, it is not surprising that they are well fitted with a power law function, similar to the average curve described in Section 5.2.1.For the slopes and standard deviations, no straightforward functional form could be fitted.Instead, we interpolated the data with cubic splines using scipy (Virtanen et al. 2020).We first binned the data into 25 equally-sized wavelength bins (i.e.every bin has approximately 110 data points), and calculated the median wavelength, median slope, and median standard deviation in every bin.These median values were then used as anchor points for the spline interpolation.The data points above 4 µm (plotted in gray in Fig. 14) are very noisy and were not included in the interpolation.The fitted intercepts, and interpolated slopes and standard deviations are shown as red dashed lines in Fig. 14, and are tabulated in Table 4 (columns 5-7) for wavelengths between 0.8 and 4 µm, and in common IR photometric bands (2MASS J, H, K S , WISE 1, L, and IRAC 1), obtained by integrating the fitted/interpolated curves over the response curves of the bands.It has to be noted here that the broad "wiggles" in the slopes are most likely a result of the telluric absorption.Our R(V )dependent model is also available as the D22 model in the dust extinction python package (Gordon et al. 2022).Fig. 15 shows the derived R(V )-dependent NIR extinction curve for 3 values of R(V ) (2.5, 3.1 and 5.5).The shaded regions represent the standard deviation about the curve (last column of Table 4).For comparison, we added the R(V )-dependent curves from Cardelli et al. (1989) and Fitzpatrick et al. (2019), which mostly fall within the shaded regions for R(V ) = 3.1 and R(V ) = 5.5.However, for R(V ) = 2.5 the literature curves deviate slightly from our curve for wavelengths 1 µm, as was also clear from Fig. 13.It has been suggested in the past that the R(V )dependence of NIR extinction curves is not real, but simply a result of the normalization of the extinction curves to A(V ), and disappears when normalizing at longer wavelengths (see e.g., Cardelli et al. 1989).We tested this by normalizing the curves to the extinction at 1 µm, and plotting the extinction as a function of 1/R(V ) − 1/3.1 as before.We still found a strong linear correlation between extinction and 1/R(V ).
We conclude that the R(V )-dependent extinction curve can reduce uncertainties when applying extinction corrections to sightlines with known R(V ) values, because it takes into account sightline variations which are correlated with R(V ).Real variations around the R(V )dependent curve are still present, but smaller than the variations around the average diffuse extinction curve derived in the previous section.We thus recommend using our R(V )-dependent extinction curve when possible.
This work completes a series of spectroscopic extinction curve studies from the far-UV to the MIR (Gordon et al. 2009;Fitzpatrick et al. 2019;Gordon et al. 2021).In future work, we will combine these results to create an R(V )-dependent Milky Way extinction curve from 912 Å to 30 µm (Gordon et al., in prep.).One focus of this future work will be to investigate the impact of underlying correlations between A(λ)/A(V ) and 1/R(V ) (and between their uncertainties) on the R(V ) relationship.Preliminary results show that such correlations mainly impact wavelengths shorter than 1 µm, confirming that the observed relationship in this work (in the NIR) is not significantly affected by underlying correlations.

Extinction features
Several NIR-MIR observations towards dense interstellar clouds have revealed features in the extinction curves that are caused by ices.Boogert et al. (2015) provide a detailed review on the detected ice features in infrared spectra.The strongest ice feature has a peak wavelength around 3 µm, and is primarily caused by the O-H stretching modes of bulk H 2 O ice.While a detailed study of this feature is beyond the scope of this work, we briefly discuss its characteristics in our sample.
The 3 µm ice feature was only clearly detected in the two dense sightlines, HD283809 and HD029647, (and very tentatively in HD183143 and HD229238).Both HD283809 and HD029647 lie toward the Taurus Dark Cloud region, close to the dense clump TMC-1, and the presence of ice features in their spectra has been previ- ously reported (e.g., Whittet et al. 1988;Smith et al. 1993).We show in Section 4.2 that the feature can be fitted reasonably well with a modified Drude profile.Fig. 7 shows that the feature has a long-wavelength wing.As explained by Boogert et al. (2015), large grains could be responsible for part but not all of the wing.The remainder of the absorption may be attributed to the O-H stretching mode in ammonia hydrates (Knacke et al. 1982;Dartois & d'Hendecourt 2001).
Structure near 3.4 µm might arise from the C-H stretching mode in hydrocarbons (Pendleton & Allamandola 2002), which have been detected along diffuse sightlines towards a few background stars (Pendleton et al. 1994).If so, it is possible that our approach to fit the ice feature masks this weaker hydrocarbon feature which could arise from the diffuse dust along these sightlines.
It has to be noted that the telluric absorption (see Fig. 1) severely complicates the measurement of any fea-tures in the NIR extinction curve.It is, for example, very difficult to constrain the underlying continuum extinction and the shape of the feature, because of the gap between 2.5 and 2.88 µm.Also, the data between 2.88 and 3.5 µm are very noisy.Fortunately, our approved cycle 1 GO program 245910 and other programs on the James Webb Space Telescope will provide continuous, well calibrated, high signal-to-noise spectra of dust extinguished Milky Way OB stars, which will enable a much more detailed study of dust extinction features in the NIR (and MIR).spectrograph between 0.8 and 5.5 µm for a sample of 25 reddened and 15 comparison stars.We were able to measure extinction curves for 15 sightlines using the pair method.This sample spans A(V ) values from 0.78 to 5.65 and R(V ) values from 2.43 to 5.33.The main conclusions of this work are: • The NIR extinction curves in our sample are well fitted by a single power law over the entire wavelength range between 0.8 and 5.2 µm, with indices and amplitudes strongly differing from sightline to sightline.
• Our average local diffuse NIR Milky Way extinction curve can be represented by a single power law with index α = 1.7 over the entire wavelength range between 0.8 and 5 µm, and agrees well with the average curve from Martin & Whittet (1990).
• The shape of any average extinction curve depends on the sample of sightlines used to measure the average, and one should thus be cautious when using an average extinction curve.
• We find that the normalized extinction curves A(λ)/A(V ) in our sample linearly depend on the selective-to-total extinction 1/R(V ), with the strength of the dependence decreasing with increasing wavelength.
• Two sightlines in our sample (HD283809 and HD029647) show a strong ice extinction feature around 3 µm, which can be approximated by a modified Drude profile.These sightlines most likely contain dense material.
• We tentatively detect (with slightly over 3σ significance) the 3 µm ice feature in the extinction curves for HD183143 and HD229238.These sightlines have the highest A(V ) values of the diffuse sample.
• We do not detect the 3 µm ice feature in the average diffuse extinction curve, with a 3σ limit of A(ice)/A(V ) = 0.0021.
• The telluric atmospheric absorption complicates the characterization of any extinction features, as well as the continuum extinction above 4 µm.
Planned cycle 1 observations with the James Webb Space Telescope will provide significant improvements to our understanding of the NIR dust extinction.
The SpeX spectra and measured extinction curves are electronically available (Decleir 2022a).The code used for the analysis and plots is available on GitHub (Decleir 2022b).

Figure 1 .
Figure1.An atmospheric transmission model, obtained from the Spextool database, and computed byCushing et al. (2004) with the atmospheric transmission tool ATRAN(Lord 1992).The red shading indicates wavelength regions where the atmospheric transmission is very low, and the spectra are significantly affected by the telluric absorption.These regions were masked in the spectra.

Figure 2 .
Figure 2. The NIR spectra of the 15 comparison stars, ordered by luminosity class (giants and supergiants at the top, and main sequence stars at the bottom), and by spectral class (O4-B8 from top to bottom).The SpeX spectra are shown as solid lines and the photometric data (JHKS) are shown as open circles.All fluxes are normalized to the flux at 1 µm, and multiplied by λ 4 to remove the strongly decreasing Rayleigh-Jeans tail.

Figure 3 .Figure 4 .
Figure 3.The NIR spectra of the 15 reddened stars that were used to measure extinction curves, ordered by A(V ).The SpeX spectra are shown as solid lines and the photometric data (JHKS) are shown as open circles.All fluxes are normalized to the flux at 1 µm, and multiplied by λ 4 to remove the strongly decreasing Rayleigh-Jeans tail.

Figure 5 .Figure 6 .
Figure 5.The normalized NIR extinction curves for the 15 suitable sightlines in our sample.The SpeX spectra are shown as solid lines and the photometric data (JHKS) are shown as open circles.The red dashed lines represent the fits to the data.

Figure 7 .
Figure7.The extinction curves of the sightlines towards HD283809 (left) and HD029647 (right) show a strong ice feature around 3 µm.The continuum extinction and the ice feature were fitted simultaneously with a linear combination of a power law and a modified Drude profile.The fit parameters are given in Table3.

Figure 8 .
Figure8.A comparison between our obtained A(V ) values from the fitting (see Tables2 and 3) to values reported in the literature for the same sightlines (taken fromCardelli et al. 1989, Valencic et al. 2004, Gordon et al. 2009 and Gordon  et al. 2021.).There is a good agreement for all sightlines.

Figure 9 .
Figure9.Scatter plots of the different parameters describing the extinction curve: the amplitude S and the index α of the power law, the V-band extinction A(V ) and the total-to-selective extinction R(V ) (see Tables2 and 3).The magenta squares indicate the dense sightlines, the black circles indicate the diffuse sightlines, and the red star shows the average diffuse extinction curve (see Section 5.2.1).The Spearman's rank correlation coefficient is shown in every plot.

Figure 13 .
Figure13.Illustration of the R(V )-dependence of the extinction curve at a handful of wavelengths.The magenta squares indicate the dense sightlines and the black circles the diffuse sightlines.The linear fits are shown as green lines.For comparison, the R(V )-dependent relationships fromCardelli et al. (1989) andFitzpatrick et al. (2019) are added as blue dashed and orange dotted lines, respectively.

Figure 14 .
Figure 14.Intercepts a (top), slopes b (middle) and standard deviations about the fit σ (bottom) of the linear 1/R(V )dependence, as a function of wavelength.The green crosses indicate the wavelengths plotted in Fig. 13.The red dashed line in the top panel shows the fitted power law to the intercepts.The red diamonds in the middle and bottom plot are the anchor points used in the spline interpolation, while the red dashed line represents the cubic spline interpolation.In these panels, data points above 4 µm are plotted in gray and were excluded from the interpolation because they are very noisy.

Figure 15
Figure 15.R(V )-dependent extinction curves for different R(V ) values: 2.5 in blue, 3.1 in orange and 5.5 in green, in comparison with the R(V )-dependent curves from Cardelli et al. (1989) in dashed lines and Fitzpatrick et al. (2019) in dotted lines.The shaded regions represent the standard deviation around the curve.

Table 2 .
MCMC fitting results for the 13 diffuse extinction curves and the average diffuse extinction curve.
a Obtained by fitting the average diffuse extinction curve (see Section 5.2.1).b Calculated as R

Table 5 .
Comparison of power law indices α from the literature.
This paper combined several extinction measurements from the literature, obtained with different methods and for different Galactic environments.bThis work used different methods based on RC stars.We report their average α value.