A publishing partnership

The following article is Free article

A Near-infrared Period–Luminosity Relation for Miras in NGC 4258, an Anchor for a New Distance Ladder

, , , , , , , , ,

Published 2018 April 16 © 2018. The American Astronomical Society. All rights reserved.
, , Citation Caroline D. Huang et al 2018 ApJ 857 67DOI 10.3847/1538-4357/aab6b3

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0004-637X/857/1/67

Abstract

We present year-long, near-infrared (NIR) Hubble Space Telescope (HST) WFC3 observations of Mira variables in the water megamaser host galaxy NGC 4258. Miras are asymptotic giant branch variables that can be divided into oxygen- (O-) and carbon- (C-) rich subclasses. Oxygen-rich Miras follow a tight (scatter ∼0.14 mag) period–luminosity relation (PLR) in the NIR and can be used to measure extragalactic distances. The water megamaser in NGC 4258 gives a geometric distance to the galaxy accurate to 2.6% that can serve to calibrate the Mira PLR. We develop criteria for detecting and classifying O-rich Miras with optical and NIR data as well as NIR data alone. In total, we discover 438 Mira candidates that we classify with high confidence as O-rich. Our most stringent criteria produce a sample of 139 Mira candidates that we use to measure a PLR. We use the OGLE-III sample of O-rich Miras in the Large Magellanic Cloud to obtain a relative distance modulus, μ4258 − μLMC = 10.95 ± 0.01 (statistical) ±0.06 (systematic) mag, that is statistically consistent with the relative distance determined using Cepheids. These results demonstrate the feasibility of discovering and characterizing Miras using the NIR with the HST and the upcoming James Webb Space Telescope and using those Miras to measure extragalactic distances and determine the Hubble constant.

Export citation and abstractBibTeXRIS

1. Introduction

The value of the Hubble constant (H0), the current expansion rate of the universe, is a source of great interest in astrophysics. The improved precision in H0 measurements (Riess et al. 2016; hereafter, R16) has revealed a 3.4σ discrepancy with the value inferred from observations of the cosmic microwave background under the assumption of a ΛCDM cosmology (Planck Collaboration et al. 2016). New parallax measurements of seven long-period Cepheids in the Milky Way (Riess et al. 2018) combined with the R16 results increases the tension with Planck to 3.7σ. Although the local results have been confirmed (Bonvin et al. 2017; Follin & Knox 2017; Dhawan et al. 2018) and the discrepancy is not dependent on any one datum (Addison et al. 2018), the standard of proof is high for new physics and additional cross-checks are warranted.

The most precise local measurement of H0 relies on Cepheid variables as distance indicators R16 to calibrate the luminosity of type Ia supernovae (hereafter, SNe Ia) in hosts at nearby distances of 10–40 Mpc. Cepheids remain the best understood and most vetted primary distance indicator (Freedman et al. 2001; Bono et al. 2010).

The next generation of space-based telescopes, the James Webb Space Telescope (JWST), will not have the optical filters equivalent to those found in the Hubble Space Telescope (HST), making it more difficult to search for Cepheids beyond ∼40 Mpc to increase the sample of SNe Ia calibrators, which is essential to improve the precision of H0. Cepheids are typically detected in the optical, where they have amplitudes ∼1 mag. Their near-infrared (NIR) amplitude variations are much smaller (∼0.3 mag) and they lose their characteristic optical saw-toothed light-curve shape, which helps in their identification. This makes them more difficult to identify as variable stars at NIR wavelengths. As an alternative, Jang & Lee (2017) have used the tip of the red giant branch (TRGB) observed with HST to check the Cepheid distances in nearby hosts, finding good agreement. However, because the TRGB is ∼2.5 mag less luminous than Cepheids in the optical and ∼0.5 mag less luminous in the NIR, this method will not be able to measure distances within the same volume as a Cepheid distance ladder.

A possible solution to this problem is the use of Mira variables to measure extragalactic distances. Miras are highly evolved asymptotic giant branch (AGB) stars. They do not follow a tight period–luminosity relation (PLR) in the optical. However, in the NIR, where they do follow a tight PLR, a 300-day Mira is roughly comparable in brightness to a 30-day Cepheid. They can provide an alternative distance indicator, allowing a simultaneous check of Cepheid distances and increasing the SNe Ia sample with JWST. Miras are long-period (P ≳ 100 days), large amplitude (ΔV > 2.5 mag, ΔI > 0.8 mag) M or later spectral-type pulsating variable stars (Kholopov et al. 1985; Soszyński et al. 2009b). They are divided into oxygen- and carbon-rich spectral classes based on surface chemistry, with oxygen-rich (O-rich) Miras having a carbon-to-oxygen C/O ratio <1 and a carbon-rich (C-rich) Miras having C/O > 1. AGB stars with C/O ∼ 1 are known as S stars. All stars are thought to enter the AGB phase as O-rich stars, but some evolve into C-rich stars due to dredge-up events (Iben & Renzini 1983). The O-rich Miras have been shown to follow a tight (σ ∼ 0.14 mag) PLR in the K-band (Whitelock et al. 2008; Yuan et al. 2017a) that is comparable to the scatter in the Cepheid PLR in that bandpass (Macri et al. 2015). The relation of Miras in the amplitude-period space relative to other classes of variables stars, including other distance indicators like Cepheids and RR Lyrae, is shown in Figure 1. As can be seen in the figure, they are distinguished from other variable stars by their large amplitudes and long, year-scale periods.

Figure 1. Refer to the following caption and surrounding text.

Figure 1. Distribution of I-band amplitude and period for LMC variable stars discovered by the Optical Gravitational Lensing Experiment (OGLE) III Survey (Soszynski et al. 2008a; Soszyński et al. 2008b, 2009a, 2009b, 2009c; Udalski et al. 2008; Poleski et al. 2010b, 2010a). Abbreviations are as follows: Cep—Classical Cepheid; DPV—Double Period Variable; DSCT—δ Scuti variable; RCB—R Coronae Borealis variable; RRLyr—RR Lyrae; T2Cep—Type II Cepheid; aCep—anomalous Cepheid; Mira—Mira; OSARG—OGLE Small Amplitude Red Giant; SRV—Semi-regular Variable. Miras and Semiregular Variables (SRVs) are separated on the basis of I band amplitude of variation, but this distinction is somewhat arbitrary (Soszynski et al. 2005). Miras, Cepheids, Type II Cepheids, RR Lyrae, and SRVs are radially pulsating variables that follow period–luminosity relations.

Standard image High-resolution image

Miras have a few advantages as distance indicators over Cepheids. Stars of a wide range of stellar masses go through the AGB phase, but Mira progenitors are typically of low-to-intermediate mass (∼1 M), while Cepheids have intermediate to high-mass progenitors (>5 M). Due to the bottom-heavy distribution of the stellar initial mass function (IMF), Mira progenitors are much more common than Cepheid progenitors. In addition, as seen in Figure 2, their NIR amplitudes are about twice as large as NIR Cepheid amplitudes, making it much easier to discover these with infrared-only observatories like JWST. Since they are older stars, Miras can also be found in galaxies without current star formation. This would allow the calibration of SNe Ia luminosities in early-type host galaxies with Miras.

Figure 2. Refer to the following caption and surrounding text.

Figure 2. Distribution of H-band amplitude and period for LMC variable stars. Cepheid (Macri et al. 2015) and Type II Cepheid (Bhardwaj et al. 2017) points use data from the LMC Near-Infrared Synoptic Survey. Abbreviations are the same as those for Figure 1. H-band amplitudes of Miras are from Yuan et al. (2017b) and are estimated based on three NIR epochs. The boxes show the target selection of O-rich Miras and the locations of SRVs.

Standard image High-resolution image

To test the efficacy of Miras as standard candles and demonstrate the feasibility of discovering and characterizing Miras with HST and JWST, we conducted a year-long, 12-epoch search for Mira variables in the megamaser-host galaxy NGC 4258 using WFC3 F160W data. Our goals are to develop criteria for identifying and classifying Miras using primarily NIR photometry, to provide an initial test of the relative distances from Cepheids, and to obtain a calibrated Mira PLR relation relative to the Large Magellanic Cloud (LMC) by using the water megamaser distance to NGC 4258.

Throughout the paper, we refer to the peak-to-trough variation in a Mira's magnitude over the course of one cycle as its "amplitude." Statistical and systematic uncertainties are represented as σr and σs or with the subscripts r and s respectively.

2. Observations, Data Reduction, and Photometry

2.1. Observations

With a roughly monthly observational cadence, we used 12 epochs of HST WFC3 data (GO 13445; PI: Bloom) collected between 2013 October 2 and 2014 August 3 to search for Miras in NGC 4258. The field was chosen to overlap with the NGC 4258 "inner" field from Macri et al. (2006), who discovered ∼50 Cepheids (see below for a description of these observations). The term "inner" references the field's proximity to the nucleus of NGC 4258. The observations were centered at R.A. = 12h18m52fs800 and decl. = +47°20'19farcs70 (J2000.0).

All 12 epochs have F160W and F125W data, but the first epoch is comprised of four 703s exposures in both the F160W and F125W (HST J-band). The other epochs each contain four 553s F160W exposures and one 353s F125W exposure.

This field was previously observed in F160W (GO 11570; PI: Riess) in 2009 December and 2010 May. These earlier images were taken to create a mosaic of NGC 4258 in F160W and to follow-up on a number of long-period Cepheids discovered in the galaxy from the ground. This provided a thirteenth epoch of observation for most objects, while a few in the overlapping regions of the mosaic were observed twice and had 14 epochs.

As already mentioned, this field was previously observed at optical wavelengths using HST (GO 9810; PI: Greenhill Macri et al. 2006). These observations were taken with the Advanced Camera for Surveys (ACS; Ford et al. 2003) in F435W, F555W, and F814W. The optical time-series consisted of 12 epochs between 2003 December 5 and 2004 January 19 with an observation spacing that followed a power law to allow for the detection of Cepheids at the largest possible range of periods. A summary of all of the observations used in the analysis is shown in Table 1 and transmission curves for every filter are shown in Figure 3.

Figure 3. Refer to the following caption and surrounding text.

Figure 3. Transmission curves for every HST filter in the images we analyzed shown with example O-rich and C-rich Mira spectra of the stars WW Sco and BH Cru, respectively, from Lançon & Wood (2000). In the bottom image, the continuum emission has been subtracted out to better show the difference in spectral features. The F435W, F555W, F814W, F125W, and F160W filters roughly correspond to B, V, I, J, and H in ground-based filter systems.

Standard image High-resolution image

Table 1.  HST Observations Used in This Work

Epoch Proposal ID Camera UT Date Exposure Time (s)
        F435W F555W F814W F125W F160W
H-01 13445 WFC3 2013 Oct 02 2812 2812
H-02 13445 WFC3 2013 Nov 17 353 2212
H-03 13445 WFC3 2013 Dec 03 353 2212
H-04 13445 WFC3 2013 Dec 22 353 2212
H-05 13445 WFC3 2014 Jan 26 353 2212
H-06 13445 WFC3 2014 Feb 11 353 2212
H-07 13445 WFC3 2014 Mar 14 353 2212
H-08 13445 WFC3 2014 Apr 10 353 2212
H-09 13445 WFC3 2014 May 11 353 2212
H-10 13445 WFC3 2014 Jun 08 353 2212
H-11 13445 WFC3 2014 Jul 05 353 2212
H-12 13445 WFC3 2014 Aug 03 353 2212
O-01 9810 ACS/WFC 2003 Dec 06 1800 1600 800
O-02 9810 ACS/WFC 2003 Dec 07 1800 1600 800
O-03 9810 ACS/WFC 2003 Dec 08 1800 1600 800
O-04 9810 ACS/WFC 2003 Dec 09 1800 1600 800
O-05 9810 ACS/WFC 2003 Dec 11 1800 1600 800
O-06 9810 ACS/WFC 2003 Dec 13 1800 1600 800
O-07 9810 ACS/WFC 2003 Dec 16 1800 1600 800
O-08 9810 ACS/WFC 2003 Dec 20 1800 1600 800
O-09 9810 ACS/WFC 2003 Dec 24 1800 1600 800
O-10 9810 ACS/WFC 2003 Dec 31 1800 1600 800
O-11 9810 ACS/WFC 2004 Jan 08 1800 1600 800
O-12 9810 ACS/WFC 2004 Jan 19 1800 1600 800
M-13 11570 WFC3 2009 Dec 17 2012
M-14 11570 WFC3 2009 Dec 17 2012
M-13b 11570 WFC3 2010 May 29 2012

Note. Exposure times are rounded to the nearest second. The epochs labeled represent the results of three campaigns to observe NGC 4258. The H-xx epochs were aimed at discovering Miras. The O-xx epochs were observed previously to search for Cepheids. The M-xx epochs were observed in order to create a mosaic image of the whole galaxy. We have no epochs with simultaneous optical and infrared observations.

Download table as:  ASCIITypeset image

2.2. Data Reduction

We use pipeline-processed images downloaded from the Canadian Astronomy Data Centre. For the F125W and F160W images taken between 2013 October and 2014 August, we generated drizzled and stacked images for each epoch and filter using v.1.1.16 of Astrodrizzle (Gonzaga et al. 2012). Each image contained four subpixel dither positions. The images were drizzled to a pixel scale of 0farcs08/pix instead of the 0farcs13/pix scale of WFC3 IR. Because the baseline of observations of the field spanned almost a year, the roll angle of the camera changed by 282° over the course of our observations. We choose the first epoch as the reference image and aligned all of the subsequent images onto it using DrizzlePac (Gonzaga et al. 2012).

We used DAOMATCH and DAOMASTER, kindly provided by P. Stetson, to match sources in common between Macri et al. (2006) and our F160W master image.

2.3. Photometry and Calibration

Given that our fields are quite crowded, we used tools specifically designed for crowded fields: DAOPHOT/ALLSTAR (Stetson 1987) and ALLFRAME (Stetson 1994). Our DAOPHOT procedure is different from the R16 NIR forced photometry because we conduct a search for Miras in WFC3 F160W and do not know their positions a priori. We created point-spread functions (PSFs) in DAOPHOT with F160W and F125W exposures of the standard star P330E. Due to the crowded nature of our fields, we were unable to find enough isolated stars to make PSFs using sources in our field. Aperture corrections, discussed in Section 2.4, were used to account for imperfections in the PSF model.

We stacked all of the F160W observations to make a deeper "master" image. Then we used the DAOPHOT routine FIND to detect sources with a greater than 3σ significance level in standard deviations from the sky background noise. Then the DAOPHOT routine PHOT was used to perform aperture photometry. We input the star list produced from the aperture photometry into ALLSTAR for PSF photometry, optimized for a 2.5 pixel full width at half maximum. We then repeated 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 were then concatenated to create a master source list of ∼1.3 × 105 entries.

We input this master source list into ALLFRAME. ALLFRAME is similar to ALLSTAR, except that it is capable of performing simultaneous fits to the profiles of all of the stars contained in all of the images of the same field. It then produces time-series PSF photometry as output. We used the master source list as the input for fitting all of the F160W epochs at the same time.

We then searched for secondary standards in the star lists by choosing bright objects that had been observed in all 12 of the last F160W epochs. We visually inspected the stellar profiles and their surroundings to choose secondary standards that were relatively isolated compared to other bright stars, removing any that showed variability or had large photometric errors. This left us with a total of 44 sources, summarized in Table 2. We calculated the celestial coordinates for all of these sources using the astrometric solutions in the FITS headers and PyAstronomy program pyasl. The mean residuals for these stars across all epochs of F160W imaging exhibited a dispersion of 0.01 mag, which is corrected for during the ALLFRAME photometry.

Table 2.  Secondary Standards

ID R.A. Decl. X Y F160W Error
  (J2000.0) (J2000.0) (Pixels) (Pixels) (mag) (mag)
49680 12 18 46.910 +47 19 58.43 987.352 1019.138 19.035 0.013
23564 12 18 48.964 +47 20 36.93 877.748 482.696 19.090 0.010
4663 12 18 50.894 +47 21 00.70 874.317 97.506 19.256 0.010
31187 12 18 51.819 +47 19 56.35 1481.836 637.801 19.298 0.011
35864 12 18 47.198 +47 20 25.75 795.737 734.066 19.309 0.012
51997 12 18 46.957 +47 19 53.17 1034.231 1065.619 19.389 0.010
44141 12 18 47.862 +47 20 02.36 1048.415 903.643 19.574 0.012
6657 12 18 50.314 +47 21 01.32 813.004 139.018 19.577 0.011
41446 12 18 47.622 +47 20 10.24 961.708 847.871 19.594 0.013
4234 12 18 51.396 +47 20 57.28 950.684 89.193 19.625 0.018
3655 12 18 51.193 +47 21 00.35 906.218 76.444 19.657 0.016
8400 12 18 50.221 +47 20 58.51 826.487 173.487 19.677 0.010
45991 12 18 46.740 +47 20 07.89 894.789 942.392 19.694 0.010
4662 12 18 50.548 +47 21 03.66 816.930 97.443 19.699 0.010
20550 12 18 48.903 +47 20 43.91 815.702 420.941 19.715 0.013
27518 12 18 48.969 +47 20 28.55 945.605 562.518 19.832 0.010
57524 12 18 48.801 +47 19 25.55 1435.620 1179.302 19.878 0.011
23751 12 18 48.816 +47 20 37.72 856.961 487.319 19.890 0.011
42544 12 18 48.425 +47 20 00.99 1114.240 870.805 19.913 0.013
7460 12 18 49.568 +47 21 06.00 702.760 155.089 19.918 0.011
3896 12 18 50.963 +47 21 01.74 872.693 81.896 19.925 0.016
10561 12 18 49.443 +47 21 00.47 735.012 218.259 19.937 0.012
4566 12 18 50.618 +47 21 03.22 827.237 95.984 19.991 0.009
22541 12 18 49.572 +47 20 33.95 960.815 461.519 20.052 0.013
14188 12 18 50.626 +47 20 42.79 992.267 290.814 20.057 0.011
695 12 18 52.032 +47 20 59.76 992.657 13.522 20.060 0.011
45068 12 18 48.285 +47 19 56.76 1134.596 922.657 20.150 0.012
8203 12 18 49.034 +47 21 09.03 626.493 169.699 20.212 0.009
42205 12 18 48.214 +47 20 03.52 1073.341 863.826 20.237 0.016
65132 12 18 45.898 +47 19 33.70 1087.567 1338.530 20.275 0.010
12889 12 18 50.029 +47 20 50.58 871.554 265.064 20.293 0.010
36017 12 18 48.265 +47 20 16.35 975.233 736.852 20.295 0.010
15965 12 18 48.742 +47 20 55.00 710.891 327.951 20.448 0.012
2890 12 18 51.222 +47 21 01.68 898.406 61.263 20.467 0.010
57706 12 18 49.091 +47 19 22.73 1486.569 1182.616 20.514 0.013
40960 12 18 44.791 +47 20 35.40 483.976 838.315 20.575 0.009
59832 12 18 44.802 +47 19 54.64 812.577 1227.611 20.609 0.009
59266 12 18 48.736 +47 19 22.40 1454.737 1214.741 20.642 0.010
4044 12 18 50.747 +47 21 03.30 839.175 84.681 20.673 0.010
48800 12 18 46.408 +47 20 04.56 889.228 1001.389 20.750 0.011
36245 12 18 45.236 +47 20 41.70 476.671 741.707 20.830 0.015
64758 12 18 46.901 +47 19 25.91 1247.840 1331.200 20.886 0.010
43584 12 18 43.617 +47 20 39.71 335.074 892.949 21.029 0.010
52996 12 18 45.086 +47 20 06.86 742.139 1087.443 21.039 0.011

Note. A list of all the secondary sources used to calibrate the F160W light curves. ID numbers are photometry IDs, X and Y positions are relative to the first epoch of the F160W image, and the errors are photometric errors as estimated by DAOPHOT.

Download table as:  ASCIITypeset image

2.4. Aperture Corrections

We use aperture corrections to account for missing flux from imperfections in our PSF model. Using the standard stars chosen for the variability search, we subtract everything but these sources from the master image using the master source list with DAOPHOT's SUB task. The standard stars are already relatively bright and isolated for our field, but this subtraction helps to ensure that we remove any additional flux from the wings of the standard stars. We then perform aperture photometry on the sources using increasing aperture radii (up to 0farcs4) and check that the growth curves look well behaved over a range of apertures.

We calculated the difference between PSF and 0farcs4 aperture magnitudes, to which we added the flux beyond this limit previously calculated by STScI. The WFC3 F160W Vegamag zeropoint was 24.5037. The overall correction from PSF to "infinite aperture" magnitudes was 0.023 ± 0.01 mag.

3. Mira Selection Criteria

RR Lyrae, Classical Cepheids, and Type II Cepheids, among other variable stars, can be identified with the use of light-curve templates, which exist for both the optical and NIR (Jones et al. 1996; Yoachim et al. 2009; Sesar et al. 2010; Inno et al. 2015; Bhardwaj et al. 2017) since the shapes of their light curves vary as a function of period in a predictable way. The templates allow them to be identified as a particular class of variable star by light-curve morphology after they have been initially identified as variable. However, Miras have irregular light-curve shapes that are not strictly a function of period. They may also exhibit variations in light-curve shape or amplitude between different cycles. These cycle-to-cycle variations can be seen in the visual light curve of the prototype Mira, oCeti, shown in Figure 4. Figures 7 and 9 of Whitelock et al. (1994) and Figure 1 of Olivier et al. (2001) show the smaller cycle-to-cycle variations present in K-band Mira light curves.

Figure 4. Refer to the following caption and surrounding text.

Figure 4. Visual light curve of oCeti (Mira) using data from the American Association of Variable Star Observers (AAVSO). The light curve shows cycle-to-cycle variations in both amplitude and mean magnitude.

Standard image High-resolution image

Miras and other long-period variables (LPVs) can also have a longer, secondary pulsational period in addition to their primary pulsational periods. These long secondary periods (LSPs; Payne-Gaposchkin 1954; Houk 1963; Nicholls et al. 2009) are typically about an order of magnitude longer than the primary periods (Wood et al. 1999) and are found in about 25%–50% of all LPVs (Wood et al. 1999; Percy et al. 2004; Soszyński 2007; Fraser et al. 2008). Several theories have been proposed to explain LSPs, but there is no general consensus (Soszyński 2007; Saio et al. 2015; Percy & Deibert 2016). Both LSPs and cycle-to-cycle variations will result in the same Mira having different magnitudes at the same phase in different cycles.

Instead of using template fitting, Mira variables are typically identified only by their large V and I amplitudes and long periods. The amplitude criterion is used to separate Miras from the more numerous, lower-amplitude, and sometimes more inconsistent semi-regular variables (SRVs), which are another type of LPV (see Figures 1 and 2). The amplitude cutoffs have traditionally been defined in the optical as ΔV > 2.5 mag or ΔI > 0.8 mag (Kholopov et al. 1985; Soszyński et al. 2009b) and the periods range from 80 to 1000 days, though there are some Miras with periods that can be significantly longer. A few previous studies also identified Miras in the infrared by using J, H, and K time-series data (Whitelock et al. 2006, 2009; Matsunaga et al. 2009) but large amplitudes in the NIR do not always correspond to large amplitudes in the visible bands. A sample of Miras selected based on large NIR amplitudes could contain objects that would not be selected based on visual criteria and vice versa.

Yuan et al. (2017b) sought to differentiate between these classes of variable stars by using a Random–Forest classifier that incorporated period, light-curve shape (O-rich Miras have more symmetric light curves), and other properties into the classification. However, our small number of epochs limits us to using simple cuts. Since we have only limited optical data (at most spanning ∼40% of a Mira cycle due to the short baseline optimized to detect Cepheids) and only one band in the NIR with time-series photometry, we need selection criteria that rely more heavily on NIR measurements.

We also note that not all Miras are good distance indicators. The two C-rich and O-rich subgroups follow different PLRs in the H-band (Ita et al. 2004; Ita & Matsunaga 2011). C-rich stars can also develop optically thick circumstellar dust shells that result in them appearing fainter even in the NIR (Yuan et al. 2017b). Therefore, we use only O-rich Miras as distance indicators and must separate them from C-rich Miras.

The photospheric differences between C- and O-rich Miras are thought to be due to the differing ratios of carbon and oxygen in their atmospheres. The carbon and oxygen in a star's atmosphere will combine to make the stable CO molecule until there is no more of the less abundant element. The excess carbon or oxygen will then be left over for dust formation and will combine to create molecules such as TiO and VO in O-rich stars or CN and C2 in C-rich stars (Cioni et al. 2001). These molecules define each spectral type and can also change a Mira's color, which has potential consequences that are discussed in greater detail in Section 3.4.3. While there are generally J − K distinctions in color between C- and O-rich stars, the cutoff in color varies in different galaxies (Cioni & Habing 2003). In addition, some O-rich stars may be very red, as in the case of OH/IR stars.

We expect to encounter a smaller ratio of C-rich Miras to O-rich Miras in NGC 4258 than in the LMC. Higher ratios of O-rich to C-rich stars are observed in galaxies with higher metallicity (Blanco & McCarthy 1983; Mouhcine & Lançon 2003; Hamren et al. 2015). The inner field of NGC 4258 is expected to be ∼0.1 dex more metal-rich than the LMC, though still approximately −0.2 dex relative to solar (Bresolin 2011).

3.1. Detection of Variability

We first used the Welch–Stetson variability index L (Stetson 1996) to identify a sample of variable objects detected in all of the F160W epochs. This is a combination of two other measures of variability (all three are defined in Equations (1)–(3) of Stetson 1996). Objects that have larger variances from their mean magnitudes, exhibit similar variability in multiple images taken at around the same time, and have non-Gaussian magnitude distribution will have a larger L value. We calculated L for every object identified in the master photometry list and then further inspected a number of objects that met our threshold of L > 1.75, about 3σ above the mean L for all objects.

Figure 5 shows the distribution of L values for candidate Miras and all sources in NGC 4258. We kept only sources with L ≥ 1.75 that were detected in all 12 epochs of the most recent F160W observations. Sources that did not show periodicity (were only continuously rising or decreasing light curves), had ΔF160W < 0.4 or periods of less than 100 days were then removed from the list of possible Miras. The ΔF160W > 0.4 cut roughly corresponds to the ΔI > 0.8 used by Soszyński et al. (2009b) to distinguish between Miras and SRVs. An additional discussion of the appropriate minimum amplitude cut can be found in Section 3.4.2. These requirements resulted in 3951 Mira candidates remaining in our sample. An example of an F160W Mira light curve and its finding charts is shown in Figure 6.

Figure 5. Refer to the following caption and surrounding text.

Figure 5. Distribution of the Welch–Stetson variability index L for all of the sources in the field (shown in white) and for the subset of objects that made it past the initial visual inspection (in blue). All objects with L ≥ 1.75 were visually examined and anything that did not pass visual inspection or had ΔF160W < 0.4 was automatically discarded.

Standard image High-resolution image
Figure 6. Refer to the following caption and surrounding text.

Figure 6. Light curve and finding charts for a 266-day Mira from NGC 4258. The Mira is the source in the blue ring in each of the finding charts. Numbers on each stamp represent the date of the observation, corresponding to MJD-50000. The vertical axis of the inset plot gives the F160W magnitude and the horizontal axis marks the time since the first epoch of observations in days. Each stamp corresponds to one point on the light curve, starting in the upper left corner and going clockwise around the image. As the Mira progresses through its light cycle, the brightness of the sources in the postage stamps noticeably changes.

Standard image High-resolution image

3.2. Estimating Cycle-to-cycle Variation

We estimated the level of cycle-to-cycle variations in the H-band mean magnitude we might measure using data from Matsunaga et al. (2009), which looked at Miras in the Galactic bulge. These Miras were observed in 12 fields in the NIR between 2001 and 2008, with observations of the same phases during various cycles for each Mira. While most of our observations took place over a span of only ∼300 days, we still need to account for possible variations in magnitude at a given phase between the main F160W/H band campaign in 2013 and 2014 and the observations taken in 2009, many oscillation cycles earlier.

The Matsunaga et al. (2009) data set covers different objects than the OGLE-III data set and is more sparsely sampled, but it is one of the largest sets of time-series observations of Miras in the NIR. Some of the OGLE-III Miras have been observed in the NIR as well as the optical bands used by OGLE, but the vast majority of the OGLE-III Miras have only a few or single epochs in the NIR. Because of this, we did not use the OGLE-III data set to estimate NIR cycle-to-cycle variations.

We binned the Galactic bulge data of each Mira candidate observed by Matsunaga et al. (2009) into 30-day bins (to simulate the frequency of our observations) and calculated the H-band mean magnitude in each bin. Next we folded the binned magnitudes by the periods measured in Matsunaga et al. (2009) to obtain their phases. We then binned the resulting points by phase to see how bright each Mira was at that particular phase over different cycles. Finally, we calculated the variance of the points in each phase bin to estimate the cycle-to-cycle variations in magnitude at a similar phase for the observations. Figure 7 shows the cycle-to-cycle variation of mean magnitudes as a function of Mira period. As expected, the shorter-period Miras have larger "cycle-to-cycle" variations because each 30-day bin averages over a larger range of phases. The median cycle-to-cycle variation overall was 0.072 mag, which we incorporated as an additional error added in quadrature when fitting periods using data from different cycles.

Figure 7. Refer to the following caption and surrounding text.

Figure 7. Distribution of cycle-to-cycle variations in Galactic bulge Miras in the H-band using data from Matsunaga et al. (2009). Blue points show the cycle-to-cycle variation for each Mira studied, as a function of the Mira's period, red points are the means of each period bin along with the standard deviations of each bin. The bins illustrate the difference in the average cycle-to-cycle variation for Miras of different periods, but only the median of the cycle-to-cycle variations for all Miras was used as an additional error when fitting periods.

Standard image High-resolution image

3.3. Determination of Periods

We used a two-step method to measure the periods of the 3951 Mira-like objects remaining after the previous cuts based on variability and preliminary ΔF160W amplitude. Unlike the previous steps, which used only the 12 recent epochs of F160W data, we incorporated the additional epochs of observation from 2009 and early 2010 into our analysis at the next stage to get more accurate periods.

First, we calculated a Lomb–Scargle periodogram for each variable source. Then the peaks of the Lomb–Scargle were each fit with a Fourier series up to the third order. Well-sampled Mira light curves from the LMC have shown that Miras can have higher-order harmonics but more could not be fit due to the limited number of available observations. We folded the light curves by each potential period and then fit every folded light curve using a Fourier series. We used a Bayesian information criterion (BIC) to determine if adding another harmonic to the fit was significant before increasing the number of harmonics. The BIC is defined as

where is the maximized value of the likelihood function for the estimated model, n is the sample size, and k is the number of free parameters to be estimated. We assumed that the errors were Gaussian and uncorrelated, and thus was equal to the error variance of the fit with k parameters.

Given the limited number of epochs, the BIC indicated that a simple sine function is most appropriate for almost all of the Miras. Finally, we used the Fourier fit parameters and period as initial guesses for a fit to the data using Levenberg–Marquardt least-squares curve-fitting. At this point, all of the parameters were fit simultaneously.

To verify that the periods were correct, we visually inspected each P < 350 -day Mira candidate light curve and checked its fit to a sinusoid using the period determined earlier. Any periods that did not produce a good fit to the data were either refit by removing outliers and overriding the original fit, or, if a good fit could not be found, flagged as a lower-quality object and not used in the analysis.

Since the baseline of our observation was only 305 days, we only have at least one cycle of good phase coverage for Mira candidates with periods less than 305 days. To calculate the period recovery rate, we used LMC Mira observations as templates and added photometric uncertainties and the observation sampling that reflected the NGC 4258 data set. We tested both the case of 13 epochs of observation (incorporating in the earlier epoch from 2009) and 12 epochs of observation. We then measured the periods of our sample Miras and considered every period measured to within 30 days of the true input period as recovered. For both, we found that the recovery rate of Mira periods less than 300 days was approximately 90%. Figure 8 shows the results of the simulation.

Figure 8. Refer to the following caption and surrounding text.

Figure 8. Fraction of Mira periods recovered as a function of period. The range of periods for the input Miras (from LMC Mira observations) ranged from 103 to 675 days. Miras were considered "recovered" if their measured periods were within 30 days of their true periods. A dashed line is shown at a recovery fraction of 0.5 to guide the eye. For the vast majority of our sources, with 13 epochs of observation, the recovery rate ∼90%.

Standard image High-resolution image

3.4. Samples and Selection Criteria

Our goal is to recover samples of the most secure Miras, rather than the most Miras, since the statistical uncertainty on the zeropoint of a PL due to our sample size will already be much smaller than the systematic error. We created three samples of Miras, which we have called Gold, Silver, and Bronze based on varying degrees of confidence in classification. Each sample contains predominantly O-rich Miras, but the Bronze sample relies only on NIR information for classification, while the Gold and Silver samples are further vetted using our short time series of optical data. This makes the Bronze sample a good test case for future NIR-only observations of Miras. The criteria for each sample are given in Table 3. We outline the consequences and motivation for each of these criteria in the following sections.

Table 3.  Mira Sample Criteria

  Bronze Silver Gold
Period Cut: P < 300 days P < 300 days P < 300 days
Amplitude Cut: 0.4 mag < ΔF160W < 0.8 mag 0.4 mag < ΔF160W < 0.8 mag 0.4 mag < ΔF160W < 0.8 mag
Color Cut: mF125W − mF160W < 1.3 mF125W − mF160W < 1.3 mF125W − mF160W < 1.3
F814W Detection: F814W detection Slope-fit to F814W data > 3σ
F814W Amplitude: ΔF814W > 0.3 mag

Download table as:  ASCIITypeset image

3.4.1. Period Cut

Because the baseline of our observations spanned only 305 days, we kept only Miras with periods less than 300 days. We also tested shorter upper period limits (down to 250 days) and found that there was no effect on the zeropoint. Choosing a more restrictive upper limit for the period cut will also make our sample less representative of the objects we would normally find in SNe Ia host galaxies. NGC 4258 is relatively close compared to most supernova hosts, so we will have an easier time finding sources at the longer end of the period range. A 300-day Mira is ∼0.5 mag more luminous than a 250-day Mira and as a result allows us to look in volume twice as large.

In addition to avoiding contamination from periods of over 300 days that may be unreliable, limiting our sample to shorter-period Miras also has additional benefits for classification. The period cut reduces the number of C-rich stars in our sample since C-rich Miras typically have longer periods than O-rich Miras as discussed in Section 3.3 (see Figure 2). A distribution of C- and O-rich Mira periods in the LMC is shown in Figure 9.

Figure 9. Refer to the following caption and surrounding text.

Figure 9. Distribution of periods for LMC O- and C-rich Mira variables. C-rich Mira periods are shown in red, the O-rich Mira periods are shown in blue. The black dashed line represents the period cut of 300 days.

Standard image High-resolution image

While this period cut will exclude longer-period O-rich Miras, some of these also do not make good distance indicators because they can be hot-bottom-burning stars (HBB). The onset of HBB depends on both mass and metallicity, but it is typically thought to occur in stars with initial masses greater than 4–5 M that are near the end of their AGB phases (Glass & Lloyd Evans 2003). Whitelock et al. (2003) showed that HBB Miras deviate from a linear PLR, making them poor candidates for distance indicators.

3.4.2. F160W Amplitude

Miras have generally been selected on the basis of their large optical (V- and I-band) amplitudes to distinguish them from SRVs. SRVs can be as consistent as Miras in their variability and have similar periods but have smaller amplitudes than Miras. SRVs can also fall on the same PLR as Miras or on various other parallel PL relations (Wood et al. 1999; Trabucchi et al. 2017), depending on their pulsation mode. In general, SRVs are brighter than Miras with the same period and thus can bias the PLR if they are not removed from the final Mira sample. Previous studies of Miras (Whitelock et al. 2008; Matsunaga et al. 2009) have suggested a minimum peak-to-trough variation of ΔJ, ΔH, ΔK ∼ 0.4 mag to classify a variable as a Mira. Thus, we have used ΔF160W > 0.4 as the cutoff for minimum change in brightness over one cycle.

In addition to removing SRVs, this minimum amplitude cut also allows us to remove constant stars and blended objects, which would not follow a PLR at all. For a variable star like a Mira or a Cepheid, the resulting blend will have a different color and amplitude from the original star in addition to being more luminous.

O- and C-rich Miras can also have different amplitude distributions. Cioni et al. (2003) found that C-rich Miras had larger optical amplitudes on average than O-rich Miras in the Small Magellanic Cloud (SMC). Yuan et al. (2017b) found that this was also the case for O-rich Miras in the LMC, especially when considering the amplitude distribution over many cycles. Over the course of a single cycle, C-rich Miras usually have larger amplitudes. This is caused by C-rich Miras having longer periods and thicker dust shells on average compared to O-rich Miras. Both longer-period Miras and heavily reddened stars are more likely to have larger amplitudes. While there is considerable overlap in the distribution of amplitudes for O- and C-rich Miras (as shown in Figure 10) the largest-amplitude objects are usually C-rich.

Figure 10. Refer to the following caption and surrounding text.

Figure 10. Left: the amplitude and period relationship for LMC Miras. Red points are C-rich Miras and blue points are O-rich Miras as classified by Soszyński et al. (2009b). They divided the Mira sample along the WI vs. WJK plane, where WI and WJK are Wesenheit indices using optical and NIR magnitudes, respectively. Amplitude information of the LMC Miras comes from Yuan et al. (2017b) with errors on the amplitude estimated to be ∼0.13 mag. The dashed vertical and horizontal lines represent the maximum period (300 days) and amplitude (0.8 mag) cuts respectively. Under our selection criterion, only objects within the bottom-left quadrant would have made it into the final Mira sample. Right: the same plot constructed for NGC 4258 Miras with the objects in the Bronze sample in blue, and all the rejected Miras in red. Period and amplitude beyond 300 days are unreliable but have been plotted with their best estimates shown for comparison. Errors on amplitudes greater than 300 days are estimated to be ∼0.2 mag. Periods below 300 days were verified by visual inspection.

Standard image High-resolution image

The left half of Figure 10 shows the distribution of LMC Miras as a function of period and ΔH using data from Yuan et al. (2017b). These Miras were classified by the OGLE team using a WIWJK diagram, where WI and WJK are the Wesenheit indicies in the optical and NIR, respectively (Madore 1982). Both types of Miras have estimated uncertainties in the amplitude of ∼0.13 mag but different distributions in amplitude. Since we are interested in obtaining a clean sample of O-rich Miras, we employ a cut of ΔF160W < 0.8 mag as a maximum amplitude cut in addition to the >0.4 mag minimum amplitude. The corresponding plot of accepted and rejected Miras in NGC 4258 is shown in the right half of Figure 10, with our "Bronze" Miras (Miras that met all of the NIR criteria) shown in blue. Some objects in the same quadrant as the Bronze sample were rejected on the basis of their uncertain periods. However, we expect that there should be a few percent overall contamination from C-rich Miras in this quadrant in the Gold sample.

In the LMC, the ratio of C-to-O Miras in this quadrant is 1/3, and the overall C-to-O ratio in the LMC is 3/1. The C-to-O ratio in the solar neighborhood, which is more similar environment to the inner field of NGC 4258 is ∼1 (Ishihara et al. 2011). This suggests an ∼10% contamination rate from C-rich Miras from using these two cuts alone. Combining these with the optical observations and the color cut, which both exclude the reddest objects, we estimate that the contamination in the Gold sample after all of the cuts should be a few percent.

3.4.3. Color

We calculated the F125W − F160W color by creating a stacked "master" image in each bandpass where each epoch was weighted evenly (despite the first epoch in both bands having a longer observation time). The F125W − F160W color was then measured from the fluxes of the objects in these two images.

Our final NIR cut uses F125W − F160W color. We see that the majority of C-rich Miras can be removed by employing a color cut of J − H < 0.9, as shown in Figure 11 (adopted from Soszyński et al. 2009b). Color cuts to separate O- and C-rich Miras are physically motivated by the differing opacity in the NIR and mid-infrared (MIR) in C- and O-rich stars. It is most effective to use both NIR and MIR colors for distinguishing between the two groups; broad NIR bands alone have not been shown to be sufficient to separate C- and O-rich Miras (Le Bertre et al. 1994). Medium-band filters on HST that target unique spectral features have been used to separate C- and O-rich stars in M31 (Boyer et al. 2013, 2017), but are not efficient for identifying Miras in more distant galaxies because they would require much longer integration times.

Figure 11. Refer to the following caption and surrounding text.

Figure 11. Top: the color–magnitude diagram for Miras in the LMC made from cross-matching the OGLE-III star catalog with near-infrared data from IRSF. The two subclasses of Miras were classified using a WI − WJK diagram. The black dashed line is at J − H = 0.9. Middle: the color–magnitude diagram for Mira candidates in the NGC 4258 made from comparing the mean F125W and F160W colors. The black dashed line represents J − H = 1.3. There are many more bluer variables in the NGC 4258 data set than in the LMC data set. The "Candidate O-rich Miras" are objects that have passed all of our NIR cuts except for color. Bottom: the color–magnitude diagram for Miras in M33, using UKIRT data from W. Y. Yuan et al. (2018, in preparation). The black dashed line is at J − H = 1.1. These Miras were observed from the ground and detected in the optical. C-rich Miras, which are redder, are less likely to have been detected.

Standard image High-resolution image

The J and H filters are also better suited for distinguishing between C- and O-rich Miras than F125W and F160W, which are much more similar in their transmission functions. As seen in Figure 11, the objects flagged as Mira-like (all objects both red and blue points) in our analysis do not follow the same distribution as the Miras in the LMC. The Mira-like objects in NGC 4258 appear to be one population rather than two seen in the LMC. This is most likely caused by differences in the HST filter system and the typical ground-based NIR filters. Stellar models of C- and O-rich AGB stars from Aringer et al. (2009, 2016) suggest that these two populations overlap more in F125WF160W color than in J and H color. Some of the reddest objects may also have been undetected.

To avoid cutting through the middle of our population, we used F125W − F160W < 1.3 mag as a color cut to remove the reddest objects, but anticipated that it would not fully remove C-rich Miras from our sample. We employ the F160W amplitude cuts from Section 3.4.2 and the period cut to remove the majority of C-rich Miras instead.

3.4.4. F555W and F814W Amplitudes

In addition to variability, amplitude, and color cuts based on NIR data, which were applied to all three samples, we looked for corroborating evidence for the variability of Mira-like objects in F555W and F814W observations for the Silver and Gold samples. Given photospheric temperatures of only ∼3000–3500 K, Miras are significantly less luminous at optical wavelengths. They also experience extremely large-amplitude variations in those bands, with ΔV ∼ 10 mag or greater in some cases. Table 4 shows a range of amplitudes in different bandpasses and magnitudes for a 200-day Mira.

Table 4.  Typical Mira Amplitudes and Absolute Magnitudes

Bandpass Amplitude Range Absolute Magnitude
  (mag) (200 day Mira)
V 2.5–>10 −1.4
I 0.8–3.5 −4.1
J 0.5–3.0 −5.8
H 0.4–3.0 −6.6
Ks 0.4–3.0 −7.0

Note. The range of amplitudes for Miras and the absolute magnitude of a 200-day Mira in various photometric bands. Upper limits on the amplitudes are approximate. V-band absolute magnitudes are not well known and have been estimated from OGLE data (Soszyński et al. 2009b). The other absolute magnitudes are calculated using LMC and M33 PLRs from Yuan et al. (2017a).

Download table as:  ASCIITypeset image

We estimated the recovery rate of Miras in the optical data. Using the F160W observations, we calculated the F160W phase at the time of the optical observations and assumed a phase lag between the optical and NIR phases that was dependent on period. The phase lag was calculated using data from Yuan et al. (2017b). We used

where ϕI is the I-band phase, ϕH is the H-band phase, and P is the period. We assumed that the I and H bands are roughly equivalent in phase to the F814W and F160W, respectively. We used the estimated differences in mean magnitude and amplitude from Table 4 to convert from F160W magnitudes and amplitudes to F555W and F814W. Any stars that had a signal-to-noise ratio greater than 3 were considered recovered in the simulation. This resulted in an expected recovery rate of ∼78% for F814W.

We used the DAOMATCH and DAOMASTER programs to match the NIR and optical master source lists and found that 296 out of 438 Miras (or 68%) of the Miras from the Bronze sample were matched with sources in F814W images. This is roughly in agreement with the simulated expected numbers. As anticipated, very few Miras in our sample could be matched with F555W light curves (a 3% recovery rate was predicted by the simulation, compared to 2% in reality). Thus, we used only information from the F814W observations in our selection criteria. Light curves for some of the few potential Miras with both F555W and F814W matches are shown alongside their F160W light curves in Figure 12. The optical observations precede the NIR epochs by about 10 years.

Figure 12. Refer to the following caption and surrounding text.

Figure 12. F160W (red), F814W (green), and F555W (blue) light curves for four Mira candidates. The magnitudes have been shifted in order to display every candidate Mira's light curve on the same plot axes and show the amplitude of the light curves. The horizontal axis marks the number of days since the start of each series of observations. The days for the optical light curves have been multiplied by a factor of six in order to better show the shape of the light curves. Numbers at the top of each subplot are photometry ID and calculated period of the object. As can be seen in Table 1, we did not have concurrent observations of near-infrared and optical data. Only 35 of the candidate Miras in the Bronze sample had F555W light curves. F555W data was not included in the analysis, but have been shown here for completeness.

Standard image High-resolution image

We interpret a Mira's chances of F814W detection as a function of both a Mira's F814W − F160W color and phase. We assumed that all of the Miras in the simulation had the same color, but C-rich or heavily dust-enshrouded O-rich Miras are known to be very red and would be difficult to detect in optical bandpasses. Similarly, Miras measured toward the trough of their light cycles would be more difficult to detect than Miras at the peak. However, we were unable to determine our Mira candidates' true phases at the time of the previous observations, so we could not test this directly through simulations.

In order to determine the significance of a change in magnitude over the observation baseline, we fit each F814W light curve fragment with a linear fit. We then kept only objects with at least a naive "3σ" significance in change of magnitude in the Gold sample. This allows us to check that the object we detected as variable in F160W is variable in F814W.

We also used the difference in magnitude between the first and last epoch of observation to estimate the I-band amplitude. Thus, we could determine if the sources that were variable in F160W were also variable in the optical observations. Given that the baseline of the optical observations were only 44 days and the shortest-period Miras have periods of about 100 days, the F555W and F814W light curves cover only about 15%–20% of an average Mira's oscillation and at most contain only about ∼60% of its total variation. Due to the short temporal baseline of observations, objects close to the peak or trough (if detected) of their observations will be rejected on account of their small overall variations.

We used the OGLE Galactic-bulge sample of Miras (Soszyński et al. 2013) to calculate the distribution of changes in I-band magnitude expected over a period of 44 days. We compared this with the distribution we obtained while looking at the changes in F814W magnitude over the same duration baseline for objects in our Silver sample. However, the two bandpasses are not completely identical and their distributions are slightly different, as shown in Figure 13. We also created some light curves of constant stars with photometric noise added in as a null case. The Silver sample from NGC 4258 and the Galactic bulge distribution both show significantly higher levels of I-band variation than the simulated "constant" stars over the same period, suggesting that they are true variable stars. The optical requirement also excludes the reddest stars, which are most likely C-rich.

Figure 13. Refer to the following caption and surrounding text.

Figure 13. Change in I-band magnitude for Miras in the Galactic bulge (detected with OGLE), Miras in the NGC 4258 Silver sample (248 objects), and constant stars. The three curves show the cumulative distribution for each class of objects as a function of change in I-band magnitude over 44 days (the baseline of our F814W observations).

Standard image High-resolution image

4. Systematics

In order to get an accurate relative distance and to compare our results with previous studies of Miras using ground-based observations, we needed to transform ground-based H and J magnitudes of LMC Miras into the HST F160W magnitude. We performed artificial star tests to account for excess background due to the density of sources in our field and loss of flux due to an imperfect PSF model. Each of these introduces a systematic error into the final result.

To estimate the systematic errors of each of our cuts, we varied each cut around a standard deviation of the values chosen or the range of values present in the literature when possible and looked at the effect that it had on the zeropoint of the PLR relation. Table 5 has a summary of each of these contributions.

Table 5.  Systematic Uncertainties

Systematic Uncertainty
F160W Amplitude Cut 0.042
F160WF125W Color Cut 0.024
Color Correction Term 0.038
Intrinsic Scatter 0.021
Total 0.06

Note. The approximate contribution to the total systematic error of the gold sample from each systematic.

Download table as:  ASCIITypeset image

4.1. Slope

Miras and other variable stars are typically fit with a linear PLR (called the Leavitt Law for Cepheids). However, there is some evidence for break at 10 days for Cepheids and previous Mira observations have suggested that the Mira PLR may have a break as well, at periods of about 400 days (Ita & Matsunaga 2011). This is likely to be caused by the onset of HBB.

Yuan et al. (2017a) used a quadratic fit for the PLR instead of fitting the sample with two linear PLRs with a break. We fit each of our subsamples of Miras using their quadratic PL relations for the H band, which is the closest match to F160W:

where P is the period in days, m is the H-band magnitude, and a0 is what we have called the zeropoint. We fit for a0 and hold the other parameters fixed to the values from Yuan et al. (2017a).

This PLR was derived using observations of about 170 LMC O-rich Miras. While the F160W filter is bluer than the H band, deriving a PLR from the NGC 4258 data alone was not possible because of the much larger scatter induced by the background brightness fluctuations, which is discussed in great detail in Section 4.2. With observations of Miras in more galaxies using F160W, we could simultaneously fit the Miras in multiple galaxies and derive a more robust fit in this bandpass.

4.2. Artificial Star Tests

The high density of sources and unresolved background objects in our images will systematically bias our magnitude measurements toward brighter values. This is caused by the superposition of several point sources. We refer to this effect as crowding.

We correct for crowding using artificial star tests. Starting with a master image created from combining all 12 epochs of F160W data, we use DAOPHOT to place fake sources in the images at the same apparent magnitude as the Mira candidates. We then compare the recovered magnitudes of the artificial stars with the input magnitudes and adjust the input magnitudes to better agree with the recovered magnitudes. On average, the artificial stars were measured to be 0.25 mag more luminous than their true magnitudes.

We then repeat the steps in the photometry process up to the creation of the master source list and then compare the recovered and input magnitudes to determine the crowding correction. Because there are ∼1700 variables in the image, we add in one fake star for half of the variable stars in the image at any given time to avoid artificially raising the background of the image. The artificial stars are dropped within a 25 pixel radius of each Mira, and at least 10 pixels away from the edges of the image. Only stars that did not fall within 3 pixels of another star up to 1.5 mag fainter were used in the analysis.

After calculating the difference in the input and recovered magnitudes, we then use a 3σ-clip about the median to remove outliers. The mean difference between input and recovered magnitude for each star is then the crowding correction we apply.

4.3. Mean Magnitude Correction

We used OGLE O-rich Miras cross-matched with J and H magnitudes from the Infrared Survey Facility (IRSF) catalog (Kato et al. 2007) to determine the LMC zeropoint. The mean magnitudes for NGC 4258 were defined as the first term in the Fourier-series fit to each object. Because Mira light curves are irregular in shape, they do not spend the same amount of time at each phase in their cycles, creating a small bias in single-epoch measurements when compared to mean magnitudes.

We used Monte Carlo simulations to calculate the difference between these two estimates of mean magnitude and found that on average, the PLRs measured using the fit mean magnitudes were 0.02 fainter than the PLRs measured using single-epoch mean magnitudes. Therefore, in order to correct between the two, we added in a mean magnitude correction of −0.02 mag to our final results.

5. Results

5.1. Color Transformation

In order to compare ground-based NIR PLRs with the F160W PLR from NGC 4258, we calculated a color correction to transform the ground-based H-band Mira observations to F160W observations. Miras have heavy molecular absorption lines in the NIR compared to M-type main-sequence stars, so we used real O-rich Mira NIR spectra observed from the ground from Lançon & Wood (2000) as input to PySynPhot to derive the transformations for Miras specifically (STScI Development Team 2013).

This resulted in an H to F160W transformation of

This is a significantly larger color term than was found for the bluer Cepheids (0.16) in R16. This difference is due to the nonlinearity of the color term, shown in Figure 14. Because Cepheids are much bluer than Miras (a typical Cepheid is an F star), and the color term is nonlinear, we must use a different color transformation for Miras. Using the same color term to transform from H to F160W for Miras as we used for Cepheids would result in a ∼0.18 mag difference for an average O-rich Mira with a J − H color of 0.8. The ground-based spectra are affected by the telluric absorption bands, as seen in Figure 15. Water bands dominate the NIR spectra of O-rich Miras, and these are difficult to separate from telluric features in ground-based observations.

Figure 14. Refer to the following caption and surrounding text.

Figure 14. Nonlinear color term for the F160WH transformation as a function of JH color. The blue points represent calculations using Castelli & Kurucz (2004) models with 3500 < T < 7000, log g = 0.1 and solar metallicity. The red points are based on observed spectra of O-rich Miras. Note that neither set of points has a constant color coefficient as a function of color. The mean of the red points is 0.39, which we adopt for our transformations.

Standard image High-resolution image
Figure 15. Refer to the following caption and surrounding text.

Figure 15. Reddest and bluest O-rich Mira spectra used to calculate the color correction (very red OH/IR stars were excluded) are shown in blue and black, respectively, with ground-based J and H filters and HST WFC3 F160W. The filters have been normalized to have a peak throughput of 1. The difference in color appears to be the result of a combination of both the continuum emission and the much stronger absorption lines in the redder spectrum.

Standard image High-resolution image

Because the O-rich Miras we used to calculate the color transformation displayed a large scatter in F160WH color, we also examined spectra of individual Miras that were particularly blue or red to check that the spectra were calibrated well enough for synthetic photometry. We found that very red Miras had much larger absorption features (the result of having more dust) than very blue Miras. Additionally, we compared the measured J − K colors of the O-rich Miras in the 2MASS catalog to the J − K colors derived synthetically (Figure 16). We found that for individual Miras, the two were in agreement and the standard deviation of the two distributions of color (0.19 from spectrophotometry and 0.14 from 2MASS measurements) was also similar. The 2MASS colors were, on average, slightly redder (by ∼0.1 mag), suggesting that the true color correction might be slightly larger. However, the distribution in Mira color appears to be real and not the result of poor calibration.

Figure 16. Refer to the following caption and surrounding text.

Figure 16. Comparison of the J − K color derived synthetically using PySynPhot and J − K colors from 2MASS. Most of the objects in the Galactic sample of Miras were observed at more than one part of their light cycle and thus had synthetic colors that varied by ∼0.15 mag. The synthetic colors had an overall standard deviation of 0.19 mag, whereas the 2MASS colors had a standard deviation of 0.14 mag.

Standard image High-resolution image

5.2. Mira Samples

We created three subsets of Mira candidates based on our estimate of their reliability as distance indicators and the possibility of contamination, applying the cuts discussed described in Table 3. A complete listing of our final sample of Miras is given in Table 6. For our Gold sample of Mira candidates, we had a total final sample size of 161 objects. Fitting these objects to the PL relation derived by Yuan et al. (2017a) for the H band gave us a zeropoint of a0 = 23.24 ±0.01 mag. The sigma clipping removed 22, ∼14% of the Gold sample Mira candidates, from the final PLR, leaving a total of 139 Mira candidates remaining.

For the larger Silver sample, we had a final sample size of 296 and determined a0 = 23.25 ± 0.01 mag for this sample. After sigma clipping, which removed 48 objects, ∼16% of the Mira candidates, leaving 248 remaining for the fit.

Table 6.  Final Sample of Miras

ID Period R.A. Decl. X Y Magnitude Amplitude Quality
  (Days) (J2000.0) (J2000.0) (Pixels) (Pixels) (F160W mag) F160W)  
60841 287.015 12 18 47.087 +47 19 32.85 1210.179 1249.445 22.029 0.776 Gold
47935 253.514 12 18 42.138 +47 20 42.89 165.583 983.280 22.918 0.684 Gold
44958 272.561 12 18 43.976 +47 20 33.73 418.105 920.864 22.790 0.712 Gold
24235 292.075 12 18 44.806 +47 21 10.89 200.187 497.320 22.455 0.716 Gold
31113 299.998 12 18 44.585 +47 20 58.19 280.781 636.955 22.539 0.796 Gold
52981 275.948 12 18 42.945 +47 20 25.17 386.517 1087.064 22.523 0.746 Bronze
64439 256.965 12 18 41.148 +47 20 15.63 288.277 1325.193 22.938 0.712 Gold
64343 298.885 12 18 40.925 +47 20 17.73 249.748 1323.194 22.792 0.758 Silver
33823 265.960 12 18 45.165 +47 20 47.61 422.330 690.870 22.681 0.654 Silver
39445 237.358 12 18 42.893 +47 20 54.83 143.229 807.290 22.897 0.798 Gold

Note. A partial list of Miras is shown here for information regarding form and content.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Finally, the Bronze sample consisted of 438 objects and had a zeropoint of a0 = 23.25 ± 0.01 mag. We sigma-clipped out 72 Mira candidates, ∼16%, comparable to the amount removed in the Silver sample. The PLR for each sample is shown in Figure 17.

Figure 17. Refer to the following caption and surrounding text.

Figure 17. Mira period–luminosity relations for the Gold, Silver, and Bronze subsamples (left, center, and right, respectively). Red points denote objects used in the final fit, while gray points represent variables that were removed through iterative 3σ clipping. The solid black curves show the best-fit relations, while the dashed lines denote the 1σ scatter (0.11, 0.13, 0.14 mag respectively). The functional forms of the PLRs are from Yuan et al. (2017a) and only the zeropoint was fit.

Standard image High-resolution image

Despite the different selection criteria, the zeropoints of all of the samples are almost the same and the fraction of their outliers is also consistent. This is especially important for the Bronze sample, which used only NIR criteria. It suggests that future studies can also be successful with only NIR HST data to select Miras.

5.3. Relative Distance to the LMC

We compared our results in NGC 4258 with a sample of O-rich Miras discovered in the LMC by the OGLE survey (Soszyński et al. 2009b) to determine the relative distance modulus between these two galaxies. We obtained random-phase J and H magnitudes for the LMC variables from the IRSF catalog of Kato et al. (2007). We used Equation (3) to transform these into the equivalent F160W magnitudes, which were fit with the PLR of Equation (3) to solve for the zeropoint. The difference between the NGC 4258 and LMC zeropoints yields the relative distance modulus (not corrected for small differences in foreground extinction between the two galaxies).

With the Gold sample defined in the previous section, we calculated a distance modulus relative to the LMC of Δμg = 10.95 ± 0.01r ± 0.06sys mag using Miras from the OGLE survey. For the Silver sample, we have Δμs =10.97 ± 0.01r ± 0.07sys mag, and for the Bronze sample, Δμb = 10.97 ± 0.01r ± 0.08sys mag. These are consistent with a previous measurement of the Cepheid relative distance modulus from R16, ΔμR16 = 10.92 ± 0.02 mag.

In order for the Cepheid scale to agree with the Planck results, it would need to be too short by ∼0.20 mag. Currently the Cepheids and Miras give consistent relative distances, but we can consider a hypothetical Mira relative distance to NGC 4258 of 10.75 ± 0.07 mag, in agreement with Planck. We find that it disagrees with the Cepheid relative distance from R16 by 2.4σ. This demonstrates that both the Cepheid and Mira distance scales have some tension with the Planck results. Because the Cepheid and Mira results are independent, we also would not expect these discrepancies with Planck to be caused by the same effect.

Finally, we also used the color transformation from Section 5.1 to derive the PLR coefficients for an F160W-band PLR using the ground J- and H-band relations from Yuan et al. (2017b). We then refit the PLRs for both the LMC and NGC 4258 and found that these two methods yield marginal differences in the results.

5.4. Absolute Calibration to NGC 4258

We use the improved megamaser distance to NGC 4258 from R16 of 29.387 ± 0.057 mag. The uncertainty in the Humphreys et al. (2013) value was reduced to 2.6% from 3% by increasing the number of Monte Carlo Markov Chain trial values in the analysis by a factor of 100. Using this distance modulus puts the absolute calibration of the PLR for the Gold sample at a0 = −6.15 ± 0.09 mag in the F160W bandpass.

5.5. Spatial Distribution

As a sanity check, we compared the spatial distributions of our Bronze Mira candidate sample with the spatial distribution of Cepheids in NGC 4258. Figure 18 shows the locations of Cepheids and Miras in the galaxy overlaid with the F160W footprints. The Cepheids trace the spiral arm of the galaxy, while the Miras are found randomly distributed in the F160W footprint. These differences in spatial distribution have a physical origin in the progenitors of Cepheids and Miras. Cepheids, with their intermediate and high-mass progenitors, are young stars that are only found in regions with active star formation. Thus, they are present in the denser spiral arms of a galaxy and are part of the disk population. Almost all Miras have progenitors of low-to-intermediate mass, which therefore have intermediate-to-old ages and can exist in areas without recent star formation. For our sample, which is limited to short-period Miras only, this is especially true since the progenitor stars will all be of low mass. Miras can additionally be found in both the disk and halo populations.

Figure 18. Refer to the following caption and surrounding text.

Figure 18. Left: the locations of the Bronze sample Miras (cyan circles) in the NGC 4258 ACS inner field. The white regions show the F160W footprint. Because our observations were taken over the course of one year, the orientation changed in each observation, leaving an approximately circular area that was observed in all epochs. We searched in this region only for Miras. Right: Cepheids (cyan circles) from Macri et al. (2006) on top of the same ACS field. The Cepheid distribution traces the spiral arms of the galaxy, while the Miras are more common and can be found evenly across the smaller F160W footprint.

Standard image High-resolution image

We calculated the autocorrelation functions for both Cepheids and Miras in the F160W footprint shown in Figure 18. The autocorrelation function for Miras discovered in this project and the Cepheids discovered by Macri et al. (2006) is shown in Figure 19. The Cepheid and Mira autocorrelation functions follow different distributions, with the Mira autocorrelation being much flatter, as expected for evenly distributed objects. The results confirm that the two distributions are different spatially, and agree with what we would expect from a physical understanding of Mira and Cepheid progenitors.

Figure 19. Refer to the following caption and surrounding text.

Figure 19. Cepheid–Cepheid and Mira–Mira standard correlation function. in the NGC 4258 F160W field. Errors are obtained through bootstrap resampling. Cepheids are blue, and the Miras are red. Due to the small sample of Cepheids (84) compared to Miras (438), the errors for the Cepheid–Cepheid autocorrelation functions are much larger.

Standard image High-resolution image

6. Discussion

The largest source of uncertainty in local measurements of H0 remains the number of SNe Ia host galaxies that have been calibrated with Cepheid distances. A 300-day Mira is roughly comparable in F160W brightness to a 30-day Cepheid, allowing them to be observed to approximately the same volume. SNe Ia used in R16 all have modern photometry, low reddening (AV < 0.5 mag), and observations prior to peak luminosity. Additionally, only late-type host galaxies were targeted to ensure the presence of Cepheids. Using Miras over Cepheids would increase the number of SN host galaxies for cross-calibration and eliminate potential biases caused by host galaxy morphology. Recent papers, such as Jones et al. (2015) and Rigault et al. (2015), have disagreed on whether host galaxy morphology can have an effect on the luminosities of SNe Ia. Rigault et al. (2015) found that SNe Ia in locally star-forming environments are dimmer than SNe Ia in locally passive environments. Jones et al. (2015) also searched for this effect but found that there was little evidence for a difference. Regardless of the effect the local environment can have, we can avoid this potential problem altogether by using Miras.

Since Miras are older population stars, they can be found in most galaxies regardless of host galaxy morphology. This can help us create a sample of cross-calibrators that is more representative of the Hubble flow SNe Ia sample. Rejkuba (2004) was able to derive a K-band Mira PLR for the giant elliptical galaxy NGC 5128, which would have been an unlikely target as an SNe Ia calibrator host. Miras are also part of the halo population so we can also potentially look for Miras in hosts that are not face-on, further increasing the number of potential targets.

Below 300 days, the proportion of O-to-C-rich Miras increases as the period decreases. The host galaxy's metallicity and IMF will also affect the relative proportions of O- to C-rich Miras, but most Miras with periods less than 300 days will be O-rich up to about Fe/H ∼ −1.0. In the SMC, with a Fe/H ∼ −1.0, even the short-period ranges (<250 days) are dominated by C-rich stars, as shown in Table 7. We would be able to find more O-rich Miras and other variables that fall on the sample PLR by searching for objects with periods below 100 days. At longer wavelengths, C-rich stars can serve as distance indicators as well. Whitelock et al. (2017) compares the O-rich and C-rich Mira PLRs in Ks. This would not work for HST, since there are no filters redder than F160W but may be possible for JWST.

Table 7.  LMC and SMC Miras

Period O-rich SMC C-rich SMC O-rich LMC C-rich LMC
≤300 days 12 80 349 187
≤250 days 7 18 302 17

Note. The numbers of O- and C-rich Miras in the SMC (Soszyński et al. 2011) and LMC (Soszyński et al. 2009b) from the OGLE survey.

Download table as:  ASCIITypeset image

In the future, we will be able to further reduce our uncertainties in the relative distance by directly comparing measurements of Miras observed in the HST NIR with our calibrated PLRs from NGC 4258. In addition, by using the same criteria and filters, we can help ensure that we selected for the same classes of objects. The consistency of the results between the Bronze, Silver, and Gold samples also demonstrates that we can conduct this search without optical data.

7. Conclusions

  • 1.  
    We discovered 438 Mira candidates in one field of NGC 4258 using only the HST F160W bandpass.
  • 2.  
    We developed criteria to reduce contamination from C-rich Miras in our sample that do not use spectroscopic information and found that amplitude cuts can help separate C- and O-rich Miras. This allows us to discover and characterize Mira candidates using solely or primarily F160W.
  • 3.  
    We determined a relative distance modulus between NGC 4258 and the LMC based on Mira variables of Δμg =10.95 ± 0.01r ± 0.06sys mag, which is consistent with the Cepheid relative distance modulus and also consistent with the relative distance modulus obtained using geometric methods.
  • 4.  
    We have calibrated the extragalactic Mira distance ladder in F160W using the geometric distance to NGC 4258.
  • 5.  
    We derived a Mira-specific color transformation from the ground-based H-band to the HST F160W.

This project was funded through HST grants associated with the following programs: GO-13445, GO-11570, and GO-9810. We acknowledge with thanks the oCeti observations from the AAVSO International Database contributed by observers worldwide and used in this research. P.A.W. acknowledges a research grant from the South African National Research Foundation. We would also like to thank Drs. Yoshifusa Ita, Noriyuki Matsunaga, Gregory Sloan, Martha Boyer, Steven McDonald, Daniel Shafer, and Sjoert van Velzen for their helpful discussions.

10.3847/1538-4357/aab6b3
undefined