The Optical/Near-infrared Extinction Law in Highly Reddened Regions

HST observations of Wd1 were obtained over an eight-year period between 2005 and 2013. The earliest observations were made in 2005 with the Advanced Camera for Surveys Wide Field Camera (ACS-WFC, GO-10172, PI: R. De Grijs) in the F814W filter. The total exposure time was 2407 s, comprised of three slightly dithered images covering a 211'' × 218'' field of view (005 pix⁻¹). A second set of observations was obtained in 2010 with the infrared channel of the Wide Field Camera 3 (WFC3-IR) in the F125W and F160W filters (GO-11708, PI: M. Andersen).⁷ The final observations were obtained in 2013 with WFC3-IR in the F160W filter to provide a third positional epoch (GO-13044, PI: J. R. Lu). These observations mimicked the 2010 F160W observations, although a different position angle was used due to new HST guide star restrictions. Since the field of view for WFC3-IR is ∼130'' (012 pix⁻¹), a 2 × 2 mosaic was used to cover the entire 2005 ACS-WFC field of view. Each pointing had seven images per filter for total exposure times of 2443 s and 2093 s in F125W and F160W, respectively.

The HST observations were reduced using the standard online HST data reduction pipeline, and the resulting "FLT" images were downloaded from the HST archive on 2011 December 14. High-quality astrometric and photometric measurements were extracted from the individual FLT images using KS2, an expansion of the software developed for the Globular Cluster Treasury Program (Anderson et al. 2008). With this software, the measurements are combined to produce a single starlist for each filter. The sources are matched across epochs and positions transformed into a common astrometric reference frame where proper motions can be calculated. A detailed explanation of this process is given in H15.

To calculate the KS2 photometric zero-points, we compared the KS2 starlists with calibrated photometry obtained using DOLPHOT, a version of HSTPHOT (Dolphin 2000), with specific modules for ACS and WFC3-IR. We used Tiny Tim (Krist et al. 2011) point-spread functions (PSFs) and the general procedure and recommendations (e.g., DOLPHOT input parameters) outlined in Williams et al. (2014). The DOLPHOT output was culled to eliminate spurious and overly crowded sources using sharpness, crowding, and cuts in signal-to-noise ratio (S/N) of <1, <0.55, and >5, respectively. The remaining sources were cross-matched with the KS2 starlists, and the KS2 zero-points calculated from the average difference between the DOLPHOT and KS2 magnitudes over a magnitude range selected to omit bright saturated stars and faint noisy stars (Figure 3). The uncertainty on the average difference is less than 0.1% for all filters, and so the KS2 zero-point uncertainty is dominated by the reported HST zero-point uncertainty of 1% (Kalirai et al. 2009). The final zero-points (in magnitudes) are 32.6783, 25.2305, 23.566, 23.088, and 24.5698 for the F814W, F125W, F127M, F153M, and F160W filters, respectively.

Zoom In Zoom Out Reset image size

The proper motions provide a reliable method to separate likely cluster members from field stars. Following H15, a Gaussian mixture model was used to describe the kinematic distributions of the cluster and field, from which a cluster membership probability is calculated for each star. The resulting proper-motion catalog contains 9922 stars with membership probabilities and is presented in detail in J. R. Lu et al. (2018, in preparation).

Using the 4 m VISTA telescope at Cerro Paranal, the VVV survey mapped 520 deg² of the Milky Way bulge and disk in the K_s band (Minniti et al. 2010). This area was observed in 1.64 deg² tiles, each composed of six individual pointings dithered such that each pixel (except for those at the extreme edges of the tile) was covered by at least four images. The VVV photometric catalog for the tile containing Wd1 was downloaded from Data Release 2 using the ESO Phase 3 Archive Interface.⁸ However, the positions in the catalog revealed that there were only limited detections near Wd1, likely due to significant stellar crowding in the field. In addition, the VVV catalog contains aperture photometry, which can be compromised in crowded regions such as Wd1. As a result, we performed PSF photometry on the VISTA tile directly to improve the measurements in the Wd1 field.

After downloading the VISTA K_s tile image using the Phase 3 Archive Interface, we trimmed the tile to a manageable region (42 × 42) that overlaps our HST observations. To perform PSF photometry we used AIROPA, an expansion of the Starfinder code (Diolaiti et al. 2000), in the legacy mode (Witzel et al. 2016). In this mode, AIROPA behaves identically to Starfinder v1.6. Briefly, an empirical PSF is derived from a subset of user-identified stars and is then cross-correlated with the image to detect stars above a user-defined threshold. We carefully selected 27 relatively isolated and non-saturated stars in the tile to define the model PSF and set the minimum correlation coefficient of 0.7. This results in the detection of ∼6000 stars in the image.

We calibrate the PSF photometry by matching sources between the PSF and VVV catalogs. A 2D first-order polynomial is used to transform the VVV positions into the PSF catalog astrometric reference frame, obtaining 1300 stars with positions matched within 1'' (the FWHM of the image). Only stars with photometric errors smaller than 0.05 mag in both catalogs and fainter than the saturation limit are considered. Furthermore, we require that each star's VVV 1'' aperture diameter magnitude ("APERMAG1") be consistent with its 2'' aperture diameter magnitude ("APERMAG3") to within 0.05 mag, in an effort to eliminate stars with close neighbors affecting the aperture photometry. After these cuts, 136 stars remain. The photometric zero-point for the PSF catalog is calculated from the median difference between the VVV magnitude and the instrumental PSF magnitude (Figure 3). The final K_s zero-point is 24.566 ± 0.012 mag, where the uncertainty is the median difference error (0.005 mag) combined in quadrature with the VVV zero-point error reported for the tile image (0.011 mag).

It has been shown that the photometric errors reported by Starfinder (and thus AIROPA in single-PSF mode) are systematically underestimated because they do not capture errors in the PSF model itself (Schödel 2010). This PSF uncertainty is largely constant with magnitude, usually dominating the error budget for bright stars where the photon noise is very low. We derive this additional error term from the same sample of matched stars that was used to derive the photometric zero-point. We assume that the standard deviation of the difference between the PSF and VVV catalog magnitudes (σ_Δ) is a combination of the VVV catalog error (σ_VVV), PSF error reported by AIROPA (σ_A), and the constant PSF error term (σ_PSF):

$\begin{eqnarray}&&{\sigma }_{{\rm{\Delta }}}^{2}={\sigma }_{\mathrm{VVV}}{(m)}^{2}+{\sigma }_{A}{(m)}^{2}+{\sigma }_{\mathrm{PSF}}^{2},\end{eqnarray} \tag{ 1 }$

where σ_VVV and σ_A are both functions of the magnitude m. We find σ_PSF = 0.058 mag, which is added in quadrature with the reported AIROPA errors to produce the final photometric errors for the PSF catalog.

The Arches cluster was observed with HST WFC3-IR using the F127M and F153M filters in 2010 (GO-11671; PI: A.M. Ghez), and then repeat observations in F153M were obtained in 2011 and 2012 for astrometric purposes (GO-12318, GO-12667; PI: A.M. Ghez). These observations, along with the photometric and astrometric measurements extracted using KS2, are presented in H15. The photometric zero-points are derived in the same manner as the Wd1 HST observations and are also shown in Figure 3. We use the same catalog as has been presented in H15, which contains ∼26,000 stars.

To model the Wd1 and RC observations, we use stellar models to represent the stellar population, apply extinction to their spectra using a generated extinction law, and then calculate synthetic photometry in the observed filters. To generate the Wd1 stars, we must first adopt a cluster age. Spectroscopic studies of the evolved star population suggest a cluster age between 4 and 6 Myr, based on the properties of the O-type supergiants (Negueruela et al. 2010) and the observed ratio of different populations of supergiants to Wolf–Rayet stars (Crowther et al. 2006). Analysis of the pre-main sequence (pre-MS) turn-on feature in the color–magnitude diagram (CMD) has found cluster ages between 3 and 5 Myr (Brandner et al. 2008; Gennaro et al. 2011), but the age is degenerate with the cluster distance (e.g., Andersen et al. 2017), which is uncertain to ±700 pc (∼18%; Kothes & Dougherty 2007). None of these studies find evidence of an age spread in the cluster, and Clark et al. (2005) point out that the large number of massive stars would be expected to remove excess gas in the cluster within a crossing time. We thus assume that the cluster is coeval and adopt an age of 5 Myr. We will perform separate analyses assuming ages of 4 and 6 Myr to assess the impact of the age uncertainty on the extinction law.

A theoretical Wd1 cluster isochrone is generated with the adopted age and assuming solar metallicity using Geneva evolution models with rotation (Ω = 0.4; Ekström et al. 2012) for the main-sequence (MS)/evolved stars and Pisa evolution models (Tognelli et al. 2011) for the pre-MS. For each star in the isochrone, the effective temperature, T_eff, and surface gravity, log(g), from the stellar evolution model are used to select a model atmosphere. The atmospheres are drawn from a combined library of ATLAS9 (Castelli & Kurucz 2004) and PHOENIX (version 16; Husser et al. 2013) synthetic spectra. ATLAS9 spectra are used for T_eff > 5500 K and PHOENIX spectra are used for T_eff < 5000 K, with a linear interpolation of the two model sets between 5000 and 5500 K. The spectra are then reddened using a custom-built extinction law (Section 3.4.1) and total extinction using Pysynphot (STScI Development Team 2013). The reddened spectra are convolved with the desired filter functions to produce synthetic photometry. The output magnitudes can then be scaled to any cluster distance (a free parameter in our model) for a direct comparison to the observations.

As discussed in Section 3.2, we limit the observed Wd1 sample to only stars that fall on the MS. These stars represent an ideal sample to examine the extinction law because their intrinsic colors are well understood and are only slightly dependent on mass in the filters used in this study. We restrict the stellar mass range of the theoretical isochrone to match the selection criterion imposed on the observed sample, selecting stars between 0.5 and 3.0 mag brighter than the pre-MS bridge in F160W (Figure 4). The pre-MS bridge, which connects the MS and pre-MS sequences in a CMD, is only dependent on the cluster age. Thus, no assumptions regarding the cluster distance, extinction law, or total extinction are required in order to appropriately set the isochrone mass range. For a 5 Myr cluster this mass range is 4.41–14.0 M_⊙. To test whether the ATLAS9 atmospheres (which assume local thermodynamic equilibrium) are appropriate for the high-mass stars in our sample, we compare the synthetic magnitudes for a 9, 12, and 15 M_⊙ MS star with those calculated using non-local thermodynamic equilibrium CMFGEN atmospheres from Fierro et al. (2015). We find that the synthetic magnitudes agree to within the photometric errors, and so we move forward with the ATLAS9 models in our analysis.

Zoom In Zoom Out Reset image size

A similar procedure is used to model the Arches field RC stars. We select the RC star model from a 10 Gyr PARSEC stellar isochrone (Bressan et al. 2012) at solar metallicity, which represents the average stellar population in the Galactic bulge (Zoccali et al. 2003; Clarkson et al. 2008). The chosen model matches the average effective temperature (T_eff) and surface gravity (log(g)) measured for solar-metallicity RC stars in the Hipparcos catalog (T_eff = 4700 K, log g = 2.40 cgs; Mishenina et al. 2006). The synthetic photometry is then scaled to the average distance of the RC population, which is a free parameter in our model.

We select a sample of high-probability Wd1 members that fall on the MS, are not saturated, and have small photometric errors. Each star in the sample has photometry in each HST filter, and a subset of the stars have VISTA K_s photometry as well. The steps used to create the sample are described below.

The sample is created directly from the HST proper-motion catalog described in Section 2.1 and J. R. Lu et al. (2018, in preparation). First, we restrict the catalog to stars with cluster membership probabilities greater than 0.6 and photometric errors smaller than 0.05 mag in each filter. MS stars are identified as those 0.5 mag brighter than the pre-MS bridge in F160W, corresponding to F160W ≤ 15.0 mag (Figure 4). This conservative criterion is adopted to minimize contamination from pre-MS stars scattered into the MS by differential reddening. At the bright end, we adopt a magnitude cut of F160W ≥ 12.5 mag in order to eliminate saturated sources in the VISTA K_s observations. After these cuts, 537 of the original 9922 stars remain in the sample.

To eliminate photometric outliers (such as field stars or binary systems with unusual colors), we apply an iterative 3-σ cut in a two-color diagram (2CD) across the HST filters: F814W–F125W versus F814W–F160W. We fit a line to the sample (via orthogonal regression) and calculate the root-mean-squared (rms) residual relative to the fit in 0.25 magnitude bins. Stars with residuals larger than three times the rms value in their magnitude bin are removed. This process is repeated until no further stars are eliminated. A total of 53 stars are rejected by this criterion.

A final cut is required to eliminate cluster stars that would otherwise be outside the adopted F160W magnitude range but have been scattered into the sample by differential extinction. For example, a high-mass star intrinsically brighter than the F160W magnitude limit can scatter into the sample if it is located in a region of higher extinction. This star would appear redward of the average MS population. Similarly, a low-mass star could scatter into the sample if it is in a region of lower extinction, placing it blueward of the average MS population. To eliminate these stars, we make a cut in F160W versus F814W–F160W CMD space, where the median color is taken to represent the MS. Two lines with a slope of 1 that intersect the bright and faint ends of the MS are used to identify and remove potentially scattered stars (Figure 5). This slope is a conservative estimate of the steepest possible reddening vector in the CMD, corresponding to A_F814W/A_F160W ∼ 2. A steeper reddening vector would require A_F814W/A_F160W < 2, which is much lower than any law reported in the literature. An additional 42 stars are removed by this criterion, resulting in a final sample size of 453 stars.

Zoom In Zoom Out Reset image size

To construct the HST–VISTA subsample, we iteratively match the HST proper-motion catalog with the VISTA catalog using a 0.25'' (∼2 HST pixels) matching radius. An additional "isolation cut" is imposed to reduce the impact of stellar crowding on the seeing-limited VISTA photometry: stars lying within 45 of another star that has a brightness within 3 mag of the star itself (as identified from the HST catalog) are rejected. All of the same cuts are applied for a final HST–VISTA subsample of 106 stars.

A summary of the adopted cuts and sample size is provided in Table 2. These stars form well-defined reddening vectors in the HST-only and HST–VISTA 2CDs, which are used in the extinction law analysis (Figure 6). The typical photometric uncertainty is ∼0.02 mag for the HST filters and ∼0.07 mag for VISTA K_s. While the magnitude cuts introduce a possible Malmquist bias, the effect is small due to the small photometric errors. Furthermore, a bias would only affect our results if the extinction law for the faint stars were different than that of the bright stars, which is highly unlikely since they are part of the same cluster.

Zoom In Zoom Out Reset image size

Red clump stars are low-mass giants that are in the core helium-burning stage of their evolution. Exhibiting a narrow range of T_eff and log(g) values, these stars form a well-defined clump in the CMD that spreads along the reddening vector in the presence of differential extinction (see Girardi 2016, for a review). While no significant RC star population is found in the Wd1 observations, H15 found a large RC field population in NIR HST observations of the Arches cluster. These stars are associated with the Galactic Bulge and have a distance distribution that peaks close to the GC along this sight-line. Since the GC distance is known to 2% (Boehle et al. 2016), the normalization of the extinction law to these stars can be constrained to a much higher precision than is possible with Wd1.

Red clump stars are visible in the Arches field CMD as a high-density "bar" (Figure 7). To identify the RC stars, we use the unsharp-masking technique described by De Marchi et al. (2016). This method increases the contrast of high-frequency features while reducing the contrast of low-frequency ones. First, we convert the F153M versus F127M–F153M CMD into a Hess diagram (i.e., a 2D histogram of stellar density), treating the position of each star as a Gaussian probability distribution with a width corresponding to the photometric error. A bin size of 0.05 mag in both color and mag space is used. Next, we create the unsharp mask by convolving the Hess diagram with a 2D Gaussian kernel having a width of 0.2 mag. The mask is then subtracted from the original Hess diagram to create the unsharp-masked diagram. We calculate a linear fit to the high-density "ridge" created by the RC population and identify RC stars as those that fall within F153M ±0.3 mag of the best-fit line (Figure 7). A total of 1119 RC stars are identified in this manner.

Zoom In Zoom Out Reset image size

Similar to the Wd1 sample, we adopt a series of cuts to ensure the quality of the RC sample for the extinction law analysis. We require each star to have a photometric error smaller than 0.05 mag in both filters and perform a 3σ iterative outlier rejection in the CMD. The final RC sample contains 819 stars with typical photometric uncertainties of 0.015 mag.

The expected distance distribution of the sample is calculated using the RC density profiles of Wegg & Gerhard (2013), which are measured for the major, intermediate, and minor axes of the Galactic bar. These profiles are rotated by 28° (the measured angle of the Galactic bar relative to the GC; Wegg & Gerhard 2013) to determine the density profile of the Arches LOS. The expected number of RC stars is then calculated as a function of distance by multiplying the solid angular area of the observations by the LOS density profile. The resulting distribution is nearly Gaussian with a peak at 96 pc beyond the GC and a width of 630 pc (Figure 8). The peak in the observed counts is not centered at the GC due to the increase in projected field area with LOS distance.

Zoom In Zoom Out Reset image size

The RC star distance distribution is robust against uncertainties in the Galactic bar rotation angle (α) and the RC density profiles. Varying α between 25° and 33° (the possible range reported by Wegg & Gerhard 2013 and Wegg et al. 2015) only shifts the peak of the distribution by 12 pc. Redrawing the RC density profiles of Wegg & Gerhard (2013) while applying the reported uncertainty of 10% on each measurement results in a typical shift in the distribution of just ±20 pc. Both sources of error are well within the uncertainty in the GC distance itself (140 pc) and are not considered in the analysis.

In the CMD, we expect the distribution of F153M residuals relative to the RC ridge to be driven by the variation in stellar distance across the population. To test this, we compare the observed residuals to those calculated for a synthetic RC population created with the predicted distance distribution in Figure 8. Details on how the synthetic population is created are provided in Appendix B. The observed F153M residual distribution is wider than the residual distribution for the synthetic data by ∼0.03 mag (Figure 9). This indicates that either the RC distance distribution is slightly wider than expected, or that there is an additional source of magnitude dispersion among the RC stars. Given that the synthetic RC population assumes a single age and metallicity, it is likely that variations in stellar age or metallicity within the bulge also affect the observed F153M residual distribution (e.g., Salaris & Girardi 2002; Chen et al. 2017).

Zoom In Zoom Out Reset image size

3.4.1. Extinction Law Model and Priors

To construct the extinction law, we define A_λ/A_Ks at five specific wavelengths and use a cubic B-spline to interpolate the law across all wavelengths. The extinction law is thus a continuous function that can be adjusted by changing the A_λ/A_Ks values at the wavelength points, which are free parameters in our model. There are 10 additional free parameters in the model: 4 for the global Wd1 and RC population parameters (Wd1 total extinction A_Ks; Wd1 distance d_wd1; average RC distance d_rc; Gaussian width of F153M residuals around the reddening vector σ_rc), and 6 for systematic offsets in the photometric zero-points of the filters (ΔZP_λ). The parameters and adopted priors are discussed below and summarized in Table 3.

Table 3.  Model Parameters and Priors

Parametera	λb (μm)	Priorc	Units	Prior Reference
${A}_{{\rm{F}}814{\rm{W}}}/{A}_{{K}_{s}}$	0.806	U(4, 14)	⋯	⋯
${A}_{y}/{A}_{{K}_{s}}$	0.962	U(4, 14)	⋯	⋯
${A}_{{\rm{F}}125{\rm{W}}}/{A}_{{K}_{s}}$	1.25	U(1, 6)	⋯	⋯
${A}_{{\rm{F}}160{\rm{W}}}/{A}_{{K}_{s}}$	1.53	U(1, 6)	⋯	⋯
${A}_{[3.6]}/{A}_{{K}_{s}}$	3.545	G(0.5, 0.05)	⋯	Nishiyama et al. (2009)
${A}_{{K}_{s}}$	2.14	U(0.3, 1.3)	mag	⋯
dwd1	⋯	G(3900, 700)	pc	Kothes & Dougherty (2007)
drc	⋯	G(7960, 140)	pc	Boehle et al. (2016), Section 3.3
σrc	⋯	G(0.17, 0.01)	mag	Section 3.3
ΔZPKs	⋯	G(0, 0.012)	mag	Section 2.2
ΔZPF160W	⋯	G(0, 0.01)	mag	Section 2.1
ΔZPF153M	⋯	G(0, 0.01)	mag	Section 2.1
ΔZPF127M	⋯	G(0, 0.01)	mag	Section 2.1
ΔZPF125W	⋯	G(0, 0.01)	mag	Section 2.1
ΔZPF814W	⋯	G(0, 0.01)	mag	Section 2.1

Notes.

^aA_λ/A_Ks: extinction law in filter; A_Ks: total extinction of Wd1; d_wd1: distance to Wd1; d_rc: average RC star distance; σ_rc: Gaussian width of the RC F153M residuals around the reddening vector; and ΔZP_λ: zerop-oint offset in filter. ^bHST + PanSTARRS filters: Pivot wavelengths of filter; IRAC [3.6] filter: isophotal wavelength from Nishiyama et al. (2009). ^cUniform distributions: U(min, max), where min and max are bounds of the distribution; Gaussian distributions: G(μ, σ), where μ is the mean and σ is the standard deviation.

Download table as:  ASCIITypeset image

In the extinction law, a wavelength point is assigned to the pivot wavelength (Tokunaga & Vacca 2005) of each Wd1 filter (F814W, F125W, and F160W). The A_λ/A_Ks values at these points are given uniform priors. Two additional wavelength points are added for the SST/IRAC [3.6] and PanSTARRS y filters, although no observations at these wavelengths are used. Our analysis revealed that the y-band point is required to capture the needed curvature in the extinction law between the F125W and F814W filters (Section 5.3). We adopt a uniform prior for A_λ/A_Ks at this point as well. The [3.6] point is included to enforce reasonable behavior through the red edge of the extinction law. We adopt a Gaussian prior using the value of Nishiyama et al. (2009) with a conservative uncertainty of 10% (A_3.5/A_Ks = 0.5 ± 0.05). The cubic B-spline interpolation is calculated using the scipy.interpolate.splrep function in python between 0.8 and 2.2 μm. The exact function call is provided in the stand-alone python code referenced in Section 4.3. This is converted into a pysynphot custom reddening law for the synthetic photometry. Since each A_λ/A_Ks point is allowed to vary independently, no assumption regarding the functional form of the extinction law is made.

For the global Wd1 population parameters, we adopt a uniform prior for A_Ks and a Gaussian prior of 3900 ± 700 pc for d_wd1. The distance constraint is derived from the kinematics of HI gas associated with the cluster (Kothes & Dougherty 2007). Although additional distance measurements exist from an eclipsing binary analysis (Koumpia & Bonanos 2012), spectrally typed evolved stars (e.g., Negueruela et al. 2010), and CMD fitting of pre-MS stars (e.g., Andersen et al. 2017), these analyses must correct for extinction and thus depend on the extinction law.

For the RC stars, we adopt a Gaussian prior of 7960 ± 140 pc for d_rc. This is 100 pc beyond the GC distance measurement of Boehle et al. (2016), matching the predicted average population distance in Section 3.3. For σ_rc, we adopt a prior of 0.17 ± 0.01 mag, corresponding to the measured width found in Figure 9.

The zero-point offsets are included in the model since errors in the zero-points would propagate through the analysis as a systematic error rather than as a random one. The offsets are assigned Gaussian priors centered at zero with a width corresponding to the uncertainty in the zero-point derivation in Section 2.1.

3.4.2. Wd1 Star Likelihood

For the Wd1 component of the likelihood, the observed sample is compared to the theoretical isochrone produced by the model in terms of the F160W magnitude and F160W–F814W, F160W–F125W, and F160W–K_s colors. The full sample is used for the HST colors and the HST–VISTA subsample is used for F160W–K_s. Since a smaller number of stars have VISTA photometry, the HST colors have greater weight in the extinction law fit. However, this is desirable as the K_s observations have larger scatter than the HST observations.

Initially, the observed photometric errors are smaller than the magnitude sampling of the isochrone. We address this by performing a cubic spline interpolation of the isochrone magnitudes as a function of stellar luminosity and resample in steps of 2.5 × 10⁻³ in log luminosity. This results in a magnitude spacing of 3.7 × 10⁻³ mag on the isochrone, which is more than twice smaller than the typical magnitude error (0.01 mag). However, there are a handful of stars with mag errors below this threshold. Since finer isochrone sampling significantly increases the computation time of the analysis, we instead add an error floor of 3.7 × 10⁻³ mag in quadrature to the observations to avoid this problem.

A single isochrone is insufficient to reproduce the data due to the differential extinction (${{dA}}_{{K}_{s}}$) in the field. Instead, we generate a grid of isochrones with a range of extinction values ±0.6 mag from the input ${A}_{{K}_{s}}$ in steps of 5 × 10⁻⁴ mag. This ensures that the color sampling between isochrones is at least twice smaller than the photometric uncertainty for each color in the model.

For each observed Wd1 MS star, we identify the nearest-neighbor synthetic star in the multidimensional magnitude-color space described above. We define the set of observed magnitudes and colors as ${\boldsymbol{m}}$ and their associated errors ${\boldsymbol{\sigma }}$. The likelihood is

$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{Wd}1}({\boldsymbol{m}},{\boldsymbol{\sigma }})=\displaystyle \prod _{j=1}^{4}{{ \mathcal L }}_{j}({{\boldsymbol{m}}}_{j},{{\boldsymbol{\sigma }}}_{j})\end{eqnarray} \tag{ 2 }$

where ${{\boldsymbol{m}}}_{j}$ and ${{\boldsymbol{\sigma }}}_{j}$ are the measurements and errors in the jth dimension, with the dimensions corresponding to F160W, F160W–F814W, F160W–F125W, and F160W–K_s. For each observed star, the nearest-neighbor synthetic star across all dimensions is found. The likelihood for each dimension is then calculated by comparing the observed sample to their corresponding nearest neighbors:

$\begin{eqnarray}&&{{ \mathcal L }}_{j}({{\boldsymbol{m}}}_{j},{{\boldsymbol{\sigma }}}_{j})=\displaystyle \prod _{i=0}^{{N}_{s,j}}\displaystyle \frac{1}{{\sigma }_{j,i}\sqrt{2\pi }}{e}^{\tfrac{{({m}_{j,i}-{{NN}}_{\mathrm{mod}}({m}_{j,i}))}^{2}}{2{\sigma }_{j,i}^{2}}},\end{eqnarray} \tag{ 3 }$

where N_s,j is the number of stars in the jth dimension, m_{j, i} and ${\sigma }_{j,i}$ are the mag/color and corresponding error of the ith star in the jth dimension, and ${{NN}}_{\mathrm{mod}}({m}_{j,i})$ is the mag/color of the ith star's nearest-neighbor synthetic star in the model. The full sample is used for the F160W, F160W–F814W, and F160W–F125W dimensions, while only the HST–VISTA subsample is used for the F160W–K_s dimension.

3.4.3. Red Clump Star Likelihood

We calculate the RC star likelihood component in a similar manner, comparing the observed sample to the reddening vector from the extinction law model in color–magnitude space (F153M and F127M–F153M). With ${{\boldsymbol{m}}}_{r}$ and ${{\boldsymbol{c}}}_{r}$ representing the set of F153M magnitudes and F127M–F153M colors with errors ${{\boldsymbol{\sigma }}}_{{{\boldsymbol{m}}}_{r}}$ and ${{\boldsymbol{\sigma }}}_{{{\boldsymbol{c}}}_{r}}$:

$\begin{eqnarray}&&{{ \mathcal L }}_{\mathrm{RC}}={{ \mathcal L }}_{{m}_{r}}({{\boldsymbol{m}}}_{r},{{\boldsymbol{\sigma }}}_{{{\boldsymbol{m}}}_{r}},{\sigma }_{\mathrm{rc}})\ast {{ \mathcal L }}_{{c}_{r}}({{\boldsymbol{c}}}_{r},{{\boldsymbol{\sigma }}}_{{{\boldsymbol{c}}}_{r}}),\end{eqnarray} \tag{ 4 }$

where σ_RC is the free parameter in the model corresponding to the Gaussian width of the F153M residuals around the reddening vector. Each component of the likelihood is defined as

$\begin{eqnarray}&&{{ \mathcal L }}_{{m}_{r}}({{\boldsymbol{m}}}_{r},{{\boldsymbol{\sigma }}}_{{{\boldsymbol{m}}}_{r}},{\sigma }_{\mathrm{rc}})=\displaystyle \prod _{i=0}^{{N}_{r}}\displaystyle \frac{1}{\sqrt{{\sigma }_{{m}_{r},i}^{2}+{\sigma }_{\mathrm{rc}}^{2}}\sqrt{2\pi }}{e}^{\tfrac{{({m}_{{m}_{r},i}-{{NN}}_{\mathrm{mod}}({c}_{r},i))}^{2}}{2({\sigma }_{{m}_{r},i}^{2}+{\sigma }_{\mathrm{rc}}^{2})}}\end{eqnarray} \tag{ 5 }$

$\begin{eqnarray}&&{{ \mathcal L }}_{{c}_{r}}({{\boldsymbol{c}}}_{r},{{\boldsymbol{\sigma }}}_{{{\boldsymbol{c}}}_{r}})=\displaystyle \prod _{i=0}^{{N}_{r}}\displaystyle \frac{1}{{\sigma }_{{c}_{r},i}^{2}\sqrt{2\pi }}{e}^{\tfrac{{({c}_{r,i}-{{NN}}_{\mathrm{mod}}({c}_{r}))}^{2}}{2{\sigma }_{{c}_{r},i}^{2}}},\end{eqnarray} \tag{ 6 }$

where ${{NN}}_{\mathrm{mod}}({c}_{r},i)$ is the nearest-neighbor model star in color space for the ith star and N_r is the total number of RC stars in the sample. Note that the nearest neighbor is only calculated for color space because of the extra dispersion in the magnitudes, as discussed in Section 3.3. The extra dispersion is captured by the σ_rc parameter in the F153M dimension of the likelihood, which manifests as an extra error term in addition to the photometric errors. Since the dispersion is primarily driven by individual RC distance variation, it is not included in the color likelihood term.

3.4.4. Final Likelihood

The final likelihood combines the Wd1 and RC likelihood components. For each to have an equal weight, we must account for the fact that the RC sample is significantly larger than the Wd1 sample, which causes ${{ \mathcal L }}_{\mathrm{rc}}\gt {{ \mathcal L }}_{\mathrm{wd}1}$ regardless of the quality of the fit. After converting into log-likelihood, we scale the RC likelihood component to the same number of stars as the Wd1 sample (N_wd1):

$\begin{eqnarray}&&\mathrm{log}{{ \mathcal L }}_{\mathrm{tot}}=\mathrm{log}{{ \mathcal L }}_{\mathrm{wd}1}+\mathrm{log}{{ \mathcal L }}_{r}\ast \displaystyle \frac{{N}_{\mathrm{wd}1}}{{N}_{\mathrm{rc}}}.\end{eqnarray} \tag{ 7 }$

With the likelihood function defined, we determine the best-fit extinction law model using Bayes' theorem:

$\begin{eqnarray}&&P({\boldsymbol{M}}| {\rm{\Theta }})=\displaystyle \frac{{{ \mathcal L }}_{\mathrm{tot}}\ast P({\boldsymbol{M}})}{P({\rm{\Theta }})},\end{eqnarray} \tag{ 8 }$

where ${\boldsymbol{M}}$ is the model with the set of free parameters $[{{\boldsymbol{A}}}_{{\boldsymbol{\lambda }}}/{{\boldsymbol{A}}}_{{Ks}},{\boldsymbol{\Delta }}{\boldsymbol{ZP}},{A}_{{Ks}},{d}_{\mathrm{wd}1},{d}_{\mathrm{rc}},{\sigma }_{\mathrm{rc}}$](with ${{\boldsymbol{A}}}_{{\boldsymbol{\lambda }}}/{{\boldsymbol{A}}}_{{Ks}},{\boldsymbol{\Delta }}{\boldsymbol{ZP}}$ representing the set values at five and six different wavelengths, respectively) and ${\rm{\Theta }}=[{\boldsymbol{m}}$, ${\boldsymbol{\sigma }}$, ${{\boldsymbol{m}}}_{r}$, ${{\boldsymbol{\sigma }}}_{r}$, ${{\boldsymbol{c}}}_{r}$, ${{\boldsymbol{\sigma }}}_{{{\boldsymbol{c}}}_{r}}]$ is the set of observations. $P({\boldsymbol{M}}| {\boldsymbol{\Theta }})$ is the posterior probability of the model given the data, and P(${\boldsymbol{M}}$) is the prior probability on ${\boldsymbol{M}}$.

The posterior probability distributions are calculated using Multinest, a nested sampling algorithm shown to be more efficient than Markov chain Monte Carlo algorithms in complex parameter spaces with multiple modes or pronounced degeneracies (Feroz & Hobson 2008; Feroz et al. 2009). This iterative technique calculates the posterior probability at a fixed number of points in the parameter space and identifies possible peaks, restricting subsequent sampling to the regions around these peaks until the change in evidence drops below a user-defined tolerance level. An evidence tolerance of 0.5 and sampling efficiency of 0.8 are adopted for 400 active points. To run the sampler, we use Pymultinest, a convenient wrapper module in python (Buchner et al. 2014).

3.4.5. Testing the Analysis

To fit the Wd1 MS sample, we set ${{ \mathcal L }}_{\mathrm{tot}}={{ \mathcal L }}_{\mathrm{wd}1}$ in Equation (8) and do not consider the model parameters related to the RC stars (d_rc, σ_rc, ZP_F127M, and ZP_F153M). A cluster age of 5 Myr is assumed. As with the simulated cluster tests in Appendix B, a large degeneracy is found between d_wd1, ${A}_{{K}_{s}}$ and the extinction law due to the uncertainty in the normalization, although the shape of the law is well constrained (Figure 10). d_wd1 is constrained to 4428 ± 309 ± 48 pc, consistent with the result from Kothes & Dougherty (2007), and ${A}_{{K}_{s}}$ is found to be 0.78 ± 0.16 ± 0.02 mag, also broadly consistent with the literature. No statistically significant zero-point offsets are obtained for any of the filters, indicating that the zero-points presented in Figure 3 are accurate to 0.006 mag for the HST filters and 0.008 mag for the VISTA K_s filter. The extinction law provides a good fit to the observations, with reduced chi-square (${\chi }_{\mathrm{red}}^{2}$) values of 0.38 for the HST CCD and 0.77 for the HST–VISTA CCD (Figure 11). That the ${\chi }_{\mathrm{red}}^{2}$ values are lower than 1.0 suggests that the photometric errors are conservative.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The shape of the extinction law is inconsistent with a single power law, which would produce a straight line in S_1/λ versus log(1/λ) plot in Figure 10. We examine this further in the final extinction law fit (Section 4.3).

We test the impact of assuming an age for Wd1 by performing identical analyses with cluster ages of 4 and 6 Myr, the range of possible ages suggested by the evolved star population (Crowther et al. 2006; Negueruela et al. 2010). Since the stellar mass of the pre-MS bridge decreases with age, changing the age affects the mass range of the MS sample. Using the selection criteria described in Section 3.1, the MS mass range becomes 4.83–16.0 M_⊙ for a 4 Myr cluster and 3.95–12.50 M_⊙ for a 6 Myr cluster. This affects the normalization of the extinction law, with the law generally becoming steeper (i.e., higher A_λ/A_Ks values) with increasing age, but the shape of the law remains almost identical (Figure 12).

Zoom In Zoom Out Reset image size

Since the distance distribution of the Arches RC population is known to significantly higher precision than Wd1, the normalization of the extinction law can be much better constrained. We perform an extinction law fit using just the RC sample, setting ${{ \mathcal L }}_{\mathrm{tot}}={{ \mathcal L }}_{\mathrm{rc}}$ in Equation (8) and removing the free parameters related to Wd1 (i.e., d_wd1, ${A}_{{K}_{s}}$). Only the extinction law ratios ${A}_{{\rm{F}}125{\rm{W}}}/{A}_{{K}_{s}}$ and ${A}_{{\rm{F}}160{\rm{W}}}/{A}_{{K}_{s}}$ are included in the model, approximately matching the wavelengths of the F127M and F153M observations. We present the ratio A_F125W/A_F160W, which combines the information from both filters.

The best fit recovers A_F125W/A_F160W = 1.527 ± 0.006 ± 0.01 and d_rc = 7938 ± 87 ± 68 pc. The statistical uncertainty in A_F125W/A_F160W is a factor of ∼8 smaller than the value obtained for the Wd1-only fit (1.468 ± 0.047 ± 0.01), eliminating much of the uncertainty in the normalization. Statistically, d_rc is constrained to ∼1.1%, which is better than the input prior and in principle places a constraint on the GC distance. However, we have adopted a simplified description of the RC (assuming an average age and metallicity for the population) and are susceptible to possible systematic uncertainties in the RC star model (Section 5.4.1). A careful analysis of the GC distance is beyond the scope of this paper.

The final extinction law derived using the combined Wd1 MS + Arches RC sample is shown in Figure 13. Once again, a Wd1 age of 5 Myr is assumed. The law is well constrained with combined uncertainties (statistical and systematic added in quadrature) better than ∼5% in ${A}_{\lambda }/{A}_{{K}_{s}}$ and ∼10% in ${{ \mathcal S }}_{1/\lambda }$ (Table 4). There is no significant change in the shape of the law relative to the Wd1-only fit, which supports the assumption that the shape is the same for both populations. Fitting a power law to the A_λ/A_Ks values results in a reduced chi-squared (${\chi }_{\mathrm{red}}^{2}$) value of 3.7, indicating that the difference between the best-fit law and a power law are statistically significant (Figure 14). This is in agreement with previous studies of the OIR extinction law (e.g., Fitzpatrick & Massa 2009). Interestingly, we find that the NIR portion of the law (K_s, F160W, F125W) is also statistically inconsistent with a power law with ${\chi }_{\mathrm{red}}^{2}=3.3$, in contrast to what is often assumed in the literature. We only consider the statistical errors on A_λ/A_Ks in the ${\chi }_{\mathrm{red}}^{2}$ calculations. Extinction law analyses with various systematics applied (Section 5.4) also show statistically significant deviations from a power law, and so the systematics do not affect this conclusion. A stand-alone python code to calculate the Wd1+RC extinction law at all wavelengths between 0.8 and 2.2 μm is available online.⁹

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

A comparison of the best-fit model and Wd1 and Arches RC observations is shown in Figure 15. The law is able to reproduce the photometry of both populations well, with ${\chi }_{\mathrm{red}}^{2}$ values for the Wd1 CCDs are nearly identical to the Wd1-only fit (0.38 and 0.77, respectively) and ${\chi }_{\mathrm{red}}^{2}=0.83$ for the RC CMD.

Zoom In Zoom Out Reset image size

Repeated analyses of the extinction law assuming Wd1 ages of 4, 6, and 7 Myr show that the extinction law is not affected by changing the cluster age. However, d_wd1 is found to systematically decrease with increasing cluster age. We explore this trend further in Section 4.4.

To demonstrate the effect of the extinction law, we show how the new extinction law changes the distance and age of the cluster. In the Wd1+RC extinction law fit, d_wd1 varies from 5222 ± 113 pc for a Wd1 age of 4 Myr to 4133 ± 66 pc for a Wd1 age of 7 Myr (statistical and systematic errors added in quadrature). Thus, an independent estimate of d_wd1 offers a constraint on the cluster age. This is especially valuable given the diverse population of evolved stars in Wd1, which include yellow hypergiants (Clark et al. 2005), red supergiants (Clark et al. 2010), WR stars (Crowther et al. 2006), luminous blue variables (Clark & Negueruela 2004; Dougherty et al. 2010), and a magnetar (Muno et al. 2006). The presence of these objects provides a strong test of stellar evolution models, which struggle to reproduce such a sizable population of cool supergiants and WR stars simultaneously (Clark et al. 2010; Ritchie et al. 2010).

We obtain an independent estimate of d_wd1 from published measurements of W13, a 9.2 day eclipsing binary within the cluster (Bonanos 2007). Koumpia & Bonanos (2012; hereafter K12) combine the optical (VRI) lightcurves from Bonanos (2007) with multi-epoch spectroscopy to derive the physical properties of the system, found to be a near-contact binary composed of a B0.5Ia+/WNVL and an O9.5-B0.5I star. From the derived effective temperatures and stellar radii, they calculate a total system luminosity of log L/L_⊙ = 5.54 ± 0.11. K12 correct for extinction by adopting A_J/A_Ks = 2.50 ± 0.15 (Indebetouw et al. 2005), ultimately reporting a distance of 3710 ± 550 pc. We repeat this calculation using the Wd1+RC extinction law, which is steeper than the Indebetouw value (A_J / A_Ks = 3.56 ± 0.15, with statistical and systematic errors added in quadrature). This does not bias the Wd1 distance calculation since the law is independent of cluster age and thus not tied to the value of d_wd1 derived in the extinction law analysis.

Following K12, we calculate a distance to W13 using its 2MASS J-band apparent magnitude, a theoretical bolometric correction (BC_λ) for O9.5I stars, and an extinction correction based on the NIR colors of nearby WR stars. The 2MASS J-band magnitude is 9.051 ± 0.16 mag, with the uncertainty set by the depth of the primary eclipse since the phase of the measurement is unknown. From Martins & Plez (2006), we adopt BC_J = −3.24 ± 0.08 mag, with the uncertainty based on the scatter in the BC_λ–T_eff relation. The total extinction (A_Ks) of W13 is calculated from the J–K color excesses of the WR stars "R" and "U" reported in Crowther et al. (2006), both of which are ∼6'' away from W13. Using our extinction law, these stars have A_Ks = 0.55 ± 0.05 mag and 0.52 ± 0.05 mag, respectively, and so we adopt A_Ks = 0.535 ± 0.035 mag for W13.

With all the pieces in place, we can calculate a distance to W13:

$\begin{eqnarray}&&\mu ={m}_{J}-{M}_{J}-{A}_{J},\end{eqnarray} \tag{ 10 }$

where μ is the distance modulus, m_J is the 2MASS J-band apparent magnitude, A_J is the total extinction at J band, and M_J is the J-band absolute magnitude. The J-band absolute magnitude is calculated from the luminosity using the bolometric correction

$\begin{eqnarray}&&{M}_{J}={M}_{\odot }^{\mathrm{bol}}-{\mathrm{BC}}_{J}-2.5\mathrm{log}\left(\displaystyle \frac{L}{{L}_{\odot }}\right),\end{eqnarray} \tag{ 11 }$

where ${M}_{\odot }^{\mathrm{bol}}=4.75$ is the bolometric magnitude of the Sun and L_⊙ is the solar luminosity. Plugging in the values and propagating the errors, we find μ = 12.958 ± 0.235 mag, or 3905 ± 422 pc. The major source of remaining uncertainty is the 2MASS photometry, since the phase of the system was unknown at the time of measurement. New multi-epoch OIR photometry of W13 would dramatically increase the precision of the derived distance. Unfortunately, W13 is saturated in the HST observations, and so our observations are not helpful in this regard.

A comparison between the W13 distance and the d_wd1 from the extinction law analysis is shown in Figure 16. The eclipsing binary distance is consistent with an older cluster age, being discrepant from the 4 and 5 Myr extinction fit distances by 3.0σ and 2.0σ, respectively. The 6 and 7 Myr distances are only discrepant by 1.2σ and 0.5σ, respectively. An older cluster age for Wd1 is not unreasonable, given that the constraint of 4–5 Myr from Crowther et al. (2006) is based stellar evolution models of the WR star population, which are not yet well understood. Negueruela et al. (2010) find that the HR diagram of OB supergiants is consistent with a cluster age greater than 5 Myr, although this analysis relies on the Rieke & Lebofsky (1985) extinction law (A_J /A_Ks = 2.35). Unfortunately, their observations are primarily shortward of I band, and so we cannot recreate their HR diagram using our extinction law. We leave a detailed analysis of the age and distance of Wd1 to a future paper.

Zoom In Zoom Out Reset image size

We compare our result with several extinction laws in the literature. These include Cardelli et al. (1989; hereafter C89), assuming R_V = 3.1 as is often adopted for the interstellar medium (however, changing R_V has little effect on the extinction law in this wavelength range); Nishiyama et al. (2009; hereafter N09), often used for the GC and Galactic bulge; Damineli et al. (2016; hereafter D16), derived specifically for Wd1; Fitzpatrick & Massa (2009; hereafter F09), assuming α = 2.5 and R_V = 3.1; and Schlafly et al. (2016; hereafter S16), assuming A_H/A_Ks ("rhk") values of 1.55 (consistent with Indebetouw et al. 2005) and 2.0. To construct the C89, D16, and F09 laws, we use the functional forms reported by the authors, and for S16, we use the python code referenced in their appendix.¹⁰ For N09, we adopt a power law with β = 2.0 between 1.25 and 2.14 μm and then use a linear interpolation in $\mathrm{log}({A}_{\lambda }/{A}_{{Ks}}$) versus log(1/λ) space to extend from A_J/A_Ks to the A_V/A_Ks value reported in Nishiyama et al. (2008). This method was adopted to have minimal impact on the shape of the N09 law shortward of J band, where no functional form is provided.

Theoretical cluster isochrones using the published extinction laws are unable to reproduce the observed colors of the Wd1 MS (Figure 17). This is especially evident in the HST2CD, where the model reddening vectors are offset between 0.29 and 0.72 mag to the blue in F814W–F125W (about 9% and 22%, respectively), and between 0.08 and 0.19 mag to the red in F125W–F160W (about 12% and 25%, respectively). These offsets are huge relative to the HST photometric errors, which are typically 0.02 mag for both colors. In the HST–VISTA 2CD, most literature extinction laws produce reddening vectors that are 0.07–0.10 mag too blue in F160W–K_s (about 12%), which is significant but not as large relative to the typical color uncertainty of ∼0.06 mag. The exception is C89, which reproduces the observed sequence in this 2CD well. These isochrones assume a cluster age of 5 Myr, although changing the age has little to no effect on the colors. Note that the discrepancies in color–color space are due to differences in the shape of the extinction law, rather than the normalization.

Zoom In Zoom Out Reset image size

In Figures 10, 13, and Table 4 we compare our extinction law results to the published laws. The Wd1+RC extinction law is generally steep (i.e., has high A_λ/A_Ks values), being most similar to F09 law with α = 2.5 and the S16 law with A_H/A_Ks = 2.0. Notably, we find A_F125W/A_Ks to be 18% larger and A_F814W/A_Ks to be 24% larger than the N09 law that is commonly used for the inner Milky Way. The shape of the law also differs from those in the literature, which is expected given the inability of the published laws to reproduce the observations in Figure 17. In particular, ${S}_{{\rm{F}}814{\rm{W}}}$/S_Ks is significantly larger than any of the published laws, indicating a steeper derivative through the F814W filter relative to the K_s filter. However, our data calls for such steepness as discussed in Section 5.3.

Our extinction law differs from the one previously derived for Wd1 by D16, who describe the law as a power law with a slope β = 2.13 ± 0.08 between 0.8 and 4.0 μm. Our Wd1+RC law is 4% larger at A_F125W/A_Ks and 16% larger at A_F814W/A_Ks. The D16 law is derived from ground-based photometry of 105 evolved stars and a sample of RC stars in the Wd1 field, using color-excess ratios and requiring a power-law functional form between the JHK_s filters. We believe that our study offers several key advantages: (1) space-based photometry with higher precision (∼0.01 mag) than is generally possible for ground-based observations; (2) a sample of kinematically selected Wd1 members that is about four times larger than the D16 sample and is composed of MS stars, where stellar models are better understood relative to evolved stars; and (3) our forward-modeling technique makes no assumption regarding the functional form of the extinction law. However, D16 does have the advantage of using RC stars in the Wd1 field, while we rely on RC stars toward the Arches cluster. Unfortunately, the HST data does not stretch as far to the red as the D16 observations, and so we do not observe a significant RC population in the Wd1 field. That said, only our law is able to reproduce the observed MS colors of the cluster, indicating its effectiveness.

Schödel et al. (2010; hereafter S10) measure the total extinction of the GC in the H and K_s filters using observations of RC stars in the region. Adopting a K_s absolute magnitude M_Ks = −1.54 mag (Groenewegen 2008), an intrinsic H–K_s color (H–K_s)₀ = 0.07 mag, and a GC distance of 8.03 kpc, they calculate absolute extinction values of A_H = 4.48 ± 0.13 mag and A_Ks = 2.54 ± 0.12 mag. This results an extinction ratio A_H/A_Ks = 1.76 ± 0.18, providing an independent test of the normalization of the extinction law.

We calculate A_H/A_Ks for the Wd1+RC extinction law as well as the N09 law, which is often adopted for the GC/Galactic Bulge. Calculating A_H and A_Ks at effective wavelengths of 1.677 μm and 2.168 μm, respectively, the Wd1+RC law produces A_H/A_Ks = 1.936 ± 0.08 (statistical and systematic combined in quadrature). This is 10% higher than the S10 result, but within 1σ given the uncertainties. The N09 law predicts a value of 1.68 ± 0.03, which is 5% lower than the S10 value (0.4σ difference). However, it is important to note that the A_H/A_Ks value in S10 and the extinction laws from this work and N09 rely on knowing the absolute magnitude, colors, and distance of RC stars in the filters of interest. This can lead to systematic errors in the extinction law, as discussed in Section 5.4. We conclude that our law is in acceptable agreement with the S10 measurement.

The Wd1+RC law is also broadly consistent with the most recent measurement of the GC NIR extinction law, which reports a power law with β = 2.31 ± 0.03 (Nogueras-Lara et al. 2017). A power-law fit to the NIR filters in our law results in β = 2.38 ± 0.15. However, we reiterate our result in Section 4.3 that the Wd1+RC law is inconsistent with a power law in the NIR regime. Non-power law behavior is indicated in Nogueras-Lara et al. (2017) as small differences between power-law exponents derived in the JH versus HK_s filters, although they conclude that a single power law is sufficient within the sensitivity of their study.

By forward-modeling the Wd1 photometry, we can account for nonlinearity in the reddening vector. This can occur in regions of high extinction, where the extinction-weighted central wavelength of a filter (i.e., the flux-weighted average wavelength of the extinguished stellar energy distribution convolved with the filter function) changes relative to another filter. Differences in the filter width (e.g., Kim et al. 2005, 2006) or slope of the extinction law between filters can cause this effect. Both processes are captured by the synthetic photometry in our extinction law analysis.

We observe reddening vector curvature in the HST 2CD, where a simple linear fit to the MS sample does not trace back to the origin at A_Ks = 0, as expected for these stars (Figure 18). Two extinction law fits to the data are shown: one with the PanStarrs y point (0.962 μm) included, and another without. Both vectors have significant curvature and are thus able to broadly match the data and trace back to the origin in the zero-extinction case. However, the model with the y point has more curvature and is able to fit the data better. This is because adding A_y/A_Ks increases the flexibility of the extinction curve between the F125W and F814W filters, allowing for a steeper slope through the F814W filter, as is called for by the data.

Zoom In Zoom Out Reset image size

Here we explore potential sources of systematic error and quantify their impact on the extinction law analysis.

5.4.1. Red Clump Star Model

5.4.2. Galactic Center Distance

5.4.3. Total Systematic Error

The Wd1+RC law can be applied to stellar populations with similar foreground dust as the Wd1 and Arches field RC stars, namely the spiral arms of the Galaxy in the Galactic Plane. Examples of future applications include studies of stellar populations near the GC, such as the Quintuplet cluster and Young Nuclear Cluster, and the structure and kinematics of the inner Bulge at low galactic latitudes. To aid future use, the Wd1+RC extinction law at the pivot wavelength of several commonly used filters is provided in Table 5. In addition, a python code to generate the law at any wavelength between 0.8 and 2.2 μm is available online (see Section 4.3). The law can also be accessed through the astropy-affiliated package PopStar, which will soon be released in a beta test.

While a detailed analysis of the dust properties is beyond the scope of this study, the measurement of significant non-power law behavior in the OIR extinction law suggests that there are subtle features that can be used to constrain dust models. The steepness of the law presents a challenge as well, as the classic silicate+graphite grain models (e.g., Mathis et al. 1977; Weingartner & Draine 2001) struggle to produce NIR laws steeper than the canonical C89 law (Moore et al. 2005; Fritz et al. 2011). Models that incorporate composite grains with organic refractory material and water ice (e.g., Zubko et al. 2004) and porous dust structures (e.g., Voshchinnikov et al. 2017) are promising in this regard.

A major caveat of this analysis is that it relies on the assumption that the extinction law is the same for Wd1 and the RC stars. While this is supported by our analysis and the literature, future studies are needed to confirm this assumption, especially at shorter wavelengths (e.g., A_F814W/A_Ks), where R_V-like variations might begin to have an effect. Currently, the law normalization is set almost entirely by the RC stars, due to the relatively large uncertainty in the Wd1 cluster distance. A reanalysis of the Wd1 data when the distance is known to higher precision would allow for a constraint on the normalization from the Wd1 stars directly. Similarly, the shape of the extinction law is dominated by the Wd1 data, since the RC observations are limited to just the F127M and F153M filters. Additional observations of the RC stars in filters across a larger wavelength range would confirm the shape of the law toward the RC stars.

In Section 5.4 we show that the systematic errors, mainly coming from uncertainties in the intrinsic RC star properties and GC distance, are roughly the size of the presented error bars in the Wd1+RC extinction law. As a result, a more precise measurement of the extinction law (at least those that use RC stars) is not possible until these uncertainties are addressed. In particular, observational calibrations of RC star models at older ages (∼10 Gyr) are needed to correctly represent the Bulge population, using multiple filters to test the stellar evolution and atmosphere models. Continued progress in understanding the Bulge age/metallicity distribution as well as the GC distance will also decrease the systematics in the extinction law analysis.