SpeX near-infrared spectroscopic extinction curves in the Milky Way

Marjorie Decleir

0000-0001-9462-5543]Marjorie Decleir
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD, 21218, USA

0000-0001-5340-6774]Karl D. Gordon
Visiting Astronomer at the Infrared Telescope Facility, which is operated by the University of Hawaii under contract 80HQTR19D0030 with the National Aeronautics and Space Administration.
Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD, 21218, USA
Sterrenkundig Observatorium, Universiteit Gent,
Gent, Belgium

0000-0003-0123-0062]Jennifer E. Andrews
Visiting Astronomer at the Infrared Telescope Facility, which is operated by the University of Hawaii under contract 80HQTR19D0030 with the National Aeronautics and Space Administration.
Gemini Observatory/NSF’s NOIRLab, 670 N. A’ohoku Place, Hilo, HI 96720, USA

0000-0002-0141-7436]Geoffrey C. Clayton
Visiting Astronomer at the Infrared Telescope Facility, which is operated by the University of Hawaii under contract 80HQTR19D0030 with the National Aeronautics and Space Administration.
Department of Physics and Astronomy, Louisiana State University, Baton Rouge, LA, 70803 USA

0000-0001-7780-3352]Michael C. Cushing
Ritter Astrophysical Research Center, Department of Physics and Astronomy, University of Toledo, 2801 W. Bancroft St., Toledo, OH 43606, USA

Karl A. Misselt
Visiting Astronomer at the Infrared Telescope Facility, which is operated by the University of Hawaii under contract 80HQTR19D0030 with the National Aeronautics and Space Administration.
Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA

0000-0001-8102-2903]Yvonne Pendleton
NASA Ames Research Center, Moffett Field, CA 94035, USA

John Rayner
Visiting Astronomer at the Infrared Telescope Facility, which is operated by the University of Hawaii under contract 80HQTR19D0030 with the National Aeronautics and Space Administration.
Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI, USA 96822

0000-0002-9123-0068]William D. Vacca
Visiting Astronomer at the Infrared Telescope Facility, which is operated by the University of Hawaii under contract 80HQTR19D0030 with the National Aeronautics and Space Administration.
SOFIA-USRA, NASA Ames Research Center, MS N232-12, Moffett Field, CA , USA

0000-0001-8539-3891]D. C. B. Whittet
Department of Physics, Applied Physics and Astronomy,
Rensselaer Polytechnic Institute, 110 Eighth Street, Troy, NY 12180, USA

Interstellar dust extinction curves provide valuable information about the dust properties, including the composition and size of the dust grains. In addition, extinction curves, and are essential to correct observations for the effects of interstellar dust. They are found to vary between different sightlines in the Milky Way. In this work, we measure a representative sample of near-infrared (NIR, 0.8–5.5 ) 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 25 reddened and 15 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.332.43 to 5.33. We confirm that the NIR extinction curves are well fit by a power law, with indices and amplitudes differing from sightline to sightline. TheOur average diffuse NIR extinction curve can be represented by a single power law with index $\alpha=1.7$, similar to previous works, but because of the sightline-to-sightline 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-selectionselective extinction $R(V)$, a rough measurement of the average dust grain sizes. Two sightlines in our sample clearly show the ice extinction feature at 3 , which can be approximatedfitted by a modified Drude profile. We find tentative ice detections with slightly over $3\sigma$ significance in two other sightlines. In theour average diffuse extinction curve, we did not detect this feature withmeasure a $3\sigma$ upper limit of 0.0021 $A(3\,\micron)/A(V)$$A(\text{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 a large amount of the starlight at ultraviolet (UV), optical,  and 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.

Extinction curves are commonly measured using theA 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 line-of-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 <cit.>, as well as in the Magellanic Clouds <cit.>. 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 <cit.>, 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, Nogueras-Lara et al., 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.

Studies of the Milky Way extinction curve at UV and optical wavelengths are omnipresent and have shown that the extinction curve (and thus the dust properties) is not universal, but varies from sightline to sightline <cit.>. However, there is a paucity of studies at longer wavelengths. For a very long timeHowever, for many years, it was believed that, while extinction curves vary significantly at UV and optical wavelengths in different regions of the Galaxy, they are rather uniform at longerNIR wavelengths $>$0.7  <cit.>. Several studies found that the IR extinction curve from $\sim0.7$$\sim0.9$ to $\sim5$  can be closely represented by a power law ($A(\lambda)\propto\lambda^{-\alpha}$), with reported values for the index $\alpha\simeq1.75$ [Draine, 1989], $\alpha\simeq1.6$ [Cardelli et al., 1989], and $\alpha\simeq1.8$ [Martin & Whittet, 1990] in the local ISM.
More recent extinction measurements toward the GC and in the Galactic disk showed steeper NIR extinction curves (between $\sim1.2$ and $\sim2.2$ ) with indices $\alpha\simeq1.99$ [Nishiyama et al., 2006], $\alpha\simeq2$ [Nishiyama et al., 2009], $\alpha\simeq2.14$ [Stead & Hoare, 2009], $\alpha\simeq2.11$ [Fritz et al., 2011], $\alpha\simeq2.47$ [Alonso-García et al., 2017], $\alpha\simeq2.30$ [Nogueras-Lara et al., 2018], $\alpha\simeq2.32$ [Nogueras-Lara et al., 2019], and $\alpha\simeq2.27$ [Maíz Apellániz et al., 2020].
Some of these GC studies reported a dependence of the power law index on the wavelength region. More specifically, [Indebetouw et al., 2005] and [Nishiyama et al., 2009] found a flattening of the extinction curve at wavelengths beyond $\sim3$ , [Fritz et al., 2011] obtained a flatter extinction curve beyond $\sim3.7$   and [Nogueras-Lara et al., 2019] found different power law indices between JH (1.27–1.65 , $\alpha\simeq2.43$) and HKS (1.65–2.16 , $\alpha\simeq2.23$).
[Nishiyama et al., 2006, Nishiyama et al., 2009], [Fritz et al., 2011], and [Alonso-García et al., 2017] found that the NIR extinction curve changes from one sightline to another, in contrast to what was previously believed. On the other hand, [Indebetouw et al., 2005], [Stead & Hoare, 2009], [Nogueras-Lara et al., 2018, Nogueras-Lara et al., 2019] and [Maíz Apellániz et al., 2020] saw no significant variation in the NIR extinction curve as a function of sightline in the Galactic regions they studied.

However, allAll of the above-mentioned IR studies (except for [Fritz et al., 2011]) are based on a limited number of broad-band photometric data points (usually a subset of the IJHKLM bands), in contrast to most of the UV/optical studies mentionedlisted above and the recent optical measurements of [Fitzpatrick et al., 2019] and [Massa et al., 2020] which use spectra. Furthermore, they report quite large uncertainties in the extinction measurements above 3 .
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 mid-infrared (MIR) extinction curve based on Spitzer photometry (3.6–24 ) and spectra (5–37 ). They showed that the average diffuse Milky Way extinction curve at these wavelengths can be represented by a power law with index $\alpha=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 , 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 <ref> describes the sample, and the processing of the SpeX spectra. In Section <ref> we explain how we measured extinction curves from these data. The fitting of the extinction curves is outlined in Section <ref>. Section <ref> 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 <ref> summarizes this work.

§ DATA PROCESSING

§.§ Sample

Our sample consists of 15 comparison and 25 reddened Milky Way OB stars. They are listed in Table <ref> 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 <ref>. To obtain these reddening values, we calculated the observed ($B-V$)-color for every comparison star from its B and V-band photometry (Table <ref>), 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 <ref>). Note that these values are not used in our analysis, but just given as a reference.[For star HD164794 (spectral type O4), no intrinsic ($B-V$)-color was available, but we did not use this star to measure extinction curves.] The sample of reddened sightlines was chosen so that it represents a large range in $R(V)$ values (2.32.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.[$R(V)$, $E(B-V)$ and $A(V)$ will be defined in Section <ref>.]

Sample stars with their spectral type and its reference., B and V-band photometry, distance, and corresponding references. For the comparison stars, also an estimate of their reddening, $E(B-V)$, is given.

star     type                type     B     V     BV     dist.     dist.     $\bm{E(B-V)}$

ref.     [mag]     [mag]     ref.     [pc]     ref.

9ccomparison stars

HD003360*          B2IV            [1]       $3.47\pm0.03$     $3.66\pm0.03$     [12]     109     [25]     0.05

HD031726*          B1V             [1]       $5.941\pm0.002$     $6.147\pm0.001$     [13]     389     [25]     0.05

HD032630          B3V             [1]       $3.00\pm0.03$     $3.18\pm0.02$     [12]     63     [25]     0.02

HD034759*          B5V             [1]       $5.08\pm0.02$     $5.22\pm0.02$     [14]     179     [25]    0.02

HD034816*          B0.5IV          [1]       $4.04\pm0.03$     $4.29\pm0.02$     [12]     270     [25]     0.03

HD036512          B0V             [1]       $4.36\pm0.03$     $4.62\pm0.02$     [12]             0.04

HD042560*          B3IV            [1]       $4.31\pm0.03$     $4.48\pm0.03$     [12]     210     [25]     0.03

HD047839          O7V((f))z var       [2]       $4.42\pm0.02$     $4.66\pm0.02$     [12]     609     [26]     0.08

HD051283*          B3II-III        [3]       $5.15\pm0.03$     $5.32\pm0.02$     [15]     736     [25]     0.01

HD078316          B8IIIp          [4]       $5.13\pm0.02$     $5.24\pm0.02$     [16]     154     [25]     -0.01

HD091316          B1Iab           [1]       $3.71\pm0.02$     $3.85\pm0.02$     [12]     334     [26]     0.05

HD164794          O4V((f))z       [2]       $5.97\pm0.03$     $5.97\pm0.02$     [12]     950     [26]

HD188209*          O9.5Ia          [1]       $5.55\pm0.03$     $5.62\pm0.02$     [12]     1112     [26]     0.20

HD204172          B0Ib            [1]       $5.85\pm0.02$     $5.93\pm0.02$     [17]             0.16

HD214680*          O9V             [1]       $4.68\pm0.03$     $4.88\pm0.02$     [12]     359     [25]     0.11

9creddened stars

[BD +56 524]BD+56d524*         B1V             [5]       $10.09\pm0.01$     $9.75\pm0.01$     [11]     2349     [25]

HD013338*          B1V             [3]       $9.27\pm0.01$     $9.03\pm0.01$     [11]     2425     [25]

HD014250          B0.5V:n         [3]       $9.28\pm0.01$     $8.96\pm0.01$     [11]      1245     [25]

HD014422          B1V:pe          [3]       $9.53\pm0.01$     $9.03\pm0.01$     [18]

HD014956*          B2Ia            [3]       $7.91\pm0.03$     $7.19\pm0.02$     [19]

HD017505*          O6.5III((f))n+O8V           [6]     $7.46\pm0.01$     $7.06\pm0.01$     [11]

HD029309*          B2V             [7]       $7.42\pm0.01$     $7.10\pm0.01$     [20]     544     [25]

HD029647*          B7IV            [8]       $9.22\pm0.03$     $8.31\pm0.02$     [21]     155     [25]

HD034921          B0IVpe          [3]       $7.65\pm0.01$     $7.51\pm0.01$     [22]      1327     [25]

HD037020          O8Vn            [7]       $6.72\pm0.01$     $6.72\pm0.01$     [11]      421     [25]

HD037022          O7Vp            [6]       $5.13\pm0.02$     $5.13\pm0.01$     [11]      430     [26]

HD037023          B0.5Vp          [5]       $6.77\pm0.01$     $6.69\pm0.01$     [11]      472     [25]

HD037061*          B1V             [9]       $7.09\pm0.01$     $6.83\pm0.01$     [11]     523     [25]

HD038087*         B3II            [10]      $8.42\pm0.01$     $8.30\pm0.01$     [11]      339     [25]

HD052721          B2Vne           [7]       $6.64\pm0.06$     $6.58\pm0.05$     [11]

HD156247*          B3V             [10]      $5.96\pm0.04$     $5.91\pm0.04$     [11]      267     [25]

HD166734          O7.5Iabf        [6]       $9.51\pm0.01$     $8.42\pm0.01$     [22]      3418     [26]

HD183143*          B7Ia            [3]       $8.08\pm0.01$     $6.86\pm0.01$     [23]

HD185418*          B0.5V           [3]       $7.67\pm0.01$     $7.45\pm0.01$     [11]     755     [25]

HD192660*          B0Ia            [3]       $8.05\pm0.04$     $7.38\pm0.03$     [15]     2003     [25]

HD204827*          B0V             [3]       $8.76\pm0.01$     $7.95\pm0.01$     [22]     1646     [26]

HD206773          B0V:pe          [3]       $7.01\pm0.02$     $6.79\pm0.02$     [24]      958     [25]

HD229238*          B0Iab           [3]       $9.78\pm0.01$     $8.88\pm0.01$     [22]     1473     [25]

HD283809*          B3V             [8]       $12.14\pm0.03$     $10.72\pm0.02$     [21]     326     [25]

HD294264          B3V             [11]      $9.83\pm0.01$     $9.47\pm0.01$     [11]      450     [25]

[1] Lesh, 1968; [2] Sota et al., 2014; [3] Morgan et al., 1955; [4] Levato, 1975; [5] Borgman, 1960; [6] Sota et al., 2011; [7] Guetter, 1968; [8] Murakawa et al., 2000; [9] Sharpless, 1952; [10] Houk & Swift, 1999; [11] Valencic et al., 2004; [12] Johnson et al., 1966; [13] Menzies et al., 1990; [14] Crawford et al., 1971; [15] Fernie, 1983; [16] Abt & Golson, 1962; [17] Lutz & Lutz, 1977; [18] Johnson & Morgan, 1955; [19] Mendoza, 1967; [20] Guetter, 1974; [21] Slutskij et al., 1980; [22] Hiltner, 1956; [23] Johnson, 1965; [24] Garrison & Kormendy, 1976; [25] Gaia Collaboration, 2018; [26] Megier et al., 2009.
Stars with an * were used to measure extinction curves, as explained in Section 3.

§.§ 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.85.5  spectrograph, in two modes: the Short wavelengths Cross-Dispersed mode (SXD, $\sim$0.82.4 ), and the Long wavelengths Cross-Dispersed mode (LXD, $\sim$1.95.5 ) [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.

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  GUI. Subsequently, multiple spectra of the same star were combined using the  GUI. Next, the stellar spectrum was corrected for the telluric absorption and the instrument throughput with the  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 , 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. <ref>, and were masked and excluded from all further analyses (such as measuring and fitting the extinction curves).

An atmospheric transmission model, obtained from the Spextool database, and computed by [Cushing 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.

All spectra were then placed on the same wavelength grid (between 0.8 and 2.45  for SXD, and between 1.9 and 5.5  for LXD) with a resolution of 2000, using the  python package [Gordon & Decleir, 2021]. Finally, with the same package, the spectra were calibrated based on photometric data points, using the  python package [Gordon & Decleir, 2021]. First, the SXD spectrum was scaled to match J, H, KS photometry from 2MASS [Skrutskie et al., 2006], obtained from the [IRSA 2MASS All-Sky Point Source Catalog]10.26131/IRSA2 [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.[The LXD part could also be scaled using IRAC or WISE photometry, but this is not available for all stars in the sample. Furthermore, the uncertainties on the available photometric data points are fairly large, and using those data did not always result in a smooth NIR spectrum.] The final spectra of all stars are electronically [available]10.5281/zenodo.5802469 [Decleir, 2022], and can be found in Figs. <ref>–<ref>. They are plotted multiplied by $\lambda^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 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 , and multiplied by $\lambda^4$ to remove the strongly decreasing Rayleigh-Jeans tail.

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 , and multiplied by $\lambda^4$ to remove the strongly decreasing Rayleigh-Jeans tail.

The NIR spectra of the 10 reddened stars that could not be used to measure extinction curves because of emission lines or other peculiarities in the spectrum, ordered by steepness. 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 , and multiplied by $\lambda^4$ to remove the strongly decreasing Rayleigh-Jeans tail.

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 suffering 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 $\sim$5.2  for most stars. The median SNRs of the stellar spectra used in this work (indicated with an * in Table <ref>) are in the range 107–141 below 2.5 , 38–137 between 2.8 and 4 , and 18–72 beyond 4.5 , 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 . 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). We matchedUsually, 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 were taken from SIMBAD [Wenger et al., 2000], and are reported in Table <ref> together with their reference.The spectral types and luminosity classes listed in Table <ref> 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 <cit.>. For these reasons, we did not use the literature spectral types in Table <ref> to match comparison and reddened stars, but merely list them as a reference.[Note that we used the literature spectral types to estimate the $E(B-V)$-values of the comparison stars, as explained in Section <ref>.] Instead, we validated empirically which comparison star is the best match. We tried severalall 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 (hydrogen) lines and jumps as well as possible.
Tables <ref> and <ref> 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. <ref>).

The extinction curves were calculated with the  package [Gordon & Decleir, 2021] using the following procedure. The absolute dust extinction $A(\lambda)$ at a wavelength $\lambda$ is given by:

\begin{equation}
A(\lambda) = -2.5 \: \text{log} \left[\frac{F(\lambda)_{\text{red}}}{F(\lambda)_{\text{comp}}} \right] + 5 \: \text{log} \left[ \frac{1/d_{\text{red}}}{1/d_{\text{comp}}} \right]
\end{equation}

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(\lambda-V)$ (also called differential extinction, dust reddening or color excess) at wavelength $\lambda$ can be calculated as:

\begin{align}
E(\lambda-V) &= A(\lambda) - A(V)\label{eq:reddening} \\
&= -2.5 \: \text{log} \left[ \frac{F(\lambda)_{\text{red}}}{F(\lambda)_{\text{comp}}} \right] + 2.5 \: \text{log} \left[ \frac{F(V)_{\text{red}}}{F(V)_{\text{comp}}} \right]
\end{align}

The V-band magnitudes of all stars are given in Table <ref>.

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(\lambda-V)$ can be converted into an absolute normalized extinction measurement $A(\lambda)/A(V)$ from Eq. <ref>:

\begin{equation}
\label{eq:alav}
\frac{A(\lambda)}{A(V)} = \frac{E(\lambda-V)}{A(V)} + 1
\end{equation}

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. IfOnce $A(V)$ is known, one can also compute the total-to-selective extinction $R(V)$:

\begin{equation}
R(V) = \frac{A(V)}{E(B-V)}
\label{eq:RV}
\end{equation}

with the color excess $E(B-V)$ calculated as

\begin{equation}
E(B-V) = (m(B)_{\text{red}} - m(B)_{\text{comp}}) - (m(V)_{\text{red}}-m(V)_{\text{comp}})
\end{equation}

with $m(B)$ the apparent B-band magnitude and $m(V)$ the apparent V-band magnitude listed in Table <ref>. The obtained values for $E(B-V)$ are given in Tables <ref> and <ref>.[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. $R(V)$ probes the dust grain size along the line-of-sight, with larger values corresponding to sightlines dominated by larger dust grains <cit.>. On average, $R(V)=3.1$ in the Milky Way diffuse ISM [Cardelli et al., 1989].

The normalized extinction curves are shown in Fig. <ref>, and are electronically [available]10.5281/zenodo.5802469 [Decleir, 2022]. We were able to measure extinction curves for 15 reddened stars of our sample (with spectra shown in Fig. <ref>). The other 10 reddened stars were not suitable to measure extinction curves, due to peculiarities in their spectrum (plotted in Fig. <ref>). 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. <ref> 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]10.26131/IRSA1 [Wright et al., 2019]. The $J-K_S$ color is a measure of the reddening and high $K_S-$WISE 4 values ($\gtrsim1$) 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 <ref> which stars we used to measure extinction curves. These are plotted as green diamonds in Fig. <ref>, 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 , 0.4–5.8 between 2.8 and 4 , and 0.3–1.5 beyond 4.5 , at a resolution of about 2000. We would like to note here that some of the 2MASS photometric data points in Fig. <ref> 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 KS) were used to scale the (SXD) spectrum (as discussed in Section <ref>), 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 JHKS photometric data points (for all stars) are within 3$\sigma$ 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.

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.

$K_S-$WISE 4 color versus $J-K_S$ color for the comparison and reddened stars in our sample. The horizontal dotted line at $K_S-\text{WISE4} = 1$ divides between sources with and without strong stellar winds. The red circles represent reddened stars that were not used to measure extinction curves (see Fig. <ref>), 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 <ref>). 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.

§ FITTING NIR EXTINCTION CURVES

§.§ Continuum extinction

In previous works, the continuum NIR extinction was usually described with a power law <cit.>, 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  with a power law. From Eq. <ref>, the measured differential extinction can be written as:

\begin{equation}
\label{eq:elv}
E(\lambda - V) = A(V) \left[ k(\lambda) - 1 \right]
\end{equation}

\begin{equation}
\label{eq:powerlaw}
k(\lambda) = \frac{A(\lambda)}{A(V)}  =  S\,\lambda^{-\alpha}
\end{equation}

where $S$ is the amplitude (by definition equal to the normalized extinction at 1 ) and $\alpha$ is the index of the power law. At infinite wavelengths, Eq. <ref> reduces to $E(\infty-V) = -A(V)$. In other words, fitting the observed extinction curve with the function in Eqs. <ref> and <ref>, gives us a direct measurement of $A(V)$. We used a combination of the  Astropy [Astropy Collaboration et al., 2013] model (for Eq. <ref>), and the  conversion function from the  python package [Gordon et al., 2022] to implement the conversion in Eq. <ref>.

The fitting was done in two steps. First, the Levenberg–Marquardt algorithm was used to obtain preliminary fit results for the three parameters $S$, $\alpha$ and $A(V)$, using the  from Astropy. These results were then used as initial guesses for the Markov chain Monte Carlo (MCMC) fitting based on the  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. <ref>.

The uncertainties on the fitted parameters given by the MCMC fitting only include random noise, computed as the difference between the 84th and 50th percentile (upper uncertainty), or between the 50th and 16th percentile (lower uncertainty) of the posterior distribution function. However, there is also a systematic uncertainty related to the fact that the extinction curves are 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  package on GitHub [Decleir, 2022].

§.§ 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 : HD283809 and HD029647, as can be seen in Fig. <ref>. This feature is caused by water ice along the line-of-sight, as will be discussed in more detail in Section <ref>. Since it is generally accepted that the detection of ice indicates the presence of a dense molecular component along the line-of-sight <cit.>, 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)\sim10$, [Whittet, 2015]). The extinction curve of HD283809 also shows a weaker bump around 3.4 , which is likely caused by hydrocarbons, as will be discussed in Section <ref>. This feature is, however, not obvious in the extinction curve towards HD029647. There also seems to be a small systematic residual around 2.3  for HD283809, which is not visible for HD029647. However, it is possible that this residual is due to the less smooth transition between the SXD and LXD spectrum around 2.4 .

The extinction curves of the sightlines towards HD283809 (left) and HD029647 (right) show a strong ice feature around 3 . 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 Table <ref>.

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 the main goalone 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. A brief discussion on the identification of these features is given in Section <ref>.

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 <cit.>. 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 introduces uncertainties on the measurement of the extinction features <cit.>. Therefore, we have opted to fit the extinction continuum and features simultaneously in the extinction curves.

The strong telluric absorption (see Fig. <ref>) significantly increases the noise between 2.5 and 3.5 , which makes it hard to obtain any meaningful constraints on the weaker feature around 3.4 . 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(\lambda)$:

\begin{equation}
\label{eq:full}
k(\lambda) = \frac{A(\lambda)}{A(V)}  =  S\,\lambda^{-\alpha} + B\, D(\lambda)
\end{equation}

\begin{equation}
D(\lambda)= \frac{(\gamma/\lambda_0)^2}{(\lambda/\lambda_0-\lambda_0/\lambda)^2+(\gamma/\lambda_0)^2}
\label{eq:drude}
\end{equation}

with $B$ the strength, $\lambda_0$ the central wavelength, and $\gamma$ the width of the profile. Because of the asymmetric appearance of the feature, we modified the standard Drude function to allow for extra asymmetry.[Note that a Drude profile is intrinsically already asymmetric, but our modified version allows for extra asymmetry by implementing an asymmetry parameter.] The modified profile is obtained by replacing the standard width in Eq. <ref> by a width $\gamma(\lambda)$ that depends on the wavelength:

\begin{equation}
\gamma(\lambda) = \frac{2 \gamma_0}{1 + \text{exp}[a(\lambda - \lambda_0)]}
\end{equation}

with $\gamma_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  silicate features in MIR extinction curves.

The fitting results (i.e. the 50th percentiles of the posterior distribution functions) and their uncertainties (including both random and systematic uncertainties, as explained in Section <ref>) for the two dense sightlines are given in Table <ref>. 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\sigma$ 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 . In sightlines HD183143 and HD229238, there is a very tentative detection of a weak ice feature with strengths $B=0.0026\pm0.0008$ ($3.3\sigma$) and $B=0.0042\pm0.0012$ ($3.5\sigma$) respectively (as compared to a $>50\sigma$ 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, theythese 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 <ref>.

MCMC fitting results for the 13 diffuse extinction curves and the average diffuse extinction curve: the amplitude $S$ and index $\alpha$ of the power law, and $A(V)$ are directly obtained from the fitting (50th percentile of the posterior distribution function), while $E(B-V)$ is obtained from the observations, and $R(V)$ is calculated as $R(V)=A(V)/E(B-V)$. The reported uncertainties include both random and systematic uncertainties, as explained in Section <ref>.
reddened     comparison     $S$     $\alpha$     $A(V)$     $E(B-V)$     $R(V)$
[BD +56 524]BD+56d524     HD034816     0.402_-0.005^+0.005     1.58_-0.02^+0.02     1.78_-0.02^+0.02     0.59_-0.04^+0.04     3.01_-0.19^+0.20

HD013338     HD031726     0.435_-0.004^+0.004     1.69_-0.02^+0.02     1.38_-0.01^+0.01     0.45_-0.01^+0.01     3.10_-0.10^+0.11

HD014956     HD214680     0.404_-0.004^+0.004     1.44_-0.01^+0.01     2.82_-0.03^+0.03     0.92_-0.05^+0.05     3.07_-0.16^+0.17

HD017505     HD214680     0.210_-0.003^+0.003     2.20_-0.03^+0.03     1.51_-0.02^+0.02     0.60_-0.04^+0.04     2.51_-0.15^+0.17

HD029309     HD042560     0.484_-0.008^+0.008     1.42_-0.01^+0.01     2.02_-0.03^+0.03     0.49_-0.05^+0.05     4.12_-0.36^+0.43

HD037061     HD034816     0.458_-0.004^+0.004     1.36_-0.01^+0.01     2.52_-0.02^+0.02     0.51_-0.04^+0.04     4.94_-0.34^+0.40

HD038087     HD051283     0.445_-0.007^+0.007     1.66_-0.01^+0.01     1.55_-0.03^+0.03     0.29_-0.04^+0.04     5.33_-0.64^+0.85

HD156247     HD042560     0.403_-0.026^+0.026     1.93_-0.06^+0.06     0.78_-0.05^+0.05     0.22_-0.07^+0.07     3.52_-0.88^+1.67

HD183143     HD188209     0.367_-0.002^+0.002     1.90_-0.01^+0.01     3.86_-0.02^+0.02     1.29_-0.04^+0.04     2.99_-0.09^+0.09

HD185418     HD034816     0.325_-0.006^+0.006     2.02_-0.03^+0.03     1.23_-0.02^+0.02     0.47_-0.04^+0.04     2.61_-0.20^+0.23

HD192660     HD214680     0.391_-0.006^+0.006     1.71_-0.01^+0.01     2.31_-0.04^+0.04     0.87_-0.06^+0.06     2.65_-0.18^+0.21

HD204827     HD003360     0.340_-0.005^+0.005     2.08_-0.01^+0.01     2.43_-0.03^+0.03     1.00_-0.05^+0.05     2.43_-0.11^+0.12

HD229238     HD214680     0.340_-0.003^+0.003     1.99_-0.01^+0.01     2.99_-0.02^+0.02     1.10_-0.04^+0.04     2.71_-0.09^+0.10

2c|average diffuse     0.386_-0.001^+0.001a     1.71_-0.01^+0.01a               3.12$\pm$0.05b

aObtained by fitting the average diffuse extinction curve (see Section <ref>).
bCalculated as $R(V)=\frac{1}{\left(\frac{A(B)}{A(V)}-1\right)}$.
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 <ref>.

MCMC fitting results for the 2 dense extinction curves: the amplitude $S$ and index $\alpha$ of the power law, the strength $B$, central wavelength $\lambda_0$, width $\gamma_0$ and asymmetry $a$ of the Drude profile, and $A(V)$ are directly obtained from the fitting (50th percentile of the posterior distribution function), while $E(B-V)$ is obtained from the observations, and $R(V)$ is calculated as $R(V)=A(V)/E(B-V)$. The reported uncertainties include both random and systematic uncertainties, as explained in Section <ref>.
reddened     comparison     $S$     $\alpha$     $B$     $\lambda_0$     $\gamma_0$     $a$     $A(V)$     $E(B-V)$     $R(V)$
HD029647     HD034759     0.369_-0.003^+0.003     1.93_-0.01^+0.01     0.0331_-0.0006^+0.0007     3.022_-0.005^+0.005     0.49_-0.03^+0.03     -1.32_-0.22^+0.19     3.52_-0.03^+0.03     1.05_-0.04^+0.04     3.35_-0.13^+0.14

HD283809     HD003360     0.389_-0.002^+0.002     1.91_-0.01^+0.01     0.0264_-0.0005^+0.0005     3.014_-0.006^+0.005     0.44_-0.02^+0.02     -4.38_-0.83^+0.64     5.65_-0.03^+0.03     1.61_-0.06^+0.06     3.51_-0.12^+0.13
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 <ref>.

§ RESULTS AND DISCUSSION

§.§ Extinction curve parameter trends

The fitting results in Tables <ref> and <ref> show a wide range in $A(V)$ (0.78–5.65) and $R(V)$ values (2.332.43–5.33).
As a cross-check, in Fig. <ref> 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(\lambda-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(\lambda-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(\lambda-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.

A comparison between our obtained $A(V)$ values from the fitting (see Tables <ref> and <ref>) to values reported in the literature for the same sightlines (taken from Cardelli et al., 1989, Valencic et al., 2004, Gordon et al., 2009 and Gordon et al., 2021.). There is a good agreement for all sightlines.

It is interesting to see if there are any correlations between the different parameters describing the extinction curve: the amplitude $S$ and the index $\alpha$ of the power law, the V-band extinction $A(V)$ and the total-to-selective extinction $R(V)$. Fig. <ref> shows scatter plots of these parameters. We found a clear anti-correlation (with Spearman's rank correlation coefficient $\rho=-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 , and that the extinction curves are normalized to $A(V)$. A steeper extinction curve (i.e. larger $\alpha$) 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  (i.e. smaller $S$). $R(V)$ correlates with the amplitude of the power law ($\rho=0.77$$\rho=0.84$), and anti-correlates with the power law index ($\rho=-0.61$$\rho=-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 $\alpha$ and larger $S$), as is also the case in the UV and optical <cit.>. The $R(V)$-dependence of the shape of the extinction curve is investigated and discussed in more detail in Section <ref>. No correlation was found between $A(V)$ and any of the other parameters ($\rho=-0.10$, $\rho=-0.09$ and $\rho=0.00$$\rho=0.03$).

Scatter plots of the different parameters describing the extinction curve: the amplitude $S$ and the index $\alpha$ of the  power law, the V-band extinction $A(V)$ and the total-to-selective extinction $R(V)$ (see Tables <ref> and <ref>). 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 <ref>). The Spearman's rank correlation coefficient is shown in every plot.

§.§ Average diffuse NIR Milky Way extinction curve

§.§.§ 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 ). The average NIR curve is plotted in the right panel of Fig. <ref>, together with severalsome 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 <ref>. Given that our average curve is close to most diffuse literature curves, we can 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]10.17909/165k-w223 <cit.>. The average UV curve for our sample is plotted in the left panel of Fig. <ref>, 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\pm0.05$ for our average curve, which is consistent with the average diffuse Milky Way value [Cardelli et al., 1989]. The peaks in the UV extinction curve show that there are small spectral mismatches between some of the reddened stars and their comparison star, which are less obvious in the NIR. A binned version of the measured average NIR extinction curveThe measured average NIR extinction curve is electronically [available]10.5281/zenodo.5802469 [Decleir, 2022], and a binned version of the curve and its uncertainty is given in Table <ref> (columns 2 and 3) for wavelengths between 0.8 and 4 , and in common IR photometric bands (2MASS J, H, KS, WISE 1, L, and IRAC 1), obtained by integrating the average extinction curve over the response curves of the bands (taken from the SVO Filter Profile Service – [Rodrigo et al., 2012, Rodrigo & Solano, 2020]).

Average diffuse Milky Way extinction curve and parameters of the linear relationship between extinction $A(\lambda)/A(V)$ and $1/R(V)$. This table is also available as part of the D22_MWAvg model in the  python package [Gordon et al., 2022].
band     3caverage extinction     3c$R(V)$-dependent ext.

or     2cmeasured     fit

$\lambda\ [\micron]$     $\frac{A(\lambda)}{A(V)}$     $\sigma\left(\frac{A(\lambda)}{A(V)}\right)$     $\frac{A(\lambda)}{A(V)}$     $a(\lambda$)     $b(\lambda$)     $\sigma(\lambda)$
J     0.269     0.018     0.268     0.258     -0.636     0.037

H     0.163     0.013     0.164     0.155     -0.458     0.034

KS     0.105     0.010     0.103     0.095     -0.349     0.028

W1     0.048     0.009     0.048     0.043     -0.290     0.020

L     0.045     0.008     0.046     0.041     -0.274     0.020

I1     0.044     0.009     0.045     0.040     -0.270     0.021

0.80     0.552     0.025     0.565     0.562     -0.874     0.068

0.85     0.523     0.023     0.509     0.504     -0.629     0.056

0.90     0.460     0.020     0.462     0.455     -0.617     0.051

0.95     0.420     0.020     0.421     0.414     -0.639     0.048

1.00     0.381     0.020     0.386     0.377     -0.647     0.045

1.05     0.349     0.019     0.355     0.346     -0.651     0.044

1.10     0.322     0.019     0.328     0.318     -0.658     0.042

1.15     0.304     0.019     0.304     0.294     -0.668     0.039

1.20     0.282     0.018     0.282     0.273     -0.635     0.038

1.25     0.264     0.017     0.263     0.254     -0.621     0.037

1.30     0.248     0.017     0.246     0.236     -0.632     0.036

1.35     0.240     0.017     0.231     0.221     -0.638     0.036

1.40     0.220     0.021     0.217     0.207     -0.617     0.034

1.45     0.210     0.019     0.204     0.195     -0.576     0.033

1.50     0.195     0.015     0.193     0.183     -0.527     0.033

1.55     0.185     0.014     0.182     0.173     -0.483     0.034

1.60     0.172     0.013     0.173     0.163     -0.461     0.034

1.65     0.162     0.013     0.164     0.155     -0.455     0.034

1.70     0.152     0.013     0.156     0.146     -0.448     0.034

1.75     0.145     0.013     0.148     0.139     -0.439     0.034

1.80     0.141     0.014     0.141     0.132     -0.431     0.033

1.85             0.135     0.126     -0.424     0.032

1.90             0.129     0.120     -0.417     0.031

1.95     0.125     0.012     0.123     0.115     -0.408     0.030

2.00     0.123     0.011     0.118     0.110     -0.395     0.029

2.05     0.116     0.011     0.113     0.105     -0.377     0.029

2.10     0.110     0.011     0.108     0.100     -0.359     0.029

2.15     0.106     0.010     0.104     0.096     -0.345     0.028

2.20     0.102     0.010     0.100     0.092     -0.335     0.028

2.25     0.097     0.010     0.096     0.089     -0.332     0.027

2.30     0.095     0.010     0.093     0.085     -0.336     0.026

2.35     0.093     0.009     0.089     0.082     -0.341     0.025

2.40     0.090     0.010     0.086     0.079     -0.347     0.025

2.45     0.081     0.010     0.083     0.076     -0.353     0.025

2.50     0.077     0.010     0.080     0.074     -0.359     0.025

2.55             0.078     0.071     -0.366     0.025

2.60             0.075     0.069     -0.372     0.025

2.65             0.073     0.066     -0.378     0.025

2.70             0.071     0.064     -0.383     0.025

2.75             0.068     0.062     -0.386     0.024

2.80             0.066     0.060     -0.388     0.024

2.85             0.064     0.058     -0.388     0.023

2.90     0.062     0.012     0.062     0.057     -0.385     0.022

2.95     0.058     0.012     0.061     0.055     -0.379     0.021

3.00     0.061     0.010     0.059     0.053     -0.370     0.020

3.05     0.060     0.010     0.057     0.052     -0.357     0.018

3.10     0.059     0.009     0.056     0.050     -0.342     0.017

3.15     0.057     0.008     0.054     0.049     -0.325     0.016

3.20     0.057     0.009     0.053     0.047     -0.307     0.016

3.25     0.055     0.008     0.051     0.046     -0.290     0.016

3.30     0.054     0.009     0.050     0.045     -0.276     0.017

3.35     0.051     0.007     0.049     0.044     -0.269     0.017

3.40     0.050     0.007     0.048     0.043     -0.266     0.018

3.45     0.047     0.007     0.046     0.041     -0.267     0.018

3.50     0.043     0.008     0.045     0.040     -0.268     0.019

3.55     0.039     0.008     0.044     0.039     -0.269     0.020

3.60     0.040     0.009     0.043     0.038     -0.269     0.021

3.65     0.039     0.008     0.042     0.037     -0.268     0.023

3.70     0.039     0.008     0.041     0.037     -0.266     0.024

3.75     0.037     0.010     0.040     0.036     -0.263     0.025

3.80     0.036     0.009     0.039     0.035     -0.259     0.025

3.85     0.033     0.009     0.038     0.034     -0.253     0.025

3.90     0.036     0.009     0.038     0.033     -0.246     0.025

3.95     0.032     0.009     0.037     0.033     -0.238     0.024

4.00     0.032     0.008     0.036     0.032     -0.229     0.022

Left: Average diffuse UV extinction curve (measured from IUE spectra [Valencic et al., 2004] and rebinned by a factor of 3 to reduce the noise), together with average diffuse curves from the literature: [Cardelli et al., 1989] and [Fitzpatrick et al., 2019]. Right: Average diffuse NIR extinction curve, together with average curves from the literature: [Rieke & Lebofsky, 1985], [Martin & Whittet, 1990], [Indebetouw et al., 2005], and [Gordon et al., 2021]. The red line shows the power law fit of the average curve.

We also fitted the average curve with a power law and obtained:[To verify that normalizing to the V-band integrated extinction, $A(V)$, does not affect the shape of the average curve, we also fitted the average curve, normalized to the monochromatic extinction at 1 , $A(1\,\micron)$, instead of $A(V)$. We find no difference in the power law index of our average extinction curve within the uncertainties.]

\begin{equation}
\label{eq:average}
\frac{A(\lambda)}{A(V)} = 0.386 \,\lambda^{-1.71}
\end{equation}

The fitted parameters are given in Table <ref>, and indicated with a red star in Fig. <ref>. The fitted model is shown in red in Fig. <ref>, and listed at specific wavelengths and in photometric bands in Table <ref> (column 4). The fitted model is also available as part of the D22_MWAvg model in the  python package [Gordon et al., 2022]. The residuals of the fitting (i.e. the data subtracted by the fit) are plotted in Fig. <ref>, 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 ), 0.0032 (1.4–1.8 ), 0.0030 (1.9–2.5 ), 0.0047 (2.8–4 ) and 0.016 (4.5–5.0 ). 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 Fig. <ref>. Some peaks (e.g. around 2 ) and dips (e.g. just below 3 ) in the residuals can mostly likely be attributed to the telluric absorption.

Bottom: Residuals (data–fit) of the average diffuse NIR extinction curve fit (see Fig. <ref>). The horizontal magenta lines indicate the standard deviations of the residuals in certain wavelength regions: 0.0050 (0.8–1.4 ), 0.0032 (1.4–1.8 ), 0.0030 (1.9–2.5 ), 0.0047 (2.8–4 ) and 0.016 (4.5–5.0 ). The blue vertical lines indicate hydrogen line residuals due to spectral mismatches. The Paschen (Pa) and Brackett (Br) jump mismatches are also indicated in blue. Top: Same atmospheric transmission model as in Fig. <ref>.

To verify whether the 3  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\pm0.0007$. This is below the $3\sigma$ detection threshold of 0.0021 $A(3\,\micron)/A(V)$$A(\text{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(\lambda)\sim\lambda^{-\alpha}$. Therefore, in this section, we use the power law index $\alpha$ to compare our result with previous measurements. Table <ref> compares the fitted power law index for our average diffuse NIR extinction curve to indices found in the literature.

Comparison of power law indices $\alpha$ from the literature.
reference     $\alpha$     method used    wavelength range     environment
this work     $1.71\pm0.01$     pair method     0.8–5.0  (spectra)     local diffuse ISM

[Draine, 1989]     $1.75$     multiplea     0.9–5  (I, J, H, K, L, M)     Galactica

[Cardelli et al., 1989]     1.61     pair method     0.9–3.3  (I, J, H, K, L)     local ISM

[Martin & Whittet, 1990]     $1.84\pm0.03$     pair method     0.9–4.8  (I, J, H, K, L, M)     local diffuse ISM

[Nishiyama et al., 2006]     $1.99\pm0.02$     RC CMD     1.2–2.2  (J, H, KS)     Galactic center

[Nishiyama et al., 2009]     2.0     RC+RGB CMD     1.2–2.2  (J, H, KS)     Galactic center

[Stead & Hoare, 2009]     $2.14^{+0.04}_{-0.05}$     RC CCD     1.2–2.2  (J, H, KS)     Galactic plane

[Fritz et al., 2011]     $2.13 \pm 0.08$     nebular H lines     1.282–2.166     Galactic center

$1.76 \pm 0.39$         2.166–2.758

$2.11 \pm 0.06$         1.282–2.758

[Alonso-García et al., 2017]     $2.47\pm0.11$     RC CMD     0.88–2.15  (Z, Y, J, H, KS)     Galactic center

[Nogueras-Lara et al., 2018]     $2.30\pm0.08$     RCb     1.27–2.16  (J, H, KS)     Galactic center

[Nogueras-Lara et al., 2019]     $2.43\pm0.03$     RCb     1.27–1.65  (J, H)     Galactic center

$2.23\pm0.03$         1.65–2.16  (H, KS)

$2.32\pm0.09$         1.27–2.16  (J, H, KS)

[Maíz Apellániz et al., 2020]     2.27     RC CCD     1.2–2.2  (J, H, KS)     Galactic plane

aThis 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 $\alpha$ value.

Our value for $\alpha$ is in reasonable agreement with that from [Draine, 1989], and not very deviant frombetween 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)\gtrsim8$) 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 <ref>), 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 line-of-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 <cit.>, which are likely only (clearly) visible in sightlines through molecular clouds and not in purely diffuse extinction regions <cit.>. On the other hand, the detection of the strong aliphatic hydrocarbon feature at 3.4  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. <ref>).
[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 <ref>), 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., 2006, Nishiyama et al., 2009], [Fritz et al., 2011], and [Alonso-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. <ref>, and as briefly discussed in Section <ref>, the variation in $\alpha$ 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 <ref>.

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 $\sim3$ . [Nogueras-Lara et al., 2019] measured a smaller index $\alpha$ in HKS than in JH (see Table <ref>). [Fritz et al., 2011] obtained a flatter extinction curve beyond $\sim3.7$ . They also found a smaller $\alpha$ between 2.2 and 2.8  than between 1.2 and 2.2  (see Table <ref>), but stated that this change in slope is not significant.
We demonstrate that our local average diffuse extinction curve between 0.8 and 4  can be represented by a single power law, since there are no long-range trends visible in the residuals in Fig. <ref>. 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. <ref>, we compare our average diffuse average NIR extinction curve to the [Draine, 2003, Draine, 2003], [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 theour average curve below 1 , while the other models correspond better beyond 3 .

Average diffuse NIR extinction curve, compared to the [Draine, 2003, Draine, 2003], [Zubko et al., 2004], [Compiègne et al., 2011] and [Jones et al., 2013] diffuse ISM dust grain models with $R(V) = 3.1$.
The [Weingartner & Draine, 2001] model is not shown in this plot as it is very similar to the [Draine, 2003, Draine, 2003] model.

§.§ 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 <ref> and in Fig. <ref>. Because of these strong sightline-to-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 . 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 JHKS measurements.

In this section, we investigate our sample of observed NIR extinction curves as a function of $R(V)$$1/R(V)$.[One could also use $R(V)$, but we found that $1/R(V)=E(B-V)/A(V)$ more naturally results in a linear relationship with the normalized extinction $A(\lambda)/A(V)$, given that both quantities have the same denominator.] Fig. <ref> shows $A(\lambda)/A(V)$ vs. $R(V)-3.1$$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 , bottom plot), the range in $R(V)$ values is limited due to the limited number of sightlines with data at wavelengths beyond 4 .

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 dashed lines. For comparison, the $R(V)$-dependent relationships from [Cardelli et al., 1989] and [Fitzpatrick et al., 2019] are added as blue dashed and orange dotted lines, respectively.

To quantify this correlation, using the  from Astropy, we fitted a linear function ( Astropy model) at every wavelength in the SpeX range where at least 5 sightlines have data ($\sim$2800 data points between 0.8 and 5.1 ):

\begin{equation}
\frac{A(\lambda)}{A(V)} = a(\lambda) + b(\lambda) \,\left[\frac{1}{R(V)}-\frac{1}{3.1}\right]
\end{equation}

By choosing $R(V)-3.1$$1/R(V)-1/3.1$ as the abscissa, the intercepts $a(\lambda)$ correspond to the average $R(V)=3.1$ Milky Way extinction curve. The slopes $b(\lambda)$ illustrate the $R(V)$-dependence. The fits are shown as green lines in Fig. <ref>. We also added the $R(V)$-dependent relationships from [Cardelli et al., 1989] and [Fitzpatrick et al., 2019] (where available) for comparison. At $\lambda=0.9\,\micron$ 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. <ref> 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(\lambda)/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(\lambda)<0$ for all NIR wavelengths (see Figs. <ref> and <ref>), and could thus not be caused by the (small) positive artificial correlation due to their common factor $A(V)$.[Note that this is different at UV wavelengths, see e.g. [Fitzpatrick] & Massa, 2007].

Intercepts $a$ (top), slopes $b$ (middle) and standard deviations about the fit $\sigma$ (bottom) of the linear $R(V)$$1/R(V)$-dependence, as a function of wavelength. The green crosses indicate the wavelengths plotted in Fig. <ref>. 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  are plotted in gray and were excluded from the interpolation because they are very noisy.

In order to provide a smooth representation of the $R(V)$-dependence for use in future extinction studies, using the  from Astropy, we fitted the intercepts with a power law ( Astropy model) and found:

\begin{equation}
a(\lambda) = 0.377 \,\lambda^{-1.78}
\end{equation}

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 <ref>. For the slopes and standard deviations, no straightforward functional form could be fitted. Instead, we interpolated the data with cubic splines using  [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  (plotted in gray in Fig. <ref>) 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. <ref>, and are tabulated in Table <ref> (columns 5–7) for wavelengths between 0.8 and 4 , and in common IR photometric bands (2MASS J, H, KS, 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  python package [Gordon et al., 2022].

Fig. <ref> 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 <ref>). 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 $\gtrsim1\,\micron$, as was also clear from Fig. <ref>.

$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.

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 <cit.>. We tested this by normalizing the curves to the extinction at 1 , and plotting the extinction as a function of $R(V)-3.1$$1/R(V)-1/3.1$ as before. We still found a strong linear correlation between extinction and $R(V)$$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  (Gordon et al., in prep.). One focus of this future work will be to investigate the impact of underlying correlations between $A(\lambda)/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 , 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 , and is primarily caused by the  stretching modes of bulk H2O 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  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 previously reported <cit.>. We show in Section <ref> that the feature can be fitted reasonably well with a modified Drude profile. Fig. <ref> 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  stretching mode in ammonia hydrates [Knacke et al., 1982, Dartois & d'Hendecourt, 2001].

Structure near 3.4  might arise from the  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. <ref>) severely complicates the measurement of any features 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 . Also, the data between 2.88 and 3.5  are very noisy. Fortunately, our approved cycle 1 GO program 2459[<https://www.stsci.edu/jwst/phase2-public/2459.pdf>] 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).

§ CONCLUSIONS

In this work, we presented the first comprehensive spectroscopic study of NIR extinction curves in the Milky Way. We obtained spectra with the IRTF/SpeX spectrograph between 0.8 and 5.5  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.332.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 , with indices and amplitudes strongly differing from sightline to sightline.

* TheOur average local diffuse NIR Milky Way extinction curve can be represented by a single power law with index $\alpha=1.7$ over the entire wavelength range between 0.8 and 5 , 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(\lambda)/A(V)$ in our sample linearly depend on the total-to-selective extinction $R(V)$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 ana strong ice extinction feature around 3 , which can be approximated by a modified Drude profile. These sightlines most likely contain dense material, and our fitting approach might mask a weak 3.4  hydrocarbon feature known to emanate from diffuse interstellar dust.

* We tentatively detect (with slightly over 3$\sigma$ significance) the 3  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  ice feature in the average diffuse extinction curve, with a $3\sigma$ limit of 0.0021 $A(3\,\micron)/A(V)$$A(\text{ice})/A(V)=0.0021$.

* The telluric atmospheric absorption complicates the characterization of any extinction features, as well as the continuum extinction above 4 . 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]10.5281/zenodo.5802469 [Decleir, 2022].
The code used for the analysis and plots is available on GitHub [Decleir, 2022].

We thank the referee for their feedback, which improved this paper.
This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.
This research has made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.
This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.
This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.
This research has made use of the SVO Filter Profile Service (<http://svo2.cab.inta-csic.es/theory/fps/>) supported from the Spanish MINECO through grant AYA2017-84089.

IRTF (SpeX), IRSA, MAST

Spextool [Cushing et al., 2004], Astropy [Astropy Collaboration et al., 2013], measure_extinction [Gordon & Decleir, 2021], dust_extinction [Gordon et al., 2022], spex_nir_extinction [Decleir, 2022]

[Abt & Golson, 1962]
Abt, H. A., & Golson, J. C. 1962, , 136, 35, 10.1086/147349

[Alonso-García et al., 2017]
Alonso-García, J., Minniti, D., Catelan, M., et al. 2017, ,
849, L13, 10.3847/2041-8213/aa92c3

[Astropy Collaboration et al., 2013]
Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013,
, 558, A33, 10.1051/0004-6361/201322068

[Boogert et al., 2015]
Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, ,
53, 541, 10.1146/annurev-astro-082214-122348

[Boogert et al., 2011]
Boogert, A. C. A., Huard, T. L., Cook, A. M., et al. 2011, , 729,
92, 10.1088/0004-637X/729/2/92

[Borgman, 1960]
Borgman, J. 1960, , 15, 255

[Cardelli et al., 1989]
Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, , 345, 245,

[Chiar & Tielens, 2006]
Chiar, J. E., & Tielens, A. G. G. M. 2006, , 637, 774,

[Clayton et al., 2015]
Clayton, G. C., Gordon, K. D., Bianchi, L. C., et al. 2015, , 815,
14, 10.1088/0004-637X/815/1/14

[Clayton et al., 2003]
Clayton, G. C., Gordon, K. D., Salama, F., et al. 2003, , 592, 947,

[Compiègne et al., 2011]
Compiègne, M., Verstraete, L., Jones, A., et al. 2011, , 525,
A103, 10.1051/0004-6361/201015292

[Crawford et al., 1971]
Crawford, D. L., Barnes, J. V., & Golson, J. C. 1971, , 76, 1058,

[Cushing et al., 2004]
Cushing, M. C., Vacca, W. D., & Rayner, J. T. 2004, , 116, 362,

[Dartois & d'Hendecourt, 2001]
Dartois, E., & d'Hendecourt, L. 2001, , 365, 144,

[Decleir, 2022]
Decleir, M. 2022a, SpeX NIR spectra and extinction curves, 2.0,
Zenodo, 10.5281/zenodo.5802469

[Decleir, 2022]
—. 2022b, mdecleir/spex_nir_extinction: Release related to
final SpeX NIR extinction paper, v1.1.0,  Zenodo,

[Draine, 1989]
Draine, B. T. 1989, in Infrared Spectroscopy in Astronomy, ed.
E. Böhm-Vitense, 93

[Draine, 2003]
Draine, B. T. 2003a, , 41, 241,

[Draine, 2003]
—. 2003b, , 598, 1017, 10.1086/379118

[Fernie, 1983]
Fernie, J. D. 1983, , 52, 7, 10.1086/190856

[Fitzgerald, 1970]
Fitzgerald, M. P. 1970, , 4, 234

[Fitzpatrick, 1999]
Fitzpatrick, E. L. 1999, , 111, 63, 10.1086/316293

[Fitzpatrick & Massa, 1986]
Fitzpatrick, E. L., & Massa, D. 1986, , 307, 286,

[Fitzpatrick & Massa, 1988]
—. 1988, , 328, 734, 10.1086/166332

[Fitzpatrick & Massa, 2005]
—. 2005, , 130, 1127, 10.1086/431900

[Fitzpatrick & Massa, 2007]
—. 2007, , 663, 320, 10.1086/518158

[Fitzpatrick & Massa, 2009]
—. 2009, , 699, 1209, 10.1088/0004-637X/699/2/1209

[Fitzpatrick et al., 2019]
Fitzpatrick, E. L., Massa, D., Gordon, K. D., Bohlin, R., & Clayton,
G. C. 2019, , 886, 108, 10.3847/1538-4357/ab4c3a

[Foreman-Mackey et al., 2013]
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, ,
125, 306, 10.1086/670067

[Fritz et al., 2011]
Fritz, T. K., Gillessen, S., Dodds-Eden, K., et al. 2011, , 737,
73, 10.1088/0004-637X/737/2/73

[Gaia Collaboration, 2018]
Gaia Collaboration. 2018, VizieR Online Data Catalog, I/345

[Garrison & Kormendy, 1976]
Garrison, R. F., & Kormendy, J. 1976, , 88, 865,

[Gordon & Decleir, 2021]
Gordon, K., & Decleir, M. 2021, karllark/measure_extinction: NIR SpeX paper
release, v1.2,  Zenodo, 10.5281/zenodo.5806655

[Gordon et al., 2022]
Gordon, K., Larson, K., McBride, A., et al. 2022, karllark/dust_extinction:
NIRSpectralExtinctionAdded, v1.1,  Zenodo, 10.5281/zenodo.6397654

[Gordon et al., 2009]
Gordon, K. D., Cartledge, S., & Clayton, G. C. 2009, , 705, 1320,

[Gordon et al., 2003]
Gordon, K. D., Clayton, G. C., Misselt, K. A., Landolt, A. U., &
Wolff, M. J. 2003, , 594, 279, 10.1086/376774

[Gordon et al., 2021]
Gordon, K. D., Misselt, K. A., Bouwman, J., et al. 2021, , 916, 33,

[Guetter, 1968]
Guetter, H. H. 1968, , 80, 197, 10.1086/128611

[Guetter, 1974]
—. 1974, , 86, 795, 10.1086/129675

[Hanson et al., 1996]
Hanson, M. M., Conti, P. S., & Rieke, M. J. 1996, , 107, 281,

[Hiltner, 1956]
Hiltner, W. A. 1956, , 2, 389, 10.1086/190029

[Houk & Swift, 1999]
Houk, N., & Swift, C. 1999, Michigan Spectral Survey, 5, 0

[Indebetouw et al., 2005]
Indebetouw, R., Mathis, J. S., Babler, B. L., et al. 2005, , 619,
931, 10.1086/426679

[Johnson, 1965]
Johnson, H. L. 1965, , 141, 923, 10.1086/148186

[Johnson et al., 1966]
Johnson, H. L., Mitchell, R. I., Iriarte, B., & Wisniewski, W. Z.
1966, Communications of the Lunar and Planetary Laboratory, 4, 99

[Johnson & Morgan, 1955]
Johnson, H. L., & Morgan, W. W. 1955, , 122, 429,

[Jones et al., 2013]
Jones, A. P., Fanciullo, L., Köhler, M., et al. 2013, , 558,
A62, 10.1051/0004-6361/201321686

[Jones & Hyland, 1980]
Jones, T. J., & Hyland, A. R. 1980, , 192, 359,

[Knacke et al., 1982]
Knacke, R. F., McCorkle, S., Puetter, R. C., Erickson, E. F., &
Kraetschmer, W. 1982, , 260, 141, 10.1086/160241

[Lanz & Hubeny, 2003]
Lanz, T., & Hubeny, I. 2003, , 146, 417, 10.1086/374373

[Lesh, 1968]
Lesh, J. R. 1968, , 17, 371, 10.1086/190179

[Levato, 1975]
Levato, H. 1975, , 19, 91

[Lord, 1992]
Lord, S. D. 1992, A new software tool for computing Earth's atmospheric
transmission of near- and far-infrared radiation, NASA Technical Memorandum

[Lutz & Lutz, 1977]
Lutz, T. E., & Lutz, J. H. 1977, , 82, 431, 10.1086/112066

[Maíz Apellániz, 2013]
Maíz Apellániz, J. 2013, in Highlights of Spanish Astrophysics VII,
ed. J. C. Guirado, L. M. Lara, V. Quilis, & J. Gorgas, 583–589.

[Maíz Apellániz & Barbá, 2018]
Maíz Apellániz, J., & Barbá, R. H. 2018, , 613, A9,

[Maíz Apellániz et al., 2020]
Maíz Apellániz, J., Pantaleoni González, M., Barbá,
R. H., García-Lario, P., & Nogueras-Lara, F. 2020, , 496,
4951, 10.1093/mnras/staa1790

[Maíz Apellániz & Rubio, 2012]
Maíz Apellániz, J., & Rubio, M. 2012, , 541, A54,

[Martin & Whittet, 1990]
Martin, P. G., & Whittet, D. C. B. 1990, , 357, 113,

[Massa et al., 2020]
Massa, D., Fitzpatrick, E. L., & Gordon, K. D. 2020, , 891, 67,

[Massa et al., 1983]
Massa, D., Savage, B. D., & Fitzpatrick, E. L. 1983, , 266, 662,

[Mathis & Cardelli, 1992]
Mathis, J. S., & Cardelli, J. A. 1992, , 398, 610,

[Megier et al., 2009]
Megier, A., Strobel, A., Galazutdinov, G. A., & Krełowski, J. 2009,
, 507, 833, 10.1051/0004-6361/20079144

[Mendoza, 1967]
Mendoza, E. E. 1967, Boletin de los Observatorios Tonantzintla y Tacubaya, 4,

[Menzies et al., 1990]
Menzies, J. W., Marang, F., & Westerhuys, J. E. 1990, South African
Astronomical Observatory Circular, 14, 33

[Morgan et al., 1955]
Morgan, W. W., Code, A. D., & Whitford, A. E. 1955, , 2, 41,

[Murakawa et al., 2000]
Murakawa, K., Tamura, M., & Nagata, T. 2000, , 128, 603,

[Nishiyama et al., 2009]
Nishiyama, S., Tamura, M., Hatano, H., et al. 2009, , 696, 1407,

[Nishiyama et al., 2006]
Nishiyama, S., Nagata, T., Kusakabe, N., et al. 2006, , 638, 839,

[Nogueras-Lara et al., 2019]
Nogueras-Lara, F., Schödel, R., Najarro, F., et al. 2019, ,
630, L3, 10.1051/0004-6361/201936322

[Nogueras-Lara et al., 2018]
Nogueras-Lara, F., Gallego-Calvente, A. T., Dong, H., et al. 2018,
, 610, A83, 10.1051/0004-6361/201732002

[Pendleton & Allamandola, 2002]
Pendleton, Y. J., & Allamandola, L. J. 2002, , 138, 75,

[Pendleton et al., 1994]
Pendleton, Y. J., Sandford, S. A., Allamandola, L. J., Tielens,
A. G. G. M., & Sellgren, K. 1994, , 437, 683, 10.1086/175031

[Potapov et al., 2021]
Potapov, A., Bouwman, J., Jäger, C., & Henning, T. 2021, Nature
Astronomy, 5, 78, 10.1038/s41550-020-01214-x

[Rayner et al., 2003]
Rayner, J. T., Toomey, D. W., Onaka, P. M., et al. 2003, , 115,
362, 10.1086/367745

[Rieke & Lebofsky, 1985]
Rieke, G. H., & Lebofsky, M. J. 1985, , 288, 618,

[Rieke et al., 1989]
Rieke, G. H., Rieke, M. J., & Paul, A. E. 1989, , 336, 752,

[Rodrigo & Solano, 2020]
Rodrigo, C., & Solano, E. 2020, in Contributions to the XIV.0 Scientific
Meeting (virtual) of the Spanish Astronomical Society, 182

[Rodrigo et al., 2012]
Rodrigo, C., Solano, E., & Bayo, A. 2012, SVO Filter Profile Service
Version 1.0, IVOA Working Draft 15 October 2012,

[Sharpless, 1952]
Sharpless, S. 1952, , 116, 251, 10.1086/145610

[Skrutskie et al., 2006]
Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, , 131,
1163, 10.1086/498708

[Skrutskie et al., 2019]
Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2019, 2MASS All-Sky
Point Source Catalog,  IPAC, 10.26131/IRSA2

[Slutskij et al., 1980]
Slutskij, V. E., Stalbovskij, O. I., & Shevchenko, V. S. 1980, Soviet
Astronomy Letters, 6, 397

[Smith, 1987]
Smith, R. G. 1987, , 227, 943, 10.1093/mnras/227.4.943

[Smith et al., 1993]
Smith, R. G., Sellgren, K., & Brooke, T. Y. 1993, , 263, 749,

[Smith Neubig & Bruhweiler, 1997]
Smith Neubig, M. M., & Bruhweiler, F. C. 1997, , 114, 1951,

[Sota et al., 2014]
Sota, A., Maíz Apellániz, J., Morrell, N. I., et al. 2014,
, 211, 10, 10.1088/0067-0049/211/1/10

[Sota et al., 2011]
Sota, A., Maíz Apellániz, J., Walborn, N. R., et al. 2011,
, 193, 24, 10.1088/0067-0049/193/2/24

[Stancik & Brauns, 2008]
Stancik, A. L., & Brauns, E. B. 2008, Vibrational Spectroscopy, 47, 66,

[Stead & Hoare, 2009]
Stead, J. J., & Hoare, M. G. 2009, , 400, 731,

[Stecher, 1965]
Stecher, T. P. 1965, , 142, 1683, 10.1086/148462

[Thi et al., 2006]
Thi, W. F., van Dishoeck, E. F., Dartois, E., et al. 2006, , 449,
251, 10.1051/0004-6361:20052931

[Vacca et al., 2003]
Vacca, W. D., Cushing, M. C., & Rayner, J. T. 2003, , 115, 389,

[Valencic et al., 2004]
Valencic, L. A., Clayton, G. C., & Gordon, K. D. 2004, , 616, 912,

[Virtanen et al., 2020]
Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17,
261, 10.1038/s41592-019-0686-2

[Weingartner & Draine, 2001]
Weingartner, J. C., & Draine, B. T. 2001, , 548, 296,

[Wenger et al., 2000]
Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, , 143, 9,

[Whittet, 2015]
Whittet, D. C. B. 2015, , 811, 110, 10.1088/0004-637X/811/2/110

[Whittet et al., 1988]
Whittet, D. C. B., Bode, M. F., Longmore, A. J., et al. 1988, ,
233, 321, 10.1093/mnras/233.2.321

[Whittet et al., 1997]
Whittet, D. C. B., Boogert, A. C. A., Gerakines, P. A., et al. 1997,
, 490, 729, 10.1086/304914

[Wright et al., 2010]
Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, ,
140, 1868, 10.1088/0004-6256/140/6/1868

[Wright et al., 2019]
Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2019, AllWISE
Source Catalog,  IPAC, 10.26131/IRSA1

[Zubko et al., 2004]
Zubko, V., Dwek, E., & Arendt, R. G. 2004, , 152, 211,