Hubble Space Telescope Observations of Mira Variables in the SN Ia Host NGC 1559: An Alternative Candle to Measure the Hubble Constant

This study uses 10 epochs of HST WFC3/IR F160W data (GO-15145; PI Riess), with an exposure time of 1006 s for each epoch. The observations are also designed to measure Cepheid variables found using epochs of optical imaging (Miras are too red to appear in the optical images). The epochs have an approximately monthly cadence between 2017 September 8 and 2018 September 11, with a field of view centered at α = 04^h17^m37^s and δ = −62°47'00'' (J2000). Table 1 contains the observational information for each epoch.

We use pipeline-calibrated images retrieved from the HST MAST Archive. We generate pixel-resampled and stacked images for each epoch using Drizzlepac 2.2.2 (Gonzaga et al. 2012). Each epoch contains two subpixel dither positions and has been resampled to a scale of 012 pixel⁻¹ (WFC3 has a native scale of 0128 pixel⁻¹), and the orientation of the images rotates by slightly over a full rotation as a consequence of HST orient constraints. We choose the first epoch as the reference image and align all of the subsequent images with it using DrizzlePac.

We begin by stacking all of the F160W observations to make a deeper "master" image. This image is used to create a master source list for the photometry and also to derive crowding used in period fitting. Given the crowded nature of our fields, we use the crowded-field photometry packages DAOPHOT/ALLSTAR (Stetson 1987) and ALLFRAME (Stetson 1994) for the data reduction and source extraction. Our procedure closely follows the steps used to search for Miras in the water megamaser host NGC 4258 (H18).

The DAOPHOT routine FIND is used to detect sources with a >3σ significance. The significance depends on sky background, readout noise, and the number of images combined. Then, the DAOPHOT routine PHOT is used to perform aperture photometry. We input the star list produced from the aperture photometry into ALLSTAR for point-spread function (PSF) photometry, using a 2.5-pixel FWHM intensity, standard input into ALLSTAR. We then repeat these steps on the star-subtracted image generated by ALLSTAR (with all of the previously discovered sources removed) to produce a second source list. The two source lists are then concatenated and used as input into another round of ALLSTAR (applied to the image without subtracted sources) to create a final master source list of ∼49,000 objects.

This master source list is input into ALLFRAME. We use the master source list to derive a consistent source list for PSF photometry of each F160W epoch. ALLFRAME is similar to ALLSTAR, except that it is capable of simultaneously fitting the profiles of all sources across the full baseline of epochs for a single field. This allows it to maintain a constant source list through multiple epochs and improves the photometry for stars closer to the detection limit. ALLFRAME produces time-series PSF photometry as output.

Secondary standards are used to correct for variations in the photometry of constant sources (detector and image quality) between different epochs. We search for secondary standards in the star lists by choosing bright objects that were observed in all 10 F160W epochs. Their stellar profiles and surroundings were visually inspected to choose secondary standards that are relatively isolated compared to other stars in the image, removing any that showed variability or had large photometric errors. This leaves us with a total of 31 sources, summarized in Table 2. We calculate the celestial coordinates for all of these sources using the astrometric solutions in the FITS headers. Photometry is conducted on the resampled images. The mean residuals for these stars across all epochs of F160W imaging exhibited a dispersion of 0.01 mag. We use the standard photometry routine TRIALP (kindly provided by Peter Stetson) to derive the epoch-to-epoch changes in zero-point using our standard stars, and this is used to correct the ALLFRAME photometry for epoch-dependent variations.

Table 2.  Secondary Standards

ID	α	δ	X	Y	F160W	Uncertainty
	(J2000)	(J2000)	(pixels)	(pixels)	(mag)	(mag)
7599	04h17m40019	−62°46'1751	843.0	203.9	18.939	0.081
10534	04 17 34.146	−62 46 41.15	461.3	280.8	19.039	0.078
4002	04 17 31.105	−62 46 26.79	335.7	111.1	19.195	0.070
25011	04 17 37.447	−62 47 23.13	525.8	673.1	19.238	0.093
10962	04 17 39.975	−62 46 28.86	809.8	292.5	19.355	0.104
8924	04 17 37.077	−62 46 28.90	653.0	238.8	19.381	0.060
25548	04 17 32.556	−62 47 36.57	225.1	687.9	19.433	0.066
16153	04 17 39.255	−62 46 48.42	717.8	433.2	19.466	0.107
4308	04 17 29.008	−62 46 32.80	206.1	119.5	19.684	0.070
31503	04 17 41.974	−62 47 34.16	740.5	844.3	19.789	0.083
25269	04 17 37.822	−62 47 23.10	546.2	679.8	19.823	0.050
11364	04 17 38.388	−62 46 34.00	710.0	303.5	19.895	0.052
5369	04 17 40.051	−62 46 10.18	864.6	146.8	19.986	0.074
24951	04 17 44.680	−62 47 05.78	963.8	671.1	20.022	0.064
7243	04 17 38.170	−62 46 20.78	734.1	195.2	20.168	0.075
29515	04 17 41.245	−62 47 29.25	714.4	791.0	20.193	0.111
26819	04 17 32.001	−62 47 42.17	179.9	721.8	20.195	0.064
25967	04 17 42.801	−62 47 13.83	840.4	699.5	20.236	0.166
31273	04 17 35.870	−62 47 47.77	373.8	837.8	20.236	0.077
1510	04 17 34.915	−62 46 09.25	589.4	43.7	20.242	0.057
35832	04 17 43.821	−62 47 43.95	813.7	955.8	20.265	0.094
23326	04 17 42.191	−62 47 06.19	828.2	627.9	20.282	0.055
29323	04 17 32.986	−62 47 48.13	216.9	787.1	20.325	0.059
30041	04 17 34.260	−62 47 47.51	287.5	805.9	20.659	0.109
27599	04 17 37.056	−62 47 32.79	478.5	741.9	20.690	0.082
27351	04 17 34.181	−62 47 38.82	306.8	735.9	20.702	0.067
28149	04 17 33.826	−62 47 42.35	278.0	757.1	20.778	0.078
35037	04 17 39.297	−62 47 51.99	547.5	934.9	20.803	0.063
9301	04 17 29.739	−62 46 47.50	205.8	248.9	20.835	0.075
5013	04 17 40.445	−62 46 08.02	891.7	137.1	20.932	0.067
6365	04 17 38.468	−62 46 17.28	759.7	173.2	21.021	0.059

Note. All of the secondary sources used to calibrate the F160W light curves. ID numbers are photometry IDs. The X and Y positions are relative to the first epoch of the F160W image, and the uncertainties are photometric errors as estimated by DAOPHOT. The magnitudes are calculated as explained in Equation (2).

Download table as:  ASCIITypeset image

To put our observations on the same zero-point as H18, which we adopt for an absolute calibration of distance, we use the same HST WFC3 IR photometric zero-point for a 04-radius aperture. In Vega magnitudes, the F160W zero-point is 24.5037 mag. The use of a consistent zero-point in a fixed aperture where we normalize the flux of the PSF compensates for small differences on the observing pattern, PSF, and resampling scale between NGC 4258 and NGC 1559. In NGC 4258, the observations of each epoch consisted of four subpixel dithers. This allowed us to drizzle images for each epoch to a finer scale in NGC 4258 (008 pixel⁻¹) than in NGC 1559 (012 pixel⁻¹), where we only have two subpixel dithers per epoch.

To define aperture corrections to the PSF, we select the brightest and most isolated stars in the field that are not saturated. These are often foreground stars, unlike secondary standards, which are chosen to be closer in expected magnitude to our variables. These stars are used to determine the correction from the PSF magnitude to the 04 aperture magnitude only, not as photometric standards.

We subtract all other stars in the master source image from the image using DAOPHOT/SUBSTAR to create a "clean" image. Aperture photometry is then performed on the bright stars in the clean image using 10 apertures ranging from 1 to 5 pixels to create a growth curve. We inspect the growth curve of each of the stars used in the aperture corrections and reject any stars that do not have smooth, monotonically increasing growth curves. The aperture correction is defined as

$$\begin{eqnarray}&&{\rm{\Delta }}{m}_{\mathrm{ac}}={m}_{\mathrm{ap}}-{m}_{\mathrm{PSF}},\end{eqnarray} \tag{ 1 }$$

where m_ap is the aperture magnitude for an aperture with a radius of 3.3 pixels, equivalent to the 04 aperture zero-point, and m_PSF is the PSF magnitude. We find an aperture correction of $0.06\pm 0.01$ mag, which we add to m_PSF to correct to the HST zero-point for a 04 aperture. PSF photometry is used throughout our analysis, and the aperture correction is included at the end to put our old and new observations on the same HST zero-point. The aperture correction may change slightly between the observations of NGC 4258 and NGC 1559 owing to "breathing" and small focus changes. Because we do not have standard stars in common between our two fields, we cannot correct for this effect, which is estimated to be ∼0.01 mag, so we include this in the uncertainty. Thus, including the aperture correction, we convert from instrumental magnitudes to calibrated magnitudes using

$$\begin{eqnarray}&&{m}_{f}={m}_{i}-2.5\mathrm{log}(t)-{a}_{\mathrm{DAO}}+{a}_{{\text{}}{HST}}+{\rm{\Delta }}{m}_{\mathrm{ac}}+{\rm{\Delta }}{m}_{b},\end{eqnarray} \tag{ 2 }$$

where m_f is the final calibrated magnitude, m_i is the instrumental magnitude output from DAOPHOT photometry, t is the exposure time, a_DAO is the DAOPHOT photometric zero-point (25.0 mag), a_HST is the HST zero-point for a 04 aperture (24.504 mag), Δm_ac is the aperture correction (0.06 mag), and Δm_b is the crowding bias correction (mean value of 0.13 mag) discussed in Section 3.4.

The first step in our Mira selection process is calculating the Welch–Stetson variability index L of all of the objects in our field (Stetson 1996). This index is only used to weed out as many nonvariable stars as possible before beginning the later parts of the search, as many objects with high variability indexes may still be non-Miras, short-period variables, or simply nonvariable stars with noisy measurements. We use a cutoff for the minimum variability index of L ≥ 0.75, consistent with the criterion H18 used for NGC 4258. With simulations, we find that this variability cut removes ∼90% of the nonvariable stars. Figure 1 shows a distribution of the L index for all of the stars we detect in our master image. Approximately 3000 variable-star candidates remain at this stage.

Zoom In Zoom Out Reset image size

Next, we calculate the most likely period (assuming that the variability is periodic) of the remaining candidates. For each, we perform a grid search on periods ranging from 100 to 1000 days, the relevant period range for Miras. At each period we fit the light curve in magnitude space with a sine function with the amplitude, mean magnitude, and phase as free parameters using a χ² minimization, as done by H18. In addition, the mean magnitude from the sine fit is used as the mean magnitude of the Mira for fitting the PLR in both this paper and H18, giving us a consistent definition of mean magnitude. The period with the lowest χ² fit to a sine is then chosen as the best-estimate period of the Mira. As our observation baseline spanned ∼370 days, and the shortest-period Miras we can effectively use in our PLR relations have P = 240 days (see Section 3.3), we would not expect to see aliasing of periods within our target range in our period estimations. Our observation timing is such that each Mira in our period range is sampled at least five phases per cycle.

Aliasing of shorter-period variables to the period range of Miras that we seek is possible. However, this should be rare since we employ an amplitude cut, which removes objects that do not have large amplitudes in F160W from our sample. The majority of variable stars with shorter periods (such as Cepheids) that could be aliased with our monthly observing cadence also have amplitudes <0.4 mag in the NIR and would not pass these amplitude cuts (see Figure 2 of H18 for a comparison of H-band amplitudes of variable stars in the LMC). Laney & Stobie (1993) found that some Cepheids with periods between 20 and 40 days may reach an amplitude of 0.4 mag in the H band. However, these Cepheids would have massive progenitors and are this expected to be to be primarily found in the star-forming regions of the galaxy, which are excluded by our surface brightness cuts (explained in Section 3.4).

Zoom In Zoom Out Reset image size

Aliasing of longer-period variables (P > 400 days) into our sample period range is also possible—particularly if we discover a Mira with a long secondary period and a bump. This bump can resemble the end of one cycle if we do not have a long-enough baseline to observe the full light curve for one of the Mira's pulsational cycle. However, these stars are typically C-rich, heavily reddened, and more massive ($\gtrapprox 3\,{M}_{\odot }$), thus making them relatively rare. A cursory examination of the OGLE data yields ∼60 stars out of a sample of >1500 (∼4% of the sample) that have these qualities for the entire LMC sample. Most of them clustered in the bar of the LMC, further suggesting that they are likely to have massive progenitors. In addition, light curves—especially those of heavily reddened stars—are more regular in the NIR than they are in optical wavelengths, so we also expect this effect to be smaller than in the optical OGLE data.

Finally, we refer to H18, where we were able to combine optical observations of Mira variables with a rapid cadence of about one visit every few days (used to search for Cepheids) and a baseline of ∼40 days with the roughly monthly visits and year-long baseline used to search for Miras. In our gold sample of Miras, which contained objects with both F814W light curves and F160W light curves, the zero-point was within 0.02 mag of that of the bronze sample Miras that had only F160W light curves to classify them, as is the case in NGC 1559. In addition, all of the optical light curves we were able to recover for our sources indicated long-period objects consistent with the periods we determined in our analysis. Additionally, the phase near the trough of the light curve at the time of the optical observations appeared to be the primary cause of Miras not having an optical counterpart. Thus, the group of gold Miras is a random sample of the bronze Miras, and not physically distinct.

The similarity of the zero-points between the three samples in NGC 4258 suggests that having additional observations with a longer baseline or more frequent visits, while useful, is not essential to accurately determining the periodicity of objects in our sample. In addition, cross-matching previously known Cepheids to sources in our NGC 4258 master frame revealed that none of the known Cepheids in this field had been falsely detected as Miras. We were also unable to recover any Cepheids as variables based solely on their NIR variability without using prior knowledge on their locations. Some contamination by objects of the wrong period is unavoidable, but we estimate the number of contaminants to be most likely less than ∼5% of our sample, and even smaller after outlier rejection.

Following H18, we consider only variables in the amplitude range 0.4 mag < A < 0.8 mag to be candidate O-rich Miras. The minimum amplitude removes semiregular variables from our sample, which typically have smaller amplitudes (Soszyński et al. 2013), or heavily blended Mira stars, which will be artificially bright and have reduced amplitudes by their blending with a constant contaminant. The maximum amplitude is used to limit the number of heavily reddened C-rich Miras in our sample, which have on average larger amplitudes than O-rich Miras with the same period range. While C-rich Miras with low reddening track the O-rich Mira PLR, heavily reddened C-rich Miras fall below the PLR and would bias our fit (Whitelock et al. 2003). From our simulations, we find that our amplitudes are well recovered, with the mean amplitude of the input and recovered Miras within a few 0.01 mag of each other.

Yuan et al. (2017b) found larger amplitudes for C-rich Miras than O-rich Miras in the LMC and Cioni et al. (2003) in SMC variables. Goldman et al. (2019) found that in IR wavelengths, both period and amplitude of pulsating AGB stars correlated with color, with redder variables having larger amplitudes. However, this effect does not appear in Soszyński et al. (2009), where the LMC C-rich stars, which are typically redder than their O-rich counterparts, do not have larger amplitudes. It is possible that part of the discrepancy may lie in inconsistencies in the definition of amplitude (since the OGLE data are longer baseline and we use only single-cycle amplitudes), or perhaps the range of pulsational amplitudes for each Mira spectral type differs between galaxies. The latter would imply that the amplitude cuts may need to be different in each galaxy. However, even in the metallicity range of the DUSTiNGS galaxies in Goldman et al. (2019) there does not appear to be an obvious trend in galaxy metallicity and amplitude. As a result, we still use an amplitude cut and include possible effects of the amplitude cut on the zero-point in our metallicity-related error budget (discussed in greater detail in Section 3.5). A more complete understanding of metallicity will be necessary to place firmer constraints on the systematic effect.

Besides using the period and amplitude of the sine fit, to remove imposters, we may also compare the quality of the fit to that of a monotonic rise or fall. We compare each object's goodness of fit to a sine and a line using the F-statistic, which is defined as

$$\begin{eqnarray}&&F={\chi }_{s}^{2}/{\chi }_{l}^{2},\end{eqnarray} \tag{ 3 }$$

where ${\chi }_{s}^{2}$ is the reduced χ² of a sine fit and ${\chi }_{l}^{2}$ is the reduced χ² for a line fit. This statistic may be used to retain only objects with light curves for which a sine fit is preferred over a straight line. We use simulations to determine the best value for the F-statistic to retain Miras while removing noisy, static stars. We find that keeping sources with F < 0.5 is a fairly optimal choice; it removes ∼92% of the nonvariable stars while keeping ∼80% of all of our simulated variables.

Examples of the straight-line and sine fits are shown in Figure 2. We further visually inspect all of the objects that pass these criteria to remove spurious fits.

We limit our sample to P < 400 days to exclude Miras undergoing HBB (although we would expect few if any HBB stars in this part of the galaxy because HBB stars typically have masses >3 M_⊙, which will be rare). An additional rationale for our period limit is that our observations span ∼370 days, making the accurate recovery of periods greater than 400 days problematic. Figure 3 illustrates the simulated period recovery rate as a function of period, and Figure 4 shows the recovered and input periods. We considered a period recovered if the measured and true values were within 15% of each other. As a statistical error, this criterion would limit a period uncertainty to less than the apparent scatter of the PLR.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also limit our sample at the short-period, faint end of the PLR, owing to magnitude incompleteness. The Mira PLR has an intrinsic width, and not all stars at the same period have the same magnitude. This means that the faintest stars at shorter periods are undetected, biasing us toward the brighter end of the PLR. Additional scatter due to the crowding and photometric errors further increases this magnitude range. We search for the minimum period above which we avoid this bias using three approaches.

First, to obtain an idea of the expected lower bound on our search for the completeness limit, we use the HST Exposure Time Calculator (ETC) to calculate the estimated signal-to-noise ratio (S/N) for each period range. We derive a preliminary fit to the candidate Mira PLR to relate Mira periods to their mean magnitudes. We then determine whether that mean magnitude would be detected at S/N > 5 at the depth of each imaging epoch, a minimum requirement for useful detection. This analysis leads us to draw the initial lower bound at ∼200 days. While this is the theoretical low-period limit to which we should discover Miras, in practice, crowding and the intrinsic width of the Mira PLR introduce additional scatter, causing some Miras near that limit to go undetected. Thus, we only use the ETC results as a starting point.

Next, we proceed to look for the expected empirical signature of incompleteness—an approximate plateau of the PLR intercept or zero-point at long periods and a brightening trend of the zero-point as we cross the period completeness limit toward the short-period direction. In order to find the location of this break, we start at the initial value of 200 days from the HST ETC analysis and use a moving 75-day window (boxcar) to reduce the noise and measure the period range average of the zero-point. The result is shown in Figure 5. At P > 240 days the trend reaches a plateau (with oscillations due to noise and irregular sampling above this).

Zoom In Zoom Out Reset image size

Finally, we use simulations (described in Section 3.6) to verify this empirical result based on the expected recovery rate.

Figure 3 displays the number of sources recovered as a function of starting period. We can see that an increasing fraction of stars is recovered as the input period increases, and at ∼240 days the recovery fraction reaches ∼90%, which is in accord with the empirical analysis. Therefore, we choose this period as our short-period end for the completeness.

Cepheid analyses (Hoffmann et al. 2016; Riess et al. 2016) in the NIR routinely use artificial-star injection and recovery to correct for the artificial brightness in detected sources due to unresolved background sources—that is, crowding. We follow this approach. Like the Cepheid analyses, we choose not to include objects in very crowded regions, thereby reducing the size of this correction (and its associated uncertainties). In the case of NGC 1559, this is most effectively addressed by including a requirement for low surface brightness background.

NGC 1559 has high surface brightness near the center, causing the image to be particularly crowded toward the nucleus of the galaxy (see Figure 6). The inclination of the galaxy, i ≈ 573 (Kassin et al. 2006), is a likely contributing factor to its high surface brightness near the center. Through our artificial-star tests (described in the following section), we found that the crowding bias corrections would generally exceed 1 mag throughout an elliptical region around the center, meaning that these sources are dominated by the background and will lose reliability compared to Miras in less crowded regions. We measured the local surface brightness in this region as the average number of counts per second in a 2''-wide box around each candidate. We identified a minimum surface brightness of ∼420 counts s^–1 arcsec^–2 as corresponding to this contiguous region and set this as the boundary for the maximum allowed surface brightness to be included in our sample. This leaves objects that are mainly distributed in the outer regions of the galaxy, as seen in Figure 6. As expected of aged populations, the Miras in NGC 1559 are evenly distributed in the regions of lower surface brightness and do not follow spiral arms or other galactic structures.

Zoom In Zoom Out Reset image size

After using surface brightness cuts to remove objects that would be in the most crowded regions, we employ artificial-star tests to determine the individual crowding bias for the remaining candidate Miras. For each remaining Mira candidate we fit an initial PLR. We then inject 100 artificial stars for every Mira, using the Mira's period and the initial PLR fit to determine a starting point for the artificial-star magnitude. The artificial star is randomly dropped within a 2'' radius of the Mira, so that both are in the same environment.

We perform aperture and PSF photometry on the artificial stars exactly as we would with the real Miras and measure the mean magnitudes at which they are recovered. The difference in magnitude between the input artificial stars and the output is used to get a first estimate of the correction for each star. The magnitude of the original Miras is then corrected before refitting the PLR to determine a new guess for the true magnitudes, and the process is repeated until the Mira magnitudes converge. The corrections apply to our use of the mean magnitude of the Mira (defined by the mean of the sine fit).

Occasionally, our artificial stars will fall close enough to a bright background star that the two will be "blended" (i.e., unresolved). If the superimposed star is sufficiently bright, this blending would suppress the amplitude of a variable star enough that it is unlikely to be recovered as a Mira based on our minimum-amplitude requirement. We approximate the blended amplitude by the expression

$$\begin{eqnarray}&&F={10}^{-0.4{\rm{\Delta }}m}+{10}^{-0.4A\phi },\end{eqnarray} \tag{ 4 }$$

where Δm is the difference in magnitude between the two stars (positive if the background star is fainter), A is the amplitude of a Mira, and ϕ is the phase of the Mira. We can likewise calculate the magnitudes at the peak and trough of the light curve and subtract the two to obtain a new blended amplitude. At an original amplitude of 0.6 mag (mean of our range), and a background star 1 mag fainter than the artificial star, we find that the blended amplitude is ∼0.4 mag, our cutoff for the low-amplitude end. Blended stars brighter than this threshold would cause the mean Mira to be removed solely on the basis of amplitude cuts. Thus, we consider the artificial star to be blended if it is within 1 pixel of a bright star (defined as <1 mag fainter than the Mira), and these (like real blended Miras) are not included in our artificial-star samples.

After a 3σ clip of the artificial-star sample (the same clip we eventually use when fitting the PLR), we take the standard deviation of these crowding corrections to get an estimate for the uncertainty in each Mira's mean magnitude. We use the mean of the clipped crowding corrections to correct the magnitude of each Mira. Figure 7 shows the distribution of crowding corrections for all of our artificial stars, after subtracting the mean correction for each Mira and dividing by the standard deviation of the corrections for each Mira. This distribution in magnitudes compares well to a Gaussian distribution (in red).

Zoom In Zoom Out Reset image size

We also conservatively eliminate five Mira candidates (3% of the sample) that have a crowding correction Δm_b > 0.5 mag. Although these may yield good measurements, this level of crowding is significantly greater than for the rest of the sample (the mean crowding correction is ∼0.13 mag, with a standard deviation of 0.14 mag). In addition, applying a crowding correction reduces the scatter of the PLR by approximately 20%. The local surface brightness criterion already excludes the majority of candidates that would have fallen into this category, so this final cut ensures that all of the objects remaining are in uncrowded regions.

Although we used amplitude cuts to reduce the number of C-rich Miras or nonvariable stars in our sample, here we quantify and correct for a potentially low level of residual contamination from C-rich Miras. C-rich Miras are typically fainter than O-rich Miras in F160W; thus, C-rich contamination in the Mira PLR can result in a fainter absolute magnitude calibration. We employ the OGLE LMC sample—the most complete sample of Mira variables with well-categorized Miras—to measure the difference in zero-point between a pure O-rich sample and a model C-rich-contaminated sample (by mixing O-rich and C-rich LMC Miras). We note that this is not a large correction because for the LMC the contaminated sample would shift the zero-point by ∼0.07 mag.

To create the LMC "contaminated sample," we convert the Two Micron All Sky Survey (2MASS) ground-based J and H data from the LMC to F160W using Equation (8) and then apply the same amplitude cuts and sigma-clipping used for NGC 4258 and NGC 1559. The single-epoch amplitudes in the H band are estimated from the H-band light-curve templates from Yuan et al. (2017b). In the LMC, at P < 400 days, the ratio of O-rich to C-rich Miras is approximately 2:3 before the amplitude cuts and 2:1 after. As expected, the amplitude cuts remove the majority of C-rich Miras in this range. Then, we calculate the zero-point (i.e., the PLR intercept) as a function of period by fitting a PLR with a fixed slope to the Mira sample binned by 50 days. As shown in Figure 8, the zero-point of the LMC contaminated sample grows fainter as the contamination by C-rich Miras increases (P ≳ 250 days), whereas a purely O-rich sample (P < 250 days) is relatively flat. Therefore, we can use the change in zero-point as a function of period of the LMC contaminated sample as a template to measure the contamination in the other hosts. After plotting the zero-point as a function of period for each galaxy, we use the amplitude and shape of the curve for each galaxy relative to the amplitude and shape LMC curve in order to estimate the fraction of C-rich contamination present in that galaxy's Mira sample.

Zoom In Zoom Out Reset image size

We take the samples of Miras in M33 (also well categorized into C-rich and O-rich subclasses) from Yuan et al. (2018) as a test of our ability to measure the contamination. The observations of M33, like NGC 1559 and NGC 4258, impose a magnitude limit on the Mira sample such that the C-rich population is relatively incomplete compared to the O-rich population. We combine the O-rich and C-rich samples from Yuan et al. (2018), which pass our amplitude cuts and measure the zero-point as a function of period. Adopting the LMC curve from Figure 8 as a model for C-rich contamination, we fit the M33 curve from Figure 8 using the form

$$\begin{eqnarray}&&{Z}_{1}=\alpha {Z}_{\mathrm{LMC}}+\beta ,\end{eqnarray} \tag{ 5 }$$

where Z₁ is the zero-point of the contaminated M33 PLR as a function of period, Z_LMC is the zero-point of the contaminated LMC PLR as a function of period, and α and β are fit parameters: β simply shifts the curve and does not contribute to a correction, while α tells us the amount to correct the zero-point as a function of the total difference in zero-point between the combined and O-rich-only values.

For M33 we find that α = 0.58 ± 0.18 is the best fit to the LMC curve. For the LMC, the difference in zero-point is −0.07 ± 0.013 mag. Therefore, in order to correct the zero-point of the combined C-rich and O-rich M33 data set to match the zero-point of the O-rich-only M33 set, m_c = −0.04 ± 0.01 mag. The known difference in zero-point of the purely O-rich and the contaminated PLR is −0.034 mag, within the uncertainties of our contamination correction. We now apply this same technique to NGC 1559 and NGC 4258.

In Figure 8, we plot the change in zero-point as a function of period for NGC 4258 and NGC 1559. Applying the LMC as a model to NGC 4258 and NGC 1559, we find that α = −0.77 ± 0.25 for NGC 4258 but α = 0.81 ± 0.35 for NGC 1559, indicating some possible C-rich contamination. We determine an overall correction to the zero-point of NGC 1559 of m_c = −0.057 ± 0.024 mag. We apply this correction to the m_f given in Equation (2),

$$\begin{eqnarray}&&{m}_{{cf}}={m}_{f}+{m}_{c},\end{eqnarray} \tag{ 6 }$$

where m_c is the C-rich correction. For NGC 4258, where $\alpha \lt 0$, we do not correct for any possible C-rich contamination, since C-rich contamination should only make the zero-point fainter, not brighter.

We would expect to find more C-rich contamination in NGC 1559 than in NGC 4258 because the range of Mira periods extends up to 400 days in NGC 1559 and only up to 300 days in NGC 4258, and there are more C-rich stars at longer periods. Depending on the amount of circumstellar extinction, C-rich Miras can be as much as 2–3 mag fainter in the H band than O-rich Miras having the same period (Ita & Matsunaga 2011). Thus, we would also expect to find less contamination in magnitude-limited samples like NGC 1559 and NGC 4258 than in the complete sample in the LMC, in good agreement with the results of our C-rich contamination corrections.

To test the efficacy of our selection criteria, we use simulated light curves of Miras and nonvariable stars. We also use the simulations to determine the amount of contamination that we might expect from nonvariable stars and the completeness limit of the sample. For each simulation, we inject 105 artificial variable stars into the star catalogs at 21 periods between 200 and 400 days. The artificial stars have random phases, similar to our Mira sample. We use a sine wave of the appropriate period and a starting amplitude of 0.7 mag, the estimated median "true" amplitude of our sources. We then include a scatter of 0.4 mag to the mean magnitude of the source to mimic the spread of our Mira PLR. Next, using the HST ETC, we estimate the S/N of each source, add in the appropriate noise, randomize the starting phase, and sample the light curve at the observation dates of our survey. We perform the rest of the variability search on these stars, just as we would for the true Miras. The same is done for nonvariable stars to see the effects of the variability cuts on both the simulated Miras and nonvariable stars. Nonvariable stars are simulated by using starting magnitudes in the same range as the Mira variables, but with zero amplitude.

To test the level of nonvariable star contamination, we simulate 10,000 nonvariable stars at the catalog level and then apply all of our Mira selection criteria (outlined in Sections 3.1–3.3). We find four stars that passed all of the variability cuts. However, all four stars also fall more than 1.5 mag below the Mira PLR and were sigma-clipped out of the final sample. Thus, we conservatively estimate that we may have 1 nonvariable star contaminant per 10,000 nonvariable stars in the master star list. Out of the observed star list of 49,000, approximately 40,000 are nonvariable stars. We expect an upper limit of four nonvariable star contaminants in the Mira PLR, and we do not expect these to meaningfully bias the PLR. A comparison of the selection criteria from this paper and H18 is shown in Table 4.

Table 4.  Mira Sample Criteria

	NGC 1559	NGC 4258 (Gold)
Period cut (days)	240 < P < 400	P < 300
Amplitude cut (mag)	0.4 < ΔF160W < 0.8	0.4 < ΔF160W < 0.8
Surface brightness cut (counts s–1)	421	⋯
F-statistic	${\chi }_{{\rm{s}}}^{2}/{\chi }_{{\rm{l}}}^{2}\lt 0.5$	⋯
Color cut (mag)	⋯	${m}_{{\rm{F}}125{\rm{W}}}-{m}_{{\rm{F}}160{\rm{W}}}\lt 1.3$
F814W detection	⋯	Slope fit to F814W data >3σ
F814W amplitude (mag)	⋯	ΔF814W > 0.3

Note. A comparison of the criteria for the final Mira samples in NGC 1559 and NGC 4258. Owing to differences in S/N and the available data, we were unable to match the criteria exactly.

Download table as:  ASCIITypeset image

Miras and other variable stars are typically fit with a linear PLR, and we apply two different choices of slopes derived from LMC Miras. Because the Mira samples in NGC 4258 and NGC 1559 have lower S/N and a shorter period range than those in the LMC, we do not attempt to use these to derive an independent slope, though this may become possible as the SN–Mira host sample grows.

The first slope is derived by Yuan et al. (2017b) using observations of about 170 LMC O-rich Miras in the H band, which is the closest match to F160W:

$$\begin{eqnarray}&&m={a}_{0}-3.64(\mathrm{log}P-2.3),\end{eqnarray} \tag{ 7 }$$

where P is the period in days, m is the H-band magnitude, and a₀ is the intercept or zero-point. For the "gold" sample of NGC 4258 Miras, we obtain ${a}_{0}=23.15\pm 0.017$ mag (statistical error only) with this slope.

Table 6.  Summary of Mira H₀ Values

Anchor	Period Range (days)	H0 (km s−1 Mpc−1)
		Slope = −3.35	Slope = −3.64
NGC 4258	240 < P < 300	$74.6\pm 5.1$	$74.7\pm 5.1$
NGC 4258	240 < P < 400	$72.7\pm 4.6$	$72.5\pm 4.6$
LMC	240 < P < 400	$73.9\pm 4.3$	$73.6\pm 4.3$
LMC+NGC 4258	240 < P < 400	73.3 ± 4.0	$73.2\pm 4.0$

Note. The H₀ with two anchors and the slope fit using F160W (−3.35) is our preferred value.

Download table as:  ASCIITypeset image

In addition to the slope determined by Yuan et al. (2017b), we use the OGLE LMC variable-star database to estimate the slope of the Mira PLR in F160W by transforming from ground-based IR data to the HST bandpass. Using 2MASS J and H data for OGLE LMC Mira variables, we adopt the color correction derived by H18,

$$\begin{eqnarray}&&{m}_{F160W}=H+0.39(J-H),\end{eqnarray} \tag{ 8 }$$

to convert their O-rich Mira sample from ground-based J and H into F160W magnitudes. This sample includes only periods below the HBB cutoff, P = 400 days, to better match our HST sample. Using this sample of 416 Miras with P < 400 days, we obtain a slope fit of

$$\begin{eqnarray}&&m={a}_{0}-3.35(\mathrm{log}P-2.3).\end{eqnarray} \tag{ 9 }$$

For this slope, we obtain ${a}_{0}=23.18\pm 0.017$ mag (statistical error only) for the NGC 4258 gold sample of Miras. In Table 5 we include a systematic uncertainty based on the differences in a₀ when using these two slopes.

Using the Schlafly & Finkbeiner (2011) dust maps, we find that the Galactic extinction in the direction of NGC 1559 ~0.022 mag in the 2MASS J band and ~0.014 mag in the 2MASS H band. However, for the purposes of calibration, we are affected by the difference in extinction between NGC 4258 and NGC 1559. The Galactic foreground extinction at the location of NGC 4258 in the 2MASS J and H bands is 0.012 and 0.007 mag, respectively, for a difference of 0.01 and 0.007 mag between the two, for which we correct as follows:

$$\begin{eqnarray}&&{m}_{\mathrm{PLR}}={m}_{{cf}}+{m}_{e},\end{eqnarray} \tag{ 10 }$$

where m_PLR is the final mean magnitude of the Mira (which we use for fitting the PLR), m_cf is the C-rich contamination-corrected mean magnitude defined in Equation (6), and m_e is the correction for the differential foreground extinction, which we estimate to be 0.01 mag. The F160W lies between the J and H bands, so we derive its extinction by interpolating between the J and H bands.

We expect the difference between the Mira interstellar extinction in NGC 1559 and NGC 4258 to be quite low. Miras are not a disk population where dust is concentrated. The positional distribution of the NGC 4258 Mira sample (see H18) is uniform and thus should mimic a lower-extinction halo distribution. Likewise for NGC 1559, most of our sources are in the outer regions of the galaxy and thus would have low extinction. Most importantly, our measurements are at 1.6 μm, which is much less sensitive to extinction than shorter wavelengths. However, there could be a small difference in the extinction between the regions where the NGC 4258 and NGC 1559 Miras are located. Riess et al. (2009) suggested a ∼0.04 mag differential extinction between Cepheids in the H band between NGC 4258 and SN hosts, and we propagate this value as a systematic uncertainty to account for a difference in interstellar extinction.

While the effects of metallicity on the Mira PLR in the NIR are not well studied, the NIR and IR wavelengths are known to be less sensitive to metallicity than the optical. In addition, there is evidence to suggest that the metallicity effect in NIR wavelengths is <0.1 mag. Goldman et al. (2019) find that at IR wavelengths the pulsation properties of AGB stars show no dependence on metallicity over a range of [Fe/H] = −1.85 (1.4% solar) to [Fe/H] = −0.38 (50% solar). Bhardwaj et al. (2019) found (without correcting for metallicity) that the mean relative distance between the LMC and the SMC measured in the K_s band using the O-rich Mira PLR is consistent with the measurements from previous studies using Miras and Cepheids.

Finally, Whitelock et al. (2008) found an upper limit of ∼0.1 mag on difference in zero-point between the LMC and Galactic O-rich Miras that have both Hipparcos parallaxes and very long baseline interferometry parallaxes and Miras in globular clusters. However, this upper limit uses a distance modulus of μ = 18.39 ± 0.05 for the LMC distance. Current best estimates of the distance to the LMC (Pietrzyński et al. 2019) list this distance as μ = 18.477 ± 0.004 (stat) ± 0.026 (sys). With this revised LMC distance, the two zero-points agree to within ∼0.01 mag. In addition, NGC 4258 and NGC 1559 are both large spiral galaxies, and we expect them to have similar metallicities, further reducing the effects of metallicity on the PLR. However, we conservatively budget a ∼0.03 mag systematic uncertainty on the zero-point due to metallicity effects to account for this and potential amplitude differences (discussed in Section 3.2).