The 3D Distribution of Long-period Mira Variables in the Galactic Disk^*

We selected bright infrared sources that satisfy the following criteria as candidates of the long-period Mira variables from the IRAS PSC (Neugebauer et al. 1984). First, sources with decl. > −25° were selected considering the latitude of our observatory site, Kagoshima, located at a latitude of 315 north. Second, targets needed to be bright in all bands; 12, 25, and 60 μm (flux density quality = 3 in the IRAS PSC). Third, targets are required to be located in the IIIa or IIIb regions on the IRAS color–color diagram (van der Veen & Habing 1988). The position of the IIIa and IIIb regions are indicated in Figure 1, together with the distribution of the 108 Mira variables examined in this paper. The IIIa and IIIb regions lie around the most evolved phase of the sequence of dust shell evolution of AGB (van der Veen & Habing 1988). Therefore, the stars located in these regions are likely to have a longer period. There are 2300 for the IIIa and 280 for IIIb, which satisfy the above three criteria.

Zoom In Zoom Out Reset image size

Figure 2 shows the l − b diagram of sources satisfying the criteria including those for which we report distances in this paper. There is a gap around l = 80°. There are few IRAS sources in this area because IRAS observation was not performed in this area. In this paper, we report on 108 sources for which accurate periods were determined so far.

Zoom In Zoom Out Reset image size

Our monitoring observations were carried out at Iriki observatory operated by Kagoshima University (Chibueze et al. 2016), Japan, since 2003. A near-IR camera with a 512 × 512 HgCdTe array detector attached to a 1 m telescope yielding 538 × 538 field of view (063 pixel scale) was used. The typical seeing measurement at the Iriki observatory is 15. Monitoring observations in the K' band were performed in the fields centered on the position of the selected IRAS sources. Each set of observations consists of five exposures at slightly dithered positions with an exposure time of 1–12 s. We adjusted the exposure time depending on the sky background level, and the brightness of the target and reference stars in the same field. Stars brighter than K' ∼ 6.5 mag would be saturated at the best focus position, and therefore we observed such bright sources with the focus off. For targets for which no reference star was found in the same field of view, we also observed the standard stars listed in Elias et al. (1982) on the same night. We searched for the counterparts of IRAS objects in the 2MASS image. In many cases, the bright near-infrared counterparts were found near the IRAS position. The coordinates of the 108 sources in Table 1 are based on the 2MASS catalog (J2000.0).

Raw images were reduced using our automatic image reduction pipeline that utilizes tasks from IRAF. The procedure has five main steps as follows. First, dark images were subtracted from the raw image to remove the effect of dark current. Dark current subtraction was made using the 3σ clipping of a combination of 10 dark frames with exposure times equal to those of the observations carried out on the same nights as the observed targets. Second, we standardized the raw image using a flat field to correct for the nonuniformity of sensitivity of our detector. Flat fields were obtained using twilight-flats from 2003 to 2012. From 2013, operations moved to using dome-flats. We obtained 10–30 images of them every week and calculated the average of difference between on-flat and off-flat. Third, we performed sky subtraction to reduce the fixed pattern of the raw images due to foreground or background radiation. The sky-frame was subtracted from the raw image after standardization to correct for the sky value. We calculated the sky value from starless area. For the on-focus images, we subtract the self-sky images, which were a median combined image without positional-shift. The blank sky fields adjacent to the target field were observed separately for out of focus images, and their median combined images were also made using the same approach as that used for the on-focus images. Fourth, we combined the slightly dithered images into a single image. We measure the (x, y) position of the same stars in the dithered images. We calculated the positional shifts among the dithered images. We combined them using the imcombine task in IRAF based on the above information. Fifth, we inserted the World Coordinate System (WCS). We obtained astrometric solutions by comparing detected stars with the 2MASS PSC (Skrutskie et al. 2006) and then inserted the WCS into the FITS file. We searched for the counterparts of the stars in our images in the 2MASS PSC. Pixel coordinates (x, y) of the detected sources were converted to equatorial coordinates (R.A., decl.) using a linear transformation. The equatorial coordinates were based on the International Celestial Reference System via the 2MASS-PSC.

We performed aperture photometry for the Mira candidates and reference stars on the same images with the APPHOT package in IRAF. For each image, reasonably bright stars with variation of less than 0.1 mag were used as reference stars. The instrumental magnitude was calibrated by comparing them to the 2MASS PSC using reference stars. For the images without any reference stars in the same field or defocused images, we calibrated the target using the standard stars in Elias et al. (1982).

We performed a period search employing a least-square method to fit sinusoidal curves to the photometric measurements of all samples in our monitoring survey. The light curve of IRAS 18511-1041 is illustrated in Figure 3 as an example and those for all the sources are given as online material. We fit our photometric data to the following equation,

$$\begin{eqnarray}&&m=B+(A/2)\sin (2\pi t/P+\phi )\end{eqnarray} \tag{ 1 }$$

where m is the observed magnitude, B is the center value of the sinusoidal, A is the peak to peak amplitude, P is the period, t is Modified Julian Date (MJD), and ϕ is the phase at t = 0.

Zoom In Zoom Out Reset image size

We searched for a period that minimizes the rms residual with a search window of 100–1500 days. Then we visually inspected the folded light curve, the fitted sinusoid and photometric data plotted against phase. As a result, we determined the period for 108 sources (Table 1). In addition, we use the average of the matching sinusoidal curve (B in Equation (1)) as the mean magnitude.

A histogram of the derived periods of our monitoring samples is given in Figure 4, which shows a peak between 500 and 600 days. Generally, populations of Mira variables peak at around 300 days, but our Mira candidate population shows a peak at around 500 days. This is as expected because our samples selected from the IIIa and IIIb regions of the IRAS color–color diagram. As explained in van der Veen & Habing (1988) the various regions of the IRAS color–color diagram host stars at different evolutionary stages and mass ejections; the stars belonging to IIIa and IIIb have more developed dust shells. Therefore, we were able to demonstrate that the stars located in these groups have longer periods.

Zoom In Zoom Out Reset image size

Another period distribution of long-period Mira variables in the Galactic bulge is given in Whitelock et al. (1991), whose sample was also selected from the IRAS point-source catalog although the selection criteria was different from ours. Their population of Mira variables had a peak at around 400 days. Therefore our sample represents a population of Mira variables that have longer periods than typical Mira variables and the sample of Whitelock et al. (1991). We found that the selection of sources in the IIIa and IIIb regions is a good method to isolate longer period Mira variables, which most likely have large initial masses.

There are two types of long-period variables; Mira variables and semi-regular (SR) variables, and our 108 Mira candidates could be contaminated by SR variables. While both classes of variables can have similar periods, their pulsation amplitudes typically differ. Mira variables typically have pulsation amplitudes of K' > 0.4 mag, while SR variables typically have smaller amplitudes of K' < 0.4 mag (e.g., Whitelock et al. 2000). Using this criterion we inspected the pulsation amplitudes of the Mira candidates, which are shown in Figure 5. The amplitude of IRAS 19494+2939 is too small to support the classification as a Mira variable, and this object is not included in the following discussions. The longer period Mira variables have larger amplitudes. To further test the validity of our Mira candidates selection we plotted K' band amplitudes against periods (Figure 6). As expected, the Mira candidate population shows a correlation; longer periods for larger pulsation amplitudes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The PLR of Mira variables was first established in the LMC and a recent evaluation of the PLR in the LMC was given in Ita & Matsunaga (2011). Since the PLR enables prediction of a star's absolute magnitude it is possible to calculate the distance modulus of Mira variables by combining measurements of the pulsation period and an extinction corrected apparent magnitude. We aim to estimate the distances of our long-period Mira variables using the PLR to reveal their distribution in our galaxy. The PLR of Mira variables at infrared wavelengths is a useful tool for distance determination within the Milky Way and at larger distances, since there is less extinction at these wavelengths.

The PLR at shorter wavelengths becomes imprecise for Mira variables with periods longer than about 300 days. C-rich Mira variables locate below the extension of the PLR along the period in the optical and JHK_s bands (Ita & Matsunaga 2011). On the other hand, O-rich Mira variables with longer periods locate above the same extended PLR in the optical and JHK_s bands. Some O-rich Mira variables with very long periods (log P > 3.0) are found below the extended PLR, probably due to circumstellar extinction, similarly to the C-rich Mira variables in the same paper. As a result of these deviations, the PLR at the K band or shorter wavelengths cannot be reliably used for the determination of distances. On the contrary, the 3.6 μm PLR in the same paper shows a tight and linear relation for both short and long-period members in the Spitzer 3.6 μm band, although it is based on single epoch data. This suggests that circumstellar extinction does not severely affect the PLR of longer period Mira variables at around 3 μm and thus can be used to more reliably determine distances.

In this paper, we combine the pulsation period of our monitoring observations we determined for the 108 long-period Mira variables with magnitudes from WISE (Wright et al. 2010) band I (3.4 μm) to estimate the distances to our sample via the 3.4 μm PLR. Our methodology is to first establish the 3.4 μm PLR. Second, we derive the interstellar extinction at 3.4 μm using a 3D reddening map. Finally, we derive the interstellar extinction corrected 3.4 μm magnitudes from the WISE data and use these to derive distance moduli for our sample.

In the present case, WISE is preferable to Spitzer since the latter does not contain data for all of our sources.

To establish the PLR at 3.4 μm we use the photometries of 1663 Mira variables in the LMC, for which the periods were determined by OGLE-III (Soszyński et al. 2009). Since these sources are located in the LMC, we assume they reside equidistantly and the distance modulus of the LMC is 18.5. Figure 7 shows the apparent magnitudes of W1 from ALLWISE (Cutri et al. 2014) against their period for 1663 Mira variables listed in the OGLE-III catalog. We can see no significant difference in the distribution between C-rich and O-rich Mira variables. We performed least-square fitting to all sources regardless of C-rich and O-rich with a linear function. A least-square fitting to these data produces the following relation.

$$\begin{eqnarray}&&{m}_{W1}=-4.48\pm 0.05\times \mathrm{log}\,P+21.39\pm 0.13\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{M}_{W1}=-4.48\pm 0.05\times \mathrm{log}\,P+2.89\pm 0.13,\end{eqnarray} \tag{ 3 }$$

where M_W1 is the absolute magnitude converted by the distance modulus of the LMC, and m_W1 is the apparent magnitude in W1. Figure 8 is a histogram of the residual from the best-fit PLR. We use the width of the Gaussian, σ = 0.29 mag, as an error source in the following analysis. The 3.4 μm PLR is likely due to the balance of extinction and thermal radiation of circumstellar dust, and the two cases; Mira variables with and without the circumstellar dust, share the same position in the PL diagram.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Riebel et al. (2010) also reported the 3.6 μm PLR of Mira variables in LMC using a similar data set, the MACHO (MAssive Compact Halo Objects survey) and the Spitzer SAGE catalog. They obtained the equations of PLR for the O-rich and C-rich Mira variables separately as well as for all sources including both the O-rich and C-rich Mira variables. The PLRs of O-rich and C-rich Mira variables give m_W1 = 9.683 mag and 9.171 mag at log P = 2.7, and they differ from our relation (Equation (2)) by 0.123 mag and 0.389 mag, respectively. The PLR of all the sources replies m_W1 = 9.156 mag at log P = 2.7 and their PLR differs from our relation (Equation (2)) by 0.138 mag. The difference for all the sources PLR is smaller than the residual of our fitting and is not significant considering the slight difference of wavelength (3.4 μm or 3.6 μm) and the difference of fitting range of period. In this paper, our purpose is to obtain distances for long-period Mira variables in Milky Way and so it is difficult to classify as O-rich or C-rich correctly for our samples and hence difficult to apply the PLR dedicated to the O-rich or C-rich Mira variables even if we know the equations of PLR for C-rich and O-rich separately. Therefore, we use our equation of PLR obtained from both C-rich and O-rich without discrimination in this paper.

We estimate distances to Mira variables using the 3.4 μm PLR. The apparent distance modulus at 3.4 μm is related to the true distance modulus μ and the extinction A_W1 at 3.4 μm as follows.

$$\begin{eqnarray}&&\mu ={m}_{W1}-{M}_{W1}-{A}_{W1}={\mu }_{0}-{A}_{W1}\end{eqnarray} \tag{ 4 }$$

Here, μ₀ (=m_W1 –M_W1) is the apparent distance modulus without correction of interstellar extinction. Our sample, consisting of longer period Mira variables, resides in the Galactic plane. Thus, corrections must be derived for interstellar extinction at 3.4 μm as described below.

We determined A_W1 using a 3D map of interstellar dust reddening, covering three-quarters of the sky (decl. > −25°) out to a distance of several kiloparsecs. (see http://argonaut.skymaps.info/; Green et al. 2018, 2015) They determined the interstellar reddening (E(B − V)) as a function of distance modulus for each line of sight using 800 million stars from Pan-STARRS 1 and 2MASS. Here, we use a universal dust extinction law with R(V) = 3.1 and ${A}_{L \mbox{-} \mathrm{band}}$/E(B − V) = 0.157 to convert E(B − V) to A_W1 (Schlafly & Finkbeiner 2011). On the other hand, interstellar extinction increases against distance modulus and A_W1 is expressed as a monotonically increasing function (f(μ)) of μ. This function is provided as the "best-fit line" from the 3D reddening map. We can estimate the interstellar extinction and the distance modulus by solving the simultaneous equations of Equation (4) and f(μ). We show an example for IRAS 01304+6211 in the Figure 9. We used an error shown in the 3D reddening map. We summarize the A_W1 values, the line-of-sight distances (D), and the height distances (z) in Table 1. We estimated the errors of μ using the error of PLR, the error of A_W1 and the WISE photometric error.

Zoom In Zoom Out Reset image size

Based on the distances derived from our monitoring survey we illustrate the Galactic face-on distribution of all of our Mira variables in Figure 10. To investigate the possibility that longer period Mira variables are associated with the spiral arm structures in the Milky Way, we compare our sample with the arm structure model derived from the distribution of HMSFRs (Reid et al. 2014). We will focus on the younger population, i.e., the subsample with 2.7 < log P < 3.0 and $| Z|$ < 300 pc, as shown in Figure 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The distribution shows that longer period Mira variables are located near the Sagittarius arm (red line) and the Scutum arm (blue line). Figure 11 shows that our sample concentrates to 40° < l < 60°. This direction corresponds to the tangential direction of the Sagittarius arm and therefore the sources tend to be on the arm regardless of their distance. To investigate this effect, we examined the distance distribution of our sample in this direction. Figure 12 presents a histogram of the distances of sources located in 40° < l < 60° in Figure 11. There are 27 sources with 2.7 < log P < 3.0 in 40 < l < 60 and $| Z|$ < 300 pc, in which 23 sources (85%) are distributed in the range of 3 kpc and 8 kpc, corresponding to the Sagittarius arm. It is suggested to extend through 2 < D < 8 kpc in the 40° < l < 60° (Reid et al. 2014). Another two stars are located in the range of 9 < D < 12 kpc, corresponding to the Perseus arm (Reid et al. 2014). The more distant sample is likely to be incomplete because of the IRAS sensitivity limit.

Zoom In Zoom Out Reset image size

The present result suggests that the distribution of longer period Mira variables probably trace the spiral arm structures of the Milky Way. In order to confirm our suggestion to the entire galaxy, it is necessary to search for more longer period Mira variables in other arms, for example, the Perseus arm and Outer arm.

The edge-on distribution of our sample is shown in Figure 13. Approximately 50% of our sample is found to exist within the Galactic plane or $| Z|$ < 0.3 kpc, corresponding to the Galactic thin disk where metallicity is high, and star formation is active. We can see that the longer period Mira variables belonging to IIIa and IIIb regions in the IRAS two-color diagram are packed in the thin Galactic disk. Almost all the remaining 50% of sources are contained within the thick disk ($| Z|$ < 1.5 kpc).

Zoom In Zoom Out Reset image size

One source (IRAS 15106-1532) with a relatively shorter period (268 days) is found at 6.35 kpc from the Galactic plane. There are three possibilities of its high $| Z|$ position; (1) it was born at the Galactic plane and ejected outward, (2) born in the stream of local dwarf spheroidal galaxies (dsph) (Huxor & Grebel 2015), and (3) born in the halo of the Milky Way. In the case of (1), the stars born in the Galactic plane are randomized according to the rotation of the Milky Way over a long period of time. AGB stars like short-period Mira variables are widely distributed to the thick disk. However, it is difficult to explain IRAS 15106-1532 as having an origin in the galactic plane although it might just be an unexpected run-away star. It is necessary to create an accurate model of this moment; however, this will require taking into account the motion of stars or sudden phenomena such as blow-off from heavy stars in the Milky Way, so it is not possible to do this at the present moment. In the case of (2) stream or dsph, stars born in the dsph galaxies are pulled off by tidal force and a stream is formed. Some Mira variables are found both in the stream and in dsph galaxies (Mauron 2008; Sakamoto et al. 2012). However, the position of IRAS 15106-1532 does not match the distribution of the already known stream; therefore, IRAS 15106-1532 does not have a stream or dsph origin. In the case of (3), Whitelock (1990) reported that some Mira variables are born in a globular cluster in the halo, which are known to have short periods. It is possible to infer the origin of IRAS 15106-1532 in globular clusters in the halo.